Determine closest point in the mesh
Hi there,
Given a point outside a (2D) mesh, is there a way to find the closest point in the mesh? The goal is to implement extrapolation by evaluating at the closest point instead, since regular extrapolation would not produce acceptable results.
Alternately, I can of course implement this myself, in which case I'd need to know/confirm a few things about how SfePy meshes work. Do I recall correctly that SfePy only supports linear elements (i.e. a secondorder element has the same shape as a firstorder elements)? And is there a clean/fast way of determining which elements contain a given mesh node?
My plan would be something like this:
 Given a point, find the n closest nodes in the mesh.
 Determine every element in the mesh which contains at least one of these nodes.
 For each element, determine the closest point in that element (which is easy if elements are linear), and then take the closest of all these points.
If I go down this road, is there also a clean way of evaluating a variable in a given element? Of course I can just use the standard PointsProbeapproach, which shouldn't be much of a speed loss, it just seems silly to do this once the correct element is already known.
Regards, Jonatan
Hi Jonathan,
On 07/03/2015 10:35 AM, Jonatan Lehtonen wrote:
Hi there,
Given a point outside a (2D) mesh, is there a way to find the closest point in the mesh? The goal is to implement extrapolation by evaluating at the closest point instead, since regular extrapolation would not produce acceptable results.
The previous (wrong) implementation of get_ref_coors() did find the closest point on the mesh to a given point outside. The current (also wrong in a general case) one does not (IIRC) have to.
In one of our applications, this function needs to be fixed, and extended for the splinebased FE fields (IGA). So I am going to focus on that in the next week(s) as time permits.
Alternately, I can of course implement this myself, in which case I'd need to know/confirm a few things about how SfePy meshes work. Do I recall correctly that SfePy only supports linear elements (i.e. a secondorder element has the same shape as a firstorder elements)? And is there a
Yes. (In FEM case, anyway...)
clean/fast way of determining which elements contain a given mesh node?
Yes.
cmesh = problem.domain.mesh.cmesh
cmesh.setup_connectivity(0, dim) conn = cmesh.get_conn(0, dim)
where dim is the topological dimension of cells (2 in your case), and the returned conn is a CSRlike object, that for an entity ii of topological dimension 0 (=vertex) contains in conn.indices[conn.offsets[ii]:conn.offsets[ii+1]] the entities of the dimension dim (cells). You can use those functions to get any kind of incidence/connectivity information. (See sfepy/discrete/fem/extmods/mesh.h for explanation what the incidence actually is for various combinations of dimensions.)
My plan would be something like this:
 Given a point, find the n closest nodes in the mesh.
 Determine every element in the mesh which contains at least one of these nodes.
 For each element, determine the closest point in that element (which is easy if elements are linear), and then take the closest of all these points.
You may want to have a look at rtree [1], that Matyas have used to workaround get_ref_coors(). It allows, unlike kdtree that indexes points, to index bounding boxes (of elements). That would solve #285, but I am reluctant to add another dependence.
[1] https://pypi.python.org/pypi/Rtree/
If I go down this road, is there also a clean way of evaluating a variable in a given element? Of course I can just use the standard PointsProbeapproach, which shouldn't be much of a speed loss, it just seems silly to do this once the correct element is already known.
Not much better, but: FieldVariable.evaluate() works on a region, so if that region contains the single cell, and the integral uses custom quadrature point
 the single point you need...
