Improve the mesh smoothing procedure

Question

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];

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

 Needs["NDSolve`FEM`"]

Clear[findInteriorNodes]
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]]
   },
  DeleteCases[
   Flatten@Select[edges, MemberQ[#, n] &],
   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 *)
ClearAll[smoothMesh]
Options[smoothMesh] = {"SmoothingLoops" -> 1};
smoothMesh[mesh_ElementMesh, opts : OptionsPattern[]] := Module[
  {originalMesh, newMesh, nLoops, interiorNodes, modCrds},
  nLoops = OptionValue["SmoothingLoops"];
  interiorNodes = findInteriorNodes[mesh];
  originalMesh = mesh;

  Do[
   modCrds = changeNodePosition[originalMesh, #] & /@ interiorNodes;
   newMesh = ToElementMesh[
     "Coordinates" -> ReplacePart[originalMesh["Coordinates"], modCrds],
     "MeshElements" -> mesh["MeshElements"]
     ];
   originalMesh = newMesh,
   {nLoops}
   ];
  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 &,
    measure
    ];
  names = "Mesh " <> ToString[#] & /@ Range[Length[meshes]];
  Print@NumberForm[
    TableForm[statistics, TableHeadings -> {measure, names}],
    3];
  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];
Row[{
  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?

Answer 1

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.

Needs["NDSolve`FEM`"]
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[
    Transpose[
      Developer`ToPackedArray[{temp2[[1]], 
        Range[Length[temp2[[1]]]]}]] -> 1];
  With[{lookup = sa, ids = temp2[[2]]},
   Function[theseNodes, 
    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"];
  Do[
   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[
   DiscretizeGraphics[
    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]

compareMeshQuality[{smoothedMesh, m}]

A second example:

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

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

compareMeshQuality[{smoothedMesh, mq1}]

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.

Update

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]];
  (*adjacencymatrix2==adjacencymatrix*)

  mass2 = DiagonalMatrix[SparseArray[Total[adjacencymatrix2, {2}]]];
  (*mass2\[Equal]mass*)
  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?)

smoothMesh2["Wireframe"]

Note that we used a second order mesh:

smoothMesh2["MeshOrder"]
2

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]
{{49}}

Now things become tricky.

mesh = ToElementMesh[
   ImplicitRegion[
    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}]

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

mesh = ToElementMesh[
   ImplicitRegion[
    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];
(*compareMeshQuality[{sm,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}]

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.

Answer 2

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

Needs["NDSolve`FEM`"]
R = ToElementMesh[DiscretizeGraphics[
    GraphicsComplex[
     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]]];
   GroupBy[
    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]

Testing the quality with the method provided by the OP yields

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.

Answer 3

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

ResourceFunction["FEMAddOnsInstall"][]

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

PacletUninstall["FEMAddOns"]

Once you have that you can load the package:

Needs["FEMAddOns`"]

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

m = ToElementMesh[
   DiscretizeGraphics[
    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}]

Alternatively you can use an iterative method:

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

compareMeshQuality[{smoothedMesh2, m}]

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