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='(set-to-None)')

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/doc-devel/examples/linear_elasticity/its2D_interactive.html [2] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html

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.

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.

[1] https://github.com/rc/sfepy/tree/evaluate-gradient

Hi,

thank you for quick implementation, I did some testing:

1. 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.

2. 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 11:10 PM, Robert Cimrman wrote:

Thanks for the debugging. I have tested it with tensor-product 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 pre-initialized).

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.