You can subscribe to this list here.
2003 
_{Jan}
(4) 
_{Feb}
(1) 
_{Mar}
(9) 
_{Apr}
(2) 
_{May}
(7) 
_{Jun}
(1) 
_{Jul}
(1) 
_{Aug}
(4) 
_{Sep}
(12) 
_{Oct}
(8) 
_{Nov}
(3) 
_{Dec}
(4) 

2004 
_{Jan}
(1) 
_{Feb}
(21) 
_{Mar}
(31) 
_{Apr}
(10) 
_{May}
(12) 
_{Jun}
(15) 
_{Jul}
(4) 
_{Aug}
(6) 
_{Sep}
(5) 
_{Oct}
(11) 
_{Nov}
(43) 
_{Dec}
(13) 
2005 
_{Jan}
(25) 
_{Feb}
(12) 
_{Mar}
(49) 
_{Apr}
(19) 
_{May}
(104) 
_{Jun}
(60) 
_{Jul}
(10) 
_{Aug}
(42) 
_{Sep}
(15) 
_{Oct}
(12) 
_{Nov}
(6) 
_{Dec}
(4) 
2006 
_{Jan}
(1) 
_{Feb}
(6) 
_{Mar}
(31) 
_{Apr}
(17) 
_{May}
(5) 
_{Jun}
(95) 
_{Jul}
(38) 
_{Aug}
(44) 
_{Sep}
(6) 
_{Oct}
(8) 
_{Nov}
(21) 
_{Dec}

2007 
_{Jan}
(5) 
_{Feb}
(46) 
_{Mar}
(9) 
_{Apr}
(23) 
_{May}
(17) 
_{Jun}
(51) 
_{Jul}
(41) 
_{Aug}
(4) 
_{Sep}
(28) 
_{Oct}
(71) 
_{Nov}
(193) 
_{Dec}
(20) 
2008 
_{Jan}
(46) 
_{Feb}
(46) 
_{Mar}
(18) 
_{Apr}
(38) 
_{May}
(14) 
_{Jun}
(107) 
_{Jul}
(50) 
_{Aug}
(115) 
_{Sep}
(84) 
_{Oct}
(96) 
_{Nov}
(105) 
_{Dec}
(34) 
2009 
_{Jan}
(89) 
_{Feb}
(93) 
_{Mar}
(119) 
_{Apr}
(73) 
_{May}
(39) 
_{Jun}
(51) 
_{Jul}
(27) 
_{Aug}
(8) 
_{Sep}
(91) 
_{Oct}
(90) 
_{Nov}
(77) 
_{Dec}
(67) 
2010 
_{Jan}
(25) 
_{Feb}
(36) 
_{Mar}
(98) 
_{Apr}
(45) 
_{May}
(25) 
_{Jun}
(60) 
_{Jul}
(17) 
_{Aug}
(36) 
_{Sep}
(48) 
_{Oct}
(45) 
_{Nov}
(65) 
_{Dec}
(39) 
2011 
_{Jan}
(26) 
_{Feb}
(48) 
_{Mar}
(151) 
_{Apr}
(108) 
_{May}
(61) 
_{Jun}
(108) 
_{Jul}
(27) 
_{Aug}
(50) 
_{Sep}
(43) 
_{Oct}
(43) 
_{Nov}
(27) 
_{Dec}
(37) 
2012 
_{Jan}
(56) 
_{Feb}
(120) 
_{Mar}
(72) 
_{Apr}
(57) 
_{May}
(82) 
_{Jun}
(66) 
_{Jul}
(51) 
_{Aug}
(75) 
_{Sep}
(166) 
_{Oct}
(232) 
_{Nov}
(284) 
_{Dec}
(105) 
2013 
_{Jan}
(168) 
_{Feb}
(151) 
_{Mar}
(30) 
_{Apr}
(145) 
_{May}
(26) 
_{Jun}
(53) 
_{Jul}
(76) 
_{Aug}
(33) 
_{Sep}
(23) 
_{Oct}
(72) 
_{Nov}
(125) 
_{Dec}
(38) 
2014 
_{Jan}
(47) 
_{Feb}
(62) 
_{Mar}
(27) 
_{Apr}
(8) 
_{May}
(12) 
_{Jun}
(2) 
_{Jul}
(22) 
_{Aug}
(22) 
_{Sep}

