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@GrayLevel@.5, 
    FaceForm@GrayLevel@.7], 
  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$
  • 1
    $\begingroup$ Please provide the code that created em. $\endgroup$ – bbgodfrey Mar 7 '16 at 0:44
  • $\begingroup$ Does it change when you make a simpler mesh (less elements, thinner shell, etc...)? $\endgroup$ – user21 Mar 7 '16 at 0:49
  • $\begingroup$ Min[em["Quality"]] seems fine. $\endgroup$ – user21 Mar 7 '16 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$ – user21 Mar 7 '16 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$ – unlikely Mar 7 '16 at 9:12

Your Answer

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

Browse other questions tagged or ask your own question.