3
$\begingroup$

The function GenerateGridMesh[a_, b_, nx_, ny_] creates a grid of 9-node finite elements mesh. I need some help to adapt this routine to create 8 node element mesh. This function has as input data the a and b which are the grid dimensions, and nx and ny representing the number of divisions of the mesh in the x and y directions. As an output this function return allcoords representing the mesh nodes organized by elements, meshnodes that are all the nodes in the mesh ordered in sequence and, the meshtopology, that contains the ids of the nodes that composes the elements. The meshtopology order is very important and is needed in finite element global vector and matrix assemblage. As a simple example consider the single 9node element illustrated bellow. The meshtopology for this example must be ordered as follows: {{1,3,9,7,2,6,8,4,5}}. This is a 9-node element.

enter image description here

I need the same for a 8-node element like the following:

enter image description here

meshtopology I need: {{1,3,14,12,2,9,13,8}}

The function GenerateGraphics[nodes_, topology_, order_] is only for plot purpose.

(*a and b are the mesh dimensions*)
(*xn and ny are the number of divisions in the x and y directions*)

GenerateGridMesh[a_, b_, nx_, ny_] := 
  Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j, 
    meshtopology = {}, allcoords, k, topolsz, l, order = 2},
   k = 0;
   meshnodes = 
    Flatten[Table[
      Table[{x, y}, {x, 0, a, a/(nx order - order)}], {y, 0, b, 
       b/(ny order - order)}], 1];
   For[i = 1, i < ny, i++,
    l = 1;
    For[j = 1, j < nx, j++,
     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;
    ];
   allcoords = 
    Table[meshnodes[[meshtopology[[i, j]]]], {i, 1, 
      Length[meshtopology]}, {j, 1, Length[meshtopology[[1]]]}];
   {allcoords, meshnodes, meshtopology}
   ];
