3
$\begingroup$

I am trying to figure out how to define and use meshes for FEM and NDSolve. I've figured out that if I want more resolution near my boundaries, I can use ToBoundaryMesh to first create the boundary mesh for my region's boundary, and then pass that boundary mesh to ToElementMesh to create the internal mesh:

mycub = Cuboid[{0, 0, 0}, {10, 10, 10}];
mysmallcube = Cuboid[{4, 4, 4}, {6, 6, 6}];
Graphics3D[{Opacity[.1], mycub, Opacity[.4], mysmallcube}]
myreg = RegionDifference[mycub, mysmallcube];

bdm = NDSolve`FEM`ToBoundaryMesh[myreg, "MaxBoundaryCellMeasure" -> .2]
bdm["Wireframe"[PlotRange -> {All, {4, 6}, All}]]
Length@First@First@bdm["BoundaryElements"]

emesh = NDSolve`FEM`ToElementMesh[bdm, 
  MaxCellMeasure -> {"Length" -> 1}]
emesh["Wireframe"["MeshElement" -> "MeshElements", 
  PlotRange -> {All, {4, 6}, All}]]
Length@First@First@emesh["MeshElements"]
Length@First@First@emesh["BoundaryElements"]

enter image description here enter image description here

Looking good. A nice fine boundary mesh and a rougher interior mesh. There are 11960 boundary elements and 24830 mesh elements.

Now, if I try and do the same exact thing, but just only with ToElementMesh, it appears to not obey the MaxBoundaryCellMeasure:

emesh2 = NDSolve`FEM`ToElementMesh[myreg, 
  MaxCellMeasure -> {"Length" -> 1}, "MaxBoundaryCellMeasure" -> .2]
emesh2["Wireframe"[PlotRange -> {All, {4, 6}, All}]]
emesh2["Wireframe"["MeshElement" -> "MeshElements", 
  PlotRange -> {All, {4, 6}, All}]]
Length@First@First@emesh2["MeshElements"]
Length@First@First@emesh2["BoundaryElements"]

enter image description here enter image description here

Clearly, the boundary mesh is much rougher.

What's the difference?

One minor related question: If we look at the first image above, but from the side, it looks like this:

enter image description here

If you can see, the corners of the center cube are rounded. Why?

edit: I think I may have a clue why. In the section "Region Approximation Quality" in the FEM documentation, it says:

For graphics primitives like Line or Polygon or a MeshRegion, a conversion to an ElementMesh is lossless. For example, an ElementMesh representation of a Rectangle is as exact or inexact as the Rectangle itself represents a region.

Which is why I was a little confused, because my cubes were both graphics primitives, like a Cuboid. However, I think the problem is that when I do RegionDifference[], that volume is no longer one, so it starts to do weird stuff (which is a little strange, because it seems like all the exact information is there).

Anyway, I think I've got around this by separately making the boundary mesh for each Cuboid, and then combining their coordinates and elements to create a new boundary mesh:

mycub = Cuboid[{0, 0, 0}, {10, 10, 10}];
mysmallcube = Cuboid[{4, 4, 4}, {6, 6, 6}];

bigcubmesh = 
  NDSolve`FEM`ToBoundaryMesh[mycub, "MaxBoundaryCellMeasure" -> 1];
smallcubmesh = 
  NDSolve`FEM`ToBoundaryMesh[mysmallcube, 
   "MaxBoundaryCellMeasure" -> .1];

allcoords = 
  Join[bigcubmesh["Coordinates"], smallcubmesh["Coordinates"]];
bigcubcoordsnum = Length@bigcubmesh["Coordinates"];

allels = QuadElement@
   Join[First@First@bigcubmesh["BoundaryElements"], 
    First@First@smallcubmesh["BoundaryElements"] + bigcubcoordsnum];
bm = ToBoundaryMesh["Coordinates" -> allcoords, 
  "BoundaryElements" -> {allels}]

em = ToElementMesh[bm, "MaxCellMeasure" -> 1];
Show[em["Wireframe"["MeshElement" -> "BoundaryElements", 
   PlotRange -> {All, {4, 6}, All}]], 
 Graphics3D[{Opacity[.1], mycub, Red, Opacity[.5], mysmallcube}]]
Show[em["Wireframe"["MeshElement" -> "MeshElements", 
   PlotRange -> {All, {4, 6}, All}]], 
 Graphics3D[{Opacity[.1], mycub, Red, Opacity[.5], mysmallcube}]]

enter image description here

Clearly, no rounding there. So I guess that's the solution for that minor question? I'm sure there's a smarter way to combine the coordinates and elements than why I did, but I couldn't find a function that does it.

$\endgroup$
4
  • 3
    $\begingroup$ I think this is the same issue as in this post $\endgroup$
    – user21
    Jun 4 '16 at 1:33
  • $\begingroup$ Concerning your minor issue, in this case the problem is that the marching cubes algorithm does not catch the interior corners well: BoundaryDiscretizeRegion[myreg, MeshCellStyle -> {{2, All} -> Opacity[0.25]}] $\endgroup$
    – user21
    Jun 4 '16 at 1:48
  • $\begingroup$ @user21 thank you, so my fix at the end for the rounding problem should work fine though right? $\endgroup$ Jun 4 '16 at 19:38
  • $\begingroup$ Yes, that should be fine. $\endgroup$
    – user21
    Jun 4 '16 at 22:16