First some comments. This same calculation works beautifully in 2D with essentially the same code. I have copied it largely unaltered into 3D as part of a migration to 3D so that I can do more complex problems there next. This particular case has an analytic solution (which I calculate by hand).

Edit Below is some code which captures the essence of what I actually want to do, and with most of the complexity I want to achieve. This runs fine, but lacks precision.

FindPsi[mesh_] := Module[{\[Psi], zMnSlv = 0, zmxOfst = 0},
  Print[Now, " Finding \[Psi]"];
  \[Psi] =
   NDSolveValue[{Inactive[Laplacian][\[Psi][x, y, z], {x, y, z}] - 
       Load[x, y, z] == 0,
     DirichletCondition[\[Psi][x, y, z] == 0, dz == 0],
     DirichletCondition[\[Psi][x, y, z] == BcPsi[x, y, z], 
      z == rng]}, \[Psi], {x, y, z} \[Element] mesh, 
    InitialSeeding -> {\[Psi][x, y, dz] == BcPsi[x, y, z]}, 
    Method -> {"FiniteElement"}, PrecisionGoal -> 6];
  Print[Now, " Finished finding \[Psi]"];
Uu[fn_] := Module[{uu},
  Print["Finding Uu"];
  uu = {D[fn[x, y, z], z], -D[fn[x, y, z], x]};
UuTh[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Module[{d, u, dz, m},
  dz = z - h[x, y];
  u = UU0 Log[(dz + z0)/z0];
z0 = 0.0001;
UU0 = 0.10857;
rng = 10.;
h[x_, y_] := 1/(1 + (x - 0.2 y)^2);
dz = z - h[x, y];
BcPsi[x_, y_, z_] := UU0 ((z + z0) Log[(z + z0)/z0] - z);
Load[x_, y_, z_] := UU0/(dz + z0);
solnRegn = 
   z >= h[x, y], {{x, -rng, rng}, {y, -rng, rng}, {z, 0, rng}}];
mesh = ToElementMesh[solnRegn, "MaxCellMeasure" -> 0.20, 
  "MaxBoundaryCellMeasure" -> 0.1]
Print["MeshOrder = ", mesh["MeshOrder"]];
Print[Magnify[mesh["Wireframe"], 1.8]];
psi = FindPsi[mesh];
psi[-3, -2, 2]
uu[x_, y_, z_] = Uu[psi];
uu[-3, -2, 2]
x = -9.;
y = 0;
LogLinearPlot[{uu[x, y, z][[1]], UuTh[x, y, z]}, {z, 0.003, 10.}, 
 PlotRange -> {{0.01, 10}, {0., 1.3}}]

And here is my output:

Complete Output

The blue curve is the output from the FEM calculation. The brown/orange curve is a good approximation to the right curve, calculated analytically. I would like better precision in the numerical solution, without it running out of memory, and without it taking hours to run :-). I don't need 6 digits of precision, but a sensible target would be a relative error of <= 3% along the whole curve, or at least down to z=0.02.

I should add that the function h[x,y] could be any single-valued, continuous surface, which may be inscribed in the box, including possibly an interpolating function.

The steps around z=0.10 would seem to be an artefact of the mesh, while the upturn at z~10 is, I think, due to the bevelled edges of the solution region.

Any advice on improving this would be very gratefully received. Thanks.

Edit Here is another output. Same code, but using "MaxBoundaryCellMeasure" -> 0.015. This uses 991000 mesh elements and takes about 20' on my machine. It is clearly better, but not yet near my desired precision. I don't think I can use more mesh elements without running out of memory.

Another example output

Further Development

I provide below one attempt to use a mesh refinement function. It does not behave as I expect, but presumably, there is a bug in it. The idea here is that for all vertices within the solution volume for which dz<1, the volume is limited to an amount which is proportional to dz itself. The idea of that is to make the elements get smaller and smaller as we reach the boundary at dz=0, where the BCs apply. First the code:

mrf = Function[{vertices, volume}, 
   Block[{x, y, z, dz}, {x, y, z} = Mean[vertices];
    dz = z - h[x, y]; If[dz < 1., volume > 0.002 dz, False]]];
mesh = ToElementMesh[solnRegn, MaxCellMeasure -> 0.25, 
  MeshRefinementFunction -> mrf]
   Show[{Plot3D[h[x, y], {x, -10, 10}, {y, -10, 10}, 
      PlotRange -> {All, {-0.3, 0.3}, {0, 10}}],
     Graphics3D[{PointSize[0.004], Point[mesh[[1]]]}, 
      AspectRatio -> Automatic, Axes -> True, 
      PlotRange -> {All, {-0.3, 0.3}, All}]}]
   , 1.8]];

