14

I am trying to create a 2d region consisting of two subregions. The inner region has several holes, where boundary conditions are applied. The figure shows the idea.

Geometry with two regions

I have tried to create this region using various region functions, but without success. The only approach that has worked so far is to replace the circles with many-sided polygons and use ToBoundaryMesh, specifying all of the points and lines that make up the polygons. This approach allows an ElementMesh to be generated, but it seems overly complex for such a simple geometry. In addition, the solution I get to Laplace's equation for this mesh looks unphysical.

Here is the mesh generated this way:

Element Mesh of region

The code to generate the mesh:

Needs["NDSolve`FEM`"]
nSides = 48;  (*number of sides for each circle*)
nCircles = 12; (*number of circles in the model *)
(*Function to create points for a single circle:*)
circlePts[{x_, y_}, r_] :=
Map[{x + r*Cos[#], y + r*Sin[#]} &, Range[0, 2 π - (2 π)/nSides, (2 π)/nSides]]

(*Generate list of coordinates and connectivity*)
cPts = Flatten[Map[circlePts[#[[1]], #[[2]]] &, circList], 1]; 
connect = Partition[Riffle[Range[nSides], RotateLeft[Range[nSides]]], 2];
nn = (Range[nCircles] - 1)*nSides;
bigConnect = Map[LineElement[connect + #] &, nn];

(*Create boundary mesh*)
bmesh = ToBoundaryMesh[
"Coordinates" -> cPts, "BoundaryElements" -> bigConnect];
bmesh["Wireframe"]

(*Create 2d mesh*)
mesh = ToElementMesh[bmesh, "RegionHoles" -> circList[[3 ;; -1, 1]]];
mesh["Wireframe"]

(*Set up boundary conditions*)
bcs = Join[{DirichletCondition[u[x, y] == 0, x^2 + y^2 >= 0.150^2]}, 
Map[DirichletCondition[
 u[x, y] == 
  20, (x - #[[1, 1]])^2 + (y - #[[1, 2]])^2 == #[[2]]^2] &, 
circList]
]

(*Solve the model*)
op=-Laplacian[u[x,y],{x,y}];
Subscript[Γ, D]=bcs;
uif=NDSolveValue[{op==0,Subscript[Γ, D]},u,{x,y}∈mesh]

ContourPlot[uif[x, y], {x, y} ∈ mesh, 
ColorFunction -> "Temperature", AspectRatio -> Automatic, 
PlotRange -> All, Contours -> 10]

Contour plot of solution

Is there a more elegant way to create subregions?

Öskå
  • 8,587
  • 4
  • 30
  • 49
Jack McInerney
  • 395
  • 1
  • 10
  • This is a nice question. I am going to suggest that creating compound regions (multi material regions) be improved. Since your PDE is not dependent on region material it's not necessary to include the interior circular region. If the PDE were dependent on the sub-region positions it were nice if NDSolve were to detect that and auto include that internal region. But, yes, I agree creating compound regions could be easier. For now you'd need to do this manually. – user21 Oct 14 '14 at 15:23
  • @user21 I think if the ability to create compound regions was simplified, Mathematica could replace a lot of modeling software, at least for simple models. There is a huge potential here. I think only one additional command would be needed. Something like SplitRegion or ImprintRegion. Of course implementing such a command is another story. – Jack McInerney Oct 14 '14 at 17:06

2 Answers2

13

Here is a clean way to do it. The idea is to specify all circular boundary regions and then to explicitly set those that are region holes, such that those are not meshed.

Needs["NDSolve`FEM`"]
\[CapitalOmega] = ImplicitRegion[(9/10)^2 <= x^2 + y^2 <= 1^2, {x, y}];
ToElementMesh[\[CapitalOmega]]["Wireframe"]

enter image description here

ToElementMesh[\[CapitalOmega], "RegionHoles" -> None]["Wireframe"]

enter image description here

disk[{x0_, y0_}, r_] := ((x + x0)^2 + (y + y0)^2 <= (r)^2)
crds = {{-1/2, 0}, {1/2, 1/2}};
sd = Or @@ (disk[#, 1/8] & /@ crds);
\[CapitalOmega]2 = 
  ImplicitRegion[Or[(9/10)^2 <= x^2 + y^2 <= 1^2, sd], {x, y}];
ToElementMesh[\[CapitalOmega]2, 
  "BoundaryMeshGenerator" -> {"Continuation"}]["Wireframe"]

enter image description here

ToElementMesh[\[CapitalOmega]2, 
  "BoundaryMeshGenerator" -> {"Continuation"}, 
  "RegionHoles" -> -crds]["Wireframe"]

enter image description here

You could also refine one of the sub regions:

(mesh = ToElementMesh[\[CapitalOmega]2, 
    "BoundaryMeshGenerator" -> {"Continuation"}, 
    "RegionHoles" -> -crds, 
    "RegionMarker" -> {{{0, 0}, 1, 0.01}, {{19/20, 0}, 2, 
       0.001}}])["Wireframe"]

enter image description here

And visualize with different colors

mesh["Wireframe"[
  "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red]}]]

enter image description here

But strictly speaking the sub regions as not necessary. as the PDE does not have a discontinuity.

user21
  • 39,710
  • 8
  • 110
  • 167
  • Thanks, this is really helpful. I hope that Wolfram clarifies the documentation on this. I find their implementation both unintuitive and quite different from other finite element solvers. – Jack McInerney Oct 26 '14 at 23:27
  • @JackMcInerney, I am one of the developers that did implement this and wrote the documentation. If you do not mind, I'd take this example and put it in the documentation. I assume you have seen the ElementMesh Generation Tutorial and the ElementMesh visualization tutorial – user21 Oct 27 '14 at 12:12
  • @JackMcInerney, if you have feedback on the documentation please let me know. I'd be happy to explain things better if needed. – user21 Oct 27 '14 at 12:15
  • please feel free to use this example in any way you like. I would suggest creating a position dependent conductivity k so the outer ring has a different value than the center. Then use an operator like -Div[k*Grad[u[x,y]...when solving the PDE. – Jack McInerney Oct 27 '14 at 20:31
  • On the documentation, I did go through the ElementMesh tutorials when I first start working on this. I found them somewhat hard to follow. I guess I am use to commercial software for PDEs and haven't yet gotten my head around the Wolfram way of tackling differential equations. I think having lots of real world examples in the documentation would help a lot. Comsol Multiphysics does this and it allows you to piece together a model from a few of their examples. – Jack McInerney Oct 27 '14 at 20:53
8

I will attempt to answer my own question. Technical support at Wolfram Research pointed out that the boundary conditions specified above are incorrect. Specifically the two outer circles have the wrong boundary conditions. The correct BCs are given by

bcs = Join[{DirichletCondition[u[x, y] == 0, x^2 + y^2 >= 0.149^2]},
Map[DirichletCondition[u[x, y] == 20, 
(x - #[[1, 1]])^2 + (y - #[[1, 2]])^2 == #[[2]]^2] &, 
circList[[3 ;;]]]
]

With that, the model runs correctly and produces this result

pressure contour

On the bigger question (is their a simple way to create a region that has both sub regions and holes) the answer appears to be no. The sense I got from Wolfram Support is that such geometries must be spelled out in a point-by-point way, with all connectivity laid out, as was done in the question above. Hopefully the next release will contain some improvements in this area.

Jack McInerney
  • 395
  • 1
  • 10