3
$\begingroup$

I have the following code which works just fine:

 Needs["NDSolve`FEM`"]
 \[CapitalDelta]z = 0.05;(*[mm]*)(*Layer thickness*)
 R = 40;(*[mm]*)
 a = 2*R;
 b = 2*R*1.2;

 cuboid= Cuboid[{0, 0, 0}, {a, b, \[CapitalDelta]z}];
 RegionPlot3D[cuboid, PlotStyle -> Red, 
              PlotLabel -> Style["Schicht", Bold, 20], Boxed -> False, 
              ImageSize -> Large, BoxRatios -> Automatic]

 meshcuboid = ToElementMesh[cuboid];
 meshcuboid["Wireframe"]

However when I describe the region as an implicit region it doesn't work! why?

 layer = ImplicitRegion[-a/2 <= x <= a/2 && -b/2 <= y <= 
 b/2 && -\[CapitalDelta]z <= z <= 0, {x, y, z}];
 plotlayer = RegionPlot3D[layer, PlotStyle -> Red,
                    PlotLabel -> Style["Schicht", Bold, 20],
                    Boxed -> True, ImageSize -> Large, BoxRatios->Automatic]

 meshlayer = ToElementMesh[layer];
 meshlayer["Wireframe"]
$\endgroup$
  • $\begingroup$ sorry was probably a misunderstanding. I was actually adding more tags to the question+ $\endgroup$ – Gustavo Meyagan May 3 at 9:07
  • $\begingroup$ Ah, I see. No problem. $\endgroup$ – Henrik Schumacher May 3 at 9:08
  • $\begingroup$ @GustavoMeyagan I have created MeshTools package with utilities for more controlled mesh generation. Check CuboidMesh function. Or alternative path through 2D triangular mesh and then using TriangleToQuadMesh -> SmoothenMesh -> ExtrudeMesh. $\endgroup$ – Pinti May 6 at 7:48
  • $\begingroup$ @Pinti Thank you very much! I will check it out. $\endgroup$ – Gustavo Meyagan May 6 at 18:53
4
$\begingroup$

The problem is that your Cuboid is very thin. When you submit a Cuboid to ToElementMesh, it will by default create a hexahedral mesh: That pretty easy: Only one layer of hexahedral elements is created.

However, when you employ ImplicitRegion, by default a tetrahedral mesh is created. For discretizations of general volumetric domains, Mathematica employs the external library TetGen. And by default TetGen tries to subdivide and subdivide further and further until all tetrahedra have good aspect ratios. For your very thin layer, this means that it has to create bazillions of tetrahedra which must be why it takes so long and won't finish.

You can observe that by increasing \[CapitalDelta]z, so that TetGen is able to finish in "finite" time:

Needs["NDSolve`FEM`"]
\[CapitalDelta]z = 0.05 10;
R = 40;
a = 2*R;
b = 2*R*1.2;

meshcuboid = ToElementMesh[Cuboid[{0, 0, 0}, {a, b, \[CapitalDelta]z}]];
layer = ImplicitRegion[-a/2 <= x <= a/2 && -b/2 <= y <= 
     b/2 && -\[CapitalDelta]z <= z <= 0, {x, y, z}];

meshLoRes = ToElementMesh[layer, "PerformanceGoal" -> "Speed", "CheckQuality" -> False];
meshHiRes = ToElementMesh[layer];

Now, let's see how many mesh elements we have in our meshes:

Length@meshcuboid["MeshElements"][[1, 1]]
Length@meshLoRes["MeshElements"][[1, 1]]
Length@meshHiRes["MeshElements"][[1, 1]]

12669

99025

138999

You see, TetGen produces about ten times as many tets (for this particular thickness) and does so with a method that has considerable higher computational complexity than just subdividing a cuboid into small once. You can try to fiddle around a bit with the discretization options of ToElementMesh, but I don't expect that this will help much for your particular choice of $\Delta z$. =/

I am not sure what you are about to do with the meshes. If you are going to solve PDEs on them, then it is not guaranteed that meshcuboid will lead to good results, because its elements have probably too bad an aspect ratio for supporting a good numerical accuracy.

That is the reason why PDE on thin layers are abstracted to PDE on surfaces. By getting rid of the third dimension, one circumvents the need for overly fine discretizations. However, one has to be sure that the 3D PDE (or rather its solutions) "converges" in a sufficiently strong sense to a well-posed PDE on the limiting 2D surfaces in the limit of $\Delta z \to 0$. When the PDE stems from a minimization problem, such is usually proven in the language of $\Gamma$-convergence.

$\endgroup$
  • $\begingroup$ Thank you very much for your help! I actaully want to simulate heat flow for the given domain with prescribed boundary conditons on the surfaces of the region.Before I impose NDSolve however, I would like to control the meshprecission ,thus the reason why I wanted the visualisation of the domain beforehand. Unfortunately the layer of the domain is fixed (the context is 3D Printing). How do I proceed further from here? $\endgroup$ – Gustavo Meyagan May 3 at 10:24
  • $\begingroup$ Hm. I hope you use Neumann boundary conditions on the large planar surfaces of the of the thick plate. Then one could model the heat flow on the thick plate by a heat flow on the thin plate with interior heat source proportional to <<Neumann boundary condition on top surface>> - <<Neumann boundary condition on bottom surface>>. Do you understand what I mean? Maybe it would help me if you wrote down the exact PDE that you want to solve. $\endgroup$ – Henrik Schumacher May 3 at 10:50
  • $\begingroup$ @Hendrik Schumacher quick question: what options can I use to enforce a regular mesh size (i.e. all mesh elements should have the same size) And how can I change the form of a mesh (i.e from triangle to quadrilateral)? $\endgroup$ – Gustavo Meyagan May 3 at 18:19
  • $\begingroup$ The element type can be selected with the option "MeshElementType". Some (but not all) possible values are "QuadElement", "TriangleElement", "HexahedralElement", and "TetrahedralElement". Currently "QuadElement" and "HexahedralElement" won't be available for must input geometries since it is actually quite hard to generate such meshes with high quality. $\endgroup$ – Henrik Schumacher May 3 at 18:58
  • $\begingroup$ For the options that control the mesh quality and size, see the documentation page of DiscretizeRegion. It works essentially the same way. In particular the option MeshRefinementFunction is o interest. $\endgroup$ – Henrik Schumacher May 3 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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