# DiscretizeRegion does not include the boundary specified in ImplicitRegion (10.1)

I am trying to do some 3D PDE solving, and I keep running into problems with boundary conditions because my meshed boundaries are not what I expected. For example, the following code

<< NDSolveFEM
l = 10.;
w = 10.;
h = 5.;
reg = ImplicitRegion[-l/2. <= x <= l/2. && -w/2. <= y <=
w/2. && -h/2. <= z <= h/2., {x, y, z}];
mesh = ToElementMesh[reg, MaxCellMeasure -> {"Area" -> .1}]


gives

ElementMesh[{{-5.00011, 5.00011}, {-5.00011, 5.00011}, {-2.50005, 2.50005}},
{TetrahedronElement["<" 120298 ">"]}]


It's strange, and unexpected to me, that the ElementMesh bounds are not even close to being within machine precision of what I expect them to be: (5, 5, 2.5). This causes problems when I then try to define a boundary condition at, for example, x=l/2. When trying to use NDSolve with such a boundary condition I get the error "No places were found on the boundary where Coordinate was True...." Replacing the boundary specification with the value (copied and pasted) from the element mesh output succeeds.

What is the best way to work around this? I will eventually have more complicated regions, and not being able to specify their boundaries algebraically is going to be a problem. Alternatively, what's the best way to find out what a region's actual boundary is?

• What version of Mathematica are you using? – user21 Apr 19 '16 at 20:54

On my mac doesn't work either.

But you can do

<< NDSolveFEM
l = 10.;
w = 10.;
h = 5.;
reg = Cuboid[{-l/2, -w/2, -h/2}, -{-l/2, -w/2, -h/2}];
mesh = ToElementMesh[reg, MaxCellMeasure -> {"Area" -> .1}]

(*ElementMesh[{{-5., 5.}, {-5., 5.}, {-2.5, 2.5}}, {HexahedronElement[
"<" 16384 ">"]}]*)


Somehow ImplicitRegion works like crap for meshing. Most probably your meshing is also not symmetric and has some ugly artifacts inside. Try to avoid its use and go for RegionUnion/RegionDifference of defined geometries instead. This is one of user21 great hints btw.

• Version 10.4 was a step in the right direction for ImplicitRegions and ToElementMesh, You might want to update and pester your IT department. – user21 Apr 19 '16 at 21:58
• I just need to download it, but I'm afraid some of the code I've written wont work anymore so I'm waiting a bit :P I do like it more to work without the implicit regions, makes the code shorter to write and easier to read. – tsuresuregusa Apr 19 '16 at 22:04
• Well, perhaps some of your code will actually work better ;-) – user21 Apr 19 '16 at 22:10
• As you pointed out, not only does this fix the boundary issue, but it cleans up the sharp edges of the mesh, resulting in a better FEM solution overall. Now I just need to think of clever ways to do this for more complicated geometries. – djphd Apr 20 '16 at 1:39
• I was thinking to make a question along the lines of "how to dynamically draw meshes with mathematica?", maybe somebody is interested on making a mouse-powered interface that lets you pick figures and move them around. What are the other geometries you want to work? I have a few myself I want to try but haven't had time yet. – tsuresuregusa Apr 20 '16 at 2:49

With version 10.4.1 on Linux I get

<< NDSolveFEM
l = 10.;
w = 10.;
h = 5.;
reg = ImplicitRegion[-l/2. <= x <= l/2. && -w/2. <= y <=
w/2. && -h/2. <= z <= h/2., {x, y, z}];
mesh = ToElementMesh[reg, MaxCellMeasure -> {"Area" -> .1}];
mesh["Bounds"]
{{-5., 5.}, {-5., 5.}, {-2.5, 2.5}}


Other alternatives are to use a small band around the boundary in the predicate like

x>=5-10^-3


As another alternative one could

ElementIncidents[mesh["BoundaryElements"]]


and then get the coordinates of the surfaces

NDSolveFEMGetElementCoordinates[mesh["Coordinates"],
Join @@ ElementIncidents[mesh["BoundaryElements"]]]


What you could also do, is add markers (see 1, 2, 3) to the boundary and use those in the boundary conditions. Then the boundary conditions will always match; but you'd need to compute how far the boundary may be off.

• Adding a buffer is nice, but it would be nicer still if I knew exactly how far off it is since that difference will affect the accuracy of the result. If I'm stuck with this weird meshing issue then the next best thing would be a way to find the actual boundary elements nearest to the desired boundary and set the condition along those elements expressly. Any suggestions? – djphd Apr 19 '16 at 21:34
• @djphd, added suggestions – user21 Apr 19 '16 at 21:56
• That GetElementCoordinates is a nice trick that's buried a little too deep in the docs for me to have hit yet. I really wish I could accept two answers on this one. Thanks for your input! – djphd Apr 20 '16 at 1:38