8
$\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$
  • $\begingroup$ DiscretizeRegion[inclusions] and DiscretizeRegion[RegionDifference[region, inclusions]]? $\endgroup$ – kglr Oct 3 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$ – Mauricio Fernández Oct 3 at 18:19
  • $\begingroup$ Show[di, dr] preserves the meshing of di and dr. $\endgroup$ – kglr Oct 3 at 18:25
  • $\begingroup$ Yes, but later on I would like do a fem computation on the discretized region. $\endgroup$ – Mauricio Fernández Oct 3 at 18:29
9
$\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$
  • $\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$ – Mauricio Fernández Oct 5 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 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 at 7:46
3
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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