Hi David,

On 08/29/2012 05:21 PM, David Libault wrote:

Hi Robert,

Enclosed are a mesh file and a problem file for the 2D electrostatic problem without volume charge of a square of 1x1 on which 1 V is applied on 2 opposite sides (Gamma_Left and Gamma_Right), and the default 0 Von Neumann boundary conditions on the 2 other sides.

The Electric field is evaluated as the gradient of V, and looks good (Ex is small, and Ey = 1 over the entire domain).

Yes, it looks good.

Flux 1 and Flux 2 evaluate the flux intensity of current with conductivity sigma=1 across Gamma_Left and Gamma_Right. They should have opposite sign and absolute value of 1.

It looks that the d_surface_flux is computed on each element of the surface so the returned value is an array, but the values in the array are all 1... Shouldn't they be 1/(element surface) so that the sum of the array gives the total flux over the all Gamma_Left or Gamma_Right surface ?

There is a glitch in using integrals I noticed when preparing the example for Alec: once you use an integral for a volume integration (volume term), you should not use it for a surface term. If you do that, you should get an error - and indeed, there is:

fmf_mulAB_nn(): ERR_BadMatch: (2 2 1) == (3 2 2) * (2 2 1)

printed just before your flux output. This should lead to raising an exception, but for some reason it does not, the exception is raised only when I step over the relevant code piece with a debugger. If I just run the code, it silently continues, and leads to wrong results you got. So this is definitely a bug, but I am not sure why it behaves this way. It used to work.

So to fix you problem: replace flux1 = pb.evaluate('d_surface_flux.i1.Gamma_Left(coef.sigma, v)', mode='el_avg')

with

flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)', mode='el_avg')

and do the same with flux2. The number is the quadrature order... Then you should get the "correct" fluxes. If fact, there is another problem - 'el_avg' mode divides the resulting integrated value in an element by the volume/area/length of that element. That's why you get 1 or -1, and the sums 9 and -9, as there are 9 edges on each of Gamma_Left, Gamma_Right.

So it seems another evaluation mode is needed (let's say called 'el') - one that just integrates a value over an element, without computing the average by dividing by volume. Do you agree?

r.

David.

Le mercredi 29 août 2012 15:38:39 UTC+2, Robert Cimrman a écrit :

Hi David,

On 08/29/2012 02:05 PM, David Libault wrote:

Hi Robert,

From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...

IMHO it's the same term as in my answer to Alec (thread "Question on the ``examples/diffusion/poisson.py''") - look at the post_process() function

- the term is d_surface_flux.
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...

OK, let us know if d_surface_flux works for you - it can take a general permeability tensor.

r.

Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :

Hi Robert,

Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to

On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a

solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!

Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.

But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?

It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.

Best regards, r.

[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html

On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:

Hi Alec,

On 08/25/2012 05:45 PM, Alec Kalinin wrote: > Dear SfePy users, > > Is it possible to evaluate a solution not only in the FEM mesh node,

but

in > any arbitrary point in the domain with the given (x, y, z) coordinates?

Yes, it is possible. Either, you could use a probe as described in

the

Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.

Cheers, r.

> For example, consider Dirichlet problem for Poisson equation. We apply > essential boundary conditions on the surface nodes and after the problem > has been solved we have the solution vector, i.e. vector of values in the > FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is > not FEM mesh node. What is the best way to obtain solution in this point v? > > Sincerely, > Alec Kalinin

[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html