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