Thanks for looking into it.  Your response makes sense, I do see higher flux near the dirchlet boundary.  I've been experimenting with denser meshes but haven't managed to get the flux error down to a satisfactory level yet.

I noticed that gmsh has the ability to generate a 2nd order mesh, which I thought was interesting as it might reduce the gradient errors.  However I haven't had any success importing this mesh into sfepy.  Do you know if it is possible?  I've been generating it in gmsh by clicking on "Set order 2", then export it from gmsh as .mesh format, but then when I load it sfepy spits out errors about 'skipping unknown entity'.  I'm attaching the 2nd order mesh to this post.

Note: this geo file is a bit different to what we were looking at earlier.  I'm now defining the sink and source to both be well inside the overall domain, the source is a sphere and the sink is a cylinder.  Both are much smaller than the outer volume.  I've found this reduces the flux error to around 20%.

Thanks for your help,

Geoff

On Monday, August 4, 2014 5:30:35 PM UTC-4, Robert Cimrman wrote:

Hi Geoff,

the explanation, AFAIK, is below. You will need the latest git sources to try
it (see the attachment), as several things were fixed or updated.

On 08/04/2014 09:11 PM, Geoff Wright wrote:
> Okay sounds good.  Keep me posted!
>
> Thank you,
>
> Geoff
>
>
> On Monday, August 4, 2014 2:21:20 PM UTC-4, Robert Cimrman wrote:
>>
>> On 08/04/2014 05:21 PM, Geoff Wright wrote:
>>> Hi Robert,
>>>
>>> When I try the extend_cell_data on the surface of a cube, I get an
>>> error:
>>> "ValueError: The weights and list don't have the same length."  I
>>> attached
>>> the file.  It uses a mesh that is provided with sfepy.. can you take a
>>> look?

I have fixed extend_cell_data() to work, see e46188f541f9bbe.

>>> This is another example where the flux doesn't make sense.  If I use
>>> only
>>> vertices on one face as the source then the flux values are good.  But
>>> if
>>> you run the file as-is, the source extends a bit into the adjacent
>>> faces,
>>> and the flux results are incorrect.  The electric field values look
>>> okay,
>>> and the voltage field looks fine.  My results are as follows:
>>>
>>> Source flux:  1.09348063065
>>> Sink flux:  -0.90821822135
>>> Block flux:  -0.040501111711
>>> Entire surf:  0.298397270469 0.0181174920822 -0.0124582155367
>>> Source surf area: 3.05555373528
>>> Sink surf area: 1.0

In [1]: 1.09348063065-0.90821822135-0.040501111711
Out[1]: 0.14476129758900003

which is not 0.298397270469 - this is because the Gamma_Left + Gamma_Right +
Blocked < EntireSurf - some faces are missing there.

To fix that, use the 'f' operators instead of 'v':

     'Gamma_Left' : ('r.EntireSurf *f vertices in (x < -0.3)', 'facet'),
     'Gamma_Right' : ('r.EntireSurf *f vertices in (x > 0.4999)', 'facet'),
     'Blocked' : ('r.EntireSurf -f (r.Gamma_Left +f r.Gamma_Right)', 'facet'),

This requires the ce7a6a73166e5d572 commit - there was a bug (or missing
feature) that made the 'r.EntireSurf *f vertices in (x < -0.3)' kind of
selection fail.

>>>
>>>
>>> I'm pretty sure there is an issue with the surface_flux function.  Is
>> this
>>> something you will be able to look at?  I'm trying to decide whether to
>>> continue investing in this or switch to another tool/software.
>>
>> I am actually up to something, stay tuned. (And the surface_flux function
>> is
>> innocent - it comes down to region definitions and gradient computation
>> accuracy).

After the above changes, the sum of partial fluxes equals to the entire surface
flux. However, the FEM is not locally conservative, so there is an inherent
error in flux computation. I modified the file to save also the partial fluxes
in the surface region parts - you may notice a relatively high error (nonzero
values) in the Blocked region near the Dirichlet boundary, where the gradients
are high. To minimalize that, you can:
- use a finer mesh (at least near the Dirichlet boundary)
- use a higher order approximation
Both will soon make the direct solver unusable. Luckily, the Laplacian can be
easily resolved by conjugate gradients, see the commented out solver_0 (PETSc
required for that).

The attached file uses approx_order 3 and gives:

Sum of partial fluxes == total: 0.0378874815651 == 0.0378874815651
2.42833097736e-18
Source flux:  1.05847824223 1.05847824223
Sink flux:  -0.884664882961 -0.884664882961
Block flux:  -0.135925877701 -0.135925877701
Entire surf:  0.0378874815651 0.0378874815651 0.0175191418749 -0.0162855737778
Source surf area: 1.5549799671
Sink surf area: 1.0
Block surf area: 3.4450200329
Entire surf area: 6.0 == 6.0

which is better than it was. Still, the mesh is too coarse.

If the above is not sufficient, another method, that is conservative, would be
needed, for example the discontinuous Galerkin (DG) method. It is not
implemented in SfePy yet, so I have talked a bit with some colleagues about a
possible implementation. You may need to look for a DG or finite volume solver
in the meantime.

Bottom line: the special case that "works" has the gradients aligned in the
x-axis direction (it is essentially a 1D problem), so the error in surface
gradient computation is not present.

That's all, let me know if it makes sense.

r.