re-hi !

I have cloned the rc/sfepy.git installed it, modified the problem file as advised, and it works : I get a current intensity of 1 and -1 on Gamma_Left and Gamma_Right which is what I expected.

Great job !

From a user point of view, the 'el' or 'el_avg' syntax, is confusing... I am also confused by the term documentation :  the term is described as an integral over a domain, but the computation is actually done over individual elements. An array is returned and we can't get the element each value of the array relates to... d_surface_flux.i1.Gamma_Left(coef.sigma, v) should return a scalar, and if we want the value over a specific element just define a domain that contains it and use d_surfacer_flux over it.

Again, that is a user point of view : I have no idea what the implementation issues are behind, nor issues related to other FE problems !


Le mercredi 29 août 2012 18:25:18 UTC+2, Robert Cimrman a écrit :

I have added the 'el' evaluation mode to d_surface_flux term, and it seems to
work - get the latest sources from [1], and use

flux1 = pb.evaluate('d_surface_flux.i1.Gamma_Left(coef.sigma, v)', mode='el')

If that's ok, I will push it into the main repository.

[1] git clone git://

On 08/29/2012 06:07 PM, Robert Cimrman wrote:
> On 08/29/2012 05:59 PM, Robert Cimrman wrote:
>> 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)
> Ha, you actually did not get this error, right? What's your platform and Python
> version?
>> 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?