FEM: how to redo a composed boundary mesh for better quality?

Suppose you have region with 2 materials: one embedding material (matrix) containing some particles (inclusions) of a second material. The particles are distributed randomly in the material. In order to solve, e.g., the heat conduction problem with FEM, I first need to create a mesh of the problem, such that the particles have to be meshed accordingly. For that I already got help in this question and @user21 was kind enough to wrap some very useful routines and provide them through the paclet FEMAddOns. I have been playing around with it again, but I run into the problem of (I think) boundary mesh quality as the start for the generation of the region mesh. For instance, in 2D I use two particles, one exceeding the region boundary, as shown below.

The mesh can be generated with the routine GenMesh[2,1] given further below. The first argument in GenMesh controls the MaxBoundaryCellMeasure in the region boundary mesh, the second the one for the inclusions. If you try to run GenMesh[4,1] the boundary meshes look as follows

but the region mesh generation crushes in my small laptop (8GB RAM are consumed completely and Mathematica crushes), even though the geometry does not seem to be too complicated here. I can only suspect that in GenMesh[4,1] the boundary mesh of the inclusions has the red highlighted point which "sadly" almost collides with the blue one of the region boundary mesh, such that the region mesh generator is forced to generate a ridiculous amount of small elements.

Is this interpretation correct? If so, then the composed boundary mesh (composition of the inclusions boundary mesh and of the region boundary mesh) needs to be redone. Is there any way to do this in an automatic manner with some sensible control over the mesh cell size?

Later on I want to put more randomly distributed inclusions in 2D and 3D (e.g., ellipsoidal inclusions in a rectangle), such that trying arbitrary combinations in mesh generation settings like MaxBoundaryCellMeasure will be very inconvenient. But first, I want to understand this in this small example. Thank you!

--

Routine GenMesh:

1. Region definition
2. Numerical region and boundary meshes
3. Create composed boundary mesh
4. Create region mesh from composed boundary mesh

Implementation

