ToElementMesh ignores Boundaries and MeshRefinemet

Question

I'm building a mesh of a pit storage in Mathematica which I then export to another FEM-software. I follow the workflow for creating a ElementMesh from the documentation, but it won't work in the way it should. I start with creating my BoundaryMesh, which works perfectly well:

boundary = ToBoundaryMesh["Coordinates" -> boundaryPoints,"BoundaryElements" enter code here-> {QuadElement[elements]}];
boundary["Wireframe"]

This is the exported BoundaryMesh, which contains the Coordinates and the Elements:

boundary mesh as .off-file

I now want to create the mesh with special focus on the center area. The center should be refined and the boundary areas inside the pit (center hole) should be preserved, because i want to apply different boundary conditions for each segment.

    refinementFunction = 
 With[{damOuterWidthX = 96, damOuterWidthY = 106, pitDepth = 10},
  Compile[{{coordinates, _Real, 2}, {vol, _Real, 0}},
   Block[{pos}, pos = Mean[coordinates];
    If[((-(damOuterWidthX/2 + 10) <= 
          pos[[1]] < (damOuterWidthX/2 + 
            10)) && (-(damOuterWidthY/2 + 10) <= 
          pos[[2]] < (damOuterWidthY/2 + 10)) && (-(pitDepth + 10) <= 
          pos[[3]])) && Volume[Tetrahedron[coordinates]] > 3, True, 
     False]
    ]]];
mesh = ToElementMesh[
   boundary,
   "MeshOrder" -> 1,
   MeshRefinementFunction -> refinementFunction,
   AccuracyGoal -> 10
   ];
mesh["Wireframe"[
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Thin], FaceForm[RGBColor[0.16, 1., 0.97]]]]]

This is the top view of the resulting mesh. As one can see, only the geometry is maintained, but the boundary areas got lost on the way and the refinement, which should go beyond the center area doesn't apply everywhere inside the refinement area. Here is a picture from my FEM software, which makes this more clear: This is a side view of the pit (cut in half), where the colors show the volume size. My volume constraint was set to 3 in the center region. Below the pit this works fine, but it fails at the pit walls. Something else one can see, in this picture is the following. The bottom and the wall of the pit are of a different material and thus I considered this geometry in the BoundaryMesh and in contrast to the divison of the pit wall into different areas, it was maintained.

I played around alot with all the options of ToElementMesh and my RefinementFunction but I can't get it right. If someone has an idea or a solution I would be very grateful.

Update: @user21: Sorry, here are the complete points and faces of my problem: Complete faces and points. I created a more simple case, that basically shows the problem with the ignored Boundaries: Simple case Here as well, an internal Boundary is dropped when the mesh is converted to the ElementMesh:

{easyPoints, easyFaces} = 
  Import[NotebookDirectory[] <> "pointsNelementsEasy.mx"];
{ToBoundaryMesh["Coordinates" -> easyPoints, 
   "BoundaryElements" -> (QuadElement[{#}] & /@ easyFaces)][
  "Wireframe"], 
 ToElementMesh[
   ToBoundaryMesh["Coordinates" -> easyPoints, 
    "BoundaryElements" -> (QuadElement[{#}] & /@ easyFaces)]][
  "Wireframe"["MeshElement" -> "MeshElements", 
   "ElementMeshDirective" -> Directive[FaceForm[Red]] ]]}

Answer 1

This is a partial answer: (The rest is in the update below)

Let's leave mesh refinement out for a second. Some of faces do not seem to be connected. (It seems you changed your logic a bit when constructed the boundary mesh.)

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[Import["~/Downloads/boundary.off"]];
bmesh2 = ToBoundaryMesh["Coordinates" -> bmesh["Coordinates"], 
   "BoundaryElements" -> {QuadElement[
      bmesh["BoundaryElements"][[1, 1]][[{(*6,*)10, 14, 18, 22, 26, 
         27, 28, 40, 47}]]]}];
bmesh2["Wireframe"]

Faces 22 and 28 do not seem to be connected. At least I was not able to find the elements they should connect to.

It would be more useful to have the actual (raw) coordinates and boundary faces than having to extract them from the off file.

I'd try to get the boundary mesh figured out before I start the refinement.

Update

If you want boundaries to not be merged they'd need different boundary element markers. If they have the same marker there is not reason to not merge them; you'd apply the same boundary condition anyways.

To keep the boundaries distinct you could assign a different marker to each face like so:

Needs["NDSolve`FEM`"]
{easyPoints, easyFaces} = Import["~/Downloads/pointsNelementsEasy.mx"];
bmesh = ToBoundaryMesh["Coordinates" -> easyPoints, 
   "BoundaryElements" -> {QuadElement[easyFaces, 
      Range[Length[easyFaces]]]}];
(*bmesh["Wireframe"]*)

bmesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementStyle" -> (FaceForm /@ 
     ColorData["SunsetColors"] /@ 
      Range[0, 1, 1/(Length[easyFaces] - 1)])]]

Each face now has it's own marker. The full mesh:

mesh = ToElementMesh[bmesh, "MaxCellMeasure" -> 100];
mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementStyle" -> (FaceForm /@ 
     ColorData["SunsetColors"] /@ 
      Range[0, 1, 1/(Length[easyFaces] - 1)])]]

Hope that helps. Two more points: if you want to use @ to ping me, you'd have to do that under my answer. I am curious to know what type of FEM analysis do you plan to do?

Answer 2

Okay, with a lot of help from user 21 this actually turned out to be a very neat example for 3d meshing. Here is the solution to my two problems: I started out with the following data: input data

{points, faces, markers, regionMarkerPoints} = 
  Import[NotebookDirectory[] <> "completeData.mx"];

As user21 pointed out, I needed to define boundary element markers if I wanted to avoid the merging of my boundary areas:

bmesh = ToBoundaryMesh[
  "Coordinates" -> points,
  "BoundaryElements" -> {QuadElement[faces, markers]}];
bmesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementStyle" -> (FaceForm /@ 
     ColorData["Rainbow"] /@ Range[0, 1, 1/13])]]

Each boundary with the same marker has the same color. The second problem was the meshing process, where I could't achieve satisfying results with a refinement function. Instead I chose a much more convenient approach now, by dividing my model into different volumes elements and assigning region markers. See the cross section of my model:

There are basically three volumes in the model, now represented by the region marker points (green=ground, red= dam, blue= pit liner). I can put a volume contraint of 80 m³ on the ground region, 1 m³ on the dam and 0.1 m³ on the pit liner:

mesh = ToElementMesh[
   bmesh,
   "MeshOrder" -> 1,
   "RegionMarker" -> {{regionMarkerPoints[[1]], 1, 
      80}, {regionMarkerPoints[[2]], 2, 1}, {regionMarkerPoints[[3]], 
      3, 0.1}}
   ];
mesh["Wireframe"["MeshElement" -> "MashElements", 
  "MeshElementStyle" -> {Directive[FaceForm[Green]], 
    Directive[EdgeForm[], FaceForm[Red]], 
    Directive[EdgeForm[], FaceForm[Blue]]}]]

I turned off the Edges of the very fine mesh in the center of the model, so one can see, that by picking a random point inside my 3 volumes I assigned RegionMarkers to them.

This is the final result after exporting the model to FEFLOW. It still needs further tuning, but fulfils my basic requirements:

Actually the process turned out to be quite simple and elegant, but in my oppinion Mathematica's documentation on 3D applications is often quite short compared to the extensive 2D cases.