From: Roy Stogner <roystgnr@ic...>  20060814 13:33:38

On Mon, 14 Aug 2006, Tim Kr=F6ger wrote: > Although obtuse angles still appear, this looks much better. > > I would like to integrate this into libMesh. One problem with this is > that it will change the behaviour of existing code. Perhaps there > should be some global option that enables the user to choose between > the two refinement procedures, with the default being the current > behaviour. On the other hand, users might forget to choose the new > behaviour in new programs although they should. > > Another point is the question how to implement this. I currently > implmented this in Tet4::embedding_matrix(). However, during the > refinement procedure, this function is called quite a lot of times so > that the decision which of the octahedron's diagonals is shortest is > computed much too often. An alternative would be some static member > variables of Tet4 that allow to remember the latest computation. Of > course, the same funcionality should also be implemented for Tet10 > elements. > > Any comments on this? It's good to worry about changing existing code and slowing the tet refinement process  but you're also fixing a misbehavior (not technically a bug, but just as bad) so heinous that I can't imagine anybody complaining. Adding a static Elem* seems simple enough  if the Elem hasn't changed, you don't need to recompute. That could break in the billiontoone chance that one tet gets deleted and another created in its place in memory, with no embedding_matrix calls on any other elements in between... It may be cleaner to add a virtual init_embedding_matrix() function; Elem::refine() and Elem::compute_children_node_keys() could call that outside their c/nc/n loops, and so putting the test there would get it called once or twice per element instead of hundreds of times.  Roy 
From: <tim@ce...>  20060814 12:12:16

Dear all, When refining a tetrahedron cell, libMesh does the following: At each vertex of the parent cell, a child of the same shape and half side length is created, and the remaining octahedron is cut into four tetraheda along one of its three diagonals. However, the four tetrahedra arising from the octahedron do in general not have the same shape as the parent cell. This is not surprising since in particular, it is impossible to fill space with regular tetrahedra only. In the current implementation, the octahedron is just cut along an arbitrary diagonal. This leads to very obtuse angled cells. I wrote a test program that starts with a single, regular tetrahedron, uniformly refines it seven times, and determines the minimal and maximal angles in the grid: #refine min max 0 70.5288 70.5288 1 54.7356 109.471 2 35.2644 125.264 3 25.2394 144.736 4 13.2627 154.761 5 9.44623 166.737 6 5.21321 170.554 7 3.67805 174.787 Appareantly, the max angle approaches 180 (degrees) and the min angle 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 be a drawback. In my opinion, a solution would be to exploit the option to choose one of the octahedron's three diagonals: If we always choose the shortest one, obtuse angles remain comparably limited. I currently use a dirty hack to do this and get the following result: 0 70.5288 70.5288 1 54.7356 109.471 2 54.7356 109.471 3 54.7356 109.471 4 54.7356 109.471 5 54.7356 109.471 6 54.7356 109.471 7 54.7356 109.471 Although obtuse angles still appear, this looks much better. I would like to integrate this into libMesh. One problem with this is that it will change the behaviour of existing code. Perhaps there should be some global option that enables the user to choose between the two refinement procedures, with the default being the current behaviour. On the other hand, users might forget to choose the new behaviour in new programs although they should. Another point is the question how to implement this. I currently implmented this in Tet4::embedding_matrix(). However, during the refinement procedure, this function is called quite a lot of times so that the decision which of the octahedron's diagonals is shortest is computed much too often. An alternative would be some static member variables of Tet4 that allow to remember the latest computation. Of course, the same funcionality should also be implemented for Tet10 elements. Any comments on this? Best Regards, Tim 
From: Roy Stogner <roystgnr@ic...>  20060814 13:33:38

On Mon, 14 Aug 2006, Tim Kr=F6ger wrote: > Although obtuse angles still appear, this looks much better. > > I would like to integrate this into libMesh. One problem with this is > that it will change the behaviour of existing code. Perhaps there > should be some global option that enables the user to choose between > the two refinement procedures, with the default being the current > behaviour. On the other hand, users might forget to choose the new > behaviour in new programs although they should. > > Another point is the question how to implement this. I currently > implmented this in Tet4::embedding_matrix(). However, during the > refinement procedure, this function is called quite a lot of times so > that the decision which of the octahedron's diagonals is shortest is > computed much too often. An alternative would be some static member > variables of Tet4 that allow to remember the latest computation. Of > course, the same funcionality should also be implemented for Tet10 > elements. > > Any comments on this? It's good to worry about changing existing code and slowing the tet refinement process  but you're also fixing a misbehavior (not technically a bug, but just as bad) so heinous that I can't imagine anybody complaining. Adding a static Elem* seems simple enough  if the Elem hasn't changed, you don't need to recompute. That could break in the billiontoone chance that one tet gets deleted and another created in its place in memory, with no embedding_matrix calls on any other elements in between... It may be cleaner to add a virtual init_embedding_matrix() function; Elem::refine() and Elem::compute_children_node_keys() could call that outside their c/nc/n loops, and so putting the test there would get it called once or twice per element instead of hundreds of times.  Roy 
From: John Peterson <peterson@cf...>  20060814 14:40:51

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. 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 diagonalselector function?=20 Could the diagonalselector function be optimized at all? > Any comments on this? >=20 > Best Regards, >=20 > Tim >=20 John 
From: Roy Stogner <roystgnr@ic...>  20060814 15:00:04

On Mon, 14 Aug 2006, John Peterson wrote: > Does a read/write static variable raise any question of thread safety? Yes. Fortunately, we already use enough such variables that libMesh is not threadsafe, so adding one more won't kill us. And actually, I'm not sure how much of an improvement multithreading would get us. A thread is just a process that shares a lot of memory, and we're already multiprocess. Right now, of course, it would be very nice to have only one Mesh in memory per motherboard rather than one Mesh per CPU core  but after we distribute the Mesh class, that won't be a serious problem.  Roy 
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 
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 