How to mesh a region with inclusions touching the boundary

Question

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.

Answer 1

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

Needs["NDSolve`FEM`"]
(* 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.

Answer 2

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