Evaluate a solution in the arbitrary point in the domain
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node, but in any arbitrary point in the domain with the given (x, y, z) coordinates?
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node, but in any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in the Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to probe a solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
Sincerely, Alexander
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node, but in any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in the Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
On 08/26/2012 12:32 PM, Alec Kalinin wrote:
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to probe a solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node, but in any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in the Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
Thank you, Robert.
Alec Kalinin.
On Monday, August 27, 2012 12:55:49 AM UTC+4, Robert Cimrman wrote:
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
Hi Robert,
From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...
David.
Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node,
in
any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in the Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the
has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z)
is
not FEM mesh node. What is the best way to obtain solution in this
but problem that point
v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
Hi David,
On 08/29/2012 02:05 PM, David Libault wrote:
Hi Robert,
From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...
IMHO it's the same term as in my answer to Alec (thread "Question on the ``examples/diffusion/poisson.py''") - look at the post_process() function - the term is d_surface_flux.
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...
OK, let us know if d_surface_flux works for you - it can take a general permeability tensor.
r.
Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node,
in
any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in the Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the
has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z)
is
not FEM mesh node. What is the best way to obtain solution in this
but problem that point
v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
Hi Robert,
Ok, sorry, I was confused by the naming : for "d_surface_flux" I was expecting a vector field as the parameter, not a scalar over which the gradient is computed before taking the actual flux. But the documentation is correct.
David.
Le mercredi 29 août 2012 15:38:39 UTC+2, Robert Cimrman a écrit :
Hi David,
On 08/29/2012 02:05 PM, David Libault wrote:
Hi Robert,
From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...
IMHO it's the same term as in my answer to Alec (thread "Question on the ``examples/diffusion/poisson.py''") - look at the post_process() function
- the term is d_surface_flux.
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...
OK, let us know if d_surface_flux works for you - it can take a general permeability tensor.
r.
Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node,
but
in
any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in
the
Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
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).
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 ?
David.
Le mercredi 29 août 2012 15:38:39 UTC+2, Robert Cimrman a écrit :
Hi David,
On 08/29/2012 02:05 PM, David Libault wrote:
Hi Robert,
From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...
IMHO it's the same term as in my answer to Alec (thread "Question on the ``examples/diffusion/poisson.py''") - look at the post_process() function
- the term is d_surface_flux.
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...
OK, let us know if d_surface_flux works for you - it can take a general permeability tensor.
r.
Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote:
Dear SfePy users,
Is it possible to evaluate a solution not only in the FEM mesh node,
but
in
any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in
the
Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
For example, consider Dirichlet problem for Poisson equation. We apply essential boundary conditions on the surface nodes and after the problem has been solved we have the solution vector, i.e. vector of values in the FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is not FEM mesh node. What is the best way to obtain solution in this point v?
Sincerely, Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
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)
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?
r.
David.
Le mercredi 29 août 2012 15:38:39 UTC+2, Robert Cimrman a écrit :
Hi David,
On 08/29/2012 02:05 PM, David Libault wrote:
Hi Robert,
From the poisson.py I would like to probe the flux of grad(T) over a surface (Gamma_Left or Gamma_right), but I am not sure which term to use...
IMHO it's the same term as in my answer to Alec (thread "Question on the ``examples/diffusion/poisson.py''") - look at the post_process() function
- the term is d_surface_flux.
I would use it in electrostatics (temperature is replaced by voltage), to compute the resistance of the volume imposing a voltage of 1 across two surfaces with dirichlet conditions and computing the electrical current. For a ohmic conductor, j = \sigma E = - \sigma grad(V), \sigma being the conductivity. j being the current density, the current intensity thru a surface is the flux of j thru that surface. For the case of homogenous \sigma, the term asked for above would do...
OK, let us know if d_surface_flux works for you - it can take a general permeability tensor.
r.
Le dimanche 26 août 2012 22:55:49 UTC+2, Robert Cimrman a écrit :
Hi Robert,
Did you mean "linear_elastic_probes.html" instead of "linear_elastic_tractions.html" example? I found the "linear_elastic_probes.html" very useful example for my purposes to
On 08/26/2012 12:32 PM, Alec Kalinin wrote: probe a
solution in the given (x, y, z) points. Also the documentation "src/sfepy/fem/probes.html" gives all necessary information to help me implement what I want to do. Thank you!
Sorry, I cut&pasted a wrong url, the correct one is [1]. But you found another one that solves the problem.
But, despite this, could you tell me more about low-level way to evaluate a variable in the given (x, y, z) point?
It's exactly how the probes do that: the key function is variable.evaluate_at() [2], where variable is an unknown or parameter variable. It takes just one compulsory parameter - the coordinates of points in which you wish to evaluate the variable. You can get the variables of a problem by problem.get_variables(), where problem is the second argument of the post_process_hook function.
Best regards, r.
[1] http://sfepy.org/doc-devel/primer.html#probing [2] http://sfepy.org/doc-devel/src/sfepy/fem/variables.html, http://sfepy.org/doc-devel/src/sfepy/fem/fields.html
On Sunday, August 26, 2012 12:47:50 AM UTC+4, Robert Cimrman wrote:
Hi Alec,
On 08/25/2012 05:45 PM, Alec Kalinin wrote: > Dear SfePy users, > > Is it possible to evaluate a solution not only in the FEM mesh node,
but
in > any arbitrary point in the domain with the given (x, y, z) coordinates?
Yes, it is possible. Either, you could use a probe as described in
the
Primer [1] - the available probes are described in [2]. Or, you could directly evaluate a variable in given points - this is a bit low-level operation, but I could provide you instructions, if the probes are not enough for you.
Cheers, r.
> For example, consider Dirichlet problem for Poisson equation. We apply > essential boundary conditions on the surface nodes and after the problem > has been solved we have the solution vector, i.e. vector of values in the > FEM mesh nodes. But I want to know the solution in point v(x, y, z) that is > not FEM mesh node. What is the best way to obtain solution in this point v? > > Sincerely, > Alec Kalinin
[1] doc-devel/examples/linear_elasticity/linear_elastic_tractions.html [2] http://sfepy.org/doc-devel/src/sfepy/fem/probes.html
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?
re-hi!
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.
Thanks, r. [1] git clone git://github.com/rc/sfepy.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?
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 !
David.
Le mercredi 29 août 2012 18:25:18 UTC+2, Robert Cimrman a écrit :
re-hi!
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.
Thanks, r. [1] git clone git://github.com/rc/sfepy.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?
On 08/30/2012 10:32 AM, David Libault wrote:
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 !
Hth!
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.
It's good to hear a feedback, thanks!
The reason for 'el_avg' is mainly to be able to visualize data defined inside elements using Mayavi, so you get a value for each element in a region. The same with 'el'. So now you can add it to the 'out' dict, and see how each border element contributes to the total flux.
If you just want a value of the integral, do:
pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)', mode='eval')
- this returns a single value (and 'eval' mode is the default - you can omit the mode argument).
Also, you can get to the elements/edges/faces etc. the values belong to, as you have the ProblemDefinition instance available - try
print pb.domain.regions['Gamma_Left'] # Edge ids. print pb.domain.regions['Gamma_Left'].edges # Edge nodes. print pb.domain.ed.facets[pb.domain.regions['Gamma_Left'].edges[0]] # (group, element, local edge number). print pb.domain.ed.indices[pb.domain.regions['Gamma_Left'].edges[0]]
etc.
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 !
It's just more general in this way - you can get what you request for. The main problem I see is that it is not properly documented in the tutorial/users' guide. My bad...
r.
David.
Le mercredi 29 août 2012 18:25:18 UTC+2, Robert Cimrman a écrit :
re-hi!
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.
Thanks, r. [1] git clone git://github.com/rc/sfepy.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?
Robert,
Understood. I was not using the correct "mode" for what I was expecting. Now the user interface looks nice.
David.
2012/8/30 Robert Cimrman <cimr...@ntc.zcu.cz>:
On 08/30/2012 10:32 AM, David Libault wrote:
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 !
Hth!
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.
It's good to hear a feedback, thanks!
The reason for 'el_avg' is mainly to be able to visualize data defined inside elements using Mayavi, so you get a value for each element in a region. The same with 'el'. So now you can add it to the 'out' dict, and see how each border element contributes to the total flux.
If you just want a value of the integral, do:
pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)', mode='eval')
- this returns a single value (and 'eval' mode is the default - you can omit the mode argument).
Also, you can get to the elements/edges/faces etc. the values belong to, as you have the ProblemDefinition instance available - try
print pb.domain.regions['Gamma_Left'] # Edge ids. print pb.domain.regions['Gamma_Left'].edges # Edge nodes. print pb.domain.ed.facets[pb.domain.regions['Gamma_Left'].edges[0]] # (group, element, local edge number). print pb.domain.ed.indices[pb.domain.regions['Gamma_Left'].edges[0]]
etc.
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 !
It's just more general in this way - you can get what you request for. The main problem I see is that it is not properly documented in the tutorial/users' guide. My bad...
r.
David.
Le mercredi 29 août 2012 18:25:18 UTC+2, Robert Cimrman a écrit :
re-hi!
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.
Thanks, r. [1] git clone git://github.com/rc/sfepy.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?
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
On 08/30/2012 11:08 AM, David Libault wrote:
Robert,
Understood. I was not using the correct "mode" for what I was expecting. Now the user interface looks nice.
David.
Nevertheless, it should be documented (-> Issue 196).
So this is my current testing function, that shows how to save the fluxes (they are displayed over the triangles that contain the boundary edges...):
def post_process(out, pb, state, extend=False): from sfepy.base.base import Struct from sfepy.fem import extend_cell_data
electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg')
out['electric_field'] = Struct(name='Electric field', mode='cell',
data=electric_field, dofs=None)
flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)',
mode='el')
print flux1
print flux1.sum()
flux1 = extend_cell_data(flux1, pb.domain, 'Gamma_Left', val=0.0,
is_surface=True)
out['flux1'] = Struct(name='output_data', mode='cell',
data=flux1, dofs=None)
flux2 = pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)',
mode='el')
print flux2
print flux2.sum()
flux2 = extend_cell_data(flux2, pb.domain, 'Gamma_Right', val=0.0,
is_surface=True)
out['flux2'] = Struct(name='output_data', mode='cell',
data=flux2, dofs=None)
print pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)')
return outr.
Robert,
I have rerun the "old" electro static poisson problem again with mac OsX port of sfepy, and the flux evaluation doesn't seem to work...
sfepy: left over: ['verbose', '__builtins__', '__file__', '__name__', 'nm', '_filename', '__package__', 'post_process', '__doc__']
sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (essai_2D.mesh)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy: warning: region Gamma_Right of cell kind has empty group indices!
sfepy: Gamma_Right
sfepy: Omega
sfepy: warning: region Gamma_Left of cell kind has empty group indices!
sfepy: Gamma_Left
sfepy: ...done in 0.01 s
sfepy: equation "Voltage":
sfepy: dw_laplace.i1.Omega( coef.val, s, v ) = 0
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: using solvers:
ts: no ts
nls: newton
ls: lssfepy: updating variables...
sfepy: ...done
sfepy: matrix shape: (144, 144)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 922 (4.45e-02% fill)
sfepy: updating materials...
sfepy: coef
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 0.000000e+00 (rel: 0.000000e+00)
sfepy: equation "tmp":
sfepy: -ev_grad.2.Omega( v )
sfepy: updating materials...
sfepy: ...done in 0.00 s
sfepy: equation "tmp":
sfepy: d_surface_flux.2.Gamma_Left(coef.sigma, v)
sfepy: updating materials...
sfepy: coef
Traceback (most recent call last):
File "/opt/local/bin/simple.py-2.7", line 155, in <module>
main()File "/opt/local/bin/simple.py-2.7", line 152, in main
app()File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/applications/application.py", line 29, in call_basic
return self.call(**kwargs)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/applications/pde_solver_app.py", line 213, in call
nls_status=nls_status)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/solvers/ts_solvers.py", line 37, in __call__
file_per_var=None)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 715, in save_state
out = post_process_hook(out, self, state, extend=extend)File "poisson_electrostatic.py", line 157, in post_process
flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)', mode='el')
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1179, in evaluate
verbose=verbose, **kwargs)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1130, in create_evaluable
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/equations.py", line 322, in time_update_materials
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 70, in time_update
mat.time_update(ts, equations, mode=mode, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 345, in time_update
self.update_data(key, ts, equations, term, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 247, in update_data
coors = qps.get_merged_values()File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/mappings.py", line 26, in get_merged_values
qps = nm.concatenate([self.values[ig] for ig in self.igs], axis=0)ValueError: need at least one array to concatenate
simple.py-2.7 --version returns simple.py-2.7 2013.4 on my system which should integrate the modifications you made two years ago.
Is your "testing function, that shows how to save the fluxes" still working ?
(I had to change 'nodes' to 'vertices' and 'elements' to 'cells' to have it work again)
Best regards,
David.
Le jeudi 30 août 2012 12:04:08 UTC+2, Robert Cimrman a écrit :
On 08/30/2012 11:08 AM, David Libault wrote:
Robert,
Understood. I was not using the correct "mode" for what I was expecting. Now the user interface looks nice.
David.
Nevertheless, it should be documented (-> Issue 196).
So this is my current testing function, that shows how to save the fluxes (they are displayed over the triangles that contain the boundary edges...):
def post_process(out, pb, state, extend=False): from sfepy.base.base import Struct from sfepy.fem import extend_cell_data
electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg') out['electric_field'] = Struct(name='Electric field', mode='cell', data=electric_field, dofs=None) flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)', mode='el') print flux1 print flux1.sum() flux1 = extend_cell_data(flux1, pb.domain, 'Gamma_Left', val=0.0, is_surface=True) out['flux1'] = Struct(name='output_data', mode='cell', data=flux1, dofs=None) flux2 = pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)', mode='el') print flux2 print flux2.sum() flux2 = extend_cell_data(flux2, pb.domain, 'Gamma_Right', val=0.0, is_surface=True) out['flux2'] = Struct(name='output_data', mode='cell', data=flux2, dofs=None) print pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)') return outr.
Hi David,
On 04/02/2014 07:09 PM, David Libault wrote:
Robert,
I have rerun the "old" electro static poisson problem again with mac OsX port of sfepy, and the flux evaluation doesn't seem to work...
<snip>
sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (essai_2D.mesh)... sfepy: equation "tmp":
sfepy: d_surface_flux.2.Gamma_Left(coef.sigma, v)
sfepy: updating materials...
sfepy: coef
Traceback (most recent call last):
File "/opt/local/bin/simple.py-2.7", line 155, in <module>
main()File "/opt/local/bin/simple.py-2.7", line 152, in main
app()File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/applications/application.py", line 29, in call_basic
return self.call(**kwargs)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/applications/pde_solver_app.py", line 213, in call
nls_status=nls_status)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/solvers/ts_solvers.py", line 37, in __call__
file_per_var=None)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 715, in save_state
out = post_process_hook(out, self, state, extend=extend)File "poisson_electrostatic.py", line 157, in post_process
flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)',mode='el')
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1179, in evaluate
verbose=verbose, **kwargs)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1130, in create_evaluable
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/equations.py", line 322, in time_update_materials
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 70, in time_update
mat.time_update(ts, equations, mode=mode, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 345, in time_update
self.update_data(key, ts, equations, term, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 247, in update_data
coors = qps.get_merged_values()File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sfepy/fem/mappings.py", line 26, in get_merged_values
qps = nm.concatenate([self.values[ig] for ig in self.igs], axis=0)ValueError: need at least one array to concatenate
simple.py-2.7 --version returns simple.py-2.7 2013.4 on my system which should integrate the modifications you made two years ago.
Is your "testing function, that shows how to save the fluxes" still working ?
(I had to change 'nodes' to 'vertices' and 'elements' to 'cells' to have it work again)
Two things are needed to make it work:
- set the kind of boundary regions to 'facet':
region_03 = { 'name' : 'Gamma_Left', 'select' : 'vertices of group 1', 'kind' : 'facet', }
region_4 = { 'name' : 'Gamma_Right', 'select' : 'vertices of group 2', 'kind' : 'facet', }
- get the latest git sources, as I have just fixed a bug in extend_cell_data() that occurred only for is_surface=True.
Then I can run your old script - let me know it the above works for you.
Best regards. r.
Best regards,
David.
Le jeudi 30 août 2012 12:04:08 UTC+2, Robert Cimrman a écrit :
On 08/30/2012 11:08 AM, David Libault wrote:
Robert,
Understood. I was not using the correct "mode" for what I was expecting. Now the user interface looks nice.
David.
Nevertheless, it should be documented (-> Issue 196).
So this is my current testing function, that shows how to save the fluxes (they are displayed over the triangles that contain the boundary edges...):
def post_process(out, pb, state, extend=False): from sfepy.base.base import Struct from sfepy.fem import extend_cell_data
electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg') out['electric_field'] = Struct(name='Electric field', mode='cell', data=electric_field, dofs=None) flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)', mode='el') print flux1 print flux1.sum() flux1 = extend_cell_data(flux1, pb.domain, 'Gamma_Left', val=0.0, is_surface=True) out['flux1'] = Struct(name='output_data', mode='cell', data=flux1, dofs=None) flux2 = pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)', mode='el') print flux2 print flux2.sum() flux2 = extend_cell_data(flux2, pb.domain, 'Gamma_Right', val=0.0, is_surface=True) out['flux2'] = Struct(name='output_data', mode='cell', data=flux2, dofs=None) print pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)') return outr.
Hi Robert,
Thank you for your kind answer. Added 'facet' to the regions did fix the error problem, but the values of the fluxes are not as expected :
flux1.sum() : -0.817686412623
flux2.sum() : 0.904600038826
It should be -1 and 1...
flux1 has 8 "terms" and flux2 has 9 "terms"...
Is it a problem in the mesh ? Or is it a "side effet" of to the "extend_cell_data()" issue ?
What values do you get ?
Best regards,
David.
2014-04-02 22:34 GMT+02:00 Robert Cimrman <cimr...@ntc.zcu.cz>:
Hi David,
On 04/02/2014 07:09 PM, David Libault wrote:
Robert,
I have rerun the "old" electro static poisson problem again with mac OsX port of sfepy, and the flux evaluation doesn't seem to work...
<snip>
sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (essai_2D.mesh)...
sfepy: equation "tmp":
sfepy: d_surface_flux.2.Gamma_Left(coef.sigma, v)
sfepy: updating materials...
sfepy: coef
Traceback (most recent call last):
File "/opt/local/bin/simple.py-2.7", line 155, in <module>
main()File "/opt/local/bin/simple.py-2.7", line 152, in main
app()File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/ applications/application.py", line 29, in call_basic
return self.call(**kwargs)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/ applications/pde_solver_app.py", line 213, in call
nls_status=nls_status)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/solvers/ts_solvers.py", line 37, in __call__
file_per_var=None)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 715, in save_state
out = post_process_hook(out, self, state, extend=extend)File "poisson_electrostatic.py", line 157, in post_process
flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)',mode='el')
File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1179, in evaluate
verbose=verbose, **kwargs)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/problemDef.py", line 1130, in create_evaluable
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/equations.py", line 322, in time_update_materials
verbose=verbose)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 70, in time_update
mat.time_update(ts, equations, mode=mode, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 345, in time_update
self.update_data(key, ts, equations, term, problem=problem)File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/materials.py", line 247, in update_data
coors = qps.get_merged_values()File "/opt/local/Library/Frameworks/Python.framework/ Versions/2.7/lib/python2.7/site-packages/sfepy/fem/mappings.py", line 26, in get_merged_values
qps = nm.concatenate([self.values[ig] for ig in self.igs], axis=0)ValueError: need at least one array to concatenate
simple.py-2.7 --version returns simple.py-2.7 2013.4 on my system which should integrate the modifications you made two years ago.
Is your "testing function, that shows how to save the fluxes" still working ?
(I had to change 'nodes' to 'vertices' and 'elements' to 'cells' to have it work again)
Two things are needed to make it work:
- set the kind of boundary regions to 'facet':
region_03 = { 'name' : 'Gamma_Left', 'select' : 'vertices of group 1', 'kind' : 'facet', }
region_4 = { 'name' : 'Gamma_Right', 'select' : 'vertices of group 2', 'kind' : 'facet', }
- get the latest git sources, as I have just fixed a bug in extend_cell_data() that occurred only for is_surface=True.
Then I can run your old script - let me know it the above works for you.
Best regards. r.
Best regards,
David.
Le jeudi 30 août 2012 12:04:08 UTC+2, Robert Cimrman a écrit :
On 08/30/2012 11:08 AM, David Libault wrote:
Robert,
Understood. I was not using the correct "mode" for what I was expecting. Now the user interface looks nice.
David.
Nevertheless, it should be documented (-> Issue 196).
So this is my current testing function, that shows how to save the fluxes (they are displayed over the triangles that contain the boundary edges...):
def post_process(out, pb, state, extend=False): from sfepy.base.base import Struct from sfepy.fem import extend_cell_data
electric_field = pb.evaluate('ev_grad.i1.Omega( v )',mode='el_avg') out['electric_field'] = Struct(name='Electric field', mode='cell', data=electric_field, dofs=None)
flux1 = pb.evaluate('d_surface_flux.2.Gamma_Left(coef.sigma, v)', mode='el') print flux1 print flux1.sum() flux1 = extend_cell_data(flux1, pb.domain, 'Gamma_Left', val=0.0, is_surface=True) out['flux1'] = Struct(name='output_data', mode='cell', data=flux1, dofs=None) flux2 = pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)', mode='el') print flux2 print flux2.sum() flux2 = extend_cell_data(flux2, pb.domain, 'Gamma_Right', val=0.0, is_surface=True) out['flux2'] = Struct(name='output_data', mode='cell', data=flux2, dofs=None) print pb.evaluate('d_surface_flux.2.Gamma_Right(coef.sigma, v)') return outr.
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sfepy-devel...@googlegroups.com. To post to this group, send email to sfepy...@googlegroups.com. Visit this group at http://groups.google.com/group/sfepy-devel.
On 04/02/2014 11:39 PM, David Libault wrote:
Hi Robert,
Thank you for your kind answer. Added 'facet' to the regions did fix the error problem, but the values of the fluxes are not as expected :
flux1.sum() : -0.817686412623
flux2.sum() : 0.904600038826
It should be -1 and 1...
flux1 has 8 "terms" and flux2 has 9 "terms"...
Is it a problem in the mesh ? Or is it a "side effet" of to the "extend_cell_data()" issue ?
What values do you get ?
-1 for Gamma_Left, 1 for Gamma_Right - with the essai_2D.mesh
r.
Yeap. That is a problem on my mesh... Sorry. So yes, it is working with the mesh sent to this list together with the problem file. I am not using the extend_cell_data() for the moment.
Thanks,
David.
2014-04-02 23:51 GMT+02:00 Robert Cimrman <cimr...@ntc.zcu.cz>:
On 04/02/2014 11:39 PM, David Libault wrote:
Hi Robert,
Thank you for your kind answer. Added 'facet' to the regions did fix the error problem, but the values of the fluxes are not as expected :
flux1.sum() : -0.817686412623
flux2.sum() : 0.904600038826
It should be -1 and 1...
flux1 has 8 "terms" and flux2 has 9 "terms"...
Is it a problem in the mesh ? Or is it a "side effet" of to the "extend_cell_data()" issue ?
What values do you get ?
-1 for Gamma_Left, 1 for Gamma_Right - with the essai_2D.mesh
r.
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sfepy-devel...@googlegroups.com. To post to this group, send email to sfepy...@googlegroups.com. Visit this group at http://groups.google.com/group/sfepy-devel.
On 04/03/2014 04:37 AM, David Libault wrote:
Yeap. That is a problem on my mesh... Sorry. So yes, it is working with the mesh sent to this list together with the problem file. I am not using the extend_cell_data() for the moment.
Thanks,
David.
Hth!
r.
2014-04-02 23:51 GMT+02:00 Robert Cimrman <cimr...@ntc.zcu.cz>:
On 04/02/2014 11:39 PM, David Libault wrote:
Hi Robert,
Thank you for your kind answer. Added 'facet' to the regions did fix the error problem, but the values of the fluxes are not as expected :
flux1.sum() : -0.817686412623
flux2.sum() : 0.904600038826
It should be -1 and 1...
flux1 has 8 "terms" and flux2 has 9 "terms"...
Is it a problem in the mesh ? Or is it a "side effet" of to the "extend_cell_data()" issue ?
What values do you get ?
-1 for Gamma_Left, 1 for Gamma_Right - with the essai_2D.mesh
r.
Hi Robert,
No I did not get this error, and I was going to tell you... My platform is Mac OSX + MacPorts. I use MacPorts python 2.7 and sfepy installed manually from git.
Le mercredi 29 août 2012 18:07:50 UTC+2, Robert Cimrman a écrit :
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?
Hi David,
thanks for the info. Anyway, I would recommend you not to use a single named integral for both the volume and surface integration, as I have written below, at least for the moment.
r.
On 08/30/2012 10:01 AM, David Libault wrote:
Hi Robert,
No I did not get this error, and I was going to tell you... My platform is Mac OSX + MacPorts. I use MacPorts python 2.7 and sfepy installed manually from git.
Le mercredi 29 août 2012 18:07:50 UTC+2, Robert Cimrman a écrit :
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?
Ok, I tried using .2 instead of .i1 and it also works.
David.
2012/8/30 Robert Cimrman <cimr...@ntc.zcu.cz>:
Hi David,
thanks for the info. Anyway, I would recommend you not to use a single named integral for both the volume and surface integration, as I have written below, at least for the moment.
r.
On 08/30/2012 10:01 AM, David Libault wrote:
Hi Robert,
No I did not get this error, and I was going to tell you... My platform is Mac OSX + MacPorts. I use MacPorts python 2.7 and sfepy installed manually from git.
Le mercredi 29 août 2012 18:07:50 UTC+2, Robert Cimrman a écrit :
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?
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
participants (3)
-
Alec Kalinin -
David Libault -
Robert Cimrman