(*This function generates a mesh graphic*)
GenerateGraphics[nodes_, topology_, order_] := 
  Block[{meshvis, nodevis, v},
   v = {1, 5, 2, 6, 3, 7, 4, 8};
   meshvis = 
    Graphics[{FaceForm[], EdgeForm[Blue], 
      GraphicsComplex[nodes, Polygon[topology[[All, v]]]]}];
   nodevis = 
    Graphics[{MapIndexed[Text[#2[[1]], #1, {-1, 1}] &, nodes], {Blue, 
       Point[nodes]}}];
   {meshvis, nodevis}
   ];
{allcoords, meshnodes, meshtopology} = GenerateGridMesh[36, 36, 4, 4];
{meshvis, nodevis} = GenerateGraphics[meshnodes, meshtopology, 2];
Show[meshvis, nodevis]

This what the code return:

enter image description here

the meshtopology i get:

{{1, 3, 17, 15, 2, 10, 16, 8, 9}, {3, 5, 19, 17, 4, 12, 18, 10, 11}, {5, 7, 21, 19, 6, 14, 20, 12, 13}, {15, 17, 31, 29, 16, 24, 30, 22, 23}, {17, 19, 33, 31, 18, 26, 32, 24, 25}, {19, 21, 35, 33, 20, 28, 34, 26, 27}, {29, 31, 45, 43, 30, 38, 44, 36, 37}, {31, 33, 47, 45, 32, 40, 46, 38, 39}, {33, 35, 49, 47, 34, 42, 48, 40, 41}}

This is what i need:

enter image description here

the meshtopology i want:

{{1, 3, 14, 12, 2, 9, 13, 8}, {3, 5, 16, 14, 4, 10, 15, 9}, {5, 7, 18, 16, 6, 11, 17, 10}, {12, 14, 25, 23, 13, 20, 24, 19}, {14, 16, 27, 25, 15, 21, 26, 20}, {16, 18, 29, 27, 17, 22, 28, 21}, {23, 25, 36, 34, 24, 31, 35, 30}, {25, 27, 38, 36, 26, 32, 37, 31}, {27, 29, 40, 38, 28, 33, 39, 32}}

$\endgroup$
  • 1
    $\begingroup$ Please explain the meaning of inputs in your code, and specify the desired outcome precisely. Add what you have tried so far, and how it has not worked. As it is, your question is a job posting, not a request for help. $\endgroup$ – MarcoB May 11 '18 at 18:59
  • $\begingroup$ @MarcoB I have modified the question trying to explain it better. $\endgroup$ – Diogo May 11 '18 at 19:38
  • $\begingroup$ Thank you for adding further details. Still, however: in your bottom specification, and referring to the original numbering scheme, why did you skip original node 16, but not nodes 18 and 20? Should node 16 not have been included? $\endgroup$ – MarcoB May 11 '18 at 19:52
  • $\begingroup$ @MarcoB I have corrected the figure. Thank you for pointing the mistake. $\endgroup$ – Diogo May 11 '18 at 20:00
6
$\begingroup$

See if this works as a start:

ClearAll[generator, meshPlot]

generator[n_Integer] := Module[
  {allpts, noBodycentered, sortedAndLabeled},
  allpts = Flatten[#, 1] &@Table[{i, j}, {i, 0, n - 1}, {j, 0, n - 1}];
  noBodycentered = DeleteCases[allpts, {i_, j_} /; OddQ[i j]];
  sortedAndLabeled = MapIndexed[#1 -> First[#2] &, SortBy[Last]@noBodycentered]
]

meshPlot[list_List] := ListPlot[Labeled[#1, #2]& @@@list, Axes -> False, AspectRatio -> 1]

You can use them as follows:

generator[7]

(* Out: 
{{0, 0} -> 1, {1, 0} -> 2, {2, 0} -> 3, {3, 0} -> 4, {4, 0} -> 5, {5, 0} -> 6, 
 {6, 0} -> 7, {0, 1} -> 8, {2, 1} -> 9, {4, 1} -> 10, {6, 1} -> 11, {0, 2} -> 12, 
 {1, 2} -> 13, {2, 2} -> 14, {3, 2} -> 15, {4, 2} -> 16, {5, 2} -> 17, {6, 2} -> 18, 
 {0, 3} -> 19, {2, 3} -> 20, {4, 3} -> 21, {6, 3} -> 22, {0, 4} -> 23, {1, 4} -> 24,
 {2, 4} -> 25, {3, 4} -> 26, {4, 4} -> 27, {5, 4} -> 28, {6, 4} -> 29, {0, 5} -> 30, 
 {2, 5} -> 31, {4, 5} -> 32, {6, 5} -> 33, {0, 6} -> 34, {1, 6} -> 35, {2, 6} -> 36,
 {3, 6} -> 37, {4, 6} -> 38, {5, 6} -> 39, {6, 6} -> 40}
*)

meshPlot[generator[7]]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thank you! This will help. $\endgroup$ – Diogo May 11 '18 at 20:31
  • $\begingroup$ @Diogo You're welcome. $\endgroup$ – MarcoB May 11 '18 at 21:19
  • 1
    $\begingroup$ I like how you detect the vertices for deletion. $\endgroup$ – Henrik Schumacher May 12 '18 at 7:35
3
$\begingroup$

You could use ToElementMesh for that:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Rectangle[], "MaxCellMeasure" -> 0.1]
Show[
 mesh["Wireframe"],
 Graphics[
  MapIndexed[Text[#2[[1]], #1, {-1, -1}] &, mesh["Coordinates"]]]
 ]

enter image description here

The element incidents can then easily be extracted with:

ElementIncidents[mesh["MeshElements"]]
{{{1, 6, 7, 2, 26, 27, 28, 29}, {2, 7, 8, 3, 28, 30, 31, 32}, {3, 8, 
   9, 4, 31, 33, 34, 35}, {4, 9, 10, 5, 34, 36, 37, 38}, {6, 11, 12, 
   7, 39, 40, 41, 27}, {7, 12, 13, 8, 41, 42, 43, 30}, {8, 13, 14, 9, 
   43, 44, 45, 33}, {9, 14, 15, 10, 45, 46, 47, 36}, {11, 16, 17, 12, 
   48, 49, 50, 40}, {12, 17, 18, 13, 50, 51, 52, 42}, {13, 18, 19, 14,
    52, 53, 54, 44}, {14, 19, 20, 15, 54, 55, 56, 46}, {16, 21, 22, 
   17, 57, 58, 59, 49}, {17, 22, 23, 18, 59, 60, 61, 51}, {18, 23, 24,
    19, 61, 62, 63, 53}, {19, 24, 25, 20, 63, 64, 65, 55}}}

Besides that you have many efficient ElementMesh utilities documented in the Scope section here.

$\endgroup$
  • 1
    $\begingroup$ Ah, "QuadElements" is the default for rectangular domains and serendipity elements are used for quad meshes. I'd have expected bilinear elements for some reason. Both is very good to know! $\endgroup$ – Henrik Schumacher May 15 '18 at 9:39
  • $\begingroup$ @HenrikSchumacher, Yes, exactly right! $\endgroup$ – user21 May 15 '18 at 9:54
2
$\begingroup$

GridGraph + VertexDelete + VertexReplace:

f[n_?OddQ, o : OptionsPattern[]] := Module[{g2, 
  g1 = GridGraph[{n, n}, VertexLabels -> "Name", ImagePadding -> 5, o, ImageSize -> 400]},
  g1 = SetProperty[g1, VertexCoordinates -> (Reverse /@ GraphEmbedding[g1])];
  g2 = VertexDelete[g1, Select[Range[n  n], 
    OddQ[#] && EvenQ[Mod[#, n]] && Positive[Mod[#, n]] &]];
  g2 = SetProperty[VertexReplace[g2, Thread[VertexList[g2] -> Range[VertexCount[g2]]]], 
   {VertexCoordinates -> GraphEmbedding[g2], VertexLabels -> "Name", 
   ImagePadding -> 5, o, ImageSize -> 400}];
  {g1, g2}]

Examples:

Row @ f @ 7

enter image description here

Row @ f @ 5

enter image description here

Row @ f @ 9

enter image description here

$\endgroup$
  • $\begingroup$ nice graphs! But what is more important to me is the element topology. In my example the elements topology are:{{1, 3, 17, 15, 2, 10, 16, 8, 9}, {3, 5, 19, 17, 4, 12, 18, 10, 11}, {5, 7, 21, 19, 6, 14, 20, 12, 13}, {15, 17, 31, 29, 16, 24, 30, 22, 23}, {17, 19, 33, 31, 18, 26, 32, 24, 25}, {19, 21, 35, 33, 20, 28, 34, 26, 27}, {29, 31, 45, 43, 30, 38, 44, 36, 37}, {31, 33, 47, 45, 32, 40, 46, 38, 39}, {33, 35, 49, 47, 34, 42, 48, 40, 41}} $\endgroup$ – Diogo May 12 '18 at 3:09
  • $\begingroup$ The topology i need for my example witout the middle nodes is:{{1, 3, 14, 12, 2, 9, 13, 8}, {3, 5, 16, 14, 4, 10, 15, 9}, {5, 7, 18, 16, 6, 11, 17, 10}, {12, 14, 25, 23, 13, 20, 24, 19}, {14, 16, 27, 25, 15, 21, 26, 20}, {16, 18, 29, 27, 17, 22, 28, 21}, {23, 25, 36, 34, 24, 31, 35, 30}, {25, 27, 38, 36, 26, 32, 37, 31}, {27, 29, 40, 38, 28, 33, 39, 32}} $\endgroup$ – Diogo May 12 '18 at 3:15
  • $\begingroup$ @Diogo, is the order of nodes important? If not SortBy[VertexList/@FindCycle[f[7][[2]] ,{8}, All] , Min] give the same list (except the ordering of nodes in each sublist). $\endgroup$ – kglr May 12 '18 at 5:31
  • $\begingroup$ .. that is, it gives {{1, 2, 3, 9, 14, 13, 12, 8}, {3, 4, 5, 10, 16, 15, 14, 9}, {5, 6, 7, 11, 18, 17, 16, 10}, {23, 24, 25, 20, 14, 13, 12, 19}, {25, 20, 14, 15, 16, 21, 27, 26}, {27, 21, 16, 17, 18, 22, 29, 28}, {23, 24, 25, 31, 36, 35, 34, 30}, {25, 31, 36, 37, 38, 32, 27, 26}, {27, 32, 38, 39, 40, 33, 29, 28}} $\endgroup$ – kglr May 12 '18 at 6:16
  • $\begingroup$ the meshtopology order is very important. But this function(Find Cycle) might help. $\endgroup$ – Diogo May 12 '18 at 12:13
2
$\begingroup$

Let's start with a quad mesh. Here are some routines I frequently use to generate grids for reactangles, cylinders, and tori:

getGridPoints = Compile[
   {{x0, _Real},
    {x1, _Real},
    {y0, _Real},
    {y1, _Real},
    {m, _Integer},
    {n, _Integer},
    {xclosed, True | False},
    {yclosed, True | False}
    },
   Block[{mm, nn, x, y, δx, δy, oo},
    mm = m - Boole[xclosed];
    nn = n - Boole[yclosed];
    δx = (x1 - x0)/(m - 1);
    δy = (y1 - y0)/(n - 1);
    x = x0 - δx;
    y = y0 - δy;
    Flatten[
     Table[
      y += δy;
      x = x0 - δx;
      Table[
       x += δx;
       {x, y},
       {i, 1, mm}],
      {j, 1, nn}],
     1]
    ],
   CompilationTarget -> "C"
   ];

getGridQuads = 
  Compile[{{m, _Integer}, {n, _Integer}, {xclosed, 
     True | False}, {yclosed, True | False}},
   Block[{a1, a2, a3, a4, b1, b2, quads, qq, mm, nn},
    b1 = Boole[xclosed];
    b2 = Boole[yclosed];
    mm = m - b1;
    nn = n - b2;

    quads = Flatten[Table[
       qq = Table[
         a1 = mm (j - 1) + i;
         a2 = mm (j - 1) + i + 1;
         a3 = mm j + i;
         a4 = mm j + i + 1;
         {a1, a2, a4, a3},
         {i, 1, mm - 1}];

       If[xclosed,
        Join[qq,
         a1 = mm (j - 1) + mm;
         a2 = mm (j - 1) + 1;
         a3 = mm (j) + 1;
         a4 = mm (j) + mm;
         {{a1, a2, a4, a3}}
         ],
        qq
        ]
       ,
       {j, 1, nn - 1}], 1];

    If[yclosed,
     qq = Table[
       a1 = mm (nn - 1) + i;
       a2 = mm (nn - 1) + i + 1;
       a3 = i;
       a4 = i + 1;
       {a1, a2, a4, a3},
       {i, 1, mm - 1}];
     If[xclosed,
      a1 = mm nn;
      a2 = mm (nn - 1) 1;
      a3 = nn;
      a4 = 1;
      qq = Join[qq, {{a1, a2, a4, a3}}]
      ];
     Join[quads, qq],
     quads
     ]
    ],
   CompilationTarget -> "C"
   ];

getEdgesFromQuads = Compile[{{q, _Integer, 1}},
   Block[{q1, q2, q3, q4},
    q1 = Compile`GetElement[q, 1];
    q2 = Compile`GetElement[q, 2];
    q3 = Compile`GetElement[q, 3];
    q4 = Compile`GetElement[q, 4];
    {
     {Min[q1, q2], Max[q1, q2]},
     {Min[q2, q3], Max[q2, q3]},
     {Min[q3, q4], Max[q3, q4]},
     {Min[q4, q1], Max[q4, q1]}
     }
    ],
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Here is our first quad mesh.

m = 8;
n = 4;
R = MeshRegion[
  getGridPoints[0, 2, 0, 1, m + 1, n + 1, False, False],
  Polygon[getGridQuads[m + 1, n + 1, False, False]]
  ]

enter image description here

The following should convert an arbitrary quad mesh into your desired format. Feel free to package that into a function.

pts = MeshCoordinates[R];
vertexcount = Length[pts];
quads = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
data = Flatten[getEdgesFromQuads[quads], 1];
edges = DeleteDuplicates[data];
edgemidpts = 0.5 Total[Partition[pts[[Flatten[edges]]], 2], {2}];

allpts = Join[pts, edgemidpts];
ordering = Ordering[Transpose[Transpose[allpts][[{2, 1}]]]];

perm = ConstantArray[0, Length[allpts]];
perm[[ordering]] = Range[Length[allpts]];
edgemidptindices = perm[[vertexcount + 1 ;;]];

edgelookup = SparseArray[
   Rule[
    Join[edges, Transpose[Transpose[edges][[{2, 1}]]]],
    Join[edgemidptindices, edgemidptindices]
    ],
   {vertexcount, vertexcount}
   ];
newtopology = Join[
   Partition[perm[[Flatten[quads]]], 4],
   Partition[Extract[edgelookup, data], 4],
   2
   ];
newpts = allpts[[ordering]];

Now, newpts contains the positions of all quad corners and all edge midpoints. newtopology contains the index lists. The first four entries point to the corners of the quad (in same order as in the original mesh). Index 5 to 8 point to the edge midpoints in order consistent with the first four indices.

A test graphic:

Graphics[{
  GraphicsComplex[
   newpts, {
    EdgeForm[Black], FaceForm[ColorData[97][1]], 
    Polygon[newtopology[[All, 1 ;; 4]]],
    EdgeForm[Black], FaceForm[ColorData[97][2]], 
    Polygon[newtopology[[All, 5 ;; 8]]]
    }
   ],
  Table[Text[i, newpts[[i]], Background -> LightBlue], {i, 1, 
    Length[newpts]}]
  }]

enter image description here

$\endgroup$
  • $\begingroup$ thank you. This is what i want. $\endgroup$ – Diogo May 13 '18 at 11:35
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher May 13 '18 at 11:45
1
$\begingroup$

After a lot of pain, I found a solution. I hope there are no bugs.

GenerateGridMesh[aa_, bb_, nx_, ny_] := 
  Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j, 
    meshtopology = {}, allcoords, k, topolsz, l, edge, vec, data, a, 
    b, c},
   vec = {};
   edge = {};
   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[vec, {x, y}];
       x += dx ;
       ];
     ,
     For[k = 1, k <= nx + 1, k++,
       AppendTo[vec, {x, y}];
       x += 2 dx ;
       ];
     ];
    x = 0;
    y += dy;
    ];
   meshtopology = {};
   a = 0;
   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[vec[[meshtopology[[i, j]]]], {i, 1, 
      Length[meshtopology]}, {j, 1, Length[meshtopology[[1]]]}];
   {allcoords, vec, meshtopology}
   ];
(*This function generates a mesh graphic*)
GenerateGraphics[nodes_, topology_, order_] := 
  Block[{meshvis, nodevis, v},
   v = {1, 5, 2, 6, 3, 7, 4, 8};
   meshvis = 
    Graphics[{FaceForm[], EdgeForm[Blue], 
      GraphicsComplex[nodes, Polygon[topology[[All, v]]]]}];
   nodevis = 
    Graphics[{PointSize[Large], 
      MapIndexed[Text[#2[[1]], #1, {-1.5, 1.5}] &, nodes], {Blue, 
       Point[nodes]}}];
   {meshvis, nodevis}];
{allcoords, meshnodes, meshtopology} = GenerateGridMesh[2, 2, 3, 2];
{meshvis, nodevis} = GenerateGraphics[meshnodes, meshtopology, 2];
Show[meshvis, nodevis]

meshtopology

{{1, 3, 14, 12, 2, 9, 13, 8}, {3, 5, 16, 14, 4, 10, 15, 9}, {5, 7, 18, 16, 6, 11, 17, 10}, {12, 14, 25, 23, 13, 20, 24, 19}, {14, 16, 27, 25, 15, 21, 26, 20}, {16, 18, 29, 27, 17, 22, 28, 21}}

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.