From: John Peterson <peterson@cf...>  20060816 05:12:26

Uhh....damn. I was using the wrong formula for computing the dihedral angles! I get the same thing as you now. Sorry, John John Peterson writes: > Hi, >=20 > 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 >=20 > http://www.cs.berkeley.edu/~jrs/papers/elemj.pdf >=20 > 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: >=20 > 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 >=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. >=20 > So the interpolation properties won't be too bad but the stiffness > matrix condition number would still be going to hell. Therefore I a= m > still intrigued by the diagonal swap trick which you implemented, > since it appears to also bound the minimum angle. >=20 > Looking forward to hearing your thoughts, > John >=20 >=20 >=20 >=20 >=20 > Tim Kr=F6ger writes: > > Dear all, > >=20 > > When refining a tetrahedron cell, libMesh does the following: At= each=20 > > vertex of the parent cell, a child of the same shape and half sid= e=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= same=20 > > shape as the parent cell. This is not surprising since in partic= ular,=20 > > it is impossible to fill space with regular tetrahedra only. > >=20 > > In the current implementation, the octahedron is just cut along a= n=20 > > arbitrary diagonal. This leads to very obtuse angled cells. I w= rote=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 a= ngle > > 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 t= o be > > a drawback. > >=20 > > In my opinion, a solution would be to exploit the option to choos= e one > > of the octahedron's three diagonals: If we always choose the shor= test > > one, obtuse angles remain comparably limited. I currently use a = dirty > > 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 th= is is > > that it will change the behaviour of existing code. Perhaps ther= e=20 > > should be some global option that enables the user to choose betw= een=20 > > the two refinement procedures, with the default being the current= =20 > > behaviour. On the other hand, users might forget to choose the n= ew=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 time= s so=20 > > that the decision which of the octahedron's diagonals is shortest= is=20 > > computed much too often. An alternative would be some static mem= ber=20 > > variables of Tet4 that allow to remember the latest computation. = Of=20 > > course, the same funcionality should also be implemented for Tet1= 0=20 > > elements. > >=20 > > Any comments on this? > >=20 > > Best Regards, > >=20 > > Tim >=20 