4
$\begingroup$

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 

<< NDSolve`FEM`
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?

Improve this question
$\endgroup$
1
  • $\begingroup$ What version of Mathematica are you using? $\endgroup$
    – user21
    Apr 19, 2016 at 20:54

2 Answers 2

Reset to default
3
$\begingroup$

On my mac doesn't work either.

But you can do

<< NDSolve`FEM`
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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Version 10.4 was a step in the right direction for ImplicitRegions and ToElementMesh, You might want to update and pester your IT department. $\endgroup$
    – user21
    Apr 19, 2016 at 21:58
  • $\begingroup$ 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. $\endgroup$
    – tsuresuregusa
    Apr 19, 2016 at 22:04
  • $\begingroup$ Well, perhaps some of your code will actually work better ;-) $\endgroup$
    – user21
    Apr 19, 2016 at 22:10
  • $\begingroup$ 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. $\endgroup$
    – djphd
    Apr 20, 2016 at 1:39
  • $\begingroup$ 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. $\endgroup$
    – tsuresuregusa
    Apr 20, 2016 at 2:49
2
$\begingroup$

With version 10.4.1 on Linux I get

<< NDSolve`FEM`
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

NDSolve`FEM`GetElementCoordinates[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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – djphd
    Apr 19, 2016 at 21:34
  • $\begingroup$ @djphd, added suggestions $\endgroup$
    – user21
    Apr 19, 2016 at 21:56
  • $\begingroup$ 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! $\endgroup$
    – djphd
    Apr 20, 2016 at 1:38

Your Answer

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

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