The rest of the code is unchanged from that posted further above. I show here the plot of part of the mesh volume, which has a dot for each vertex. The brown surface is defined by dz=0.

mesh volume

The thing that surprises me is that despite the fact that there are many more mesh vertices in the region around 0.3 <~ dz < 1, there is a marked gap for dz <~ 0.4, containing just a single "surface" of points not the dense cloud I expected. Why the gap?

Finally, the plot of the quantity I am hoping to extract:

enter image description here

Which is very similar to before.

After Adopting Provided Answer

Just to close the loop, I provide here my test plot after adopting the solution provided by @Tim Laska and adapting the step size and ratio.

enter image description here

It is clearly very much improved. I believe the small residual disagreement in the unstepped part of the curve is "physics", ie. a small mismatch such as this is anticipated between the theoretical solution and the "true" solution, including the cross-overhear z=4 (z is shown on the "x-axis").

  • $\begingroup$ How many elements does this mesh have and how much RAM does your computer have? $\endgroup$
    – user21
    Commented Apr 20, 2020 at 13:07
  • $\begingroup$ @user21 HexahedronElement["<" 1431644 ">"]. 8GB of RAM (2018 MacBook air). I am running MMA 12.1. I think it would work with fewer elements in the 5 "upper" faces of the box. I need (if anything) more on the face at z=0, but I don't know how to control the MaxCellBoundaryMeasure on just one face. $\endgroup$ Commented Apr 20, 2020 at 13:15
  • 1
    $\begingroup$ You are running out of memory. Why don't you post what you really want to do. There is not much point in answering this question as this is not what you actually want to do. $\endgroup$
    – user21
    Commented Apr 20, 2020 at 13:41
  • $\begingroup$ @user21 I have posted a much more complete description of what I am trying to achieve. $\endgroup$ Commented Apr 20, 2020 at 15:14
  • $\begingroup$ I have edited the later post to include a precision target. $\endgroup$ Commented Apr 20, 2020 at 15:25

1 Answer 1


Currently, I would not use Mathematica to mesh complex 3D CAD object consisting of many parts. I would, however, consider importing the node and connectivity information from a commercial mesher and map them ToElementMesh following the Element Mesh Genertation Tutotrial. I would also borrow the many of the techniques of commercial meshers to acheive accuracy and efficiency, such as extrusions and boundary layer meshing, and apply them to smaller problems Mathematica.

Although I am unfamiliar with the OP problem, I would expect that one would desire refinement around the peak to capture curvature and also a "boundary layer" mesh $z_{min}$ I will present a possible workflow to map a refinement plan in $i,j,k$ space to $x,y,z$ cooridnates.

3D Anisotropic Meshing Workflow

Initial Parameters and Mesh Helper Functions

