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

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:
- Is it possible to get a flux on the $\Gamma_{D1}$ and $\Gamma_{D2}$
surfaces after the problem has been solved?
- 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?
- 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

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:
- 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):
- 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
- Add the function to the options:
options = { 'post_process_hook' : 'post_process', }
- 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.
- 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_tractio...

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:
- 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):
- 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
- Add the function to the options:
options = { 'post_process_hook' : 'post_process', }
- 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.
- 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_tractio...

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:
- 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):
- 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
- Add the function to the options:
options = { 'post_process_hook' : 'post_process', }
- 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.
- 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.
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tractio...

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.
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:
- 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):
- 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
- Add the function to the options:
options = { 'post_process_hook' : 'post_process', }
- 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.
- 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.
http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tractio...
participants (2)
-
Alec Kalinin
-
Robert Cimrman