# How to mesh a region with inclusions touching the boundary

Suppose you have a region composed of two materials. You define the inclusions of material 1 through, e.g., disks, ellipsoids, rectangle and so on. Material 2 is the region minus the inclusions. The inclusions may touch the boundary of the region, see following example.

region = Rectangle[{0, 0}, {10, 5}];
inclusions =
RegionUnion[Disk[{3, 2}, 1.5], Rectangle[{9, 4}, {10, 5}]];
Show[{RegionPlot@region, RegionPlot[inclusions, PlotStyle -> Orange]},
AspectRatio -> Automatic] How do you mesh such a region? By that I mean creating a mesh for the whole region with a "nice" boundary mesh for the inclusions of material 1 (inclusions). I have been trying to read through the meshing tutorials but I do not understand which functions are to be used. Later on I want to compute the solution of a PDE with FEM with varying coefficients, e.g.,

c[x_, y_] :=
If[Element[{x, y}, inclusions], DiagonalMatrix@{100, 20},
DiagonalMatrix@{3, 2}]
pde = Inactive[Div][
c[x, y].Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = {
DirichletCondition[u[x, y] == 0, x == 0]
, DirichletCondition[u[x, y] == 100, x == 10]
};
usol = NDSolveValue[{pde, bc}, u, Element[{x, y}, region]];
ContourPlot[usol[x, y], Element[{x, y}, region],
AspectRatio -> Automatic] I know I can immediately put the coefficient function c[x,y] into the FEM formulation in Mathematica (as in the code above), but I can surely imagine, that the results will improve if a better, geometry coherent mesh is provided. Thank you!

edit:

sorry, I should have pointed out, that I tried DiscretizeRegion but I do not understand how to combine the individual discretized region to a good mesh. For example

di = DiscretizeRegion[inclusions]
dr = DiscretizeRegion[RegionDifference[region, inclusions]]
RegionUnion[di, dr] is a viable solution but I think that the mesh quality of the region union is not a good one.

• DiscretizeRegion[inclusions] and DiscretizeRegion[RegionDifference[region, inclusions]]? – kglr Oct 3 '19 at 18:09
• sorry, I forgot to mention that I tried that but I do not know how to properly combine these (see edit). – Mauricio Fernández Oct 3 '19 at 18:19
• Show[di, dr] preserves the meshing of di and dr. – kglr Oct 3 '19 at 18:25
• Yes, but later on I would like do a fem computation on the discretized region. – Mauricio Fernández Oct 3 '19 at 18:29

## 2 Answers

I cannot take credit for the following code, but it allows one to join several boundary meshes together.

Needs["NDSolveFEM"]
(* Code to join multiple boundary meshes *)
ClearAll[validInputQ]
validInputQ[bm1_, bm2_] :=
BoundaryElementMeshQ[bm1] &&
BoundaryElementMeshQ[
bm2] && (bm1["EmbeddingDimension"] ===
bm2["EmbeddingDimension"]) && (bm1["MeshOrder"] ===
bm2["MeshOrder"] === 1)
BoundaryElementMeshJoin[bm1_, bm2_,
opts : OptionsPattern[ToBoundaryMesh]] /; validInputQ[bm1, bm2] :=
Module[{c1, c2, nc1, newBCEle, newPEle, eleTypes, markers},
c1 = bm1["Coordinates"];
c2 = bm2["Coordinates"];
nc1 = Length[c1];
newBCEle = bm2["BoundaryElements"];
eleTypes = Head /@ newBCEle;
If[ElementMarkersQ[newBCEle], markers = ElementMarkers[newBCEle],
markers = Sequence[]];
newBCEle =
MapThread[#1[##2] &, {eleTypes, ElementIncidents[newBCEle] + nc1,
markers}];
newPEle = bm2["PointElements"];
eleTypes = Head /@ newPEle;
If[ElementMarkersQ[newPEle], markers = ElementMarkers[newPEle],
markers = Sequence[]];
newPEle =
MapThread[#1[##2] &, {eleTypes, ElementIncidents[newPEle] + nc1,
markers}];
ToBoundaryMesh["Coordinates" -> Join[c1, c2],
"BoundaryElements" -> Flatten[{bm1["BoundaryElements"], newBCEle}],
"PointElements" -> Flatten[{bm1["PointElements"], newPEle}], opts]]
BoundaryElementMeshJoin[r1_, r2_, r3__] :=
BoundaryElementMeshJoin[BoundaryElementMeshJoin[r1, r2], r3];


Now, we can create a domain region and inclusion region boundary meshes and join them together to create a multiregion element mesh.

(* Code to create and join several distinct boundary meshes *)
(* Create Regions *)
srdom = Rectangle[{0, 0}, {10, 5}];
nrdom = ToNumericalRegion[srdom];
srdisk = Disk[{3, 2}, 3/2];
srrect = Rectangle[{9, 4}, {10, 5}];
(* Estimate Region Bounds of Full Domain *)
symbolicBounds = RegionBounds[srdom];
(* Create Boundary Meshes for Each Region *)
(bm1 = ToBoundaryMesh[srdom, symbolicBounds])["Wireframe"];
(bm2 = ToBoundaryMesh[srdisk, symbolicBounds])["Wireframe"];
(bm3 = ToBoundaryMesh[srrect, symbolicBounds])["Wireframe"];
(* Join Boundary Meshes *)
(bm = BoundaryElementMeshJoin[bm1, bm2])["Wireframe"];
(bm = BoundaryElementMeshJoin[bm, bm3])["Wireframe"];
SetNumericalRegionElementMesh[nrdom, bm];
meshTriangle =
ToElementMesh[nrdom,
"RegionMarker" -> {{{1, 1}, 1, 0.5}, {{3, 2}, 2, 0.5}, {{9.5, 4.5},
3, 0.5}}];
meshTriangle[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Red], FaceForm[Green],
FaceForm[Yellow]}, ImageSize -> Medium]] Now, we can solve your PDE using a piecewise function and the region markers for the inclusions and domain defined above.

