From: John Peterson <peterson@cf...>  20060816 01:58:19

Hi, I've had a chance to look into this issue in more detail. I'm wondering precisely which angles you are measuring...? According to JR Shewchuk's paper http://www.cs.berkeley.edu/~jrs/papers/elemj.pdf the relevant angles which control the stiffness matrix conditioning and the interpolation properties of the grid are the *dihedral* angles, i.e. the angles formed between any two planar faces of the tetrahedron. If you compute the global min/max dihedral angles in the test you mentioned, you should get the following: 0 70.5288 70.5288 (1 equilateral tet) 1 54.7356 90. 2 35.2644 90. 3 25.2394 90. 4 13.2627 90. 5 9.44623 90. 6 5.21321 90. 7 3.67805 90.=20 In other words, the min dihedral angle is going to zero but the max dihedral angle remains bounded by pi/2. I think this makes sense since the octahedron split will generate 4 tets meeting at a corner, thus the 90deg angle. So the interpolation properties won't be too bad but the stiffness matrix condition number would still be going to hell. Therefore I am still intrigued by the diagonal swap trick which you implemented, since it appears to also bound the minimum angle. Looking forward to hearing your thoughts, John 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 Mmatrix any more (and hence the= > discrete maximum principle might not hold any more), this seems to b= e > a drawback. >=20 > 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. >=20 > 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. >=20 > Any comments on this? >=20 > Best Regards, >=20 > Tim 