_{Oct}
(17) 
_{Nov}
(20) 
_{Dec}
(12) 
2015 
_{Jan}
(25) 
_{Feb}
(2) 
_{Mar}
(16) 
_{Apr}
(13) 
_{May}
(21) 
_{Jun}
(5) 
_{Jul}
(1) 
_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 



1

2

3
(3) 
4
(6) 
5

6

7
(4) 
8
(12) 
9
(5) 
10

11
(1) 
12

13

14
(5) 
15

16
(2) 
17
(1) 
18

19

20

21

22
(1) 
23
(1) 
24

25
(1) 
26

27

28

29

30
(2) 
31



From: Roy Stogner <roystgnr@ic...>  20060830 18:43:37

On Wed, 30 Aug 2006, Tim Kr=F6ger wrote: > Please find attached a patch that makes all versions of inverse_map() u= se a=20 > default of TOLERANCE and accept an optional parameter for overriding th= at.=20 > Also, they all now accept the `secure' parameter. This patch wouldn't compile on my system (with infinite element support turned on); adding the additional parameters to the InfFE versions of those functions fixed it, so it's all committed now.  Roy 
From: <tim@ce...>  20060830 15:27:23

Dear libMesh developers, In FEInterface::inverse_map(), the default value for tolerance is 1.e4, while in FE::inverse_map(), it is the macro TOLERANCE, which is defined in libmesh_common.h and depends on the compiler. I find this somehow irritating. In particular, if I use the variant of FEInterface::inverse_map() that takes a vector of points (rather than a single point), the default tolerance of FE::inverse_map() (not of FEInterface::inverse_map()) applies, i.e. I get a different default tolerance depending on whether I use single points or a vector. Also, I think, the user should be able to specify a different than the default tolerance also for the vector version if inverse_map(). (Currently, this is only possible for the point version.) The reason why I found out this is that I face a situation where, in a KellyErrorEstimator, the inverse_map() fails to find the result up to the default tolerance (which is TOLERANCE, which is 1.e6 in my case), although the cells have affine maps. I think it's just because of the very small cells that I use. Please find attached a patch that makes all versions of inverse_map() use a default of TOLERANCE and accept an optional parameter for overriding that. Also, they all now accept the `secure' parameter. Best Regards, Tim 
From: <tim@ce...>  20060825 13:35:35

Dear all, Another patch about the TetGen interface. Without this, the code will not compile in debug mode. Sorry. Best Regards, Tim 
From: <tim@ce...>  20060823 13:44:50

Dear libMesh developers, PointLocatorTree::init() creates an OctTree (or a QuadTree) with level=100, which means that each TreeNode can hold up to 100 items before it refines. However, I just faced a situation in which 100 is not enough. Note that for tetrahedra, it is quite usual that up to 36 elements can meet at one node (I discussed that with a colleague of mine who is an expert in such geometric issues). In my situation, there are 66 tetrehedra meeting at one node (seemingly due to irregularities of the initial grid). These 66 elements will always belong to the same TreeNode, no matter how many times it is refined. Since TreeNode::insert() now decides by virtue of the element's bounding box, there are some more elements that will always be in the same TreeNode, although they do not contain this point. Hence, 100 can easily be reached. In my case, this led to a segmentation fault due to stack overflow. A number of 200 instead of 100 (see attached patch) solved the problem for me. It seems to be difficult so find a secure upper bound, but 100 is obviously not enough. Comments welcome. Best Regards, Tim 
From: <tim@ce...>  20060822 12:58:37

Dear libMesh developers, After submitting my patch about the TetGenWrapper class, I noticed that using region attributes is only useful if you can access the cell attributes later. For this reason, please find attached another patch that implements a function get_element_attribute(). Best Regards, Tim 
From: <tim@ce...>  20060817 15:16:07

Dear all, The TetGenWrapper class did until now not support region attributes. Since I needed them, I extended the class appropriately; see the attached patch. Could please someone commit this change to cvs? Best Regards, Tim  Dr. Tim Kroeger CeVis  Center of Complex Systems and Visualization University of Bremen Universitaetsallee 29 (Office 3.13) tim@... D28359 Bremen Phone +494212187710 Germany Fax +494212184236 
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 
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: 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...>  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: Kirk, Benjamin \(JSCEG\) <benjamin.kirk1@na...>  20060814 14:17:44

