How to evaluate gradient at arbitrary point?
Hello,
after much trying I figured I need your help. I managed to calculate the scalar field variable using sfepy
field = Field.from_args('fu', np.float64, 'scalar', omega, space='H1', poly_space_base='lagrange', approx_order=2) u = FieldVariable('u', 'unknown', field) ...
I can evaluate it at arbitrary points using probes or directly:
p = zeros((500,2)) p[:,0] = linspace(10,10,500) u.evaluate_at(p)
Now I need the gradient, I don't want the average cell value as I am using it to integrate particle path through the field.
This is what I tried:
vecfield = Field.from_args('fgradu', field.dtype, 'vector', field.region, space=field.space, poly_space_base=field.poly_space_base, approx_order=field.approx_order) gradu = FieldVariable('gradu', 'parameter', vecfield, primary_var_name='(settoNone)')
from sfepy.discrete.fem.geometry_element import geometry_data gdata = geometry_data['2_3'] nc = len(gdata.coors) ivn = Integral('ivn', 'custom', gdata.coors, [gdata.volume / nc] * nc) d = u.evaluate('grad', integral=i) gradu.set_data_from_qp(d, i)
Now I can use gradu.evaluate_at(p), however it is not precise. Even (u.evaluate_at(x+dx)u.evaluate_at(x))/dx is better.
I suspect some interpolation is still going on. I thought the special Integral I got from one of the examples was supposed to eliminate the interpolation. Where is the catch?
Thanks for help
Jozef
Hi Jozef,
On 07/15/2015 05:27 PM, Jozef Vesely wrote:
Hello,
after much trying I figured I need your help. I managed to calculate the scalar field variable using sfepy
field = Field.from_args('fu', np.float64, 'scalar', omega, space='H1', poly_space_base='lagrange', approx_order=2) u = FieldVariable('u', 'unknown', field) ...
I can evaluate it at arbitrary points using probes or directly:
p = zeros((500,2)) p[:,0] = linspace(10,10,500) u.evaluate_at(p)
Now I need the gradient, I don't want the average cell value as I am using it to integrate particle path through the field.
Where do you need the gradient? If you need it in the quadrature points, you can use the ev_grad term in 'qp' mode.
This is what I tried:
vecfield = Field.from_args('fgradu', field.dtype, 'vector', field.region, space=field.space, poly_space_base=field.poly_space_base, approx_order=field.approx_order)
gradu = FieldVariable('gradu', 'parameter', vecfield, primary_var_name='(settoNone)')
from sfepy.discrete.fem.geometry_element import geometry_data gdata = geometry_data['2_3'] nc = len(gdata.coors) ivn = Integral('ivn', 'custom', gdata.coors, [gdata.volume / nc] * nc) d = u.evaluate('grad', integral=i) gradu.set_data_from_qp(d, i)
Now I can use gradu.evaluate_at(p), however it is not precise. Even (u.evaluate_at(x+dx)u.evaluate_at(x))/dx is better.
I suspect some interpolation is still going on. I thought the special Integral I got from one of the examples was supposed to eliminate the interpolation. Where is the catch?
Yes, .set_data_from_qp() averages (it calls VolumeField.average_qp_to_vertices()).
Your gradient is linear over each element and not continuous over element boundaries  you get different values in a vertex for each element it is in, right? So some averaging is needed, but you could bypass .set_data_from_qp(), take a simple arithmetic average, have the values in the right order (using the field DOF connectivity) and use Variable.set_data() directly.
Or, the most general way is to use a projection, see [1] (line 261) or [2].
r.
[1] http://sfepy.org/docdevel/examples/linear_elasticity/its2D_interactive.html [2] http://sfepy.org/docdevel/examples/diffusion/time_poisson_interactive.html
On 07/15/2015 05:55 PM, Robert Cimrman wrote:
Hi Jozef,
Hi,
thanks for quick answer,
Where do you need the gradient? If you need it in the quadrature points, you can use the ev_grad term in 'qp' mode. I need it at arbitrary point, unrelated to mesh, much like evaluate_at().
Your gradient is linear over each element and not continuous over element boundaries  you get different values in a vertex for each element it is in, right? So some averaging is needed, but you could bypass .set_data_from_qp(), take a simple arithmetic average, have the values in the right order (using the field DOF connectivity) and use Variable.set_data() directly.
I see that it is difficult to treat gradient discontinuity at element boundaries. Therefore I am not really looking into generating new field representing the gradient (which would need averaging at vertices).
I suppose the proposed method using projection does just that, but more in a more clever way than simple average.
Rather than that I would like to only evaluate the gradient of my field at given point. Something along the lines:
To evaluate the the variable at given point you first evaluate basis functions at given point and then take their linear combination depending on the values on vertices (quarature points) of given element. Isn't it possible to evaluate gradients of basis functions at given point weight them by vertex values of original field and get the result?
Arbitrary point has little chance being exactly on the vertex (or facet) where gradient would be discontinuous.
J.
On 07/15/2015 07:07 PM, Jozef Vesely wrote:
On 07/15/2015 05:55 PM, Robert Cimrman wrote:
Hi Jozef,
Hi,
thanks for quick answer,
Where do you need the gradient? If you need it in the quadrature points, you can use the ev_grad term in 'qp' mode. I need it at arbitrary point, unrelated to mesh, much like evaluate_at().
OK, I missed the subject :)
Your gradient is linear over each element and not continuous over element boundaries  you get different values in a vertex for each element it is in, right? So some averaging is needed, but you could bypass .set_data_from_qp(), take a simple arithmetic average, have the values in the right order (using the field DOF connectivity) and use Variable.set_data() directly.
I see that it is difficult to treat gradient discontinuity at element boundaries. Therefore I am not really looking into generating new field representing the gradient (which would need averaging at vertices).
I suppose the proposed method using projection does just that, but more in a more clever way than simple average.
Yes.
Rather than that I would like to only evaluate the gradient of my field at given point. Something along the lines:
To evaluate the the variable at given point you first evaluate basis functions at given point and then take their linear combination depending on the values on vertices (quarature points) of given element. Isn't it possible to evaluate gradients of basis functions at given point weight them by vertex values of original field and get the result?
No, but I will add that functionality to Variable.evaluate_at() tomorrow  it should be pretty easy. I will let you know when it's ready.
Arbitrary point has little chance being exactly on the vertex (or facet) where gradient would be discontinuous.
Yes, and on a good mesh the gradients should not differ too much between neighbouring elements.
r.
Hi Jozef,
On 07/15/2015 10:14 PM, Robert Cimrman wrote:
On 07/15/2015 07:07 PM, Jozef Vesely wrote:
To evaluate the the variable at given point you first evaluate basis functions at given point and then take their linear combination depending on the values on vertices (quarature points) of given element. Isn't it possible to evaluate gradients of basis functions at given point weight them by vertex values of original field and get the result?
No, but I will add that functionality to Variable.evaluate_at() tomorrow  it should be pretty easy. I will let you know when it's ready.
you can try [1] (use u.evaluate_at(p, mode='grad') or probe call with the same argument). Let me know if that works for you and I will then merge it into master. It passes the updated test test_invariance_qp() in in tests/test_mesh_interp.py.
r.
On 07/16/2015 02:47 PM, Robert Cimrman wrote:
you can try [1] (use u.evaluate_at(p, mode='grad') or probe call with the same argument). Let me know if that works for you and I will then merge it into master. It passes the updated test test_invariance_qp() in in tests/test_mesh_interp.py.
Hi,
thank you for quick implementation, I did some testing:
shape of data returned from evaluate_at changed. There is extra second axis (even without mode='grad'), this would be compatibility problem. However I might have had old version and the change happened already earlier.
check attached test.py: For first order elements finite difference derivative and evaluate_at(...,mode='grad') match. For higher order elements there are discrepancies. I suppose this should not happen.
Again thanks for good work and help.
J.
On 07/16/2015 05:14 PM, Jozef Vesely wrote:
On 07/16/2015 02:47 PM, Robert Cimrman wrote:
you can try [1] (use u.evaluate_at(p, mode='grad') or probe call with the same argument). Let me know if that works for you and I will then merge it into master. It passes the updated test test_invariance_qp() in in tests/test_mesh_interp.py.
Hi,
thank you for quick implementation, I did some testing:
 shape of data returned from evaluate_at changed. There is extra second axis (even without mode='grad'), this would be compatibility problem. However I might have had old version and the change happened already earlier.
