# ToElementMesh ignores Boundaries and MeshRefinemet

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]] ]]}


• I added an update with the simple data. – user21 Aug 24 '18 at 13:29

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["NDSolveFEM"]
bmesh2 = ToBoundaryMesh["Coordinates" -> bmesh["Coordinates"],
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["NDSolveFEM"]
bmesh = ToBoundaryMesh["Coordinates" -> easyPoints,
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?

• I uploaded my initial coordinates here: link {boundaryPoints, elements} = Import[NotebookDirectory[] <> "pointsNelements.mx"]; With my original data I have trouble recreating the problem with the missing faces. They seem to be all connected. However, and I don't know if this might be the problem, they are areas connected on the boundary. Does ToElementMesh ignore subdivisions of areas in 3D? Do i have to create multiple 3D volumes as boundary? – jufo Aug 24 '18 at 8:40
• @jufo, unfortunately, the data you send is not complete. Please include the link to the complete data in the actual questions. Try to start with a simple geometry and add more subdivisions to your pit later. – user21 Aug 24 '18 at 9:46
• Thank you very much for the help. This solves the first part of my problem. Maybe the second part will be much easier now. I'm doing research on the impact of pit thermal energy stroages (large water basins to store solar thermal energy seasonaly) on the groundwater. The software I use for the simulation (FEFLOW) can't generate unstructured meshes by itself. I need to maintain the subdivision of the boundary inside the pit, because this is where I want to apply my boundary conditions from measured operation data or co-simulations. – jufo Aug 24 '18 at 13:57
• @jufo, interesting. What is the PDE model for that? Do you perhaps have a paper/book to recommend I could look at? – user21 Aug 24 '18 at 14:02
• Here is a description of PTES technology: link This is the FEFLOW book: link If you don't have access to it, there is a comprehensive collection of white papers on the FEFLOW methodology. Vol. 1 probably being the one you look for: link I personally use mostly heat and sometimes fluid flow. Mass transport is covered as well. – jufo Aug 24 '18 at 16:08

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,
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.

• I can add an example that shows that one needs to define boundary element markers if one wants to avoid the merging of boundary facets. Is there anything else that should be there? – user21 Aug 27 '18 at 5:30
• I added the example to the ElementMesh Generation Tutorial. Pending a review it will be available in a next version of Mathematica (>V11.3). If there is anything else you think should be added to the documentation let me know. – user21 Aug 27 '18 at 8:38
• Thank you, I guess this will be very valuable. Another thing would be some documentation on ElementMarkers, which right now comes kind of like out of the blue. Another inconvenience I encountered was the extraction of the Points of the different BoundaryRegions. Extracting the Elements with ElementMarkers by the defined "RegionsMarkers" was straight forward, but I couldn't extract the indices of points by "BoundaryRegionMarkers" easily since I had to do something nested like pointElementList[Position[ElementMarkers[pointElementList], markerID]]. – jufo Aug 27 '18 at 9:18
• ...but after I figured out I couldn't use ElementMarkers[mesh["PointElements"]]straight away it was quite convenient. If I couldn't make myself clear I could deliver a example as well. The result turned out to be very nice. Keep up the good work THUMBSUP – jufo Aug 27 '18 at 10:17
• Not sure if you are aware there is a section on Markers in the ElementMesh Generation Tutorial. Just in case you had not seen that. If have an example that you could send me then I can have a look what I can add to the documentation. – user21 Aug 27 '18 at 14:32