My question is inspired by this answer about meshing with quadrilateral elements.

The following functions are my (naive) try of implementing the Laplacian smoothing. I find interior nodes in the mesh and move them to the average coordinates of all other connected nodes.

This is the simplest example of the procedure.

<< NDSolve`FEM`
mesh = ToElementMesh[
   "Coordinates" -> {{0., 0.}, {0.5, 0.}, {1., 0.}, {0., 0.5}, {0.4,0.4}, {1., 0.5}, {0., 1.}, {0.5, 1.}, {1., 1.}},
   "MeshElements" -> {QuadElement[{{1, 2, 5, 4}, {2, 3, 6, 5}, {4, 5, 8, 7}, {5, 6, 9, 8}}]}
smoothed = smoothMesh[mesh];

  Show[mesh["Wireframe"], ImageSize -> 100],
  Show[smoothed["Wireframe"], ImageSize -> 100]



findInteriorNodes[mesh_ElementMesh] := With[{
   allNodes = Range@Length[mesh["Coordinates"]],
   boundaryNodes = Flatten@ElementIncidents[mesh["BoundaryElements"]]
  Complement[allNodes, boundaryNodes]

(* Find nodes of the "mesh" connected to node "n". *)
findConnectedNodes[mesh_ElementMesh, n_Integer] := With[{
   edges = MeshCells[MeshRegion[mesh], 1][[All, 1]]
   Flatten@Select[edges, MemberQ[#, n] &],

(* Returns the rule to be later used in ReplacePart *)
changeNodePosition[mesh_ElementMesh, n_Integer] := With[
  {connectedNodes = findConnectedNodes[mesh, n]},
  n -> Mean[mesh["Coordinates"][[connectedNodes]]]

(* The main function *)
Options[smoothMesh] = {"SmoothingLoops" -> 1};
smoothMesh[mesh_ElementMesh, opts : OptionsPattern[]] := Module[
  {originalMesh, newMesh, nLoops, interiorNodes, modCrds},
  nLoops = OptionValue["SmoothingLoops"];
  interiorNodes = findInteriorNodes[mesh];
  originalMesh = mesh;

   modCrds = changeNodePosition[originalMesh, #] & /@ interiorNodes;
   newMesh = ToElementMesh[
     "Coordinates" -> ReplacePart[originalMesh["Coordinates"], modCrds],
     "MeshElements" -> mesh["MeshElements"]
   originalMesh = newMesh,

(* a helper function to compare mesh quality *)
compareMeshQuality[meshes : {__}] := Module[
  {qualityList, measure, statistics, names},
  qualityList = Flatten[Through[meshes["Quality"]], 1];
  measure = {Min, Mean, Median};
  statistics = Map[
    Map[#]@qualityList &,
  names = "Mesh " <> ToString[#] & /@ Range[Length[meshes]];
    TableForm[statistics, TableHeadings -> {measure, names}],
  Histogram[qualityList, {0.05}, ImageSize -> 200, ChartLegends -> names]

Results of the example given in the linked question for L shaped domain where quadrilateral elements are created by subdivision of triangular elements.

smoothed = smoothMesh[mq1, "SmoothingLoops" -> 3];
  Show[mq1["Wireframe"], ImageSize -> 300],
  Show[smoothed["Wireframe"], ImageSize -> 300]


compareMeshQuality[{mq1, smoothed}]


We can see from the last figure that the quality of the smoothed mesh is (slightly) improved.

My questions are:

1) How can we make this code faster, because it is pretty slow for larger meshes and multiple smoothing loops?

2) Is there any better algorithm to improve the smoothing procedure?


3 Answers 3


Here is an answer to question 1. I'll comment on question 2 at the end.

The thing to realize is that the connectivity never changes during the smoothing algorithm. So to speed things up we compute the connectivity once and then use that information repetitively.

We use the "VertexElementConnectivity" from the ElementMesh (See documentation for specifics) and take a dot product of that. We then extract the position where we have two connections of the nonzero values. We construct a SparseArray that we use as a lookup table for a node to connected nodes lookup. This is packaged into a function and returned.

mkFindConnectedNodes[mesh_ElementMesh] := 
 Block[{vec, vecTvec, nzv, nzp, neigbours, temp2, sa},
  vec = mesh["VertexElementConnectivity"];
  vecTvec = vec.Transpose[vec];
  nzv = vecTvec["NonzeroValues"];
  nzp = vecTvec["NonzeroPositions"];
  neigbours = Join @@ Position[nzv, 2];
  temp2 = Transpose[nzp[[neigbours]]];
  sa = SparseArray[
        Range[Length[temp2[[1]]]]}]] -> 1];
  With[{lookup = sa, ids = temp2[[2]]},
    ids[[Flatten[lookup[[#]]["NonzeroPositions"]]]] & /@ theseNodes]

A slightly more efficient version of your findInteriorNodes:

findInteriorNodes[mesh_ElementMesh] := 
 With[{allNodes = Range[Length[mesh["Coordinates"]]], 
   boundaryNodes = 
    Join @@ (Flatten /@ ElementIncidents[mesh["BoundaryElements"]])}, 
  Complement[allNodes, boundaryNodes]]

And the smoothing function:

smoothMesh[m_ElementMesh, iter_] := 
 Block[{inodes, findConnectedNodes2, temp, movedCoords, oldCoords},
  inodes = findInteriorNodes[m];
  findConnectedNodes2 = mkFindConnectedNodes[m];
  temp = findConnectedNodes2[inodes];
  oldCoords = m["Coordinates"];
   movedCoords = 
    Developer`ToPackedArray[Mean[oldCoords[[#]]] & /@ temp];
   oldCoords[[inodes]] = movedCoords;
   , {iter}];
  ToElementMesh["Coordinates" -> oldCoords, 
   "MeshElements" -> m["MeshElements"], 
   "BoundaryElements" -> m["BoundaryElements"], 
   "PointElements" -> m["PointElements"]]

Note, the updates of the coordinates is now the only thing in the loop.

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

For 1000 iterations:

AbsoluteTiming[smoothedMesh = smoothMesh[m, 1000];]
{1.60202`, Null}

GraphicsGrid[{{smoothedMesh["Wireframe"], m["Wireframe"]}}, 
 ImageSize -> Large]

enter image description here

compareMeshQuality[{smoothedMesh, m}]

enter image description here enter image description here

A second example:

AbsoluteTiming[smoothedMesh = smoothMesh[mq1, 1000];]
{2.454281`, Null}

GraphicsGrid[{{smoothedMesh["Wireframe"], mq1["Wireframe"]}}, 
 ImageSize -> Large]

enter image description here

compareMeshQuality[{smoothedMesh, mq1}]

enter image description here enter image description here

As you can see the results are OK-ish but they do not rock you out of your socks.

A different approach that could be tried is to consider the connectivity as a set of springs and have these smooth out the mesh. But this is nonlinear and needs connectivity changes, so it's going to be slow.


Here is a new version. This is based on Henrik's nice answer which is very efficient in time. I put that in a function and extended a bit. This now works for any element type, any mesh order and any dimension. Internal boundary nodes (for example specified via "IncludePoints") are preserved.

Here is the function:

LaplacianElementMeshSmoothing[mesh_] := 
 Block[{n, vec, mat, adjacencymatrix2, mass2, laplacian2, 
   bndvertices2, interiorvertices, stiffness, load, newCoords},
  n = Length[mesh["Coordinates"]];
  vec = mesh["VertexElementConnectivity"];
  mat = Unitize[vec.Transpose[vec]];
  vec = Null;
  adjacencymatrix2 = mat - DiagonalMatrix[Diagonal[mat]];

  mass2 = DiagonalMatrix[SparseArray[Total[adjacencymatrix2, {2}]]];
  stiffness = N[mass2 - adjacencymatrix2];
  adjacencymatrix2 = Null;
  mass2 = Null;
  bndvertices2 = 
   Flatten[Join @@ ElementIncidents[mesh["PointElements"]]];
  interiorvertices = Complement[Range[1, n], bndvertices2];
  stiffness[[bndvertices2]] = 
   IdentityMatrix[n, SparseArray][[bndvertices2]];
  load = ConstantArray[0., {n, mesh["EmbeddingDimension"]}]; 
  load[[bndvertices2]] = mesh["Coordinates"][[bndvertices2]]; 
  newCoords = 
   LinearSolve[stiffness, load(*,Method\[Rule]"Pardiso"*)];

  ToElementMesh["Coordinates" -> newCoords, 
   "MeshElements" -> mesh["MeshElements"], 
   "BoundaryElements" -> mesh["BoundaryElements"], 
   "PointElements" -> mesh["PointElements"], 
   "CheckIncidentsCompletness" -> False, 
   "CheckIntersections" -> False, 
   "DeleteDuplicateCoordinates" -> False]

A few comments: I use the vertex element connectivity to construct the adjacency matrix. This will make it work for any element type. To compute the boundary points I use the PointElements - that will take every thing that is on a boundary as specified by the mesh, regardless if it's an internal boundary point or a point on the surface. If you want more speed and less memory usage you could try to use "Method"->"Pardiso". Then I use ToElementMesh to construct the new mesh. Setting the coordinates with Part is not a good idea for various reasons. I am not a 100% sure that ToElementMesh should not be checking for intersections, duplicates etc. We will get to that later.

Here are a few examples:

coords = Sort[RandomReal[{0, 1}, {10, 1}]];
m1D = ToElementMesh["Coordinates" -> coords, 
   "MeshElements" -> {LineElement[Partition[Range[10], 2, 1]]}];

sm = LaplacianElementMeshSmoothing[m1D]

Subtract @@@ Partition[Flatten[m1D["Coordinates"]], 2, 1]
{-0.0969347, -0.0755736, -0.0836963, -0.0438433, -0.0958689, \
-0.00254611, -0.0298039, -0.180217, -0.145562}

Subtract @@@ Partition[Flatten[sm["Coordinates"]], 2, 1]

{-0.0837829, -0.0837829, -0.0837829, -0.0837829, -0.0837829, \
-0.0837829, -0.0837829, -0.0837829, -0.0837829}

A quad element mesh example:

smoothMesh2 = 
 LaplacianElementMeshSmoothing[MeshOrderAlteration[mq1, 2]]

(* ElementMesh[{{-1.18329*10^-29, 8.}, {-9.86076*10^-31, 
   8.}}, {QuadElement["<" 1338 ">"]}] *)

(Note the 10^-29 - I am not sure if one needs to worry about that. Comments?)


enter image description here

Note that we used a second order mesh:


Here is an example with an interior node that is then not moved around:

mesh = ToElementMesh[Disk[], "IncludePoints" -> {{1/2, 1/2}}];
smoothMesh2 = LaplacianElementMeshSmoothing[mesh];
Position[smoothMesh2["Coordinates"], {1/2, 1/2} // N]

Now things become tricky.

mesh = ToElementMesh[
    x^2 + y^2 + z^2 <= 
     1, {{x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1}}], 
   MaxCellMeasure -> 1];
sm = LaplacianElementMeshSmoothing[mesh];
compareMeshQuality[{sm, mesh}]

enter image description here

enter image description here

We see that mean and median improve but the min deteriorated! This can lead to a negative quality measure:

mesh = ToElementMesh[
    x^2 + y^2 + z^2 <= 
     1, {{x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1}}], 
   "MeshOrder" -> 1];
sm = LaplacianElementMeshSmoothing[mesh];

ToElementMesh::femimq: The element mesh has insufficient quality of -0.148995. A quality estimate below 0. may be caused by a wrong ordering of element incidents or self-intersecting elements.

Even if we extract the coordinates and have ToElementMesh regenerate a delaunay tetrahedralization we get a very small min mesh element quality:

sm2 = ToElementMesh[sm["Coordinates"]];
compareMeshQuality[{sm2, mesh}]

enter image description here

enter image description here

This means that even if ToElementMesh could auto fix the element's node ordering we'd still get a close to zero min quality.

I can not quite say if this is expected; if someone has thoughts about this I'd be interested in hearing them.

  • $\begingroup$ Thank you for showing me an efficient approach to find connected nodes. Why do you reset some symbols in that function (e.g. vec = Null;)? Does it have to do with saving the memory? Doesn't Block automatically take care of that? $\endgroup$
    – Pinti
    Sep 28, 2017 at 9:09
  • 1
    $\begingroup$ @Pinti, Block will 'free' it's memory on exit. If one want to allow for that earlier I usually just overwrite that with Null. In this case I am not sure if it helps; it;s a remainder from when I developed that function. I think it's OK to remove it in this case. $\endgroup$
    – user21
    Sep 28, 2017 at 22:11
  • 1
    $\begingroup$ Nice! I am not too familiar with ElementMesh and I learnt quite a lot from your reworked code! Towards the bad behavior for tets: Even for triangle mesh in the plane, the method is not bulletproof, I guess. If the star of a vertex is nonconvex, the average of the neighbors may lie outside of the star, so the vertex might get pushed outwards. This might not happen too often in the plane, but it might be quite common for tet meshes... But I have to admit that I did not expect that at all, in particular with such a simple tet mesh. $\endgroup$ Sep 30, 2017 at 8:58
  • 2
    $\begingroup$ @HenrikSchumacher, one possibly could go ahead and hold a second layer of nodes next to the boundary fixed but then things get more complicated. I doubt that smoothing of the mesh will very much improve FEM results. My feeling is that mesh smoothing is not worth the trouble if the original mesh is somewhat decent. For the original problem it might be best to implement a quad/hex mesh generator in a different manner - though that's non trivial. $\endgroup$
    – user21
    Sep 30, 2017 at 12:01
  • 1
    $\begingroup$ @user21 If I really wanted to write a mesh smoother, I would try optimization methods based on elastic energies that blow up as the element qualities deteriorate, using Cauchy-Green-tensors from regularly shapes elements as reference.This would also prevent selfintersections. I would solve this nonlinear optimization problem iteratively. This would need quite some effort and it might prove not worth it, in the end. $\endgroup$ Oct 3, 2017 at 10:36

The conditions for a Tutte embedding can be formulated as a single (sparse) linear system that involves the graph Laplacian of the mesh. The graph Laplacian can be obtained from the adjacency matrix (see also Laplacian smoothing). In your case, we have to fix the boundary vertices such that this problem becomes the Dirichlet problem for the graph Laplacian.

Your approach resembles the analogous heat flow (for the vertex positions) which has the Tutte embedding as steady state (limit for the number of iterations going to $\infty$). So you might save a lot of time by skipping the heat flow and solve the steady state case directly by solving a single linear and sparse system. And even if you prefer your iterative approach (which gives you more control over the distortion), computing the graph Laplacian once and utilizing it for the iterations should speed up the process. With implicit time discretizations of the heat flow (such as Crank-Nicholson), you may also gain a lot since implicit methods propagate local changes much more efficiently through the whole mesh.

Alas, I am not sure wether the interior vertices will stay within the domain (if the domain is convex, I am quite confident about that).

Setting up the graph Laplacian

R = ToElementMesh[DiscretizeGraphics[
     N@{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}}, 
     Polygon[{1, 2, 3, 4, 5, 6}]]],
   "MeshOrder" -> 1
n = Length[R["Coordinates"]];
edgesbyvalence = Module[{f1, f2, f3},
   {f1, f2, f3} = Transpose[R["MeshElements"][[1, 1]]];
    Tally[Sort /@ ArrayReshape[Transpose[{f1, f2, f2, f3, f3, f1}],{Length[f1] 3, 2}]],
    Last -> First]
adjacencymatrix = With[{a = SparseArray[Join @@ Values[edgesbyvalence] -> 1., {n, n}, 0.]}, a + a\[Transpose]];
mass = DiagonalMatrix[SparseArray[Total[adjacencymatrix, {2}]]];
laplacian = N[mass - adjacencymatrix];
bndvertices = Union @@ edgesbyvalence[1];
interiorvertices = Complement[Range[1, n], bndvertices];

Solving for the steady state (Tutte embedding)

This is the somewhat ideomatic way to solve the Dirichlet problem for the graph Laplacian.

A = laplacian;
A[[bndvertices]] = N@IdentityMatrix[n, SparseArray][[bndvertices]];
b = ConstantArray[0., {n, 2}];
b[[bndvertices]] = R["Coordinates"][[bndvertices]];
x = LinearSolve[A, b];
S = R;
S[[1]] = x;

Even cleverer is the following which employs a smaller, symmetric, and positive definite linear system:

S = R; 
S[[1]][[interiorvertices]] =  LinearSolve[
  laplacian[[interiorvertices, interiorvertices]], 
  -laplacian[[interiorvertices, bndvertices]].R["Coordinates"][[bndvertices]],
  Method -> "Cholesky"

GraphicsRow[{S["Wireframe"], R["Wireframe"]}, ImageSize -> Large]

enter image description here

Testing the quality with the method provided by the OP yields

enter image description here

Doing that with the mesh with 69695 triangles that is generated like R but with MaxCellMeasure -> 0.001 takes about 0.6 s on my machine.

  • 1
    $\begingroup$ Yeah, I'd like to see an implementation of that. Sounds interesting. $\endgroup$
    – user21
    Sep 29, 2017 at 7:09
  • $\begingroup$ I'm also interested to see your solution. $\endgroup$
    – Pinti
    Sep 29, 2017 at 7:37
  • 1
    $\begingroup$ @user21 There it is for the steady state. If you are interested, I can also do the flow approach tonight. $\endgroup$ Sep 29, 2017 at 9:00
  • $\begingroup$ @HenrikSchumacher Impressive! Do you think you could make it work also for quadrilateral elements? $\endgroup$
    – Pinti
    Sep 29, 2017 at 12:42
  • 1
    $\begingroup$ I am not sure. But splitting quadrilaterals into triangles, applying the smoothing and reversing the splitting afterwards might work well. (I do not know if it is guaranteed that the quads stay convex in this process but it somewhat feels like that.) Or if you obtain your quad mesh from a triangle mesh via splitting (e.g. your mesh in the question), then perform the smoothing first on the triangle mesh and then split. $\endgroup$ Sep 29, 2017 at 13:02

There is a paclet called FEMAddOns on the WRI github acount that contains, among other things, a function ElementMeshSmoothing. This has both the iterative and the direct approach of this answer implemented and documented.

You can install the paclet via


This will install the current paclet version. You can uninstall it with:


Once you have that you can load the package:


This will auto load NDSolve`FEM too. The examples above then work in the following way:

m = ToElementMesh[
    GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}}, 
     Polygon[{1, 2, 3, 4, 5, 6}]]]];
AbsoluteTiming[smoothedMesh = ElementMeshSmoothing[m];]
{0.027129, Null}

This will use the default method which is based on Henrik's method.

compareMeshQuality[{smoothedMesh, m}]

enter image description here

Alternatively you can use an iterative method:

 smoothedMesh2 = ElementMeshSmoothing[m, Method -> "Iterative"];]
{0.746984, Null}

compareMeshQuality[{smoothedMesh2, m}]

enter image description here

The paclet contains some additional documentation. I hope it is useful.


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.