I have changed it now in order to have the same number of dimensions in the output in both modes. But on second thought it might not be a good idea. What do you think? I am now inclined to revert that change.
 check attached test.py: For first order elements finite difference derivative and evaluate_at(...,mode='grad') match. For higher order elements there are discrepancies. I suppose this should not happen.
Thanks for the debugging. I have tested it with tensorproduct element mesh (meshes/2d/square_quad.mesh), where your test worked. The problem was in wrong initialization of the basis gradient in simplex element case (not manifesting elsewhere, as the function had been called with arguments preinitialized).
Your test file works now.
r.
On 07/16/2015 11:10 PM, Robert Cimrman wrote:
Thanks for the debugging. I have tested it with tensorproduct element mesh (meshes/2d/square_quad.mesh), where your test worked. The problem was in wrong initialization of the basis gradient in simplex element case (not manifesting elsewhere, as the function had been called with arguments preinitialized).
Your test file works now.
Thank you very much, my magnetic lens is now finally focusing :) When I polish the code and verify cylindrical transformations I'll post the code.
Thanks for help, bye.
J.
On 07/17/2015 12:23 PM, Jozef Vesely wrote:
On 07/16/2015 11:10 PM, Robert Cimrman wrote:
Thanks for the debugging. I have tested it with tensorproduct element mesh (meshes/2d/square_quad.mesh), where your test worked. The problem was in wrong initialization of the basis gradient in simplex element case (not manifesting elsewhere, as the function had been called with arguments preinitialized).
Your test file works now.
Thank you very much, my magnetic lens is now finally focusing :)
Glad to hear that :)
When I polish the code and verify cylindrical transformations I'll post the code.
Great, nice examples are welcome!
Thanks for help, bye.
Hth. I will revert the change of shape and merge it then.
r.
participants (2)

Jozef Vesely

Robert Cimrman