Help to create a 2D mesh generator (FEM)
Mathematica Asked on August 22, 2021
I’m working in a finite element mesh generator. I built this function which generates an 8 node mesh (polynomials of order 2) without any interior node:
(*Generate Grid Mesh of dimensions axb with nx divisions in x and ny
divisions in y*)
GenerateGridMesh[aa_, bb_, nx_, ny_, order_] :=
Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j,
meshtopology = {}, allcoords, k, topolsz, l, data, c, a, b},
k = 0;
meshnodes = {};
dx = aa/(2 nx);
dy = bb/(2 ny);
For[i = 1, i <= 2 ny + 1, i++,
If[OddQ[i] == True,
For[j = 1, j <= 2 nx + 1, j++,
AppendTo[meshnodes, {x, y}];
x += dx ;
];
,
For[k = 1, k <= nx + 1, k++,
AppendTo[meshnodes, {x, y}];
x += 2 dx ;
];
];
x = 0;
y += dy;
];
meshtopology = {};
b = 0;
a = 1;
l = 0;
c = 3 nx + 2;
For[i = 1, i <= ny, i++,
For[j = 1, j <= nx, j++,
data = {a, a + 2, 3 nx + 4 + a, 3 nx + 3 + b, a + 1,
2 nx + 3 + l, 3 nx + 4 + b, 2 nx + 2 + l};
AppendTo[meshtopology, data];
a += 2;
b += 2;
l += 1;
];
l = 3 nx + 2 + c (i - 1);
a = 3 nx + 3 + c (i - 1);
b = 3 nx + 2 + c (i - 1);
];
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, Length[meshtopology[[1]]]}];
{allcoords, meshnodes, meshtopology}
];
(*Generates graphics to visualize mesh and nodes*)
GenerateGraphics[nodes_, topology_, order_] :=
Block[{meshvis, nodevis, v},
If[order == 1, v = {1, 2, 3, 4}, v = {1, 5, 2, 6, 3, 7, 4, 8}];
meshvis =
Graphics[{FaceForm[], EdgeForm[Black],
GraphicsComplex[nodes, Polygon[topology[[All, v]]]]}];
(*nodevis=Graphics[{MapIndexed[Text[#2[[1]],#1,{-1,1}]&,
nodes],{Blue,Point[nodes]}}];*)
nodevis =
Graphics[{MapIndexed[
Style[Text[#2[[1]], #1, {-1.8, 1.8}], FontSize -> 9] &,
nodes], {PointSize[Large], Black, Point[nodes]}}];
{meshvis, nodevis}
];
L = 5;
h = 5;
nx = 2;
ny = 2;
order = 2;
{allcoords, meshnodes, meshtopology} =
GenerateGridMesh[L, h, nx, ny,
order];(*Generate finite element mesh*)
{meshvis, nodevis} =
GenerateGraphics[meshnodes, meshtopology,
order];(*Generates graphics to visualize mesh*)
Show[meshvis, nodevis, AspectRatio -> Automatic, ImageSize -> Large]
which results in the following mesh:
I want to build a generic mesh generator for any polynomial order. Here is an example of what I need:
L = 5;
h = 5;
x = 0;
y = 0;
nx = 2;
ny = 2;
order = 3;
meshnodes = {};
dx = L/(nx order);
dy = h/(ny order);
For[irow = 1, irow <= order nx + 1, irow++,
For[icol = 1, icol <= order ny + 1, icol++,
AppendTo[meshnodes, {x, y}];
If[OddQ[Mod[irow, 3]] == True,
x += dx ;
,
x += 3 dx ;
icol += 2;
];
];
y += dy;
x = 0;
];
meshtopology = {{1, 4, 17, 14, 2, 9, 16, 11, 3, 12, 15, 8}, {4, 7, 20,
17, 5, 10, 19, 12, 6, 13, 18, 9}, {14, 17, 30, 27, 15, 22, 29,
24, 16, 25, 28, 21}, {17, 20, 33, 30, 18, 23, 32, 25, 19, 26, 31,
22}};
{meshvis, nodevis} =
GenerateGraphics[meshnodes, meshtopology,
order];(*Generates graphics to visualize mesh*)
Show[meshvis, nodevis, AspectRatio -> Automatic, ImageSize -> Large]
I need this to be created automatically for any dimensions of L and h, and for any node quantity.
Below is an example of a code that generates a curved mesh (nine noded elements, not serendipity).
GenerateGridMesh[R0_, RE_, nx_, ny_, order_] :=
Block[{x = 0., y = 0., dx, dy, meshnodes, i, j, meshtopology = {},
allcoords, k, l},
meshnodes = {};
k = 0;
(*meshnodes=Flatten[Table[Table[{R Cos[[Theta]],
R Sin[[Theta]]},{R,R0,RE,(RE-R0)/(nx order-2)}],{[Theta],0,Pi/2,
Pi/2 /(ny order-2)}],1]//N;*)
r = (RE/R0)^(1/(-2 + nx order));
meshnodes =
Flatten[Table[
Table[{ R0 r^(n - 1) Cos[[Theta]],
R0 r^(n - 1) Sin[[Theta]]}, {n, 1.,
nx order - 1}], {[Theta], 0, Pi/2, Pi/2 /(ny order - 2)}],
1] // N;
k = 0;
For[i = 1 , i < ny, i++,
l = 1;
For[j = 1, j < nx, j++,
(*AppendTo[meshtopology,{j+k,j+2+k,4 nx+j+k,4 nx-2+j+k,j+1+k,j+1+
nx 2+k,j+nx 4-1+k,2 nx+ j-1+k,2 nx+ j+k}];*)
AppendTo[
meshtopology, {l + k, l + 2 + k, 4 nx + l + k, 4 nx - 2 + l + k,
l + 1 + k, l + 1 + nx 2 + k, l + nx 4 - 1 + k, 2 nx + l - 1 + k,
2 nx + l + k}];
l += 2;
];
k += 4 nx - 2;
];
If[order == 2,
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, 9}];
,
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, 4}];
];
{allcoords, meshnodes, meshtopology}
]
GenerateGraphics[nodes_, topology_, order_] :=
Block[{meshvis, nodevis},
If[order == 2,
meshvis =
Graphics[{FaceForm[], EdgeForm[Blue],
GraphicsComplex[nodes,
Polygon[topology[[All, {1, 5, 2, 6, 3, 7, 4, 8}]]]]}];
,
meshvis =
Graphics[{FaceForm[], EdgeForm[Blue],
GraphicsComplex[nodes,
Polygon[topology[[All, {1, 2, 3, 4}]]]]}];
];
nodevis =
Graphics[{MapIndexed[Text[#2[[1]], #1, {-1, 1}] &, nodes], {Blue,
Point[nodes]}}];
{meshvis, nodevis}
];
interpolatingQuadBezierCurve[pts_List] /; Length[pts] == 3 :=
BezierCurve[{pts[[1]], 1/2 (-pts[[1]] + 4 pts[[2]] - pts[[3]]),
pts[[3]]}];
interpolatingQuadBezierCurve[ptslist_List] :=
interpolatingQuadBezierCurve /@ ptslist;
interpolatingQuadBezierCurveComplex[coords_, indices_] :=
interpolatingQuadBezierCurve[Map[coords[[#]] &, indices]]
GenerateGraphics[nodes_, topology_] := Block[{meshvis, nodevis},
nodevis =
Graphics[{MapIndexed[
Style[Text[#2[[1]], #1, {-1.8, 1.8}], FontSize -> 12] &,
nodes], {PointSize[Large], Black, Point[nodes]}}];
Show[nodevis]]
order = 2;
serendipity = False;
{allcoords, nnodes, topol} = GenerateGridMesh[100, 200, 5, 4, order];
linestopology = Flatten[Table[
{{topol[[i]][[1]], topol[[i]][[5]], topol[[i]][[2]]},
{topol[[i]][[2]], topol[[i]][[6]], topol[[i]][[3]]},
{topol[[i]][[3]], topol[[i]][[7]], topol[[i]][[4]]},
{topol[[i]][[4]], topol[[i]][[8]], topol[[i]][[1]]}
}, {i, 1, Length[topol]}], 1];
Show[GenerateGraphics[nnodes, topol],
Graphics[interpolatingQuadBezierCurveComplex[nnodes, linestopology]],
ImageSize -> Automatic]
One Answer
I am not sure, if this fully answers your question, but you should be able to work from here. I tried to explain my code with the comments above each for-loop and I maintained your overall structure.
(*Generate Grid Mesh of dimensions axb with nx divisions in x and ny
divisions in y*)
GenerateGridMesh[aa_, bb_, nx_, ny_, p_] :=
Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j,
meshtopology = {}, allcoords, k, topolsz, l, data, c, a, b}, k = 0;
meshnodes = {};
(*determine the distance between each node*)
dx = aa/(p nx);
dy = bb/(p ny);
(*Generate node coordinates,
meshnodes should contain (p nx+1)(ny+1)+(p
ny+1)(nx+1)-(nx+1)(ny+1) nodes*)
For[hl = 0, hl < ny, hl++, (*loop over ny horizontal lines,
the last one is done below, after the loop *)
For[hln = 0, hln < p nx + 1,
hln++, (*loop over p nx+1 nodes on the horizontal line*)
AppendTo[meshnodes, {hln dx, dy p hl}];
];
For[vl = 1, vl <= p - 1,
vl++, (*loop over the p-1 horizontal "lines" that are not part
of the mesh so we can assign the coordinates to the nodes on the
vertical lines*)
For[vln = 0, vln < nx + 1,
vln++,(*loop over the nx+1 nodes on the vertical lines*)
AppendTo[meshnodes, {vln p dx , (p hl + vl) dy}];
];
];
];
(*Now for the last horizontal line,
note that we need to set hl to ny,
since we started counting from 0:*)
hl = ny;
For[hln = 0, hln < p nx + 1,
hln++, (*loop over p nx+1 nodes on the horizontal line*)
AppendTo[meshnodes, {hln dx, dy p hl}];
];
(*generate the list of cells - each cell is a list of node-
IDs that are on its border.
We have nx ny cells, each cell has 4p nodes*)
meshtopology = Table[{}, nx ny];
(*label the cells (cx,cy), i.e. (0,0),(1,0,1),...,(nx-1,0),(0,
1),...,(nx-1,ny-1) etc.*)
For[cy = 0, cy < ny, cy++,
For[cx = 0, cx < nx, cx++,
(*bottom edge*)
For[i = 0, i < p + 1, i++,
AppendTo[meshtopology[[cx + cy nx + 1]],
i + cx p + cy ((nx + 1) (p - 1) + p nx + 1) + 1];
];
(*right edge*)
For[i = 0, i < p - 1, i++,
AppendTo[meshtopology[[cx + cy nx + 1]],
cx + 1 + i (nx + 1) + (cy + 1) (p nx + 1) +
cy (p - 1) (nx + 1) + 1];
];
(*top edge, from right to left*)
For[i = p, i >= 0, i--,
AppendTo[meshtopology[[cx + cy nx + 1]],
i + cx p + (cy + 1) ((nx + 1) (p - 1) + p nx + 1) + 1];
];
(*left edge,from top to bottom*)
For[i = p - 2, i >= 0, i--,
AppendTo[meshtopology[[cx + cy nx + 1]],
cx + i (nx + 1) + (cy + 1) (p nx + 1) +
cy (p - 1) (nx + 1) + 1];
];
];
];
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, Length[meshtopology[[1]]]}];
{allcoords, meshnodes, meshtopology}
];
(*Generates graphics to visualize mesh and nodes*)
GenerateGraphics[nodes_, topology_, p_] := Block[{meshvis, nodevis, v},
If[order == 1,
v = {1, 2, 3, 4},
v = Table[i, {i, 1, 4 p}];
];
meshvis =
Graphics[{FaceForm[], EdgeForm[Black],
GraphicsComplex[nodes, Polygon[topology[[All, v]]]]}];
(*nodevis=Graphics[{MapIndexed[Text[#2[[1]],#1,{-1,1}]&,
nodes],{Blue,Point[nodes]}}];*)
nodevis =
Graphics[{MapIndexed[
Style[Text[#2[[1]], #1, {-1.8, 1.8}], FontSize -> 9] &,
nodes], {PointSize[Large], Black, Point[nodes]}}];
{meshvis, nodevis}
];
There is really not that much complicated stuff in my solution, I just used some basic for loop indexing wizardry.
When trying to understand my solution, you should note, that I started counting from 0 for all my indices and added 1 at the very end where necessary. This is in part because I am very experienced in C++ (which starts counting from 0 while Mathematica counts from 1), but also because we have many multiplications where it is usefull to have the zeroth row. I recommend making a general sketch of the mesh you want of order $p$ and count the nodes for each horizontal line, for each cell, etc. That way you should arrive at the same equations that I found.
Note, that I changed the ordering of the nodes in each cell. I am going counter-clockwise around the cell, starting at the bottom left corner. I achieved this in part by letting the For-loops for the top and left edge "run backwards".
I am sure that one could work more with Table and similar commands, but since you used For-loops, I did too.
Using the functions like so:
L = 12;
h = 10;
nx = 4;
ny = 5;
order = 6;
{allcoords, meshnodes, meshtopology} =
GenerateGridMesh[L, h, nx, ny,
order];(*Generate finite element mesh*){meshvis, nodevis} =
GenerateGraphics[meshnodes, meshtopology,
order];(*Generates graphics to visualize mesh*)Show[meshvis,
nodevis, AspectRatio -> Automatic, ImageSize -> Large]
yields the following image
I am also not sure, what you meant with the dimensionality of L and h, but since your examples are all 2-D meshes, I assumed that they are only lengths (and hence have basically no influence on the mesh apart from the distance between the nodes).
Correct answer by Andrea on August 22, 2021
Add your own answers!
Ask a Question
Get help from others!
Recent Questions
- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?
Recent Answers
- Jon Church on Why fry rice before boiling?
- Peter Machado on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?
- haakon.io on Why fry rice before boiling?