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.