3
$\begingroup$

This question is about ToBoundaryMesh. I want to add additional border line to close the hole in region on picture below and then fill it with use of RegionMarkeroption, to get 2 different domains and mesh it with ToElementMesh. Nodes on the border between these two domains needs to connect.

 << NDSolve`FEM`
 reg = ImplicitRegion[(x - 1)^2 + (y - 1)^2 >= 2, {{x, -2, 2}, {y, -2, 2}}]
 RegionPlot[reg]

enter image description here

    bm = ToBoundaryMesh[reg, "RegionHoles" -> None,
    "RegionMarker" -> {{{-1,-1},1}, {{1, 1},2}},
    "MaxBoundaryCellMeasure" -> 1];
    bm2 = ToElementMesh[bm, "RegionHoles" -> None, 
    "RegionMarker" -> {{{-1, -1}, 1}, {{1, 1},2}},
    "MaxBoundaryCellMeasure" -> 1];
    bm2["Wireframe"]

enter image description here

I managed to solve problem, with adding corner points, searching neighbouring points on border and then adding the connections, then I did ToBoundaryMesh again with newcoordinates and connections. My solution is complicated and maybe someone knows about simpler option, maybe it already exists inside ToBoundaryMesh.

ToBoundaryMesh["Coordinates" -> newcoordinates, 
  "BoundaryElements" -> connections, "RegionHoles" -> None];

I also need to transfer my problem to 3D, and there searching for this points and connecting them gets really complicated. So basically I would just like simple way to enforce master boundary or something, in this case rectangle, or in 3D Hexahedron.

$\endgroup$
5
$\begingroup$

I assume you want to do some FEA with. NDSolve does what you want automatically; it detects the discontinuity. As an example:

uif = NDSolveValue[{Laplacian[u[x, y], {x, y}] + 
      If[(x - 1)^2 + (y - 1)^2 >= 2, 1, 2]*u[x, y] == 0, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} \[Element] 
    ImplicitRegion[-2 <= x <= 2 && -2 <= y <= 2, {x, y}]];

If we look at the mesh it generated we see:

(mesh = uif["ElementMesh"])["Wireframe"]

enter image description here

It's a bit hard to see so we look at the boundary mesh:

ToBoundaryMesh[mesh]["Wireframe"]

enter image description here

If you are not using NDSolve it's a bit harder and hackier:

Needs["NDSolve`FEM`"]
boundaryMeshGenerator[region_, opts : OptionsPattern[]] := 
 Module[{pvars, const},
  pvars = region["PredicateVariables"];
  const = {NDSolve`DiscontinuitySurface[1, 
      2 - (-1 + x)^2 - (-1 + y)^2 == 0, True]} /. 
    Thread[Rule[{x, y}, pvars]];
  ContinuationBoundaryMeshGenerator[region, 
   Flatten[{{opts}, "RegionConstraints" -> const}]]
  ]

ToBoundaryMesh[Rectangle[{-2, -2}, {2, 2}], 
  "BoundaryMeshGenerator" -> boundaryMeshGenerator]["Wireframe"]

The advance of this one is that now even the internal curve is second order accurate if you use

ToElementMesh[Rectangle[{-2, -2}, {2, 2}], 
  "BoundaryMeshGenerator" -> boundaryMeshGenerator]["Wireframe"]

Since there is on documented way of specifying region discontinuities, we intercept the boundary mesh generation and add an undocumented

NDSolve`DiscontinuitySurface[1, 2 - (-1 + x)^2 - (-1 + y)^2 == 0, True]

1 - means the dimension (not the embedding dimension), then comes the constraint and the third argument - I don't recall.

Using Alexei's code:

em = ToElementMesh[Rectangle[{-2, -2}, {2, 2}], 
   "RegionMarker" -> {{{-1, -1}, 1}, {{1, 1}, 2}}, 
   "BoundaryMeshGenerator" -> boundaryMeshGenerator];
Show[{em["Wireframe"[
    "MeshElementStyle" -> {FaceForm[LightGray], FaceForm[LightYellow],
       FaceForm[LightRed], FaceForm[LightGreen], 
      FaceForm[LightMagenta]}]], 
  em["Wireframe"["MeshElementMarkerStyle" -> Black]]}]

enter image description here

$\endgroup$
4
$\begingroup$

You can do the following. These are the boundary coordinates including the points where the circle crosses the square:

squareCoords = {{-2, -2}, {-2, 2}, {0, 2}, {2, 2}, {2, 0}, {2, -2}};

Here are the points of the circle:

pts = Table[{1 + Sqrt[2]*Cos[\[CurlyPhi]], 
    1 + Sqrt[2]*Sin[\[CurlyPhi]]}, {\[CurlyPhi], 3 \[Pi]/4, 
    7 \[Pi]/4, \[Pi]/40}];

One can make sure that this yields the desired domain:

Graphics[{Line[Join[squareCoords, {{-2, -2}}]], Line[pts]}]

enter image description here

Now let us create the list of coordinates by joining these two and the list of incidents:

coordinates = Join[coords, pts1];
incidents = 
  Join[Join[Table[{i, i + 1}, {i, 1, 5}], {{6, 1}}], 
   Join[{{3, 7}}, 
    Table[{i, i + 1}, {i, 7, 6 + Length[pts1] - 2}], {{46, 5}}]];

Now one can make a boundary mesh:

<< NDSolve`FEM`

bm = ToBoundaryMesh["Coordinates" -> coordinates, 
  "BoundaryElements" -> {LineElement[incidents]}]
bm["Wireframe"]

yielding this boundary mesh:

enter image description here

which looks exactly as the previous picture, as it should be. After that we can make the element mesh:

em = ToElementMesh[bm, "RegionMarker" -> {{{-1, -1}, 1}, {{1, 1}, 2}},
   MaxCellMeasure -> 0.1]

and let us now have a look at it:

Show[{
  em["Wireframe"[
    "MeshElementStyle" -> {FaceForm[LightGray], FaceForm[LightYellow],
       FaceForm[LightRed], FaceForm[LightGreen], 
      FaceForm[LightMagenta]}]],
  em["Wireframe"["MeshElementMarkerStyle" -> Black]]
  }]

enter image description here

Done. Have fun!

$\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.