Menu
▾
▴
libmesh-devel
|
From: Kirk, B. \(JSC-EG\) <ben...@na...> - 2006-06-21 17:18:16
|
Which type of elements are these? For Quads the zero point corresponds to the centroid of the element, which should be as good a guess as any. For triangles, the reference element is an isosceles triangle with the right angle at the origin, hence the zero point is at a node. Probably not the best place, especially if that angle on the physical element approaches (or exceeds) 180 degrees. -Ben =20 -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Derek Gaston Sent: Wednesday, June 21, 2006 12:05 PM To: lib...@li... Subject: [Libmesh-devel] FE::inverse_map() not converging... Hey guys... thought I would get this out there before our meeting today so you have a chance to think about it. The newton iteration in FE::inverse_map() is not converging for me after I move nodes with my smoothing stuff. Here's what I'm getting: ############ WARNING: Newton scheme has not converged in 31 iterations: physical_point=3D(x,y,z)=3D( 0, 0.459733, 0) physical_guess=3D(x,y,z)=3D(0.115107, 0.343355, 0) dp=3D(x,y,z)=3D( 15.3103, -19.1814, 0) p=3D(x,y,z)=3D(-0.862915, -0.356698, 0) error=3D24.5424 WARNING: diff is 0.233629 point=3D(x,y,z)=3D( 0, 0.459733, 0) local=3D(x,y,z)=3D(0.164291, 0.293627, 0) lref=3D (x,y,z)=3D(-0.862915, -0.356698, 0) ########## I get a _whole bunch_ of these. I have tried to up the number of newton iterations (from 10 to 30) but that didn't seem to help... so it is definitely not converging. What's weird is that it doesn't seem to be really affecting my solutions... to the eye they look good! After looking through the inverse_map() code a bit it looks like we are taking an initial guess of the zero point... it is also mentioned that a better choice might be the centroid. Do you think this could have some effect on me? Obviously my elements must not be the "good" elements it mentions in the comments ;-) We can talk about this more this afternoon.... I just wanted to get the juices flowing. Later, Derek _______________________________________________ Libmesh-devel mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-devel |
|
From: Roy S. <roy...@ic...> - 2006-06-21 21:19:13
|
On Wed, 21 Jun 2006, Kirk, Benjamin (JSC-EG) wrote: > Which point locator class is he using? In my experience the Tree based > ones with a fairly small number of elements per bin help avoid this > problem. > > Furthermore, it is not clear to me why this is getting tripped at all, > since Elem::contains_point() calls FEInterface::inverse_map() with the > secure flag set to false? That's a very good question. I didn't notice that contains_point() turned off the secure flag. Derek, would you run your code in gdb? Stack traces of the places where convergence is failing would help a lot. --- Roy |
|
From: Derek G. <fri...@gm...> - 2006-06-22 03:08:11
|
Here's the stack: FE<2u, (libMeshEnums::FEFamily)0>::inverse_map (elem=0xed4cf0, physical_point=@0xf68920, tolerance=0.0001, secure=false) at fe_map.C:1220 FEInterface::inverse_map(dim=2, fe_t=@0x7fffff8c05a0, elem=0xed4cf0,p=@0xf68920, tolerance=0.0001, secure=false) at fe_interface.C1186 Elem::contains_point(this=0xed4cf0, p=@0xf68920) at elem.C:629 PointLocatorTree::operator()(this=0xf6d8b0, p=@0xf68920) at point_locator_tree.C:162 MeshFunction::operator()(this=0x7fffff8c0a40, p=@0xf68920, output=@0x7fffff8c07f0) at mesh_function.C:171 MeshFunction::operator()(this=0x7fffff8c0a40, p=@0xf68920, time=0) at mesh_function.C:157 main(argc=5, argv=0x7fffff8c11d8) at ex10.C:461 gdb is having a tough time actually stopping at the point I want it too... so this isnt a call that actually produces one of the error messages... but it should be close enough. The weird part to me is that secure is definitely false... but somehow the messages are still coming out?! Looking over the code the "if(secure)" part should definitely never be true.... weird. I'll keep looking. Derek On 6/21/06, Roy Stogner <roy...@ic...> wrote: > On Wed, 21 Jun 2006, Kirk, Benjamin (JSC-EG) wrote: > > > Which point locator class is he using? In my experience the Tree based > > ones with a fairly small number of elements per bin help avoid this > > problem. > > > > Furthermore, it is not clear to me why this is getting tripped at all, > > since Elem::contains_point() calls FEInterface::inverse_map() with the > > secure flag set to false? > > That's a very good question. I didn't notice that contains_point() > turned off the secure flag. > > Derek, would you run your code in gdb? Stack traces of the places > where convergence is failing would help a lot. > --- > Roy > |
|
From: Derek G. <fri...@gm...> - 2006-06-22 03:24:40
|
This is driving me nuts. GDB absolutely refuses to automatically stop on breakpoints set in fe_map. If I "Step Into" a bunch I can dig down into it manually, but I have no idea why it won't stop automatically. This is making it pretty hard to get a trace on one of the actual failing calls. One weird thing though... I changed the "Warning" print statement in inverse_map so that it would print out "secure" as well... and here's what I get: WARNING: Newton scheme has not converged in 31 iterations: physical_point=(x,y,z)=(0.158938, 0.422823, 0) physical_guess=(x,y,z)=(0.497654, 0.17708, 0) dp=(x,y,z)=(-9.77611, 7.05057, 0) p=(x,y,z)=(0.685321, 0.746925, 0) secure=1 error=12.0533 How in the _HELL_ is secure=1?!???!???!?!! Continuing to dig.... Derek On 6/21/06, Derek Gaston <fri...@gm...> wrote: > Here's the stack: > > FE<2u, (libMeshEnums::FEFamily)0>::inverse_map (elem=0xed4cf0, > physical_point=@0xf68920, tolerance=0.0001, secure=false) at > fe_map.C:1220 > FEInterface::inverse_map(dim=2, fe_t=@0x7fffff8c05a0, > elem=0xed4cf0,p=@0xf68920, tolerance=0.0001, secure=false) at > fe_interface.C1186 > Elem::contains_point(this=0xed4cf0, p=@0xf68920) at elem.C:629 > PointLocatorTree::operator()(this=0xf6d8b0, p=@0xf68920) at > point_locator_tree.C:162 > MeshFunction::operator()(this=0x7fffff8c0a40, p=@0xf68920, > output=@0x7fffff8c07f0) at mesh_function.C:171 > MeshFunction::operator()(this=0x7fffff8c0a40, p=@0xf68920, time=0) at > mesh_function.C:157 > main(argc=5, argv=0x7fffff8c11d8) at ex10.C:461 > > gdb is having a tough time actually stopping at the point I want it > too... so this isnt a call that actually produces one of the error > messages... but it should be close enough. > > The weird part to me is that secure is definitely false... but somehow > the messages are still coming out?! Looking over the code the > "if(secure)" part should definitely never be true.... weird. > > I'll keep looking. > > Derek > > On 6/21/06, Roy Stogner <roy...@ic...> wrote: > > On Wed, 21 Jun 2006, Kirk, Benjamin (JSC-EG) wrote: > > > > > Which point locator class is he using? In my experience the Tree based > > > ones with a fairly small number of elements per bin help avoid this > > > problem. > > > > > > Furthermore, it is not clear to me why this is getting tripped at all, > > > since Elem::contains_point() calls FEInterface::inverse_map() with the > > > secure flag set to false? > > > > That's a very good question. I didn't notice that contains_point() > > turned off the secure flag. > > > > Derek, would you run your code in gdb? Stack traces of the places > > where convergence is failing would help a lot. > > --- > > Roy > > > |
|
From: Roy S. <roy...@ic...> - 2006-06-22 03:38:09
|
On Wed, 21 Jun 2006, Derek Gaston wrote: > what I get: > > WARNING: Newton scheme has not converged in 31 iterations: > physical_point=(x,y,z)=(0.158938, 0.422823, 0) > physical_guess=(x,y,z)=(0.497654, 0.17708, 0) > dp=(x,y,z)=(-9.77611, 7.05057, 0) > p=(x,y,z)=(0.685321, 0.746925, 0) > secure=1 error=12.0533 > > How in the _HELL_ is secure=1?!???!???!?!! Keep in mind that secure==true is the default argument, and inverse_map gets called from more than just the point locator. Even when you're not using adaptivity, your code may be calling inverse_map from FE::nodal_soln, FE::edge_reinit, or MeshFunction::operator(). However, none of those functions should ever give you an inverse_map that doesn't converge, so it's important that we find what's wrong. At the very least, it would be a good idea to print out the node locations of the element for which inverse_map is failing. If the physical point is inside the element and the element has a positive lower bound on it's jacobian, then we'll know it's a serious problem and we'll probably have enough information to replicate it. --- Roy |
|
From: Derek G. <fri...@gm...> - 2006-06-22 19:55:28
|
Ok got it... here's a stack trace on a failing call: ########### #0 dumb_ass () at fe_map.C:34 #1 0x00002afe1d47fd5a in FE<2u, (libMeshEnums::FEFamily)0>::inverse_map (elem=0xed52b0, physical_point=@0xfe4db0, tolerance=0.0001, secure=true) at fe_map.C:1364 #2 0x00002afe1d45b565 in FEInterface::inverse_map (dim=2, fe_t=@0xeec470, elem=0xed52b0, p=@0xfe4db0, tolerance=0.0001, secure=true) at fe_interface.C:1186 #3 0x00002afe1d5f7d6c in MeshFunction::operator() (this=0x7fffffc8b2d0, p=@0xfe4db0, output=@0x7fffffc8b080) at mesh_function.C:190 #4 0x00002afe1d5f7c2c in MeshFunction::operator() (this=0x7fffffc8b2d0, p=@0xfe4db0, time=0) at mesh_function.C:157 #5 0x000000000044a155 in main (argc=5, argv=0x7fffffc8ba68) at ex10.C:461 ############ Note that I had to put an un-templated function into fe_map.C... then I called it just before the WARNING was printed... I was able to put a break point on that function therefore finally recovering the stack at one that fails. The problem was that gdb couldn't figure out which templated version of inverse_map to put a break on... and it was always guess wrong and therefore never stopping... very weird. But I worked around it ;-) The "problem" appears to be happening when mesh_function is looking for the mapped_point. I don't know enough about this stuff to interpret though. Let me know if you guys need me to do anything else. Derek On 6/21/06, Roy Stogner <roy...@ic...> wrote: > On Wed, 21 Jun 2006, Derek Gaston wrote: > > > what I get: > > > > WARNING: Newton scheme has not converged in 31 iterations: > > physical_point=(x,y,z)=(0.158938, 0.422823, 0) > > physical_guess=(x,y,z)=(0.497654, 0.17708, 0) > > dp=(x,y,z)=(-9.77611, 7.05057, 0) > > p=(x,y,z)=(0.685321, 0.746925, 0) > > secure=1 error=12.0533 > > > > How in the _HELL_ is secure=1?!???!???!?!! > > Keep in mind that secure==true is the default argument, and > inverse_map gets called from more than just the point locator. Even > when you're not using adaptivity, your code may be calling inverse_map > from FE::nodal_soln, FE::edge_reinit, or MeshFunction::operator(). > However, none of those functions should ever give you an inverse_map > that doesn't converge, so it's important that we find what's wrong. > > At the very least, it would be a good idea to print out the node > locations of the element for which inverse_map is failing. If the > physical point is inside the element and the element has a positive > lower bound on it's jacobian, then we'll know it's a serious problem > and we'll probably have enough information to replicate it. > --- > Roy > |
|
From: Roy S. <roy...@ic...> - 2006-06-22 20:05:31
|
On Thu, 22 Jun 2006, Derek Gaston wrote: > #1 0x00002afe1d47fd5a in FE<2u, > (libMeshEnums::FEFamily)0>::inverse_map (elem=0xed52b0, > physical_point=@0xfe4db0, tolerance=0.0001, secure=true) at > fe_map.C:1364 > #2 0x00002afe1d45b565 in FEInterface::inverse_map (dim=2, > fe_t=@0xeec470, elem=0xed52b0, p=@0xfe4db0, tolerance=0.0001, > secure=true) at fe_interface.C:1186 > #3 0x00002afe1d5f7d6c in MeshFunction::operator() > (this=0x7fffffc8b2d0, p=@0xfe4db0, output=@0x7fffffc8b080) at > mesh_function.C:190 > #4 0x00002afe1d5f7c2c in MeshFunction::operator() > (this=0x7fffffc8b2d0, p=@0xfe4db0, time=0) at mesh_function.C:157 > #5 0x000000000044a155 in main (argc=5, argv=0x7fffffc8ba68) at ex10.C:461 > ############ > > The "problem" appears to be happening when mesh_function is looking > for the mapped_point. I don't know enough about this stuff to > interpret though. It looks like the failing call to inverse_map occurs after mesh_function thinks it's already found the element containing p. So, the first thing we should do is check and make sure that's true. Would you give us the x/y coordinates of p and of every node (in local order) of elem? --- Roy |
|
From: Derek G. <fri...@gm...> - 2006-06-22 21:55:27
|
It does appear that it is finding the wrong element... here is an element and a point for which the newton iterations fail (as they very well probably should since the point is outside the element...) ################ Node 0:0.131279 0 Node 1:0.17075 0 Node 2:0.167717 0.032511 Node 3:0.131771 0.0454457 Point: 0.184449 0.121701 ##################### They are x y pairs. The nodes are the nodes of the "element" pointer at around line 192 of mesh_function.C the point is the x y pairs of "p" in the same place. It would appear that "p" is _not_ in "element".... Derek On 6/22/06, Roy Stogner <roy...@ic...> wrote: > On Thu, 22 Jun 2006, Derek Gaston wrote: > > > #1 0x00002afe1d47fd5a in FE<2u, > > (libMeshEnums::FEFamily)0>::inverse_map (elem=0xed52b0, > > physical_point=@0xfe4db0, tolerance=0.0001, secure=true) at > > fe_map.C:1364 > > #2 0x00002afe1d45b565 in FEInterface::inverse_map (dim=2, > > fe_t=@0xeec470, elem=0xed52b0, p=@0xfe4db0, tolerance=0.0001, > > secure=true) at fe_interface.C:1186 > > #3 0x00002afe1d5f7d6c in MeshFunction::operator() > > (this=0x7fffffc8b2d0, p=@0xfe4db0, output=@0x7fffffc8b080) at > > mesh_function.C:190 > > #4 0x00002afe1d5f7c2c in MeshFunction::operator() > > (this=0x7fffffc8b2d0, p=@0xfe4db0, time=0) at mesh_function.C:157 > > #5 0x000000000044a155 in main (argc=5, argv=0x7fffffc8ba68) at ex10.C:461 > > ############ > > > > The "problem" appears to be happening when mesh_function is looking > > for the mapped_point. I don't know enough about this stuff to > > interpret though. > > It looks like the failing call to inverse_map occurs after > mesh_function thinks it's already found the element containing p. So, > the first thing we should do is check and make sure that's true. > Would you give us the x/y coordinates of p and of every node (in local > order) of elem? > --- > Roy > |
|
From: Roy S. <roy...@ic...> - 2006-06-22 22:08:05
|
On Thu, 22 Jun 2006, Derek Gaston wrote: > It does appear that it is finding the wrong element... here is an > element and a point for which the newton iterations fail (as they very > well probably should since the point is outside the element...) > > ################ > Node 0:0.131279 0 > Node 1:0.17075 0 > Node 2:0.167717 0.032511 > Node 3:0.131771 0.0454457 > Point: 0.184449 0.121701 > ##################### > > They are x y pairs. The nodes are the nodes of the "element" pointer > at around line 192 of mesh_function.C the point is the x y pairs of > "p" in the same place. > > It would appear that "p" is _not_ in "element".... Okay, then our next step is figuring out why the code thinks otherwise. The MeshFunction shouldn't be trying to use an element for which element->contains_point(p) hasn't returned true. If you run Elem::contains_point on that element with that p, does it return true or false? --- Roy |
|
From: Derek G. <fri...@gm...> - 2006-06-23 04:25:15
|
It would appear that for that element and that point element->contains_point(p) _is_ returning true! There's definitely something going on there! Derek On 6/22/06, Roy Stogner <roy...@ic...> wrote: > On Thu, 22 Jun 2006, Derek Gaston wrote: > > > It does appear that it is finding the wrong element... here is an > > element and a point for which the newton iterations fail (as they very > > well probably should since the point is outside the element...) > > > > ################ > > Node 0:0.131279 0 > > Node 1:0.17075 0 > > Node 2:0.167717 0.032511 > > Node 3:0.131771 0.0454457 > > Point: 0.184449 0.121701 > > ##################### > > > > They are x y pairs. The nodes are the nodes of the "element" pointer > > at around line 192 of mesh_function.C the point is the x y pairs of > > "p" in the same place. > > > > It would appear that "p" is _not_ in "element".... > > Okay, then our next step is figuring out why the code thinks > otherwise. The MeshFunction shouldn't be trying to use an element for > which element->contains_point(p) hasn't returned true. If you run > Elem::contains_point on that element with that p, does it return true > or false? > --- > Roy > |
|
From: Derek G. <fri...@gm...> - 2006-06-23 04:38:08
|
Ok... so here's the deal. elem::contains_point() is using the same flawed logic to decide whether or not the point is contained in the element (meaning the newton iterations). What's happening is that the newton iterations are failing to converge so severely that the final guess as where the point is... actually lies within the element! Thus giving a false positive. So here's the deal... inverse_map is _NOT_ a good idea to use for finding out if the point lies within the element! My question is this: contains_point() is a virtual function right? Can we not just override it in the concrete classes so that we don't have to do this newton thing at all? I mean a quad4 should be able to figure out pretty quickly whether or not a point lies in it.... without all the newton iteration stuff. Sure we can leave the newton stuff to fall back on in case we don't want to implement contains_point in a concrete_class... but for elements for which it would be easy (maybe even trivial?) to do it for manually... why not? Open for comments Derek PS: Yes, this is me "volunteering" to go write contains_point() functions for all our dang elements ;-) On 6/22/06, Derek Gaston <fri...@gm...> wrote: > It would appear that for that element and that point > element->contains_point(p) _is_ returning true! There's definitely > something going on there! > > Derek > > On 6/22/06, Roy Stogner <roy...@ic...> wrote: > > On Thu, 22 Jun 2006, Derek Gaston wrote: > > > > > It does appear that it is finding the wrong element... here is an > > > element and a point for which the newton iterations fail (as they very > > > well probably should since the point is outside the element...) > > > > > > ################ > > > Node 0:0.131279 0 > > > Node 1:0.17075 0 > > > Node 2:0.167717 0.032511 > > > Node 3:0.131771 0.0454457 > > > Point: 0.184449 0.121701 > > > ##################### > > > > > > They are x y pairs. The nodes are the nodes of the "element" pointer > > > at around line 192 of mesh_function.C the point is the x y pairs of > > > "p" in the same place. > > > > > > It would appear that "p" is _not_ in "element".... > > > > Okay, then our next step is figuring out why the code thinks > > otherwise. The MeshFunction shouldn't be trying to use an element for > > which element->contains_point(p) hasn't returned true. If you run > > Elem::contains_point on that element with that p, does it return true > > or false? > > --- > > Roy > > > |
|
From: Roy S. <roy...@ic...> - 2006-06-23 15:38:20
|
On Thu, 22 Jun 2006, Derek Gaston wrote: > elem::contains_point() is using the same flawed logic to decide > whether or not the point is contained in the element (meaning the > newton iterations). What's happening is that the newton iterations > are failing to converge so severely that the final guess as where the > point is... actually lies within the element! Thus giving a false > positive. This is what I was afraid of - the fix is as I suggested the other day: if the Newton iteration doesn't converge, return a point like (1.e6,1.e6,1.e6) so we don't get a false positive from contains_point(). > So here's the deal... inverse_map is _NOT_ a good idea to use for > finding out if the point lies within the element! I think I disagree. These Newton failures should only happen if the element is badly twisted or if the point is outside the element, right? So in that case, if we just return a junk point when Newton fails, we should still get correct results. > My question is this: contains_point() is a virtual function right? > Can we not just override it in the concrete classes so that we don't > have to do this newton thing at all? I mean a quad4 should be able to > figure out pretty quickly whether or not a point lies in it.... > without all the newton iteration stuff. Sure we can leave the newton > stuff to fall back on in case we don't want to implement > contains_point in a concrete_class... but for elements for which it > would be easy (maybe even trivial?) to do it for manually... why not? > PS: Yes, this is me "volunteering" to go write contains_point() > functions for all our dang elements ;-) Works for me! I'm not sure how much of an improvement you can make in general, but at the very least it might be a little optimization for TRI3/QUAD4. For the more complicated geometric elements you might try testing has_affine_map() first - if the element transformation is affine (and I don't know about your weird smoothed meshes but in my meshes most of the elements are affine) then you can implement contains_point with n_sides inequality tests. --- Roy |
|
From: Kirk, B. \(JSC-EG\) <ben...@na...> - 2006-06-23 17:03:39
|
Nice effort, y'all. I agree with both (heh, how political...) of you. inverse_map() is probably as good as anything to robustly determine if an elem->contains_point(), but I must confess when I wrote it I envisioned that the point *would* be in the element. It was only later I thought about using it for elem->contains_point(). Clearly inverse_map() needs to return a bad point when iteration fails, right now it is luck of the draw... Similarly, I envisoned elem->contains_point() as generally being called when the element actually contains the point. A simple improvement could be to check the bounding box of the element's nodes before beginning Newton iteration. Clearly if a point is not contained in the Cartesian bounding box of an element's nodes (with some allowance for delta_z~0 in a 2D mesh) we have no business expecting the inverse_map() to work. Certainly optimizations are possible for certain elements, especially linear simplices, but you might want to profile the code to see if it is worth the effort. If you are volunteering to write some code I'd love to have you help me update the XDR mesh format... ;-) -Ben -----Original Message----- From: lib...@li... [mailto:lib...@li...] On Behalf Of Roy Stogner Sent: Friday, June 23, 2006 10:38 AM To: Derek Gaston Cc: lib...@li... Subject: Re: [Libmesh-devel] FE::inverse_map() not converging... On Thu, 22 Jun 2006, Derek Gaston wrote: > elem::contains_point() is using the same flawed logic to decide=20 > whether or not the point is contained in the element (meaning the=20 > newton iterations). What's happening is that the newton iterations=20 > are failing to converge so severely that the final guess as where the=20 > point is... actually lies within the element! Thus giving a false=20 > positive. This is what I was afraid of - the fix is as I suggested the other day: if the Newton iteration doesn't converge, return a point like (1.e6,1.e6,1.e6) so we don't get a false positive from contains_point(). > So here's the deal... inverse_map is _NOT_ a good idea to use for=20 > finding out if the point lies within the element! I think I disagree. These Newton failures should only happen if the element is badly twisted or if the point is outside the element, right? So in that case, if we just return a junk point when Newton fails, we should still get correct results. > My question is this: contains_point() is a virtual function right? > Can we not just override it in the concrete classes so that we don't=20 > have to do this newton thing at all? I mean a quad4 should be able to > figure out pretty quickly whether or not a point lies in it.... > without all the newton iteration stuff. Sure we can leave the newton=20 > stuff to fall back on in case we don't want to implement=20 > contains_point in a concrete_class... but for elements for which it=20 > would be easy (maybe even trivial?) to do it for manually... why not? > PS: Yes, this is me "volunteering" to go write contains_point()=20 > functions for all our dang elements ;-) Works for me! I'm not sure how much of an improvement you can make in general, but at the very least it might be a little optimization for TRI3/QUAD4. For the more complicated geometric elements you might try testing has_affine_map() first - if the element transformation is affine (and I don't know about your weird smoothed meshes but in my meshes most of the elements are affine) then you can implement contains_point with n_sides inequality tests. --- Roy Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D120709&bid=3D263057&dat=3D= 121642 _______________________________________________ Libmesh-devel mailing list Lib...@li... https://lists.sourceforge.net/lists/listinfo/libmesh-devel |
|
From: Roy S. <roy...@ic...> - 2006-06-23 17:11:40
|
On Fri, 23 Jun 2006, Kirk, Benjamin (JSC-EG) wrote:
> Similarly, I envisoned elem->contains_point() as generally being called
> when the element actually contains the point.
If that was the case, wouldn't contains_point(){return true;} have
been much easier to code? ;-)
> A simple improvement could be to check the bounding box of the
> element's nodes before beginning Newton iteration.
> Clearly if a point is not contained in the Cartesian bounding box of
> an element's nodes (with some allowance for delta_z~0 in a 2D mesh)
> we have no business expecting the inverse_map() to work.
Yes, but the opposite isn't true. Just picture a QUAD4 with one
corner gradually being brought toward the center. There will
eventually be a point where the Jacobian inside the element is
positively bounded below (and so it's still an acceptable element) but
where there are points in the bounding box with inverted Jacobians
that might trash Newton.
> Certainly optimizations are possible for certain elements, especially
> linear simplices, but you might want to profile the code to see if it is
> worth the effort. If you are volunteering to write some code I'd love
> to have you help me update the XDR mesh format... ;-)
Hey! If we've got a volunteer on our hands, I want to see the FE
shape function initialization better optimized and the build process
moved to automake.
---
Roy
|
|
From: Roy S. <roy...@ic...> - 2006-06-23 18:43:11
|
On Fri, 23 Jun 2006, Roy Stogner wrote: > Hey! If we've got a volunteer on our hands, I want to see the FE > shape function initialization better optimized and the build process > moved to automake. You know, as long as I'm writing my wishlist: For quite a while now (ever since I wrote that "subactive elements" hack), we've had the capability to get multigrid linear systems out of any System class with a side effect free assembly() function. We've also got most of the code necessary to write transfer matrices in System::project_vector(). If anyone's looking for a thankless coding task, it would be nice to figure out how to hook all this stuff into PETSc and LASPACK. --- Roy |
|
From: Derek G. <fri...@gm...> - 2006-06-23 19:29:34
|
Wow... now I know better than to volunteer my time out loud ;-) After thinking more about this while laying in bed last night.... I think just returning a point that is obviously outside of the reference element will be fine for now. I looked at some of the inverse_map calls that converged and they usually do it very quickly... and really it's not that many ops (surely not enough to pull me off my other tasks). So for now I think this will do. Also... I think it would actually be kinda difficult to manually code up contains_point() for some of the elements... BTW... do we _not_ have this inverse map coded somewhere? I understand that inverse_map is used in other circumstances as well... but for contains_point() don't we have a jacobian sitting around somewhere that we could use to compute where the point is in the reference domain? Or maybe I'm not understanding things completely.... Roy... moving the build process to Automake would be a great project for me... and I would love to do it. Let me get a little more off my plate and I'll start looking into it. One of the reasons I would love to see this happen is so that IDE's (such as Kdevelop) would work a lot better with libmesh. Currently they are just glorified editors.... Thanks to everyone who helped me with debugging this! I really appreciate it! Derek On 6/23/06, Roy Stogner <roy...@ic...> wrote: > On Fri, 23 Jun 2006, Roy Stogner wrote: > > > Hey! If we've got a volunteer on our hands, I want to see the FE > > shape function initialization better optimized and the build process > > moved to automake. > > You know, as long as I'm writing my wishlist: > > For quite a while now (ever since I wrote that "subactive elements" > hack), we've had the capability to get multigrid linear systems out of > any System class with a side effect free assembly() function. We've > also got most of the code necessary to write transfer matrices in > System::project_vector(). If anyone's looking for a thankless coding > task, it would be nice to figure out how to hook all this stuff into > PETSc and LASPACK. > --- > Roy > > Using Tomcat but need to do more? Need to support web services, security? > Get stuff done quickly with pre-integrated technology to make your job easier > Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo > http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642 > _______________________________________________ > Libmesh-devel mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > |
|
From: Roy S. <roy...@ic...> - 2006-06-21 19:40:51
|
On Wed, 21 Jun 2006, Kirk, Benjamin (JSC-EG) wrote: > Which type of elements are these? They're QUAD4 elements. > For Quads the zero point corresponds to the centroid of the element, > which should be as good a guess as any. For triangles, the reference > element is an isosceles triangle with the right angle at the origin, > hence the zero point is at a node. Probably not the best place, > especially if that angle on the physical element approaches (or exceeds) > 180 degrees. Starting at element centroids (or approximate element centroids) would probably be more robust, but we've looked at things in detail and found two more potential problems: 1. Elem::contains_point passes a tolerance of 1e-4 to inverse_map() but then uses the default tolerance of 1e-6 for on_reference_element() This means that it's slightly more likely for contains_point to return a false positive (not a big problem), but it also means it's possible for contains_point to return a false negative (which can break programs). 2. Because Elem::contains_point uses inverse_map() and because it gets called even for points which are far away from the physical elements (in the point locator functions is where Derek's code tripped it), if the physical elements are badly distorted then the Newton iteration won't converge. For now Derek's just going to silence the Newton warning messages and return a grossly invalid point when Newton fails to converge. Is that safe to do in general, though? It seems like there are many situations where inverse_map gets called and we expect to get a valid answer, and it would be nice to keep asserting convergence in DEBUG mode for such calls. We should avoid changing the declaration of inverse_map() if possible, so there's no way for inverse_map() itself to distinguish between the "must converge" and "may not converge" cases, but on the other hand inverse_map() get called dosens of times by the library and it would be tedious to add asserts() after every such call. --- Roy |
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X