r.
Hi Robert,
Thank you for the swift reply. That information was just what I needed. The rtree does also seem appealing, but unfortunately I have the same issue in that I don't want to add another dependence. I know that the solution I described doesn't work in general (as we've discussed before), but it should definitely be good enough for what I'm working on.
I do have some other questions though, related to higherorder elements. Suppose I have a secondorder triangle mesh. conn.indices only gives me the three corner nodes for each element; how do I find the indices (and coordinates) of the other three nodes for the secondorder element? And more generally, how do I find the indices and coordinates of every node in the mesh?
For example, I have a secondorder FEM problem where problem.domain.mesh.coors has shape (188381, 2), but the solution vector from state.get_parts()['u'] is 751637 elements long (where u is a scalar variable). The first 188381 elements in the solution vector correspond to the corner nodes of the elements, but how do I find the coordinates of the nodes corresponding to the remaining values in the solution vector?
Regards, Jonatan
Hi Jonatan,
On 07/06/2015 05:48 PM, Jonatan Lehtonen wrote:
Hi Robert,
Thank you for the swift reply. That information was just what I needed. The rtree does also seem appealing, but unfortunately I have the same issue in that I don't want to add another dependence. I know that the solution I described doesn't work in general (as we've discussed before), but it should definitely be good enough for what I'm working on.
OK
I do have some other questions though, related to higherorder elements. Suppose I have a secondorder triangle mesh. conn.indices only gives me the three corner nodes for each element; how do I find the indices (and coordinates) of the other three nodes for the secondorder element? And more generally, how do I find the indices and coordinates of every node in the mesh?
Yes, the cmesh functions work only with the underlying (linear) mesh.
To get a full connectivity of a field, use Field.get_econn()  this function works for all kinds of fields (FE with Lagrange basis, FE with hierarchic basis, IG with NURBS basis). The extended connectivity in the FE case has always this structure in each row: vertex DOFs, edge DOFs, face DOFs, bubble DOFs. See field.node_desc attribute after a field is constructed.
The DOFs corresponding to each basis function have nodes with defined coordinates only in the FE Lagrange case, though  for the other fields, the higher order nodes have no coordinates  their DOFs do not correspond to a solution at some point. For the FE field with Lagrange basis, you can use field.get_coor().
What do you need that information for? There might be other functions that do what you need.
For example, I have a secondorder FEM problem where problem.domain.mesh.coors has shape (188381, 2), but the solution vector from state.get_parts()['u'] is 751637 elements long (where u is a scalar variable). The first 188381 elements in the solution vector correspond to the corner nodes of the elements, but how do I find the coordinates of the nodes corresponding to the remaining values in the solution vector?
Yes, the initial items correspond to vertex DOFs, the remaining to the edge ones. There are no face (it is 2D) or bubble DOFs in your case.
r.
Hi Robert,
There's two reasons I asked. Firstly, I wanted to understand the data contained in State (mainly to sate my curiosity); secondly, I've been toying with the idea of determining the inverse of a FEM solution directly using FEM, and in order to try that, I need both the full mesh connectivity and the value of the solution at each mesh node. I might not end up doing that though, since it seems like an errorprone idea even though analytically the FEMsolution should define an invertible function (it's the conformal mapping that I've talked about earlier).
Anyway, thank you again for your help, that should be everything I needed to know. =)
Regards, Jonatan
On Tuesday, July 7, 2015 at 5:07:37 PM UTC+3, Robert Cimrman wrote:
Hi Jonatan,
On 07/06/2015 05:48 PM, Jonatan Lehtonen wrote:
Hi Robert,
Thank you for the swift reply. That information was just what I needed. The rtree does also seem appealing, but unfortunately I have the same issue in that I don't want to add another dependence. I know that the solution I described doesn't work in general (as we've discussed before), but it should definitely be good enough for what I'm working on.
OK
I do have some other questions though, related to higherorder elements. Suppose I have a secondorder triangle mesh. conn.indices only gives me the three corner nodes for each element; how do I find the indices (and coordinates) of the other three nodes for the secondorder element? And more generally, how do I find the indices and coordinates of every node in the mesh?
Yes, the cmesh functions work only with the underlying (linear) mesh.
To get a full connectivity of a field, use Field.get_econn()  this function works for all kinds of fields (FE with Lagrange basis, FE with hierarchic basis, IG with NURBS basis). The extended connectivity in the FE case has always this structure in each row: vertex DOFs, edge DOFs, face DOFs, bubble DOFs. See field.node_desc attribute after a field is constructed.
The DOFs corresponding to each basis function have nodes with defined coordinates only in the FE Lagrange case, though  for the other fields, the higher order nodes have no coordinates  their DOFs do not correspond to a solution at some point. For the FE field with Lagrange basis, you can use field.get_coor().
What do you need that information for? There might be other functions that do what you need.
For example, I have a secondorder FEM problem where problem.domain.mesh.coors has shape (188381, 2), but the solution vector from state.get_parts()['u'] is 751637 elements long (where u is a scalar variable). The first 188381 elements in the solution vector correspond to the corner nodes of the elements, but how do I find the coordinates of the nodes corresponding to the remaining values in the solution vector?
Yes, the initial items correspond to vertex DOFs, the remaining to the edge ones. There are no face (it is 2D) or bubble DOFs in your case.
r.
Hi Jonatan,
On 07/09/2015 04:44 PM, Jonatan Lehtonen wrote:
Hi Robert,
There's two reasons I asked. Firstly, I wanted to understand the data contained in State (mainly to sate my curiosity); secondly, I've been toying with the idea of determining the inverse of a FEM solution directly using FEM, and in order to try that, I need both the full mesh connectivity and the value of the solution at each mesh node. I might not end up doing that though, since it seems like an errorprone idea even though analytically the FEMsolution should define an invertible function (it's the conformal mapping that I've talked about earlier).
I do not understand what you mean by "the inverse of a FEM solution", but I guess you know what you are doing :)
Anyway, thank you again for your help, that should be everything I needed to know. =)
Hth!
BTW. could you try [1]? The new code (strategy='general') might be slower than the current one (strategy='convex'), but it should work for any mesh. (It seems to work with the awful.mesh you sent me.)
r.
[1] https://github.com/rc/sfepy/tree/fixfindrefcoors
Regards, Jonatan
On Tuesday, July 7, 2015 at 5:07:37 PM UTC+3, Robert Cimrman wrote:
Hi Jonatan,
On 07/06/2015 05:48 PM, Jonatan Lehtonen wrote:
Hi Robert,
Thank you for the swift reply. That information was just what I needed. The rtree does also seem appealing, but unfortunately I have the same issue in that I don't want to add another dependence. I know that the solution I described doesn't work in general (as we've discussed before), but it should definitely be good enough for what I'm working on.
OK
I do have some other questions though, related to higherorder elements. Suppose I have a secondorder triangle mesh. conn.indices only gives me the three corner nodes for each element; how do I find the indices (and coordinates) of the other three nodes for the secondorder element? And more generally, how do I find the indices and coordinates of every node in the mesh?
Yes, the cmesh functions work only with the underlying (linear) mesh.
To get a full connectivity of a field, use Field.get_econn()  this function works for all kinds of fields (FE with Lagrange basis, FE with hierarchic basis, IG with NURBS basis). The extended connectivity in the FE case has always this structure in each row: vertex DOFs, edge DOFs, face DOFs, bubble DOFs. See field.node_desc attribute after a field is constructed.
The DOFs corresponding to each basis function have nodes with defined coordinates only in the FE Lagrange case, though  for the other fields, the higher order nodes have no coordinates  their DOFs do not correspond to a solution at some point. For the FE field with Lagrange basis, you can use field.get_coor().
What do you need that information for? There might be other functions that do what you need.
For example, I have a secondorder FEM problem where problem.domain.mesh.coors has shape (188381, 2), but the solution vector from state.get_parts()['u'] is 751637 elements long (where u is a scalar variable). The first 188381 elements in the solution vector correspond to the corner nodes of the elements, but how do I find the coordinates of the nodes corresponding to the remaining values in the solution vector?
Yes, the initial items correspond to vertex DOFs, the remaining to the edge ones. There are no face (it is 2D) or bubble DOFs in your case.
r.
Hi Robert,
I do not understand what you mean by "the inverse of a FEM solution", but I guess you know what you are doing :)
I guess I could've been clearer. =) What I meant was that I have a FEM solution mapping coordinates (x, y) in a given domain to coordinates (u, v) in the unit square. The inverse would then be the mapping from (u, v) to (x, y), which could basically be obtained by just swapping the coordinates and the FEM solution around. It would be much faster than inverting the mapping pointwise (which is what I'm currently doing), but there's some obvious problems with it since there's really no guarantee that the new mesh created by using (u, v) as node coordinates isn't degenerate (bad/overlapping elements et cetera). It's just on my list of things that would be cool to try if I find the time, since it probably won't work but it would be great if it did.
Anyway, thank you again for your help, that should be everything I needed
to know. =)
Hth!
BTW. could you try [1]? The new code (strategy='general') might be slower than the current one (strategy='convex'), but it should work for any mesh. (It seems to work with the awful.mesh you sent me.)
r.
I'd love to help, but I'm kind of swamped right now. I'll have a look if I find the time, but I can't guarantee anything.
Regards, Jonatan
On 07/10/2015 12:05 PM, Jonatan Lehtonen wrote:
Hi Robert,
I do not understand what you mean by "the inverse of a FEM solution", but I guess you know what you are doing :)
I guess I could've been clearer. =) What I meant was that I have a FEM solution mapping coordinates (x, y) in a given domain to coordinates (u, v) in the unit square. The inverse would then be the mapping from (u, v) to (x, y), which could basically be obtained by just swapping the coordinates and the FEM solution around. It would be much faster than inverting the mapping pointwise (which is what I'm currently doing), but there's some obvious problems with it since there's really no guarantee that the new mesh created by using (u, v) as node coordinates isn't degenerate (bad/overlapping elements et cetera). It's just on my list of things that would be cool to try if I find the time, since it probably won't work but it would be great if it did.
OK, thanks for the explanation.
Anyway, thank you again for your help, that should be everything I needed
to know. =)
Hth!
BTW. could you try [1]? The new code (strategy='general') might be slower than the current one (strategy='convex'), but it should work for any mesh. (It seems to work with the awful.mesh you sent me.)
r.
I'd love to help, but I'm kind of swamped right now. I'll have a look if I find the time, but I can't guarantee anything.
No problem :)
r.
participants (2)

Jonatan Lehtonen

Robert Cimrman