How does variables and integrals are evaluated?

[Alec Kalinin]

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:
1. 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.

[Robert Cimrman]

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.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:
> 1. 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...

> 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?

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.

[Alec Kalinin]

Hello,

> 1. 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

[Robert Cimrman]

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.

[Alec Kalinin]

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.

[Robert Cimrman]

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.

[Alec Kalinin]

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

[Robert Cimrman]

Hi Alec,

I meant volume integral as an Integral instance=quadrature points, not the volume term :). I will try your file when I get to my laptop.

r.

--
You received this message because you are subscribed to the Google Groups "sfepy-devel" group.
To view this discussion on the web visit https://groups.google.com/d/msg/sfepy-devel/-/hMvweknCgkkJ.
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.

[Robert Cimrman]

The following runs for me:

def post_process(out, pb, state, extend = False):

import numpy as nm
from sfepy.fem import Integral

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)')

coors = nm.array([[0.25, 0.25, 0.25]], dtype=nm.float64)
weights = nm.array([1.0], dtype=nm.float64)
integral = Integral('ii', coors=coors, weights=weights)

flux_var.data_from_qp(flux, integral)

val = flux_var()
print val.min(), val.max()

out.update(flux_var.create_output())

return out

As i1 has 1 udrature point, I have manually create a volume integral
with a single quadrature point, so that the array dimensions match. But
it is not the way the correct values should be obtained :)

The correct way would be to average from surface faces to surface
vertices, and set the other vertex values to zeros.

r.