[Libmesh-devel] Tetrehedron refinement

[John P. <pet...@cf...>]

Tim Kr=F6ger writes:
 > Dear all,
 >=20
 > When refining a tetrahedron cell, libMesh does the following:  At ea=
ch=20
 > vertex of the parent cell, a child of the same shape and half side=20=

 > length is created, and the remaining octahedron is cut into four=20
 > tetraheda along one of its three diagonals.  However, the four=20
 > tetrahedra arising from the octahedron do in general not have the sa=
me=20
 > shape as the parent cell.  This is not surprising since in particula=
r,=20
 > it is impossible to fill space with regular tetrahedra only.
 >=20
 > In the current implementation, the octahedron is just cut along an=20=

 > arbitrary diagonal.  This leads to very obtuse angled cells.  I wrot=
e=20
 > a test program that starts with a single, regular tetrahedron,=20
 > uniformly refines it seven times, and determines the minimal and=20
 > maximal angles in the grid:
 >=20
 > #refine=09min=09max
 > 0=0970.5288=0970.5288
 > 1=0954.7356=09109.471
 > 2=0935.2644=09125.264
 > 3=0925.2394=09144.736
 > 4=0913.2627=09154.761
 > 5=099.44623=09166.737
 > 6=095.21321=09170.554
 > 7=093.67805=09174.787
 >=20
 > Appareantly, the max angle approaches 180 (degrees) and the min angl=
e
 > approaches 0.  Since it is known that obtuse angles are likely to
 > cause the stiffness matrix to be no M-matrix any more (and hence the=

 > discrete maximum principle might not hold any more), this seems to b=
e
 > a drawback.

To say the least!!

 > In my opinion, a solution would be to exploit the option to choose o=
ne
 > of the octahedron's three diagonals: If we always choose the shortes=
t
 > one, obtuse angles remain comparably limited.  I currently use a dir=
ty
 > hack to do this and get the following result:
 >=20
 > 0=0970.5288=0970.5288
 > 1=0954.7356=09109.471
 > 2=0954.7356=09109.471
 > 3=0954.7356=09109.471
 > 4=0954.7356=09109.471
 > 5=0954.7356=09109.471
 > 6=0954.7356=09109.471
 > 7=0954.7356=09109.471
 >=20
 > Although obtuse angles still appear, this looks much better.
 >=20
 > I would like to integrate this into libMesh.  One problem with this =
is
 > that it will change the behaviour of existing code.  Perhaps there=20=

 > should be some global option that enables the user to choose between=
=20
 > the two refinement procedures, with the default being the current=20=

 > behaviour.  On the other hand, users might forget to choose the new=20=

 > behaviour in new programs although they should.

You're right, the new refinement pattern would produce different meshes=

and mess up regression tests etc.  Luckily we are unencumbered by those=
 :)
This might actually explain a problem I've recently faced in an applica=
tion
code (without knowing what it was).  Good find.

 > Another point is the question how to implement this.  I currently=20=

 > implmented this in Tet4::embedding_matrix().  However, during the=20=

 > refinement procedure, this function is called quite a lot of times s=
o=20
 > that the decision which of the octahedron's diagonals is shortest is=
=20
 > computed much too often.  An alternative would be some static member=
=20
 > variables of Tet4 that allow to remember the latest computation.  Of=
=20
 > course, the same funcionality should also be implemented for Tet10=20=

 > elements.

Does a read/write static variable raise any question of thread safety?
(I admit to not being knowledgeable in this area.)  How much slower is
the overall code when calling the diagonal-selector function?=20
Could the diagonal-selector function be optimized at all?


 > Any comments on this?
 >=20
 > Best Regards,
 >=20
 > Tim
 >=20



-John