2
$\begingroup$

The appearance of the following ElementMesh is bad as you can see for example on the inside surface of the shell. The same problem happens converting to a MeshRegion (conversion does not warn about problems).

Needs["NDSolve`FEM`"] em = Import[ "https://www.dropbox.com/s/wc9nd0kyf55t9bg/ElementMeshRenderingIssue.mx?dl=1"] em["Wireframe"["MeshElement" -> "MeshElements", "ElementMeshDirective" -> Directive[Specularity[.2], EdgeForm@[email protected], FaceForm@[email protected]], Lighting -> {{"Ambient", GrayLevel[0.45]}, {"Directional", GrayLevel[0.3], ImageScaled[{2, 0, 2}]}, {"Directional", GrayLevel[0.33], ImageScaled[{2, 2, 2}]}, {"Directional", GrayLevel[0.3], ImageScaled[{0, 2, 2}]}}]] 

Mathematica graphics

I don't understand if this is just a rendering problem or there is also a problem in the underlying mesh.

I have built this mesh by myself. As I stated in another question the vertices of the QuadElement are not exactly coplanar, and I don't know how to adjust the coordinates so that they eventually become coplanar.

I noticed some strange behavior using this mesh to solve a PDE with FEM: this can be another clue of some problem with the mesh.

To summarize I want to try to understand

  1. if the visualization problem is caused by some problem with the mesh construction
  2. if the problem is related with the non-exactly-planar faces, if there is some way to adjust

Any idea on how to proceed? Thanks


After more investigation, and thanks to @user21, I discovered that the problem of the failed ToBoundaryMesh/ToElementMesh roundtrip is probably related to the fact that the mesh uploaded and shown here is no a "complete" spherical shell (this is what I want for a figure to illustrate the process of mesh creation).

The mesh I use as a domain for my PDE is a complete spherical shell, so this should not be the reason why I get strange result when solving over this mesh.

I made another simpler mesh with the same algorith, starting with just 1 face of the cubed sphere (i.e. to get 1/6 of the spherical shell), and with a coarse grid.

mesh = Import[ "https://www.dropbox.com/s/46e93mbacvgn2i6/\ ElementMeshRenderingIssue2.mx?dl=1"] 

I also made a function to extract only part of this mesh.

ElementMeshExtract[mesh_, mkform : {memkform_, bemkform_, pemkform_}, idform : {meidform_, beidform_, peidform_}] := Module[{elems}, elems = MapThread[ Cases[#1, type_[el_, mk_] :> With[{sel = MatchQ[#2] /@ mk}, If[! Or @@ sel, Nothing, type[Pick[el, sel], Pick[mk, sel]]]] , {1}] &, {mesh /@ {"MeshElements", "BoundaryElements", "PointElements"}, mkform}]; elems = MapThread[ Cases[#1, type_[el_, mk_] :> With[{sel = MatchQ[#2] /@ Range@Length@el}, If[! Or @@ sel, Nothing, type[Pick[el, sel], Pick[mk, sel]]]] , {1}] &, {elems, idform}]; incidents = Union@Flatten@DeleteCases[elems, {}][[All, 1, 1]]; assoc = AssociationThread[incidents, Range@Length@incidents]; elems[[All, All, 1]] = Map[Lookup[assoc, #] &, elems[[All, All, 1]], {-2}]; ToElementMesh[ "Coordinates" -> mesh["Coordinates"][[Keys@assoc]], "MeshElements" -> elems[[1]], "BoundaryElements" -> elems[[2]], "PointElements" -> elems[[3]] ] ] ElementMeshExtract[mesh_, mkform : {_, _, _}, idform_] := ElementMeshExtract[mesh, mkform, {idform, idform, idform}] ElementMeshExtract[mesh_, mkform_, idform : {_, _, _}] := ElementMeshExtract[mesh, {mkform, mkform, mkform}, idform] ElementMeshExtract[mesh_, mkform_, idform_] := ElementMeshExtract[ mesh, {mkform, mkform, mkform}, {idform, idform, idform}] ElementMeshExtract[mesh_, mkform_] := ElementMeshExtract[mesh, mkform, _] 

Putting all together in a Manipulate I still cannot find anything wrong with that mesh that could explain the bad rendering or the strange results od NDSolve.

Manipulate[( lmesh = ElementMeshExtract[ mesh, {Alternatives @@ Range@layers, Except@_, Except@_}]; lc = Total[Length@*First /@ lmesh["MeshElements"]]; pmesh = ElementMeshExtract[lmesh, _, {i_ /; 1 <= i <= progress, _, _}]; Column[{ Show[ MeshRegion[pmesh, Method -> {"CoplanarityTolerance" -> -100}], ImageSize -> Medium, PlotRange -> lmesh["Bounds"], Method -> {"ShrinkWrap" -> False}], BoxWhiskerChart[ElementMeshQuality[pmesh], BarOrigin -> Left, BarSpacing -> $MachineEpsilon, PlotLabel -> "Quality", AspectRatio -> 1/4, PlotRange -> {{0.5, 1}, Automatic}, Frame -> {{False, False}, {True, True}}, GridLines -> {Automatic, None}, ImageSize -> Small] }, Center] ), {{lmesh, lmesh}, None}, {{pmesh, pmesh}, None}, {{lc, lc}, None}, {{layers, 2}, 1, 5, 1, Appearance -> "Labeled"}, {{progress, Floor[lc/2]}, 1, lc, 1, Appearance -> "Labeled"}, ControlPlacement -> Bottom ] 

Mathematica graphics

$\endgroup$
16
  • 1
    $\begingroup$ Please provide the code that created em. $\endgroup$ Commented Mar 7, 2016 at 0:44
  • $\begingroup$ Does it change when you make a simpler mesh (less elements, thinner shell, etc...)? $\endgroup$ Commented Mar 7, 2016 at 0:49
  • $\begingroup$ Min[em["Quality"]] seems fine. $\endgroup$ Commented Mar 7, 2016 at 0:53
  • $\begingroup$ This reveals that there is a problem in the boundary of the mesh: bmesh = ToBoundaryMesh[em]; ToElementMesh[bmesh] as it fails to generate a tetrahedralized mesh. $\endgroup$ Commented Mar 7, 2016 at 1:15
  • $\begingroup$ @bbgodfrey In the linked question I provided a simplified version of the code I use to create the spferical surface. The code I use to create the spherical shells is very long and complex because of the "wedge" elements and of many edge cases. I'm not sure if I can extract and post a simplified enough version. $\endgroup$ Commented Mar 7, 2016 at 9:12

0

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.