From: Roy S. <roy...@ic...> - 2006-07-24 17:07:50
|
On Mon, 24 Jul 2006, Tim Kr=F6ger wrote: > On Thu, 20 Jul 2006, John Peterson wrote: > >> I got your message. Was there another question you had? I didn't >> see one, hence the lack of response. > > Well, the question is still whether you would accept such a port in you= r=20 > repository. I think we're all willing to tolerate moving variable declarations out of for loops. That's a little confusing to me as a C++ programmer (I'm then tempted to hunt through the code after the loop and see where the index variable is going to be reused), but it's worth it for wider compiler support. You're never going to be able to get other people to write new loops that way, but that just means you'll occasionally need to submit another small patch to keep the CVS head VC6-compatible. I'm just doubtful that this variable scoping problem is going to be the only problem Visual C++ 6.0 has. We use a lot of templates, and I find it hard to imagine a C++ compiler getting even simple template code right if it still has troubles with the basics of the language. If you're asking whether we'd accept a patch reverting all the declarations inside for loops, the answer is yes (unless Ben strongly objects), but if you're asking whether we'd accept a patch which makes every change that might be required for VC6 compatibility, nobody can answer that until we've seen it. > However, I just face a different question: It seems to me that libMesh = has=20 > got a memory leak. Please find attached a very simple test program tha= t=20 > creates a grid and initializes an EquationSystems object on that grid. = Both=20 > then go out of scope, and the whole procedure is repeated infinitely ma= ny=20 > times. If I watch the memory usage of the program (e.g. using "top"), = I see=20 > that it increases roughly by about 10MB in each cycle. After some minu= tes,=20 > my computer's memory has been eaten up completely. Can you confirm thi= s=20 > behaviour? Using the CVS libMesh with --enable-everything and METHOD=3Dopt, the memory use seems to oscillate between 97MB and 101MB here, and stays that way after a couple hundred cycles. Is this a memory leak you're seeing with 0.5.0, or can you duplicate it on the CVS head too? > By the way, is there an official bug reporting procedure? If you're certain it's a bug in the library (and if you're certain you can replicate it in CVS - we're not backporting bugfixes to version 0.x code), send it to the libmesh-devel mailing list; if you're not srue where the bug is then send it to libmesh-users. There's probably some bug tracking system on sourceforge, but I don't know if it gets as much attention. Actually, this whole discussion ought to be moved to libmesh-devel so there's a public archive of it.. and I'm not just saying that because someone typoed my email address and left me out of half of it. ;-) --- Roy Stogner |
From: Roy S. <roy...@ic...> - 2006-08-03 14:54:27
|
On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > Dear libMesh developers, You ought to send stuff like this to libmesh-devel, unless it's really time critical and Sourceforge's mailing lists are acting up again. > Could you please commit the attached patch to the libMesh cvs repositor= y? It=20 > consists of the following changes: > > 1.) Three changes in the MeshFunction class (see discussion in libmesh-= users=20 > dated July 31, 2006); You might want look over the previous recent discussion on the mailing lists before you come to rely on MeshFunction too heavily. There's a lot of suspicious code in there that makes me believe it won't work on an adaptively h-refined mesh. > 2.) EquationSystems::read() now also recognizes the type ExplicitSystem. It looked good to me, so I'll commit it just as soon as I'm sure gcc agrees. I don't use MeshFunction myself, however - could Derek or someone who does doublecheck that part of the patch? Bearing in mind that if you've already had to work around MeshFunction's bugs yourself, Tim's fix is likely to break your workarounds. --- Roy |
From: <ti...@ce...> - 2006-08-03 15:19:48
|
Dear Roy, On Thu, 3 Aug 2006, Roy Stogner wrote: > On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > You might want look over the previous recent discussion on the mailing > lists before you come to rely on MeshFunction too heavily. There's a > lot of suspicious code in there that makes me believe it won't work > on an adaptively h-refined mesh. Well, it seems to work for my meshes (which are h-refined but not=20 p-refined) with the mentioned modifications. Without these=20 modifications, however, it was impossible to use MeshFunction at all=20 (unless using some workarounds). Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-08-03 15:31:39
|
On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: >> You might want look over the previous recent discussion on the mailing >> lists before you come to rely on MeshFunction too heavily. There's a >> lot of suspicious code in there that makes me believe it won't work >> on an adaptively h-refined mesh. > > Well, it seems to work for my meshes (which are h-refined but not p-ref= ined)=20 > with the mentioned modifications. You might want to set up some validation problems and make sure that "seems to work" is really "works exactly". I don't use MeshFunction myself, so it's possible that I'm just being paranoid, but it looks to me like the PointLocator classes might return a coarse element when they should be returning one of it's descendant elements. If that happens, then as long as you're using Lagrange elements MeshFunction may still give you a solution without crashing, but the solution will be slightly inaccurate. --- Roy |
From: <ti...@ce...> - 2006-08-04 07:16:41
|
Dear Roy and all, On Thu, 3 Aug 2006, Roy Stogner wrote: > On Thu, 3 Aug 2006, Tim Kr=F6ger wrote: > >>> You might want look over the previous recent discussion on the mailing >>> lists before you come to rely on MeshFunction too heavily. There's a >>> lot of suspicious code in there that makes me believe it won't work >>> on an adaptively h-refined mesh. >>=20 >> Well, it seems to work for my meshes (which are h-refined but not=20 >> p-refined) with the mentioned modifications. > > You might want to set up some validation problems and make sure that > "seems to work" is really "works exactly". I don't use MeshFunction > myself, so it's possible that I'm just being paranoid, but it looks to > me like the PointLocator classes might return a coarse element when they > should be returning one of it's descendant elements. If that happens, > then as long as you're using Lagrange elements MeshFunction may still > give you a solution without crashing, but the solution will be > slightly inaccurate. Thank you for your hint; I just looked at the code more closely and=20 now agree that there might be two problems in the code: 1.) Unless I missed something, not only active elements are stored in=20 the Tree object, so that, as you mentioned, the PointLocator class=20 might return a coarse element instead of a descendant. This could be fixed by letting Tree::Tree() loop over active elements=20 only. 2.) Elements are arranged by the location of their centroid only,=20 although they in fact might happen to have non-empty intersection with=20 the bounding boxes of more than one of the children of a TreeNode=20 object. After fixing Problem 1, this will cause Tree::find_element()=20 to occasionally return NULL -- i.e. if the point to look for lies in a=20 different TreeNode's bounding box than the centroid of the=20 corresponding element. In this case, PointLocatorTree::operator()=20 will fall back to a linear search and thus still find the correct=20 element but be slow. I think that this problem also can be solved quite easily: One would=20 just have to insert such elements in all the TreeNode's children that=20 they have non-empty intersection with. I would not refuse to make these changes to the code myself, but=20 before, I would like to wait for your comments. Another point is that I would like MeshFunction to be evaluable=20 outside the mesh domain and return a used-defined constant there=20 (instead of crashing). I think I would be able to do these changes as=20 well (using optional parameters so that existing code will not change=20 its behaviour). Any comments on this? Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-08-04 11:46:01
|
On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > On Thu, 3 Aug 2006, Roy Stogner wrote: > > 1.) Unless I missed something, not only active elements are stored in t= he=20 > Tree object, so that, as you mentioned, the PointLocator class might re= turn a=20 > coarse element instead of a descendant. > > This could be fixed by letting Tree::Tree() loop over active elements o= nly. I don't think that's a complete fix. The problem with mesh refinement / coarsening isn't just that the mesh is more complicated, it's that the mesh keeps *changing*. We need to make sure that new elements are inserted into the Tree after every refinement step and removed after every coarsening - and just as importantly for the modularity of the code, we ought to have that happen without messing too badly with the Mesh and Elem classes, and it ought to work even if the user never calls anything but MeshFunction::operator() after the AMR/C. > 2.) Elements are arranged by the location of their centroid only, altho= ugh=20 > they in fact might happen to have non-empty intersection with the bound= ing=20 > boxes of more than one of the children of a TreeNode object. After fix= ing=20 > Problem 1, this will cause Tree::find_element() to occasionally return = NULL=20 > -- i.e. if the point to look for lies in a different TreeNode's boundin= g box=20 > than the centroid of the corresponding element. In this case,=20 > PointLocatorTree::operator() will fall back to a linear search and thus= still=20 > find the correct element but be slow. > > I think that this problem also can be solved quite easily: One would ju= st=20 > have to insert such elements in all the TreeNode's children that they h= ave=20 > non-empty intersection with. That sounds like the best way to do things. > I would not refuse to make these changes to the code myself, but before= , I=20 > would like to wait for your comments. You'd be welcome (and thanked!) for the fixes. I'd like it if we could figure out how to make fix 1 a little more robust, but even just looping over active elements would be better than nothing. I'd rather place a warning in the MeshFunction documentation than let this class hold up a 0.6.0 release, and it would be better for that warning to say "rebuild this object after every mesh change" and not "don't use this object on a refined mesh". > Another point is that I would like MeshFunction to be evaluable outside= the=20 > mesh domain and return a used-defined constant there (instead of crashi= ng).=20 > I think I would be able to do these changes as well (using optional=20 > parameters so that existing code will not change its behaviour). Any=20 > comments on this? Adding an optional parameter to the constructor would be fine. Adding a method that sets that parameter might be better - although that's kind of a matter of taste. --- Roy |
From: <ti...@ce...> - 2006-08-04 20:50:46
Attachments:
patch
|
Dear Roy, On Fri, 4 Aug 2006, Roy Stogner wrote: > On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > >> This could be fixed by letting Tree::Tree() loop over active elements on= ly. > > I don't think that's a complete fix. The problem with mesh refinement > / coarsening isn't just that the mesh is more complicated, it's that > the mesh keeps *changing*. Okay, I did not think about this. I don't use MeshFunction with=20 changing grids -- i.e., I rebuild the MeshFunction anyway each time=20 the mesh might change. > You'd be welcome (and thanked!) for the fixes. I'd like it if we > could figure out how to make fix 1 a little more robust, but even just > looping over active elements would be better than nothing. I'd rather > place a warning in the MeshFunction documentation than let this class > hold up a 0.6.0 release, and it would be better for that warning to > say "rebuild this object after every mesh change" and not "don't use > this object on a refined mesh". I did these changes now; see the attached patch. Please commit it to=20 the repository (provided that you are satisfied with it). I did set up a test scenario in which the old implementation gave=20 wrong results and the new one got it right. I also removed the fall-back linear search since in my opinion it is=20 not necessary any more. Anyway, if the linear search includes=20 non-active elements (as it did), it would again give wrong results,=20 and if not, it is very likely to crash, because the involved inverse=20 map seems to be very instable in the case that the point is very far=20 away from the element. At least, this is true (and not surprising)=20 for distortet HEX8 elements. Additionally, I added a function get_point_locator() to MeshFunction,=20 because in my application, I want to reuse the MeshFunction's=20 PointLocator for additional purposes. >> Another point is that I would like MeshFunction to be evaluable outside = the=20 >> mesh domain and return a used-defined constant there (instead of crashin= g).=20 >> I think I would be able to do these changes as well (using optional=20 >> parameters so that existing code will not change its behaviour). Any=20 >> comments on this? > Adding an optional parameter to the constructor would be fine. Adding > a method that sets that parameter might be better - although that's > kind of a matter of taste. I did not implement this yet, but I will probably do this later=20 (unless someone else does) because I am likely to need it. Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-08-04 21:27:21
|
On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > I did these changes now; see the attached patch. Please commit it to t= he=20 > repository (provided that you are satisfied with it). One problem I see: in the case of quadratic mappings (or any more complicated mappings, if the University of Auckland folks decide to code them!) it's possible for an element's nodes to all be inside a bounding box even though part of the element interior extends outside the box. I don't know of any fast way to get the real bounding box for even a quadratic element. This will be a very uncommon problem, though - perhaps an adequate solution for now is just to leave in the fallback linear search. I'll add that search back in (but changed to loop over only active elements) and commit the patch to CVS this afternoon. > I also removed the fall-back linear search since in my opinion it is no= t=20 > necessary any more. Anyway, if the linear search includes non-active=20 > elements (as it did), it would again give wrong results, and if not, it= is=20 > very likely to crash, because the involved inverse map seems to be very= =20 > instable in the case that the point is very far away from the element. = At=20 > least, this is true (and not surprising) for distortet HEX8 elements. The inverse map iteration definitely won't converge for some elements if the physical point is outside the element, but that shouldn't cause a problem for Elem::constains_point() anymore; if its call to inverse_map() doesn't converge, the result is a distant point that will make sure on_reference_element() is false. > Additionally, I added a function get_point_locator() to MeshFunction, b= ecause=20 > in my application, I want to reuse the MeshFunction's PointLocator for=20 > additional purposes. Looks good. --- Roy |
From: <ti...@ce...> - 2006-08-07 09:43:38
|
Dear Roy, On Fri, 4 Aug 2006, Roy Stogner wrote: > On Fri, 4 Aug 2006, Tim Kr=F6ger wrote: > >> I did these changes now; see the attached patch. Please commit it to th= e=20 >> repository (provided that you are satisfied with it). > > One problem I see: in the case of quadratic mappings (or any more > complicated mappings, if the University of Auckland folks decide to > code them!) it's possible for an element's nodes to all be inside a > bounding box even though part of the element interior extends outside > the box. You are right, I did not think about this. > I don't know of any fast way to get the real bounding box for even a > quadratic element. This will be a very uncommon problem, though - > perhaps an adequate solution for now is just to leave in the fallback > linear search. > The inverse map iteration definitely won't converge for some elements > if the physical point is outside the element, but that shouldn't cause > a problem for Elem::constains_point() anymore; if its call to > inverse_map() doesn't converge, the result is a distant point that > will make sure on_reference_element() is false. No: If the nonlinear solver doesn't converge, it will call error() and=20 thus crash. At least, this call to error() should be removed --=20 perhaps by giving the solver class an optional argument. The other problem that I see when leaving the linear search in is that=20 when I will later implement the possibiliy to evaluate the=20 MeshFunction outside the domain covered by the grid, the linear search=20 will slow this down considerabely. What about adding some heuristic safety margin around the bounding box=20 in the case of quadratic mappings? Well, not really satisfying ... Perhaps better: Enable the linear search only if quadratic mappings=20 are involved. Any suggestions? Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-08-07 20:18:39
|
On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > On Fri, 4 Aug 2006, Roy Stogner wrote: > >> The inverse map iteration definitely won't converge for some elements >> if the physical point is outside the element, but that shouldn't cause >> a problem for Elem::constains_point() anymore; if its call to >> inverse_map() doesn't converge, the result is a distant point that >> will make sure on_reference_element() is false. > > No: If the nonlinear solver doesn't converge, it will call error() and = thus=20 > crash. At least, this call to error() should be removed -- perhaps by=20 > giving the solver class an optional argument. 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. Of course, inverse_map() is *supposed* to exit by returning an out-of-bounds point, not just whatever its last iterate was. Derek, didn't we fix this? Or did you send me a patch I forgot to commit to CVS? > The other problem that I see when leaving the linear search in is that = when=20 > I will later implement the possibiliy to evaluate the MeshFunction outs= ide=20 > the domain covered by the grid, the linear search will slow this down=20 > considerabely. That is undesirable. > What about adding some heuristic safety margin around the bounding box = in=20 > the case of quadratic mappings? Well, not really satisfying ... As long as the margin is a proveable upper bound, I'll be happy with it. Trying to figure out what would be an upper bound on a quadratic Hex makes my head spin, though. > Perhaps better: Enable the linear search only if quadratic mappings are= =20 > involved. That sounds good too, but deciding whether or not quadratic mappings are involved is easier said than done. I'd go so far as to say that second-order geometric elements with linear mappings is one of the most common usage patterns in libMesh, since we need the extra nodes in the second-order elements to store degree of freedom numbers for any p>1 finite element. > Any suggestions? No good ones. All I can think to do is warn about the problem. When you add an interface for out-of-mesh 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. --- Roy |
From: <ti...@ce...> - 2006-08-08 06:45:10
|
Dear Roy, On Mon, 7 Aug 2006, Roy Stogner wrote: > On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > >> On Fri, 4 Aug 2006, Roy Stogner wrote: >>=20 >>> The inverse map iteration definitely won't converge for some elements >>> if the physical point is outside the element, but that shouldn't cause >>> a problem for Elem::constains_point() anymore; if its call to >>> inverse_map() doesn't converge, the result is a distant point that >>> will make sure on_reference_element() is false. >>=20 >> No: If the nonlinear solver doesn't converge, it will call error() and t= hus=20 >> crash. At least, this call to error() should be removed -- perhaps by= =20 >> giving the solver class an optional argument. > > 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=20 code such that only active elements were inserted into the tree but=20 they were still handled only by their centroid: The code crashed=20 because Newton didn't converge within 20 iterations (with warnings=20 from the 11th iteration on). I have no explanation for that because=20 when I look at the code, it seems to me that you are right. >> Any suggestions? > > No good ones. All I can think to do is warn about the problem. When > you add an interface for out-of-mesh 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=20 out-of-mesh-extension of MeshFunction is enabled, it will check=20 Elem::has_affine_map() for all elements (at least in debug mode), and=20 if any element returns true, an error will be generated. But,=20 however, I wonder whether Elem::has_affine_map() could occasionally=20 return false due to rounding errors. What do you think? Should the=20 equality checks be replaced by some fabs()<eps stuff? Best Regards, Tim |
From: John P. <pet...@cf...> - 2006-08-08 13:45:02
|
Tim Kr=F6ger writes: > Dear Roy, >=20 > On Mon, 7 Aug 2006, Roy Stogner wrote: >=20 > > On Mon, 7 Aug 2006, Tim Kr=F6ger wrote: > > > >> On Fri, 4 Aug 2006, Roy Stogner wrote: > >>=20 > >>> The inverse map iteration definitely won't converge for some ele= ments > >>> if the physical point is outside the element, but that shouldn't= cause > >>> a problem for Elem::constains_point() anymore; if its call to > >>> inverse_map() doesn't converge, the result is a distant point th= at > >>> will make sure on_reference_element() is false. > >>=20 > >> No: If the nonlinear solver doesn't converge, it will call error(= ) and thus=20 > >> crash. At least, this call to error() should be removed -- perha= ps by=20 > >> giving the solver class an optional argument. > > > > That's exactly why the inverse_map function has the "bool secure" > > argument. If secure is false (as it is when contains_point() call= s > > 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=20= > code such that only active elements were inserted into the tree but=20= > they were still handled only by their centroid: The code crashed=20 > because Newton didn't converge within 20 iterations (with warnings=20= > from the 11th iteration on). I have no explanation for that because= =20 > when I look at the code, it seems to me that you are right. >=20 > >> Any suggestions? > > > > No good ones. All I can think to do is warn about the problem. W= hen > > you add an interface for out-of-mesh 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 th= at > > it doesn't see any quadratic mappings. >=20 > Okay, so I think the best thing to do is: Whenever the=20 > out-of-mesh-extension of MeshFunction is enabled, it will check=20 > Elem::has_affine_map() for all elements (at least in debug mode), an= d=20 > if any element returns true, an error will be generated. But,=20 > however, I wonder whether Elem::has_affine_map() could occasionally=20= > return false due to rounding errors. What do you think? Should the= =20 > equality checks be replaced by some fabs()<eps stuff? I think the point=3D=3D tests do have fabs() tests in them. -J |
From: <ti...@ce...> - 2006-08-08 14:07:59
|
Dear John and Roy, On Tue, 8 Aug 2006, John Peterson wrote: > Tim Kr=F6ger writes: > > Okay, so I think the best thing to do is: Whenever the > > out-of-mesh-extension of MeshFunction is enabled, it will check > > Elem::has_affine_map() for all 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 think? Should the > > equality checks be replaced by some fabs()<eps stuff? > > I think the point=3D=3D tests do have fabs() tests in them. You are right, thank you (except that the doxygen docu tells something=20 wrong about Point::operator=3D=3D and Point::operator!=3D). Hence,=20 everything is clear now unless someone objects to the aforementioned=20 suggestion. Best Regards, Tim |
From: Roy S. <roy...@ic...> - 2006-08-08 15:26:45
|
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=20 > that only active elements were inserted into the tree but they were sti= ll=20 > handled only by their centroid: The code crashed because Newton didn't=20 > converge within 20 iterations (with warnings from the 11th iteration on= ). I=20 > have no explanation for that because when I look at the code, it seems = to me=20 > 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 double-check 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 out-of-mesh 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 out-of-mesh-exte= nsion=20 > of MeshFunction is enabled, it will check Elem::has_affine_map() for al= l=20 > elements (at least in debug mode), and if any element returns true, an = error=20 > will be generated. But, however, I wonder whether Elem::has_affine_map= ()=20 > could occasionally return false due to rounding errors. What do you th= ink?=20 > 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 |
From: Derek G. <fri...@gm...> - 2006-08-08 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 <roy...@ic...> 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 double-check 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 out-of-mesh 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 out-of-mesh-exte= 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 pre-integrated technology to make your job ea= sier > Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronim= o > http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D120709&bid=3D263057&dat= =3D121642 > > _______________________________________________ > Libmesh-devel mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > > > |
From: Roy S. <roy...@ic...> - 2006-08-08 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 user-defined 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 G. <fri...@gm...> - 2006-08-08 17:45:48
|
Roy, Thanks for the summary of the discussion! I will do a cvs-update 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 <roy...@ic...> 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 user-defined 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: <ti...@ce...> - 2006-08-09 09:00:14
Attachments:
patch
|
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 >> out-of-mesh-extension 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: Roy S. <roy...@ic...> - 2006-08-09 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 >>> out-of-mesh-extension 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: John P. <pet...@cf...> - 2006-07-24 17:38:16
|
Hi, > Well, the question is still whether you would accept such a port in your > repository. I think the answer is "it depends". I can't speak for the entire libmesh development team, and I can't give a blanket acceptance to code I have not yet seen. > However, I just face a different question: It seems to me that > libMesh has got a memory leak. Please find attached a very simple > test program that creates a grid and initializes an > EquationSystems object on that grid. I can confirm: no memory leak in METHOD=devel CVS libmesh for your sample code after ~ 200 iterations. If you have valgrind you might be able to find it... What OS/compiler are you using? Is this with or without petsc? -J |
From: <ti...@ce...> - 2006-07-25 08:56:34
|
Dear Roy and all, On Mon, 24 Jul 2006, Roy Stogner wrote: > If you're asking whether we'd accept a patch reverting all the > declarations inside for loops, the answer is yes (unless Ben strongly > objects), but if you're asking whether we'd accept a patch which makes > every change that might be required for VC6 compatibility, nobody can > answer that until we've seen it. This is true, of course. So, we will contact you when we can estimate what amount of changes would be necessary. For now, many thanks for your help and cooperation. > Is this a memory leak you're seeing with 0.5.0, or can you duplicate > it on the CVS head too? It was in 0.5.0. I just tested it with the CVS head and found that is does not appear there any more. I will only use the CVS head version from now on for my work. > Actually, this whole discussion ought to be moved to libmesh-devel so > there's a public archive of it.. and I'm not just saying that because > someone typoed my email address and left me out of half of it. ;-) That must have been my fault, sorry for that. Best Regards, Tim |