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

BoundaryMesh of my pit model

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, 
mesh = ToElementMesh[
   "MeshOrder" -> 1,
   MeshRefinementFunction -> refinementFunction,
   AccuracyGoal -> 10
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Thin], FaceForm[RGBColor[0.16, 1., 0.97]]]]]

top view of the final mesh

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: enter image description here 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)][
   ToBoundaryMesh["Coordinates" -> easyPoints, 
    "BoundaryElements" -> (QuadElement[{#}] & /@ easyFaces)]][
  "Wireframe"["MeshElement" -> "MeshElements", 
   "ElementMeshDirective" -> Directive[FaceForm[Red]] ]]}

Simple BoundaryMesh and ElementMEsh

  • $\begingroup$ I added an update with the simple data. $\endgroup$
    – user21
    Aug 24, 2018 at 13:29

2 Answers 2


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

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

enter image description here

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.


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:

{easyPoints, easyFaces} = Import["~/Downloads/pointsNelementsEasy.mx"];
bmesh = ToBoundaryMesh["Coordinates" -> easyPoints, 
   "BoundaryElements" -> {QuadElement[easyFaces, 

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

enter image description here

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

enter image description here

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?

  • $\begingroup$ 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? $\endgroup$
    – jufo
    Aug 24, 2018 at 8:40
  • $\begingroup$ @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. $\endgroup$
    – user21
    Aug 24, 2018 at 9:46
  • $\begingroup$ 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. $\endgroup$
    – jufo
    Aug 24, 2018 at 13:57
  • $\begingroup$ @jufo, interesting. What is the PDE model for that? Do you perhaps have a paper/book to recommend I could look at? $\endgroup$
    – user21
    Aug 24, 2018 at 14:02
  • $\begingroup$ 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. $\endgroup$
    – jufo
    Aug 24, 2018 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,
  "BoundaryElements" -> {QuadElement[faces, markers]}];
bmesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementStyle" -> (FaceForm /@ 
     ColorData["Rainbow"] /@ Range[0, 1, 1/13])]]

boundary mesh with markers 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: cross section

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

Final Mesh

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:

enter image description here

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.

  • $\begingroup$ 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? $\endgroup$
    – user21
    Aug 27, 2018 at 5:30
  • $\begingroup$ 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. $\endgroup$
    – user21
    Aug 27, 2018 at 8:38
  • $\begingroup$ 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]]. $\endgroup$
    – jufo
    Aug 27, 2018 at 9:18
  • $\begingroup$ ...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 $\endgroup$
    – jufo
    Aug 27, 2018 at 10:17
  • $\begingroup$ 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. $\endgroup$
    – user21
    Aug 27, 2018 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.