How does variables and integrals are evaluated?
Hello SfePy devels,
Could you, please, explain me some SfePy internals?
(1) Question 1. Surface processing.
As I see in examples the 3D meshes contains only volume elements (tetrahedra) and no surface elements (triangles). But in problem solution pipeline we need to apply essential boundary conditions on surface nodes and need to evaluate some surface integrals, for example, for Neumann boundary conditions. How does SfePy working with surfaces?
For essential boundary conditions, yes, we can take only the far nodes in tetrahedra. But for the surface integrals we need the surface elements -- triangles. Does SfePy extract boundary triangles from volume tetrahedra?
(2) Question 2. Dimensions of the SfePy arrays.
Consider the postprocesing hook: def post_process(out, pb, state, extend = False): electric_field = pb.evaluate('ev_grad.2.Omega(u)', mode='el_avg') print electric_field.shape
flux = pb.evaluate('d_surface_flux.2.Gamma(m.K, u)', mode='el_avg')
print flux.shape
In my example shapes are: electric_field: (5235, 1, 3, 1) flux: (3898, 1, 1, 1)
Could you explain me why the variables are 4-dimensional arrays. My understanding is following:
- 1st dimension is the number of the tetrahedra element. 5235 is the total
number of the volume element. By what is 3898? Is it a number of the tetrahedra which are contains in boundary or is it a number of triangles on the surface? 2. The 3d dimension is the dimension of variable. The grad is a 3 dimensional vector, the flux is a scalar. But what are the 2nd and 4th dimensions?
Sincerely, Alexander.
On 09/06/2012 10:11 AM, Alec Kalinin wrote:
Hello SfePy devels,
Could you, please, explain me some SfePy internals?
Sure!
(1) Question 1. Surface processing.
As I see in examples the 3D meshes contains only volume elements (tetrahedra) and no surface elements (triangles). But in problem solution pipeline we need to apply essential boundary conditions on surface nodes and need to evaluate some surface integrals, for example, for Neumann boundary conditions. How does SfePy working with surfaces?
Every integration occurs over a region, and every term "knows" whether it corresponds to a volume or surface integral or something else. So defining the regions is crucial here.
For essential boundary conditions, yes, we can take only the far nodes in tetrahedra. But for the surface integrals we need the surface elements -- triangles. Does SfePy extract boundary triangles from volume tetrahedra?
Yes - each region contains a set of vertices, and also the edges, face, and cells that are completely defined by the vertices in the region. So if a region contains just a single vertex of a face, that face is not a part of the region, if it contains all vertices of a face, that face belongs to the region.
Then the term chooses, depending on its integration kind (volume, surface, ...) the correct entities from the region.
(2) Question 2. Dimensions of the SfePy arrays.
Consider the postprocesing hook: def post_process(out, pb, state, extend = False): electric_field = pb.evaluate('ev_grad.2.Omega(u)', mode='el_avg') print electric_field.shape
flux = pb.evaluate('d_surface_flux.2.Gamma(m.K, u)', mode='el_avg') print flux.shapeIn my example shapes are: electric_field: (5235, 1, 3, 1) flux: (3898, 1, 1, 1)
Could you explain me why the variables are 4-dimensional arrays. My understanding is following:
- 1st dimension is the number of the tetrahedra element. 5235 is the total
number of the volume element. By what is 3898? Is it a number of the tetrahedra which are contains in boundary or is it a number of triangles on the surface?
It is the number of the entities of the region. If you integrate over a whole volume domain, it's the number of elements. If you integrate over a surface region, it's the number of faces.
Then, to save it for visualization, one usually wants to extend the data over the whole mesh...
- The 3d dimension is the dimension of variable. The grad is a 3
dimensional vector, the flux is a scalar. But what are the 2nd and 4th dimensions?
The scheme is (n_entities, n_qp, n_row, n_col) - n_qp = number of quadrature points = 1 if the value is integrated, >= 1 for 'qp' term evaluation mode. So it is an array of 2d tensors/matrices of shape (n_row, n_col) over all elements/faces/edges and quadrature points in a region. The same data structure (4d array) is used internally in the C term evaluation functions in sfepy/terms/extmods/
Does it help? Feel free to ask more!
r.
Hello,
- 1st dimension is the number of the tetrahedra element. 5235 is the
total
number of the volume element. By what is 3898? Is it a number of the tetrahedra which are contains in boundary or is it a number of triangles
on
the surface?
It is the number of the entities of the region. If you integrate over a whole volume domain, it's the number of elements. If you integrate over a surface region, it's the number of faces.
Then, to save it for visualization, one usually wants to extend the data over the whole mesh...
OK, I see. So, for example, we integrate over a whole volume domain and obtain the data in each volume element (tetrahedra). Is it possible to convert this element data to the data in each volume vertex (node). I think, each vertex (node) is contained in several elements and we need to do some average to obtain data in the vertex. Does extend data procedure over the whole mesh do such thing?
Alec
Alec
On 09/07/2012 08:07 AM, Alec Kalinin wrote:
Hello,
- 1st dimension is the number of the tetrahedra element. 5235 is the
total
number of the volume element. By what is 3898? Is it a number of the tetrahedra which are contains in boundary or is it a number of triangles
on
the surface?
It is the number of the entities of the region. If you integrate over a whole volume domain, it's the number of elements. If you integrate over a surface region, it's the number of faces.
Then, to save it for visualization, one usually wants to extend the data over the whole mesh...
OK, I see. So, for example, we integrate over a whole volume domain and obtain the data in each volume element (tetrahedra). Is it possible to convert this element data to the data in each volume vertex (node). I think, each vertex (node) is contained in several elements and we need to do some average to obtain data in the vertex. Does extend data procedure over the whole mesh do such thing?
Yes, it is possible, Check the nodal_stress() function in examples/linear_elasticity/its2D_3.py - it evaluates the stress in quadrature points and then averages them into the mesh vertices. The key function is FieldVariable.data_from_qp(). After that, the dofs can be obtained by calling the variable instance.
To extend element-averaged data to all elements, use extend_cell_data().
Cheers, r.
Hello!
Yes, it is possible, Check the nodal_stress() function in
examples/linear_elasticity/its2D_3.py - it evaluates the stress in quadrature points and then averages them into the mesh vertices. The key function is FieldVariable.data_from_qp(). After that, the dofs can be obtained by calling the variable instance.
To extend element-averaged data to all elements, use extend_cell_data().
I tried to evaluate flux using data_from_qp, but I got a dimensional error:
# Integrals
integral_1 = {
'name' : 'i1',
'kind' : 's',
'order' : 2
}
...
def post_process(out, pb, state, extend = False):
flux = pb.evaluate('d_surface_flux.i1.Gamma(m.K, u)', mode = 'qp')
flux = extend_cell_data(flux, pb.domain, 'Gamma', is_surface = True)
flux_field = Field('flux', np.float64, (1,), pb.domain.regions['Omega'])
flux_var = FieldVariable('flux', 'parameter', flux_field, 1, primary_var_name='(set-to-None)')
flux_var.data_from_qp(flux, pb.integrals['i1'])
The error was: ValueError: incompatible shapes! ((5235, 4, 3, 4), out: (5235, 1, 1, 1), arr: (5235, 1, 1, 1))
What is wrong in my code?
Sincerely, Alec.
On 09/10/2012 02:09 PM, Alec Kalinin wrote:
Hello!
Yes, it is possible, Check the nodal_stress() function in
examples/linear_elasticity/its2D_3.py - it evaluates the stress in quadrature points and then averages them into the mesh vertices. The key function is FieldVariable.data_from_qp(). After that, the dofs can be obtained by calling the variable instance.
To extend element-averaged data to all elements, use extend_cell_data().
I tried to evaluate flux using data_from_qp, but I got a dimensional error:
# Integrals
integral_1 = {
'name' : 'i1',
'kind' : 's',
'order' : 2
}
...
def post_process(out, pb, state, extend = False):
flux = pb.evaluate('d_surface_flux.i1.Gamma(m.K, u)', mode = 'qp')
flux = extend_cell_data(flux, pb.domain, 'Gamma', is_surface = True)
flux_field = Field('flux', np.float64, (1,), pb.domain.regions['Omega'])
flux_var = FieldVariable('flux', 'parameter', flux_field, 1, primary_var_name='(set-to-None)')
flux_var.data_from_qp(flux, pb.integrals['i1'])
The error was: ValueError: incompatible shapes! ((5235, 4, 3, 4), out: (5235, 1, 1, 1), arr: (5235, 1, 1, 1))
It seems that the flux shape should be (5235, 4, 1, 1) (4 quadrature points). So you should use a volume integral in data_from_qp(), that you use in solving the problem, as extend_cell_data() returns the flux in the volume elements (averaged from the surface). This is probably not optimal for surface terms (a better way would be to directly average the face data to face vertices), but it is how things are now. If it does not help, send me your code so that I can debug it.
Cheers, r.
Hello!
It seems that the flux shape should be (5235, 4, 1, 1) (4 quadrature
points). So you should use a volume integral in data_from_qp(), that you use in solving the problem, as extend_cell_data() returns the flux in the volume elements (averaged from the surface).
Hmm... But the volume integral that is used in problem definition has the "dw_" prefix and cannot be used in evaluation mode. I think only "dq_" and "ev_" terms can be used in "qp" mode. I try to used "ev_grad" term but this term returns gradient (vector). It is also very nice, but I need scalar flux. :) To be more clear I attached my standalone problem definition python script and the mesh
This is probably not optimal for surface terms (a better way would be to directly average the face data to face vertices), but it is how things are now.
Yes, yes. This is not the problem to directly average the face elements data to the face vertices, it is about several lines of code. But I want to understand more about sfepy internals. :)
Sincerely, Alexander
participants (2)
-
Alec Kalinin -
Robert Cimrman