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:
- Region definition
- Numerical region and boundary meshes
- Create composed boundary mesh
- Create region mesh from composed boundary mesh
Implementation
Needs["NDSolve`FEM`"];
(*Install FEMAddOns if not already done*)
(*ResourceFunction["FEMAddOnsInstall"][];*)
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 \
inclusions - thanks to FEMAddOns*)
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
, "BoundaryMeshGenerator" -> "OpenCascade"
, "RegionHoles" -> None
, "RegionMarker" -> Join[{#, 2, 5} & /@ inclusioncenters]
];
Print["Generated mesh with RegionMarkers"];
Print@mesh[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Blue], FaceForm[Orange]}]];
];
