## Re: [Libmesh-devel] MeshFunction without interpolation

 Re: [Libmesh-devel] MeshFunction without interpolation From: John Peterson - 2010-10-15 15:03:41 ```On Fri, Oct 15, 2010 at 8:08 AM, Derek Gaston wrote: > Tim, we recently did something similar... and I think you can do this > much more efficiently another way.  Instead of looking for what > element the point is in then finding the closest node.... just look > for the closest node directly then access the dof of the variable out > of the node and use that to index into the solution vector. > > Writing a loop over all nodes (or all nodes in a sideset) that checks > the distance to a point is a trivial exercise... and MUCH faster than > calling is_in_elem().  The only time I could see the search algorithm > in meshfunction winning is if you started off with a single element > mesh and uniformly refined it 10 or more times (because then the built > in search will be able to traverse the refinement tree and quickly > hone in on the position).  But this is rarely the case. I agree the code is simpler, and will definitely be faster if you only have to check a subset of nodes, but you still have to check *all* the nodes in that list, whereas you can stop once you find the element which contains the point. In a 2D mesh of triangles, for example, there are about half as many elements as there are nodes... Perhaps this doesn't change the theoretical complexity of the algorithm from O(N) but it seems like it has the potential to be faster. Maybe we should look at how our test for whether a point is inside an element is actually implemented... perhaps it could be made faster. In terms of a general library function, we should decide what happens when the search point is equi-distant, to within some tolerance, to two or more nodes of the containing element. Always pick the node of lower ID? -- John ```

 [Libmesh-devel] MeshFunction without interpolation From: Tim Kroeger - 2010-10-15 09:01:20 Attachments: patch ```Dear all, While I'm still on the subset stuff (I had to do a couple of other things in the last weeks, but I'm on it again, and my application does not yet work as expected), I'd suddenly like to be able to disable the interpolation in the MeshFunction. That is, what I want is this: I pass a point to the MeshFunction, and the MeshFunction finds the element in which this point is and then finds the node of this element which is nearest to the point and then returns the dof value of that node. I have implemented this, see attached patch. What I do is, I loop over all dofs, but instead of multiplying the dof value with the value of the shape function, I find out which shape function takes the largest value, and then I return the corresponding dof's value. This does what I want in my situation. However, I would like to hear your comments, in particular: * Would you like this functionality in the library? * Do you like the API? * Will this always do what I mean? (Actually, for me it is not really important that it's definitely the nearest node's value what I get. Rather, it suffices that in some reasonable neighborhood a node, I get that node's value. Possible, the name of the mode should be changed then.) * What should this mode do to the gradient and the hessian? Best Regards, Tim PS: I won't be here on Monday. -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236```
 Re: [Libmesh-devel] MeshFunction without interpolation From: Derek Gaston - 2010-10-15 13:08:20 ```Tim, we recently did something similar... and I think you can do this much more efficiently another way. Instead of looking for what element the point is in then finding the closest node.... just look for the closest node directly then access the dof of the variable out of the node and use that to index into the solution vector. Writing a loop over all nodes (or all nodes in a sideset) that checks the distance to a point is a trivial exercise... and MUCH faster than calling is_in_elem(). The only time I could see the search algorithm in meshfunction winning is if you started off with a single element mesh and uniformly refined it 10 or more times (because then the built in search will be able to traverse the refinement tree and quickly hone in on the position). But this is rarely the case. At the very least... The idea of checking the largest shape function might not be right. There might be some shape functions where the closest dof's shape function might not be the largest. I would just loop through the nodes connected to the element and check their distance to the point. Derek Sent from my iPad On Oct 15, 2010, at 3:01 AM, Tim Kroeger wrote: > Dear all, > > While I'm still on the subset stuff (I had to do a couple of other things in the last weeks, but I'm on it again, and my application does not yet work as expected), I'd suddenly like to be able to disable the interpolation in the MeshFunction. That is, what I want is this: I pass a point to the MeshFunction, and the MeshFunction finds the element in which this point is and then finds the node of this element which is nearest to the point and then returns the dof value of that node. > > I have implemented this, see attached patch. What I do is, I loop over all dofs, but instead of multiplying the dof value with the value of the shape function, I find out which shape function takes the largest value, and then I return the corresponding dof's value. > > This does what I want in my situation. However, I would like to hear your comments, in particular: > > * Would you like this functionality in the library? > > * Do you like the API? > > * Will this always do what I mean? (Actually, for me it is not really important that it's definitely the nearest node's value what I get. Rather, it suffices that in some reasonable neighborhood a node, I get that node's value. Possible, the name of the mode should be changed then.) > > * What should this mode do to the gradient and the hessian? > > Best Regards, > > Tim > > PS: I won't be here on Monday. > > -- > Dr. Tim Kroeger > CeVis -- Center of Complex Systems and Visualization > University of Bremen tim.kroeger@... > Universitaetsallee 29 tim.kroeger@... > D-28359 Bremen Phone +49-421-218-7710 > Germany Fax +49-421-218-4236 > > ------------------------------------------------------------------------------ > Download new Adobe(R) Flash(R) Builder(TM) 4 > The new Adobe(R) Flex(R) 4 and Flash(R) Builder(TM) 4 (formerly > Flex(R) Builder(TM)) enable the development of rich applications that run > across multiple browsers and platforms. Download your free trials today! > http://p.sf.net/sfu/adobe-dev2dev > _______________________________________________ > Libmesh-devel mailing list > Libmesh-devel@... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel ```
 Re: [Libmesh-devel] MeshFunction without interpolation From: John Peterson - 2010-10-15 15:03:41 ```On Fri, Oct 15, 2010 at 8:08 AM, Derek Gaston wrote: > Tim, we recently did something similar... and I think you can do this > much more efficiently another way.  Instead of looking for what > element the point is in then finding the closest node.... just look > for the closest node directly then access the dof of the variable out > of the node and use that to index into the solution vector. > > Writing a loop over all nodes (or all nodes in a sideset) that checks > the distance to a point is a trivial exercise... and MUCH faster than > calling is_in_elem().  The only time I could see the search algorithm > in meshfunction winning is if you started off with a single element > mesh and uniformly refined it 10 or more times (because then the built > in search will be able to traverse the refinement tree and quickly > hone in on the position).  But this is rarely the case. I agree the code is simpler, and will definitely be faster if you only have to check a subset of nodes, but you still have to check *all* the nodes in that list, whereas you can stop once you find the element which contains the point. In a 2D mesh of triangles, for example, there are about half as many elements as there are nodes... Perhaps this doesn't change the theoretical complexity of the algorithm from O(N) but it seems like it has the potential to be faster. Maybe we should look at how our test for whether a point is inside an element is actually implemented... perhaps it could be made faster. In terms of a general library function, we should decide what happens when the search point is equi-distant, to within some tolerance, to two or more nodes of the containing element. Always pick the node of lower ID? -- John ```
 Re: [Libmesh-devel] MeshFunction without interpolation From: Tim Kroeger - 2010-10-15 13:33:08 ```Dear Derek, On Fri, 15 Oct 2010, Derek Gaston wrote: > Tim, we recently did something similar... and I think you can do this > much more efficiently another way. Instead of looking for what > element the point is in then finding the closest node.... just look > for the closest node directly then access the dof of the variable out > of the node and use that to index into the solution vector. I see, sounds good. However, it's not quite the same as what I meant. Consider the following situation (use a monospaced font to display this): A - - - - - - - B - - - C - - - D | | | | | | | | | * | | | | E - - - F - - - G | | | | | | | | | | | | H - - - - - - - I - - - J - - - K At the position indicated by "*", your algorithm would return the value at node E, but mine would return the value at node B because the element in which "*" is contained does not know about node E. On the other hand, for the purpose I need this for, this is actually not really too much of importance, so I might be satisfied with your algorithm. Anyway, if you have already implemented such a thing, why don't you just send it to me? Or, if it is clean enough, why don't you just check it in? > At the very least... The idea of checking the largest shape function > might not be right. There might be some shape functions where the > closest dof's shape function might not be the largest. I would just > loop through the nodes connected to the element and check their > distance to the point. Yes, I certainly should do this. This makes it even slower. I really think, I should use your method. Anyway, I think it would be nice to have such a thing in the library. What would be the correct API? I would propose Mesh::get_nearest_node_to(const Point& point); System::get_value_at_nearest_node_to(const Point& point, std::vector& values, unsigned int var=invalid_int); What do you think? BTW: I will be calling this function extremly often, and in succesive calls the points will often be close to each other, so that I think that the method should do some fast check to see whether the node might be the same as before. One possibility to do that would be that Mesh::get_nearest_node_to() will remember the Point and the node that was returend, as well as the difference between the distance to this node and the distance to the second-nearest node. Then, by the triangle inequality, the next call can be short-circuited if the distance of the new point and the remembered point is smaller than the remembered difference. The only thing is that I don't know how to figure out whether the mesh has changed since the last call, because then of course this fast check cannot be used. Anyway, I can implement this if we agree on the API. Also, if you send me your implementation, this would simplify things for me. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236 ```
 Re: [Libmesh-devel] MeshFunction without interpolation From: Derek Gaston - 2010-10-15 14:11:46 ```On Oct 15, 2010, at 7:32 AM, Tim Kroeger wrote: > I see, sounds good. However, it's not quite the same as what I meant. Consider the following situation (use a monospaced font to display this): > > A - - - - - - - B - - - C - - - D > | | | | > | | | | > | * | | | > | E - - - F - - - G > | | | | > | | | | > | | | | > H - - - - - - - I - - - J - - - K > > At the position indicated by "*", your algorithm would return the value at node E, but mine would return the value at node B because the element in which "*" is contained does not know about node E. I think this is a vote _for_ my algorithm, not against! > On the other hand, for the purpose I need this for, this is actually not really too much of importance, so I might be satisfied with your algorithm. > > Anyway, if you have already implemented such a thing, why don't you just send it to me? Or, if it is clean enough, why don't you just check it in? Unfortunately, my algorithm isn't exactly this... it's just related. > Yes, I certainly should do this. This makes it even slower. I really think, I should use your method. Anyway, I think it would be nice to have such a thing in the library. What would be the correct API? I would propose > > Mesh::get_nearest_node_to(const Point& point); > System::get_value_at_nearest_node_to(const Point& point, std::vector& values, unsigned int var=invalid_int); > > What do you think? > > BTW: I will be calling this function extremly often, and in succesive calls the points will often be close to each other, so that I think that the method should do some fast check to see whether the node might be the same as before. One possibility to do that would be that Mesh::get_nearest_node_to() will remember the Point and the node that was returend, as well as the difference between the distance to this node and the distance to the second-nearest node. Then, by the triangle inequality, the next call can be short-circuited if the distance of the new point and the remembered point is smaller than the remembered difference. The only thing is that I don't know how to figure out whether the mesh has changed since the last call, because then of course this fast check cannot be used. I like the general idea of these two functions... however I'm not sure how useful they would actually be. Like you start to get into here in this last paragraph... if you are doing this kind of thing then you want it to be highly optimized (for instance we've hand tuned these algorithms in our code to perfectly match our situation). Optimizing them for the general case can be tricky. My point is that when you get into this spatial search game and you are depending on it being fast... you are almost always going to need to write something yourself to cater to your particular needs. That said... again, I don't have anything against adding these functions... I just wanted to inject a bit of reality. For us, while they are somewhat related to what we are doing... we wouldn't actually be able to use them. But maybe they could help people get started. As for helping you out with the code... here is some pseudocode: get_nearest_node_to_point(point) { Node * closest_node; Real closest_distance = RealMax; for(cur_node in active_nodes) { Real cur_distance = (cur_node - point).size(); if(cur_distance < closest_distance) { closest_node = cur_node; closest_distance = cur_distance; } } return closest_node; } There are a couple more things to think about here as well. This kind of algorithm will only work with SerialMesh. If you are using ParallelMesh then the closest node might be on another processor. In that case the whole algorithm gets a TON more complicated (you need to look for batches of closest nodes to keep down parallel communication). Also, like you mentioned you may want to somehow "cache" the closest node... however I would caution against doing that _in_ this function. Instead I would do it in your own code on the outside of this function. The reason for that is that the mesh might change and this function won't know it. Finally, after you have a working function I would rewrite it using NodeRange and parallel_reduce() so that it actually runs threaded. Once you get to that point I can help you out more with that. Hope that points you in the right direction! Derek ```
 Re: [Libmesh-devel] MeshFunction without interpolation From: Tim Kroeger - 2010-10-15 14:25:41 ```On Fri, 15 Oct 2010, Derek Gaston wrote: > On Oct 15, 2010, at 7:32 AM, Tim Kroeger wrote: > >> A - - - - - - - B - - - C - - - D >> | | | | >> | | | | >> | * | | | >> | E - - - F - - - G >> | | | | >> | | | | >> | | | | >> H - - - - - - - I - - - J - - - K >> >> At the position indicated by "*", your algorithm would return the value at node E, but mine would return the value at node B because the element in which "*" is contained does not know about node E. > > I think this is a vote _for_ my algorithm, not against! Well, let's say, it's just a different purpose. I anyway only use it for some visualization that I use for debugging (e.g. in the case where the hanging node value is wrong and I want to analyze why that happend). That is, it doesn't matter what the method does exactly, as long as I know what it does. But intuitively, I would find it more sensible to show the (possibly wrong) value of the hanging node only on the right side of the separating line here. > I like the general idea of these two functions... however I'm not > sure how useful they would actually be. Like you start to get into > here in this last paragraph... if you are doing this kind of thing > then you want it to be highly optimized (for instance we've hand > tuned these algorithms in our code to perfectly match our > situation). Optimizing them for the general case can be tricky. > My point is that when you get into this spatial search game and you > are depending on it being fast... you are almost always going to > need to write something yourself to cater to your particular needs. I see... > As for helping you out with the code... here is some pseudocode: > > get_nearest_node_to_point(point) > { > Node * closest_node; > Real closest_distance = RealMax; > for(cur_node in active_nodes) > { > Real cur_distance = (cur_node - point).size(); > if(cur_distance < closest_distance) > { > closest_node = cur_node; > closest_distance = cur_distance; > } > } > return closest_node; > } Thank you! > There are a couple more things to think about here as well. This > kind of algorithm will only work with SerialMesh. If you are using > ParallelMesh then the closest node might be on another processor. > In that case the whole algorithm gets a TON more complicated (you > need to look for batches of closest nodes to keep down parallel > communication). Also, like you mentioned you may want to somehow > "cache" the closest node... however I would caution against doing > that _in_ this function. Instead I would do it in your own code on > the outside of this function. The reason for that is that the mesh > might change and this function won't know it. I see, I didn't think of ParallelMesh. It really seems easier to me now not to put that into the library. I'll wait whether Roy says something about this and then do something on Tuesday, where "something" will most probably mean that I just use your algorithm in my application code and put nothing into the library. Anyway, thank you very much for your help! Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim.kroeger@... Universitaetsallee 29 tim.kroeger@... D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236 ```
 Re: [Libmesh-devel] MeshFunction without interpolation From: Roy Stogner - 2010-10-15 15:19:26 ```On Fri, 15 Oct 2010, Tim Kroeger wrote: > It really seems easier to me now not to put that into the library. > I'll wait whether Roy says something about this I'm afraid I just skimmed the discussion after seeing what looked like a Lagrange-only algorithm. I generally use the DiscontinuityMeasure to highlight possible hanging node issues on C0 elements or the KellyErrorEstimator on C1 elements. --- Roy ```