(* Reduced Parameter Set *)
z0 = 0.0001;
UU0 = 0.10857;
rng = 10.;
h[x_, y_] := 1/(1 + (x - 0.2 y)^2);
dz = z - h[x, y];
BcPsi[x_, y_, z_] := UU0 ((z + z0) Log[(z + z0)/z0] - z);
Load[x_, y_, z_] := UU0/(dz + z0);
(* RegionProduct Helper Functions *)
pointsToMesh[data_] :=
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
meshGrowth[x0_, xf_, n_, ratio_] := Module[{k, fac, delta},
  k = Log[ratio]/(n - 1);
  fac = Exp[k];
  delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
  N[{x0}~Join~(x0 + 
      delta Rest@
        FoldList[(#1 + #2) &, 0, 
         PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]
unitMeshGrowth[n_, ratio_] := meshGrowth[1, 0, n, ratio]
unitMeshGrowth2Sided [nhalf_, 
  ratio_] := (1 + Union[-Reverse@#, #])/2 &@
  unitMeshGrowth[nhalf, ratio]

Map Meshing Plan on a Unit Cube

I will make use of the Tensor Product Grid example in the RegionProduct documentation to create 3D regions.

(* Create a non-Uniform i,j,k Hex Mesh *)
rx = pointsToMesh@unitMeshGrowth2Sided[40, 1/10]
ry = pointsToMesh[Subdivide[0, 1, 10]]
rz = pointsToMesh@unitMeshGrowth[40, 1/10]
RegionProduct[rx, rz];
RegionProduct[rx, ry];
rpcube = RegionProduct[rx, ry, rz];
HighlightMesh[rpcube, Style[1, Orange]]

RegionProduct Image

We have a nice non-uniform mesh that is refined at the floor and around the peak.

Determine $i,j,k\rightarrow x,y,z$ Mapping

I used the following code to help me determine the mapping.

(* Use RescalingTransform to Map from i,j,k\[Rule]x,y,z *)
bndm = {{0, 1/2}, {0, 1}, {0, 1}};
bndp = {{1/2, 1}, {0, 1}, {0, 1}};
bndms = {{xmin, xmid}, {ymin, ymax}, {zmin, zmax}};
bndps = {{xmid, xmax}, {ymin, ymax}, {zmin, zmax}};
RescalingTransform[bndm, bndms][{i, j, k}] /. {xmid -> xmidfn, 
  zmin -> zminfn}
RescalingTransform[bndp, bndps][{i, j, k}] /. {xmid -> xmidfn, 
  zmin -> zminfn}
(* Inspection of RescalingTransforms suggested following functions *)

ycj[j_, ymin_ , ymax_] := j (ymax - ymin) + ymin
xcij[i_ /; i <= 1/2, xmin_, xmax_, xmidfn_] := 
 2 i (xmidfn - xmin) + xmin
xcij[i_ /; i > 1/2, xmin_, xmax_, xmidfn_] := -xmax + 
  2 i (xmax - xmidfn) + 2 xmidfn
zcijk[k_, zmin_, zmax_, zminfn_] := k (zmax - zminfn) + zminfn
xmid[y_] := y/5
zminijk[x_, y_] := h[x, y]
trans[xmin_, xmax_, ymin_, ymax_, zmax_][{i_, j_, k_}] := {xcij[i, 
   xmin, xmax, 1/5 (j (ymax - ymin) + ymin)], 
  j (ymax - ymin) + 
   ymin, (k zmax + (1 - k) zminijk[
      xcij[i, xmin, xmax, xmid[ycj[j, ymin, ymax]]], 
      ycj[j, ymin, ymax]])}

Convert RegionProduct $i,j,k$ Mesh to $x,y,z$ ElementMesh

The following workflow show how to convert and visualize the RegionProduct mesh into an ElementMesh.

(* Create non-uniform mesh in x,y,z space *)
(*Use associations for clearer bc assignment later*)
bnd = <|"t" -> 1, "n" -> 2, "e" -> 3, "s" -> 4, "w" -> 5, "b" -> 0, 
   "default" -> 0|>;
eps = 0.000001;
ex = {1, 0, 0};
ey = {0, 1, 0};
ez = {0, 0, 1};
(* Convert RegionProduct mesh to x,y,z mesh *)
marker = 0;
crd = MeshCoordinates[rpcube];
(* transform coordinates *)
crd = trans[-rng, rng, -rng, rng, rng][#] & /@ crd;
(* grab hexa element incidents from RegionProduct mesh *)
inc = Delete[0] /@ MeshCells[rpcube, 3];
(* incidents have negative Jacobian and need re-ordering *)
inc = RotateLeft[#, 4] & /@ inc;
(* setup markers if desired *)
mrkrs = ConstantArray[marker, First@Dimensions@inc];
(* create element mesh *)
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {HexahedronElement[inc, mrkrs]}];
(* This step is for visualization *)
Short[bn = Flatten[mesh["BoundaryNormals"], 1]];
nFaces = First@Dimensions@bn;
faceMarkers = ConstantArray[bnd["b"], nFaces];
quads = Flatten[ ElementIncidents[mesh["BoundaryElements"]], 1];
posOpN = Position[{x_, y_, z_} /; Abs[x + 1] <= eps];
posOpE = Position[{x_, y_, z_} /; Abs[y - 1] <= eps];
posOpS = Position[{x_, y_, z_} /; Abs[x - 1] <= eps];
posOpW = Position[{x_, y_, z_} /; Abs[y + 1] <= eps];
posOpT = Position[{x_, y_, z_} /; Abs[z - 1] <= eps];
(* Reassign faceMarkers based on positions *)
faceMarkers[[Flatten@posOpN@bn]] = bnd["n"];
faceMarkers[[Flatten@posOpE@bn]] = bnd["e"];
faceMarkers[[Flatten@posOpS@bn]] = bnd["s"];
faceMarkers[[Flatten@posOpW@bn]] = bnd["w"];
faceMarkers[[Flatten@posOpT@bn]] = bnd["t"];
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {HexahedronElement[inc, mrkrs]}, 
   "BoundaryElements" -> {QuadElement[quads, faceMarkers]}];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]

XYZ Mesh

The mesh looks pretty good and it only took 30,420 hexahedron elements or about 30x fewer elements than stated in the OP.

Apply the Mesh to NDSolve

Now, apply the mesh to NDSolve using defaults.

sol = NDSolveValue[{Inactive[Laplacian][\[Psi][x, y, z], {x, y, z}] - 
      Load[x, y, z] == 0, 
    DirichletCondition[\[Psi][x, y, z] == 0, dz == 0], 
    DirichletCondition[\[Psi][x, y, z] == BcPsi[x, y, z], 
     z == rng]}, \[Psi], {x, y, z} \[Element] mesh];
 sol[x, y, z], {"YStackedPlanes", 3}, {x, y, z} \[Element] mesh, 
 PlotRange -> {All, All, {0, 5}}, PlotPoints -> 50, Contours -> 40, 
 AxesLabel -> Automatic, ColorFunction -> "BrightBands", 
 PlotLegends -> Automatic]


The solution does not look too bad. Of course, one needs to play with mesh parameters (e.g., boundary layer refinement in $r_z$ RegionProduct) to discover the "mesh insensitive" solution.

  • 1
    $\begingroup$ Thank-you, I am very grateful. This looks very promising indeed. I will study it carefully. $\endgroup$ Commented Apr 22, 2020 at 7:33
  • $\begingroup$ This is a very elegant solution and produces some very nice results. I have adapted the step size and ratio and the results now compare very favourably with the 2D version. Again, thank-you very much for inputting your time and expertise to help. $\endgroup$ Commented Apr 22, 2020 at 8:55
  • 1
    $\begingroup$ (+1) Excellent as always. $\endgroup$
    – user21
    Commented Apr 22, 2020 at 11:05
  • $\begingroup$ @PaulHarrison You are welcome! In my opinion, if you want to do FEM right, then you need precise control over your meshing (especially in 3D). That is why it may appear that I am harping over Element Mesh Tutorial. $\endgroup$
    – Tim Laska
    Commented Apr 22, 2020 at 16:31
  • 1
    $\begingroup$ @user21 Thank you as always. Initially, I entertained the idea of sweeping an "Acoustic Cloak" type mesh, but I abandoned it because I don't believe prism elements are supported. Going from 2D to 3D, it would be useful to support additional volumetric element types. Does the underlying mesher support additional types that have not been exposed yet or would it be a development effort? Just curious. $\endgroup$
    – Tim Laska
    Commented Apr 23, 2020 at 3:20

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.