c = Evaluate[Piecewise[{{DiagonalMatrix@{100, 20}, ElementMarker == 2},
{DiagonalMatrix@{100, 20}, ElementMarker == 3},
{DiagonalMatrix@{3, 2}, True}}]];
pde = Inactive[Div][c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = {DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[u[x, y] == 100, x == 10]};
usol = NDSolveValue[{pde, bc}, u, {x, y} \[Element] meshTriangle];
ContourPlot[usol[x, y], Element[{x, y}, meshTriangle],
AspectRatio -> Automatic] Let's make the disk a resistor.

c = Evaluate[
Piecewise[{{DiagonalMatrix@{0.0003, 0.0002}, ElementMarker == 2},
{DiagonalMatrix@{100, 20}, ElementMarker == 3},
{DiagonalMatrix@{3, 2}, True}}]];
pde = Inactive[Div][c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = {DirichletCondition[u[x, y] == 0, x == 0],
DirichletCondition[u[x, y] == 100, x == 10]};
usol = NDSolveValue[{pde, bc}, u, {x, y} \[Element] meshTriangle];
ContourPlot[usol[x, y], Element[{x, y}, meshTriangle],
PlotPoints -> All, AspectRatio -> Automatic] It seems to be behaving as expected.

• Oh wow, ok I would have definitely not figured that one out, that is way longer from what I expected. Thanks for the code, I will go through it and see if I can simplify it for my needs. – Mauricio Fernández Oct 5 '19 at 14:07
• Mathematica is great for setting up quick and dirty simple models. As your model builds in complexity, you generally need more control over the meshing and preprocessing steps and I will move quickly to the methods found in the Element Mesh Generation Tutorial. Thank you for the accept. – Tim Laska Oct 5 '19 at 14:31
• @MauricioFernández, I have wrapped BoundaryElementMeshJoin up and put it in the FEMAddOns paclet. For details see my answer below. – user21 Nov 13 '19 at 7:46

This is an extentsion to Tim Laska's answer. I have added the BoundaryElementMeshJoin (and a few other Boolean operations) for boundary element meshes into the FEMAddOns paclet. The installation of the paclet is now very easy since the installation can be done via the FEMAddOnsInstall resource function.

Install and load the paclet:

ResourceFunction["FEMAddOnsInstall"][]
Needs["FEMAddOns"]
(* Paclet[FEMAddOns, 1.3.2] *)


The run the above code as usual:

(*Code to create and join several distinct boundary meshes*)(*Create \
Regions*)srdom = Rectangle[{0, 0}, {10, 5}];
nrdom = ToNumericalRegion[srdom];
srdisk = Disk[{3, 2}, 3/2];
srrect = Rectangle[{9, 4}, {10, 5}];
(*Estimate Region Bounds of Full Domain*)

symbolicBounds = RegionBounds[srdom];
(*Create Boundary Meshes for Each Region*)
(bm1 =
ToBoundaryMesh[srdom, symbolicBounds])["Wireframe"];
(bm2 = ToBoundaryMesh[srdisk, symbolicBounds])["Wireframe"];
(bm3 = ToBoundaryMesh[srrect, symbolicBounds])["Wireframe"];
(*Join Boundary Meshes*)
(bm =
BoundaryElementMeshJoin[bm1, bm2, bm3])["Wireframe"];
SetNumericalRegionElementMesh[nrdom, bm];
meshTriangle =
ToElementMesh[nrdom,
"RegionMarker" -> {{{1, 1}, 1, 0.5}, {{3, 2}, 2, 0.5}, {{9.5, 4.5},
3, 0.5}}];
meshTriangle[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Red], FaceForm[Green],
FaceForm[Yellow]}, ImageSize -> Medium]] • Thanks for the wonderful addon! One question: if you change the size of the rectangle to Rectangle[{9, 4}, {9.8, 4.8}], then it does not touch the boundary and the mesh generation assumes the inclusions are to be removed from the region (see final mesh for new rectangle). Any way to keep the nicely meshed inclusions in the region? – Mauricio Fernández Dec 29 '19 at 17:13
• Ignore my last comment, just found out that ToElementMesh[regiondiff,"RegionHoles"->None]` does the job for me when the inclusions do not touch the boundary based on the region difference of the whole region minus the inclusions. – Mauricio Fernández Dec 29 '19 at 17:34
• @user21, nice example! if we want to refine the circle boundary (not this rectangle), which "RegionMarker" commands should we use? – ABCDEMMM Jun 2 at 20:56
• @ABCDEMMM, I can tell you how to do that. What I need to know though, is what trouble did you encounter when you tried this? How do I need to improve the documentation such that it would work better for you? – user21 Jun 3 at 3:45