Needs["NDSolveFEM"];
Needs["FEMAddOns"];
GenMesh[bmregionMax_, bminclusionsMax_] := Block[
{},

(*[1] Region definitions*)
L1 = 10;
L2 = 5;
region = Rectangle[{0, 0}, {L1, L2}];
inclusioncenters = {{2, 2}, {7, 4.5}};
inclusions =
RegionIntersection[RegionUnion[Disk[#, 1] & /@ inclusioncenters],
region];
matrix = RegionDifference[region, inclusions];

(*[2] Numerical region and boundary meshes*)
nregion = ToNumericalRegion@region;
regionbounds = RegionBounds@region;
bm1 = ToBoundaryMesh[region, regionbounds,
"MaxBoundaryCellMeasure" -> bmregionMax];
bm2 = ToBoundaryMesh[inclusions, regionbounds,
"MaxBoundaryCellMeasure" -> bminclusionsMax];
Print["Generated boundary mesh for inclusions"];
Print@Show[
bm2["Wireframe"]
, Graphics@Point@bm2["Coordinates"]
];

(*[3] Create composed boundary mesh by joining bm1 of region and bm2 of \
bm = BoundaryElementMeshJoin[bm1, bm2];
Print["Joined boundary mesh"];
Print@Show[
bm["Wireframe"]
, Graphics@Point@bm["Coordinates"]
];

(*[4] Create region mesh from composed boundary mesh*)
SetNumericalRegionElementMesh[nregion, bm];
mesh = ToElementMesh[
nregion
, "RegionHoles" -> None
, "RegionMarker" -> Join[{#, 2, 5} & /@ inclusioncenters]
];
Print["Generated mesh with RegionMarkers"];
Print@mesh[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Blue], FaceForm[Orange]}]];
];

• If we run this script, then Mma will run out of memory and break. What is the reason for that. Unfortunately, I was not able to find a workaround. Jun 4, 2021 at 13:51

I believe your hypothesis is correct that you have nodes that are shared or almost shared between your inclusions and your matrix causing the code to fault (I have 128 GB and the kernel crashes long before that limit is reached).

I modified the code to remove the matrix nodes that were contained within the inclusion regions. Now, your test case as well as a number of other cases (although I was not exhaustive) run without crashing.

Needs["NDSolveFEM"];
Clear[GenMesh2]
GenMesh2[bmregionMax_, bminclusionsMax_] :=
Block[{},(*[1] Region definitions*)L1 = 10;
L2 = 5;
region = Rectangle[{0, 0}, {L1, L2}];
inclusioncenters = {{2, 2}, {7, 4.5}};
inclusions =
RegionIntersection[RegionUnion[Disk[#, 1] & /@ inclusioncenters],
region];
RegionUnion[Disk[#, 1.005] & /@ inclusioncenters];
matrix = RegionDifference[region, inclusions];
(*[2] Numerical region and boundary meshes*)
nregion = ToNumericalRegion@region;
regionbounds = RegionBounds@region;
bm1 = ToBoundaryMesh[region, regionbounds,
"MaxBoundaryCellMeasure" -> bmregionMax];
bm2 = ToBoundaryMesh[inclusions, regionbounds,
"MaxBoundaryCellMeasure" -> bminclusionsMax];
(*Remove nodes from bm1 that are shared with bm2*)
newCrd = Pick[#, Thread[! (rmf /@ #)]] &@bm1["Coordinates"];
inc = Partition[FindShortestTour[newCrd][[2]], 2, 1];
bm3 = ToBoundaryMesh["Coordinates" -> newCrd,
"BoundaryElements" -> {LineElement[inc]}];
Print["Generated boundary mesh for inclusions"];
Print@Show[bm2["Wireframe"], Graphics@Point@bm2["Coordinates"]];
(*[3] Create composed boundary mesh by joining bm1 of region and \
bm = BoundaryElementMeshJoin[bm2, bm3];
Print["Joined boundary mesh"];
Print@Show[bm["Wireframe"], Graphics@Point@bm["Coordinates"]];
(*[4] Create region mesh from composed boundary mesh*)
SetNumericalRegionElementMesh[nregion, bm];
mesh =
ToElementMesh[nregion, "MaxBoundaryCellMeasure" -> bmregionMax,
"RegionHoles" -> None,
"RegionMarker" -> Join[{#, 2, 5} & /@ inclusioncenters]];
Print["Generated mesh with RegionMarkers"];
Print@mesh[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Blue], FaceForm[Orange]}]];];
GenMesh2[4, 1]


Inclusions in 3D

In response Mauricio's question in the comments about extending the workflow to 3D, I indicated that it would be nontrivial. I suggested that the OpenCascadeLink might be a good starting point to begin such an endeavor. Here is a workflow to get you started with an OpenCascade approach.

The following model has two ellipsoids in a cube matrix. One of the ellipsoids intersects with three faces of the cube. Here is an example workflow:

(*Import necessary OpenCascadeLink package*)
Needs["OpenCascadeLink"]
(*Create two ellipsoids and a cube*)
shapeell1 =
el1 = Ellipsoid[{0, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}]];
shapeell2 =
el2 = Ellipsoid[{4, 4, 4}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}]];
shapecube = OpenCascadeShape[cu1 = Cuboid[{0, 0, 0}, {8, 8, 8}]];
(*Find Boolean intersection*)
(*Create boundary meshes*)
(*Joined boundary meshes*)
bmeshj = BoundaryElementMeshJoin[bmeshint, bmeshell2, bmeshcube];
mesh = ToElementMesh[bmeshj, "RegionHoles" -> None,
"RegionMarker" -> {{{0.1, 0.1, 0.1}, 1}, {{4, 4, 4},
2}, {{7.9, 7.9, 7.9}, 3}}];
mesh["Wireframe"]


The following will display each volumetric mesh region.

(*Display each region*)
parts = Map[
mesh["Wireframe"[ElementMarker == #[[1]],
"MeshElement" -> "MeshElements",
"ElementMeshDirective" ->
Directive[EdgeForm[], FaceForm[#[[2]]]]]] &, {{1, Gray}, {2,
Pink}, {3, Orange}}]


Print["Internal ellipsoid and matrix"];
Show[parts, PlotRange -> {{4, 8}, {0, 8}, All}]
Print["Ellipsoid that intersects domain boundary"];
Show[parts, PlotRange -> {All, {0, 8}, All}]


To take it to the next level of being able to add refinements to the different surfaces and regions still could be challenging, but this is a start.

• Thank you very much! One small question: you gather the data for the new line boundary elements in your inc which you then use in "BoundaryElements" -> {LineElement[inc]}. Just in principle, for 3D cases (e.g., ellipsoid in a rectangle), the approach should be extendable, right? Of course only using TriangleElement or QuadElement. Or do you thing I should look out for additional things? Apr 25, 2021 at 9:57
• @MauricioFernández You are welcome! Just looking at the BoundaryElementMeshJoin documentation (the 3D cylinder inclusion in a cube), I do not think the extension to 3D will be trivial. You get many high aspect ratio triangles at the cube surface that can extend far from the cylinder surface. I would investigate the capabilities of OpenCascadeLink. It does a pretty good job of snapping to feature edges with Boolean operations. Good feature detection could aid in extending the approach. Apr 25, 2021 at 17:18
• A small suggestion: I think you want to use ToElementMesh[nregion... in place of ToElementMesh[bm... in the first part of your code. Apr 26, 2021 at 5:43
• @TimLaska super nice! Thanks a lot for the extension, that is already a huge help! Apr 26, 2021 at 11:46
• @ABCDEMMM Does GenMesh2[\[Pi]/12, \[Pi]/48]` give you what you need? Jun 3, 2021 at 0:36