That should work, but it should be more efficient to check if the number of elements is 1 at the beginning of the partitioning and, if so, assign the element to processor 0 and return. Sorry, I clearly did not think enough about calling the partitioner with one element. Ben=20 Original Message From: libmeshdevelbounces@... [mailto:libmeshdevelbounces@...] On Behalf Of Roy Stogner Sent: Friday, August 11, 2006 4:44 PM To: libmeshdevel@... Subject: [Libmeshdevel] Bug in MetisPartitioner? When there's only one element in the mesh (as happens to me a lot: I build a single cube and then uniformly refine), the vector adjncy is empty, and adjncy[0] is then undefined. Adding: if (adjncy.empty()) adjncy.push_back(0); Seems to fix the problem, so I'll commit that to CVS, but I'd appreciate it if someone more familiar than me with Metis and our interface code would take a look at this.  Roy   Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with preintegrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.asus.falkag.net/sel?cmd=3Dlnk&kid=3D120709&bid=3D263057&dat=3D= 121642 _______________________________________________ Libmeshdevel mailing list Libmeshdevel@... https://lists.sourceforge.net/lists/listinfo/libmeshdevel 
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...>  20060811 21:44:15

When there's only one element in the mesh (as happens to me a lot: I build a single cube and then uniformly refine), the vector adjncy is empty, and adjncy[0] is then undefined. Adding: if (adjncy.empty()) adjncy.push_back(0); Seems to fix the problem, so I'll commit that to CVS, but I'd appreciate it if someone more familiar than me with Metis and our interface code would take a look at this.  Roy 
From: Karl Tomlinson <k.tomlinson@au...>  20060809 20:59:49

Roy Stogner writes: > On Wed, 9 Aug 2006, Karl Tomlinson wrote: > >>> Do you have any references to macroelementspecific quadrature rules >>> in the literature? I'd been thinking about deriving something more >>> efficient for the cubic CT, but I'd hate to reinvent the wheel. >> >> I was merely thinking that implementing the three Gaussian copies >> as a one set of weights and abscissae would keep loops over >> quadrature points simple. I haven't looked at libMesh's >> implementation. Maybe that is what you already do. > > No, we duplicate the weights and pretransform the point locations  > that only wastes a tiny bit of RAM, and it keeps the code simpler. Sounds sensible. (I think we might be trying to describe the same thing.) > > What I was thinking when talking about macroelementspecific rules is > that it's very unlikely that the duplicatedGaussian is the most > efficient quadrature rule on a macroelement. For example it seems > like you could shave a few points off of the rule by coming up with > something that would exactly integrate all of the pairs of > CloughTocher functions (and/or their derivatives) but wouldn't > necessarily give an exact answer on any product of piecewise cubic > polynomials (or their derivatives), since most of the piecewise cubics > aren't C1 and so aren't in the CT space. Sounds possible. I don't know if has been done before. Would be an interesting study. But if the Jacobian is not constant, or the integral is not just a product of shape functions/derivatives, it could mess things up a bit. I'd feel safer with something that at least handled the cubic(squared)order variations. If your goal is just to make the code faster, I can't help wondering if the same time could be spent making bigger gains elsewhere. 
From: Roy Stogner <roystgnr@ic...>  20060809 13:57:31

On Wed, 9 Aug 2006, Tim Kr=F6ger wrote: > Dear Roy, > > On Tue, 8 Aug 2006, Roy Stogner wrote: > >> Is it possible that the crash was when inverse_map() was being called >> from some other library function? > > It must have been like this although I have no idea from where inverse_= map()=20 > would be called with a point far away from the cell. If the problem comes up again and you can reproduce it, let us know. >>> Okay, so I think the best thing to do is: Whenever the=20 >>> outofmeshextension of MeshFunction is enabled, it will check=20 >>> Elem::has_affine_map() for all elements (at least in debug mode), and= if=20 >>> any element returns true, an error will be generated. > > I impmeneted this now; please find the patch attached. Looks good; I've committed to CVS. > Note that I had to include the new functions also in PointLocatorBase, = so=20 > that they also had to be included in PointLocaterList, where they now j= ust=20 > call error(). I roughly looked at the code of PointLocatorList and fou= nd it=20 > not really reliable anyway. I wonder whether anyone will be using it. My guess is that Daniel wrote the simpler List implementation first, then everybody ignored it once Tree became the default. The Tree should be more efficient for most meshes; we ought to just mark List's constructor deprecated() before releasing 0.6.0, and get rid of the whole class later unless somebody wants to go through it and bugfix.  Roy 
From: Roy Stogner <roystgnr@ic...>  20060809 12:50:58

