# How to generate a mesh with quadrilateral elements?

I have the following code that generates a finite element mesh:

Needs["NDSolveFEM"]
order = 2;
mesh = ToElementMesh[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]], "MeshElementType" -> QuadElement,
"MeshOrder" -> order, "NodeReordering" -> True]
mesh["Wireframe"]


I'm indicating that the element type i want is a quadrilateral type ("MeshElementType" -> QuadElement), but it still creates a mesh of triangles. Can i change this? • it would be easy to hand code for such simple geometry. Feb 10, 2017 at 14:16
• @george2079 i need a very fine mesh. For simplicity i have chosen to shown a coarse mesh in this question. It is very dificult to create more than 1000 nodes by hand. Feb 10, 2017 at 14:22
• i did not literally mean by hand, you would use a loop of course. The density doesn't matter, but the exact shape of the boundary does. Feb 10, 2017 at 14:49

Based on the tutorial here

QuadElement meshes behave exactly the same as TriangleElement meshes, with the exception that, for linear quad elements, four incidents per element are needed, and, for quadratic elements, eight incidents per element are needed.

Needs["NDSolveFEM"]

rectangleA = {{0, 4}, {8, 8}};
rectangleB = {{5, 0}, {8, 4}};

pitch = 0.25;

rectAx = 1/pitch (rectangleA[[2, 1]] - rectangleA[[1, 1]]) + 1;
rectAy = 1/pitch (rectangleA[[2, 2]] - rectangleA[[1, 2]]) + 1;
rectBx = 1/pitch (rectangleB[[2, 1]] - rectangleB[[1, 1]]) + 1;
rectBy = 1/pitch (rectangleB[[2, 2]] - rectangleB[[1, 2]]) + 1;

coordinatesA =
Flatten[Table[{x, y}, {x, rectangleA[[1, 1]], rectangleA[[2, 1]], pitch},
{y, rectangleA[[1, 2]], rectangleA[[2, 2]], pitch}], 1];
coordinatesB =
Flatten[Table[{x, y}, {x, rectangleB[[1, 1]], rectangleB[[2, 1]], pitch},
{y, rectangleB[[1, 2]], rectangleB[[2, 2]], pitch}], 1];

coordinates = Join[coordinatesA, coordinatesB];

incidents = Join[
Flatten[
Table[{j*rectAy + i, j*rectAy + i + 1, (j - 1)*rectAy + i + 1, (j - 1)*rectAy + i},
{i, 1, rectAy - 1},
{j, 1, rectAx - 1}], 1],
Flatten[
Table[{j*rectBy + i, j*rectBy + i + 1, (j - 1)*rectBy + i + 1, (j - 1)*rectBy + i},
{i, 1 + Length[coordinatesA], Length[coordinatesA] + rectBy - 1},
{j, 1, rectBx - 1}], 1]
];

mesh1 = ToElementMesh["Coordinates" -> coordinates,
"NodeReordering" -> True];
mesh2 = MeshOrderAlteration[mesh1, 2];
mesh2["MeshOrder"] Show[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]],
mesh2["Wireframe"],
ListPlot[coordinates]
]


Sample Problem:

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100, y == 4 && x <= 5 || y <= 4 && x == 5],
DirichletCondition[u[x, y] == 0, y == 8 || x == 8]},
u, {x, y} ∈ mesh2]

sol["ElementMesh"]
ContourPlot[sol[x, y], {x, y} ∈ mesh2, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic] Here is a different approach. The idea is to first generate a second order triangle mesh. Next a center coordinate is added in every triangle. Then we split every triangle into three first order quads making use of the newly added center coordinate and the initial second order triangle incidents. So we add a node in the center of the element and create the quad elements {c,6,1,4}, {c,4,2,5} and {c,5,3,6}.

Needs["NDSolveFEM"]
m = ToElementMesh[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]]];


Two helper functions that work for second order triangle mesh elements. If one want to extend this to 3D then you'd need to add a tetToHex function.

getMidEleCoords = Compile[{{eleCoords, _Real, 2}},
Mean[eleCoords], RuntimeAttributes -> {Listable}];

triToQuad = Compile[{{triInci, _Integer, 1}, {midInci, _Integer, 0}}, {
{midInci, triInci[], triInci[], triInci[]},
{midInci, triInci[], triInci[], triInci[]},
{midInci, triInci[], triInci[], triInci[]}
}
, RuntimeAttributes -> {Listable}
];

getNumOfResEle[TriangleElement] := 3


First, we compute the new coords within the triangle/tet elements:

coords = m["Coordinates"];
numOfCoords = Length[coords];
meshEle = m["MeshElements"];
inci = ElementIncidents[meshEle];
markers = ElementMarkers[meshEle];
eleCoords = GetElementCoordinates[coords, #] & /@ inci;
midEleCoords = getMidEleCoords[eleCoords];
numOfMidCoords = Length /@ midEleCoords;
temp = Join[{numOfCoords}, numOfMidCoords];
temp = FoldList[Plus, temp] + 1;
midEleInci = MapThread[Range, {Most[temp], Rest[temp - 1]}];
newCoords = Join[coords, Join @@ midEleCoords];
(*Max[midEleInci]===Length[newCoords]*)


The we generate the new elements:

newEle = Table[
Block[{thisMeshEle, thisType, thisInci, thisMarker, converter,
newEle, newEleMarkers},
thisMeshEle = meshEle[[i]];
thisInci = ElementIncidents[thisMeshEle];
thisMarker = ElementMarkers[thisMeshEle];
converter = getConverterCode[thisType];
newEle = Join @@ converter[thisInci, midEleInci[[i]]];
newEleMarkers =
Join @@ Transpose[
ConstantArray[thisMarker, getNumOfResEle[thisType]]];
getConvertedEleType[thisType][newEle, newEleMarkers]
]
, {i, Length[meshEle]}];


mq1 = ToElementMesh["Coordinates" -> newCoords,
"MeshElements" -> newEle]
(* ElementMesh[{{0., 8.}, {0., 8.}}, {QuadElement["<" 1338 ">"]}] *)


Look at the mesh:

mq1["Wireframe"] This is now a first order mesh:

mq1["MeshOrder"]
1


One can use MeshOrderAlteration[mq1, 2] to then subsequently generate a second order mesh.

Histogram[Join @@ mq1["Quality"]] Quality looks OK, but a smoothing of some kind (e.g. Laplacian smoothing) might be well worth trying out.

(Update: There is a follow up question that discusses the idea of Laplacian smoothing. Spoiler: Does not look like it's worth the trouble)

Solve a model problem over the mesh:

sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1,
DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] mq1];
Plot3D[sol[x, y], {x, y} \[Element] sol["ElementMesh"]] Since this will generate three quad elements for every triangle the mesh cell measure is lower than what it needs to be, one could adjust for that by making use of "MaxCellMeasure" in the initial triangle mesh.

Hope this helps.

That's not possible in general in V11.0. The only primitive from which a quad mesh can be generated is

ToElementMesh[Rectangle[]]


or

ToElementMesh[FullRegion, {{0, 1}, {0, 1}}]


So to get a quad mesh for a more complicated shape you'd need to program a bit and use

ToElementMesh["Coordinates"->..., "MeshElements"->{QuadElement[....]}]


As another approach you could merge 2 triangles to a quad and thus get a mostly quad mesh. If you do that. I'd be interested in seeing the result.

Here is another piece of information: If you have a first order mesh you can make that second order.

