Thank you very much, Robert! With the latest git version the given example works and I got the flux.
Alec
On Tuesday, August 21, 2012 2:50:30 PM UTC+4, Robert Cimrman wrote: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.
On 08/21/2012 12:23 PM, Robert Cimrman wrote:
> Hi Alec,
>
> On 08/21/2012 10:50 AM, Alec Kalinin wrote:
>> 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?
>
> 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
>
>