On Wed, 9 Aug 2006, Karl Tomlinson wrote: >> Do you have any references to macroelementspecific quadrature rules >> in the literature? I'd been thinking about deriving something more >> efficient for the cubic CT, but I'd hate to reinvent the wheel. > > I was merely thinking that implementing the three Gaussian copies > as a one set of weights and abscissae would keep loops over > quadrature points simple. I haven't looked at libMesh's > implementation. Maybe that is what you already do. No, we duplicate the weights and pretransform the point locations  that only wastes a tiny bit of RAM, and it keeps the code simpler. What I was thinking when talking about macroelementspecific rules is that it's very unlikely that the duplicatedGaussian is the most efficient quadrature rule on a macroelement. For example it seems like you could shave a few points off of the rule by coming up with something that would exactly integrate all of the pairs of CloughTocher functions (and/or their derivatives) but wouldn't necessarily give an exact answer on any product of piecewise cubic polynomials (or their derivatives), since most of the piecewise cubics aren't C1 and so aren't in the CT space.  Roy 
From: Karl Tomlinson <k.tomlinson@au...>  20060809 09:46:03

Thanks for the explanations, Roy. I now have a better idea of some of the issues involved. I'll have to have a deeper look into solution and geometry implementations to work out how they can be made more similar and provide extra dofs for the geometry (while still keeping [or making] the API abstract enough to allow other representations in the future). Although, maybe I'll spend a bit more time getting used to libMesh as is, before diving into changing core functionality. Roy Stogner writes: > On Tue, 8 Aug 2006, Karl Tomlinson wrote: > >>> ... I'd be worried about what CloughTocher mappings >>> might do to your quadrature rules. Granted, any kind of nonaffine >>> map can mess up your nice exact Gaussian quadrature, but the >>> CloughTocher basis functions have internal subelement boundaries >>> which might be even worse. >> >> The quadrature rules should be applied over each of the >> subelements individually (or through a macroquadrature rule). >> Isn't this also an issue when integrating solution variables or >> their shape functions? > > Yes; how we do things now is to just copy a Gaussian rule over each > subelement. I guess it's not a serious problem, but tripling the cost > of FEM calculations isn't something to be happy about. That is what I thought was appropriate, with the weights and abscissae the same in each subelement (apart from a simple transformation from macroelement to subelement xispace). > > Do you have any references to macroelementspecific quadrature rules > in the literature? I'd been thinking about deriving something more > efficient for the cubic CT, but I'd hate to reinvent the wheel. I was merely thinking that implementing the three Gaussian copies as a one set of weights and abscissae would keep loops over quadrature points simple. I haven't looked at libMesh's implementation. Maybe that is what you already do. 
From: <tim@ce...>  20060809 09:00:14

Dear Roy, On Tue, 8 Aug 2006, Roy Stogner wrote: > On Tue, 8 Aug 2006, Tim Kr=F6ger wrote: > >>> That's exactly why the inverse_map function has the "bool secure" >>> argument. If secure is false (as it is when contains_point() calls >>> inverse_map() ), then Newton isn't expected to converge and >>> inverse_map just exits after 10 iterations. >>=20 >> Strange enough, this didn't happen last week when I had changed the code= =20 >> such that only active elements were inserted into the tree but they were= =20 >> still handled only by their centroid: The code crashed because Newton=20 >> didn't converge within 20 iterations (with warnings from the 11th iterat= ion=20 >> on). I have no explanation for that because when I look at the code, it= =20 >> seems to me that you are right. > > Is it possible that the crash was when inverse_map() was being called > from some other library function? It must have been like this although I have no idea from where=20 inverse_map() would be called with a point far away from the cell. >> Okay, so I think the best thing to do is: Whenever the=20 >> outofmeshextension of MeshFunction is enabled, it will check=20 >> Elem::has_affine_map() for all elements (at least in debug mode), and if= =20 >> any element returns true, an error will be generated. I impmeneted this now; please find the patch attached. Note that I had to include the new functions also in PointLocatorBase,=20 so that they also had to be included in PointLocaterList, where they=20 now just call error(). I roughly looked at the code of=20 PointLocatorList and found it not really reliable anyway. I wonder=20 whether anyone will be using it. Best Regards, Tim 
From: Steffen Petersen <steffen.petersen@tu...>  20060808 19:23:15

