Question on the ``examples/diffusion/poisson.py''

[Alec Kalinin]

Dear SfePy users,

I am new to the SfePy and first of all I want to say thank you to the developers of this FEM engine. Now I am trying to use SfePy for my tasks and I got several questions.

Consider problem in the example ``examples/diffusion/poisson.py'' (I will use LaTeX notations further for the math). We need to find unknown function $t(x)$ such that
$$
  c \Delta t = 0, x \in \Omega,
$$
$$
  t(x) = T_1, x \in \Gamma_{D1},
$$
$$
  t(x) = T_2, x \in \Gamma_{D2},
$$
$$
  \frac{\partial t(x)}{\partial n} = 0, x \in \Gamma_{N},
$$
where $T_1$ and $T_2$ is known functions, $\Gamma_{D1} = \{x \mid x < 0.00001\}$, $\Gamma_{D2} = \{x \mid x > 0.099999\}$, $\Gamma_{N} = \partial \Omega \setminus (\Gamma_{D1} \cup \Gamma_{D2})$. This is the mixed boundary value problem with both Dirichlet and Neumann boundary conditions.

Could you, please, answer me on following questions:

1. Is it possible to get a flux on the $\Gamma_{D1}$ and $\Gamma_{D2}$ surfaces after the problem has been solved?

2. Suppose we solve pure Dirichlet problem (no Neumann conditions). We need to specify essential boundary conditions on the all surface nodes? What is the best way to select all surface nodes?

3. Suppose the Neumann boundary conditions is not zero. In this case we need to add surface integral with the flux. Do we need to add new term for the flux in the problem description? And this term will be the surface integral but in the problem description we have only tetrahedral domain 3D mesh. What is the best way to select only surface part of the mesh for the surface integration?

Thanks,
Alec Kalinin

[Robert Cimrman]

Hi Alec,

Not yet - we actually never saved a solution on a surface. I will try to fix
this issue now. Then, something like the following code should work (see [1]
for definition of d_surface_flux term, K is a permeability matrix):

1. Define a postprocessing function.

def post_process(out, pb, state, extend=False):
"""
Calculate fluxes.
"""
from sfepy.base.base import Struct
from sfepy.fem import extend_cell_data

ev = pb.evaluate
flux = ev('d_surface_flux.2.Gamma(m.K, t)', mode='el_avg')
flux = extend_cell_data(flux, pb.domain, 'Gamma') # Extend to whole domain
- this does not work now.
out['flux'] = Struct(name='output_data', mode='cell',
data=flux, dofs=None)
return out

2. Add the function to the options:

options = {
'post_process_hook' : 'post_process',

}


> 2. Suppose we solve pure Dirichlet problem (no Neumann conditions). We need
> to specify essential boundary conditions on the all surface nodes? What is
> the best way to select all surface nodes?

You can use the following region definition:

regions = {
'Gamma' : ('nodes of surface', {}),
}

to select all the surface nodes.

> 3. Suppose the Neumann boundary conditions is not zero. In this case we
> need to add surface integral with the flux. Do we need to add new term for
> the flux in the problem description? And this term will be the surface
> integral but in the problem description we have only tetrahedral domain 3D
> mesh. What is the best way to select only surface part of the mesh for the
> surface integration?

Yes, you need to add the flux term "dw_surface_ndot.2.Gamma(m.c, s)" to the
equation, where m.c. is the given flux (a material), and s is the test
function, again check [1], c is the flux vector. As for the surface selection -
see above.

See also example [2] - linear elasticity with surface tractions.

Cheers,
r.

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

[Robert Cimrman]

Hi again,

This is the modified poisson.py showing how to do the post-processing and how
to add flux term to the equations. You need the latest git version for it to
run, though... Let us know whether it works - it's completely untested, so any
bug-report is welcome!

r.

[Alec Kalinin]

Thank you very much, Robert! With the latest git version the given example works and I got the flux.

Alec

[Robert Cimrman]

On 08/21/2012 05:23 PM, Alec Kalinin wrote:
> Thank you very much, Robert! With the latest git version the given example
> works and I got the flux.

Good, hth!
r.