m = ToElementMesh["Coordinates" -> coordinates,
m["MeshOrder"]
1

m2 = MeshOrderAlteration[m, 2];
m2["MeshOrder"]
2


Since your case does not involve curved boundaries no mid side nodes (second order nodes) need to be moved to the proper position.

• don't even try merging triangles, you will end up with horribly unacceptable shape elements. Feb 10, 2017 at 14:57
• @george2079, probably not easy, but it has been done Feb 10, 2017 at 15:22
• Are quad elements constrained to single order? Feb 10, 2017 at 17:30
• @Young, no. They can be quadratic. Feb 10, 2017 at 17:46
• @Young, see update on getting a second order mesh. Feb 10, 2017 at 18:00

this is a bit of a hack but it works for this example:

generate separate rectangle meshes -- note the connected edge nodes must be exactly coincident.

<< NDSolveFEM
m1 = ToElementMesh[Rectangle[{-10, 0}, {0, 10}],
"MeshElementType" -> QuadElement, MaxCellMeasure -> 4,
"MeshOrder" -> 1];
m2 = ToElementMesh[Rectangle[{-20, 10}, {0, 30}],
"MeshElementType" -> QuadElement, MaxCellMeasure -> 4,
"MeshOrder" -> 1];
Show[{m1["Wireframe"] /. EdgeForm[GrayLevel] -> EdgeForm[Red],
m2["Wireframe"]}]


all of this works fine for "MeshOrder"->2 btw join the node lists and shift the second mesh to the new list index:

allnodes = Join[m1["Coordinates"], m2["Coordinates"]];
common = Intersection[m1["Coordinates"], m2["Coordinates"]];
m2el = Map[Length@m1["Coordinates"] + # &, m2["MeshElements"] , {2}];


now point the duplicate nodes in the second mesh to the nodes in the first mesh:

m2el = m2el /. ((Rule @@ Reverse@Flatten[Position[allnodes, #]]) & /@
common);


Edit: This step is not needed for NDSolve NDsolve evidently merges coincident nodes (beware other solvers might not do that).

create a new mesh:

allelements = Join[m1["MeshElements"], m2el];
ToElementMesh["Coordinates" -> allnodes,
"MeshElements" -> allelements];
Show[%["Wireframe"]] Note we left stray unconnected nodes in the allnodes list, I don't think that should be a problem ..

Edit: test problem: ( with MaxCellMeasure -> 1/4 )

eqns = {
Laplacian[u[x, y], {x, y}] == 0,
DirichletCondition[u[x, y] == 1, y == 30],
DirichletCondition[u[x, y] == 0, y == 0]};
uif = NDSolveValue[eqns, u, {x, y} \[Element] mesh]
Point[#1]} &,
{uif["Grid"], uif["ValuesOnGrid"]}]] a better way to plot..

vertvals = ColorData["SunsetColors"]@Extract[uif["ValuesOnGrid"],
First@Position[   uif["Grid"] , #] ] & /@ allnodes;
(* note the "Grid" list and allnodes are not the same because of the unused nodes in allnodes *)
Graphics[GraphicsComplex[ allnodes,
Cases[ allelements , {n1_, n2_, n3_, n4_,___} :>
Polygon[{n1, n2, n3, n4}], Infinity], VertexColors -> vertvals ]] Edit:

This method is now implemented as StructuredMesh function (2D and 3D) in FEMAddOns package. With minor enhancements it is also implemented in MeshTools package.

If the shape of geometric region can be somehow described as a rectangular array of points then structured quadrilateral mesh can be generated on it. I hope the next example will more clearly illustrate my point.

First we create a raster of 2D points and then pass it to structuredMesh function which uses 2D interpolation of each coordinate. It returns ElementMesh object according to specified number of elements in each direction.

 raster = With[{end = 4 Pi, n = 20},
N@{
Table[{x, Sin[x] - 2}, {x, 0, end, end/n}],
Table[{x, Sin[x]/4}, {x, 0, end, end/n}],
Table[{x, Cos[x] + 2}, {x, 0, end, end/n}]
}
];

mesh = structuredMesh[raster, {50, 10}, InterpolationOrder -> 2];
Show[
mesh["Wireframe"],
Graphics[{Red, PointSize[0.01], Point[Flatten[raster, 1]]}]
] Definition of structuredMesh and helper function(s) is given below.

getElementConnectivity[nx_,ny_]:=Flatten[
Table[
{(i-1)*(nx+1)+j,i*(nx+1)+j,i*(nx+1)+j+1,(i-1)*(nx+1)+j+1},
{i,1,ny},
{j,1,nx}
],
1
]

structuredMesh::array="Raster of input points must be full array of numbers with depth of 3.";

structuredMesh//Options={InterpolationOrder->1};
structuredMesh[raster_,{nx_,ny_},opts:OptionsPattern[]]:=Module[
{order,dim,restructured,xInt,yInt,zInt,nodes,connectivity},
If[Not@ArrayQ[raster,3,NumericQ],Message[structuredMesh::array];Return[\$Failed]];

order=OptionValue[InterpolationOrder]/.Automatic->1;
dim=Last@Dimensions[raster];

restructured=Transpose[raster,{2,3,1}];
xInt=ListInterpolation[restructured[],{{0,1},{0,1}},InterpolationOrder->order];
yInt=ListInterpolation[restructured[],{{0,1},{0,1}},InterpolationOrder->order];

nodes=Flatten[#,1]&@If[dim==3,
zInt=ListInterpolation[restructured[],{{0,1},{0,1}},InterpolationOrder->order];
Table[{xInt[i,j],yInt[i,j],zInt[i,j]},{i,0,1,1./ny},{j,0,1,1./nx}]
,
Table[{xInt[i,j],yInt[i,j]},{i,0,1,1./ny},{j,0,1,1./nx}]
];

connectivity=Reverse/@getElementConnectivity[nx,ny];

If[dim==3,
]
]


And this is how the function can be used for OP's shape. In general the more complicated shapes can be manually split to "topological quadrialterals", meshed separately with structured meshes and then merged in one mesh.

raster2 = {{{5, 0}, {8, 0}}, {{5, 4}, {8, 8}}, {{0, 4}, {0, 8}}};
mesh = structuredMesh[raster2, 10*{1, 2}];
mesh["Wireframe"] • This is a very cool way to do it. Feb 26, 2018 at 10:04

Edit:

This method is now also implemented as ToQuadMesh function in FEMAddOns package and as TriangleToQuadMesh in MeshTools package.

Once the open source package is downloaded and installed generating a quad element mesh then works like this:

Needs["FEMAddOns"]
triMesh =
ToElementMesh[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]]];

mesh["Wireframe"] Note that this also returns a second order mesh:

mesh["MeshOrder"]
2


AceFEM package contains a useful function SMTTriangularToQuad which converts triangular mesh to quadrilaterals. It works only for 1st order mesh.

<< AceFEM
triMesh = ToElementMesh[
DiscretizeGraphics[
GraphicsComplex[{{0,4},{5,4},{5,0},{8,0},{8,8},{0,8}},Polygon[{1,2,3,4,5,6}]]],
"MeshOrder" -> 1
]

mesh["Wireframe"] • ToQuadMesh is probably the most general method of the solutions presented for this question. Jun 18, 2018 at 15:06

For integer valued coordinates like you have here, you can make a grid of rational step size and use RegionMember + ArrayMesh.

Something like:

mr = MeshRegion[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}}, Polygon[{1, 2, 3, 4, 5, 6}]];

{{x1, x2}, {y1, y2}} = bds = RegionBounds[mr];

δ = 1/10;
xrng = Range[x1, x2 - δ, δ];

onoff = Boole[Table[ MeshCellCount[quad]

{4561, 8960, 4400}