No, that is not the problem. I had noticed that change a while ago. The last time I had access to the repository was probably about a month ago. Has anything changed since? Steffen > > The CVS/Root recently changed to > > user@...:/cvsroot/libmesh > > instead of > > user@...:/cvsroot/libmesh > > Maybe that's it? > J > > >  >  > Using Tomcat but need to do more? Need to support web > services, security? > Get stuff done quickly with preintegrated technology to make > your job easier > Download IBM WebSphere Application Server v.1.0.1 based on > Apache Geronimo > http://sel.asus.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&; dat=121642 > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel > 
From: John Peterson <peterson@cf...>  20060808 18:26:15

The CVS/Root recently changed to user@...:/cvsroot/libmesh instead of user@...:/cvsroot/libmesh Maybe that's it? J 
From: Steffen Petersen <steffen.petersen@tu...>  20060808 18:01:05

I have some problems with the libMesh cvs. I currently cannot acces the actual libMesh version with my login (no updates, no commits, no diffs, ...). Perhaps I was just missing some sorceforge news. Has anyone else noticed similar problems? Steffen 
From: Derek Gaston <friedmud@gm...>  20060808 17:45:48

Roy, Thanks for the summary of the discussion! I will do a cvsupdate and try to use the new stuff... the only problem being that I haven't updated since the beginning of the summer.... so I might have a bunch of stuff to manually merge... sigh (it's my own damn fault!) Also... the chunk of code I wrote that uses MeshFuncion... has simply _not_ been working lately. I don't know what the problem is, but it's just not doing what it should be... so I'm not sure I can even reliably test MeshFunction right now.... fixing that code is on my todo list though (it's my projection code... it just doesn't seem to be doing the projection anymore....). Derek On 8/8/06, Roy Stogner <roystgnr@...> wrote: > On Tue, 8 Aug 2006, Derek Gaston wrote: > > > Firstly about the newton iterations for inverse_map(). I ran into > > this same bit of code a while ago... and through discussions with Roy > > and John decided the best thing to do was to change the code so that > > it returns a point that is definitely wrong when the newton iterations > > aren't converging. Therefore I have something like this in my code: > > > > if (cnt > 10) > > { > > //Return a point that is definitely wrong > > for(unsigned int i=0;i<Dim;i++) > > p(i)=1e6; > > return p; > > } > > > > Now whether or not that's right or good... it does seem to work. I > > just wanted to report what I was doing. > > You reported what you were doing weeks ago; but either you never sent > me a patch or I never got it committed! > > Anyway, I'm changing fe_map.C now to put that behavior into CVS; > double check it and make sure it works the way your current version > does. > > > I am definitely up for trying out new point_locator stuff! Should I > > try to apply the patches from 4 days ago give feedback? > > They should all be in CVS now; just doing a "cvs update" should get > everything ready for you to test. > > > Since I am late to the discussion I'm still trying to figure out > > what the latest state of the patches is. > > The PointLocatorTree (the default PointLocator implementation that > MeshFunction uses) should be significantly more efficient now and > should work correctly on refined meshes; that's all in CVS. > > There will be an new MeshFunction API (not yet in CVS) that will > return a userdefined constant value when called on points outside the > mesh. > > There's still a chance of the PointLocatorTree failing (and falling > back on a slow linear search) if there are mesh elements with curved > sides. Nobody knows how to fix that efficiently; we'll probably put > in some hack to try and minimize the problem. >  > Roy > 
From: Roy Stogner <roystgnr@ic...>  20060808 17:38:07

On Tue, 8 Aug 2006, Derek Gaston wrote: > Firstly about the newton iterations for inverse_map(). I ran into > this same bit of code a while ago... and through discussions with Roy > and John decided the best thing to do was to change the code so that > it returns a point that is definitely wrong when the newton iterations > aren't converging. Therefore I have something like this in my code: > > if (cnt > 10) > { > //Return a point that is definitely wrong > for(unsigned int i=0;i<Dim;i++) > p(i)=1e6; > return p; > } > > Now whether or not that's right or good... it does seem to work. I > just wanted to report what I was doing. You reported what you were doing weeks ago; but either you never sent me a patch or I never got it committed! Anyway, I'm changing fe_map.C now to put that behavior into CVS; double check it and make sure it works the way your current version does. > I am definitely up for trying out new point_locator stuff! Should I > try to apply the patches from 4 days ago give feedback? They should all be in CVS now; just doing a "cvs update" should get everything ready for you to test. > Since I am late to the discussion I'm still trying to figure out > what the latest state of the patches is. The PointLocatorTree (the default PointLocator implementation that MeshFunction uses) should be significantly more efficient now and should work correctly on refined meshes; that's all in CVS. There will be an new MeshFunction API (not yet in CVS) that will return a userdefined constant value when called on points outside the mesh. There's still a chance of the PointLocatorTree failing (and falling back on a slow linear search) if there are mesh elements with curved sides. Nobody knows how to fix that efficiently; we'll probably put in some hack to try and minimize the problem.  Roy 
From: Derek Gaston <friedmud@gm...>  20060808 17:15:54

Wow, I apologize for being out of town all last week and missing this whole discussion! Firstly about the newton iterations for inverse_map(). I ran into this same bit of code a while ago... and through discussions with Roy and John decided the best thing to do was to change the code so that it returns a point that is definitely wrong when the newton iterations aren't converging. Therefore I have something like this in my code: if (cnt > 10) { //Return a point that is definitely wrong for(unsigned int i=3D0;i<Dim;i++) p(i)=3D1e6; return p; } Now whether or not that's right or good... it does seem to work. I just wanted to report what I was doing. I am definitely up for trying out new point_locator stuff! Should I try to apply the patches from 4 days ago give feedback? Since I am late to the discussion I'm still trying to figure out what the latest state of the patches is. Derek On 8/8/06, Roy Stogner <roystgnr@...> wrote: > On Tue, 8 Aug 2006, Tim Kr=F6ger wrote: > > >> That's exactly why the inverse_map function has the "bool secure" > >> argument. If secure is false (as it is when contains_point() calls > >> inverse_map() ), then Newton isn't expected to converge and > >> inverse_map just exits after 10 iterations. > > > > Strange enough, this didn't happen last week when I had changed the cod= e such > > that only active elements were inserted into the tree but they were sti= ll > > handled only by their centroid: The code crashed because Newton didn't > > converge within 20 iterations (with warnings from the 11th iteration on= ). I > > have no explanation for that because when I look at the code, it seems = to me > > that you are right. > > Is it possible that the crash was when inverse_map() was being called > from some other library function? We turn secure on in several places > where we expect Newton to converge, just to doublecheck us in case > one of those expectations is false. > > >> No good ones. All I can think to do is warn about the problem. When > >> you add an interface for outofmesh MeshFunction extensions, let that > >> interface turn off (well, replace) the linear search, and add an > >> assert() or #ifdef DEBUG block in that code path which verifies that > >> it doesn't see any quadratic mappings. > > > > Okay, so I think the best thing to do is: Whenever the outofmeshexte= nsion > > of MeshFunction is enabled, it will check Elem::has_affine_map() for al= l > > elements (at least in debug mode), and if any element returns true, an = error > > will be generated. But, however, I wonder whether Elem::has_affine_map= () > > could occasionally return false due to rounding errors. What do you th= ink? > > Should the equality checks be replaced by some fabs()<eps stuff? > > As John pointed out there is an epsilon in the the point operator=3D=3D > that handles that for us. > > That actually may cause some problems elsewhere in the library, now > that I think about it. If we perturb mesh nodes by something like > TOLERANCE/10, that should be enough to make a noticeable difference in > a sufficiently fine finite element solution, but it won't be enough to > convince the FE class that the elements are all distinct  we'd end up > erroneously caching shape function & derivative calculations from > element to element, as well as erroneously caching mapping derivative > calculations from quadrature point to quadrature point. >  > Roy > >  > Using Tomcat but need to do more? Need to support web services, security? > Get stuff done quickly with preintegrated technology to make your job ea= sier > Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronim= o > http://sel.asus.falkag.net/sel?cmd=3Dlnk&kid=3D120709&bid=3D263057&dat= =3D121642 > > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel > > > 