10
$\begingroup$

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]

region with inclusions

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]

fem solution

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]

enter image description here

is a viable solution but I think that the mesh quality of the region union is not a good one.

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

2 Answers 2

10
$\begingroup$

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

Multiregion mesh

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]

Contour Plot

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]

Resistor Contour Plot

It seems to be behaving as expected.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Oct 5, 2019 at 14:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – Tim Laska
    Oct 5, 2019 at 14:31
  • $\begingroup$ @MauricioFernández, I have wrapped BoundaryElementMeshJoin up and put it in the FEMAddOns paclet. For details see my answer below. $\endgroup$
    – user21
    Nov 13, 2019 at 7:46
7
$\begingroup$

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

enter image description here

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$ Dec 29, 2019 at 17:13
  • $\begingroup$ 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. $\endgroup$ Dec 29, 2019 at 17:34
  • $\begingroup$ @user21, nice example! if we want to refine the circle boundary (not this rectangle), which "RegionMarker" commands should we use? $\endgroup$
    – ABCDEMMM
    Jun 2, 2021 at 20:56
  • $\begingroup$ @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? $\endgroup$
    – user21
    Jun 3, 2021 at 3:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.