72 views

Skip to first unread message

Sep 6, 2012, 4:11:45 AM9/6/12

to sfepy...@googlegroups.com

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.

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.

Sep 6, 2012, 5:57:48 AM9/6/12

to sfepy...@googlegroups.com

On 09/06/2012 10:11 AM, Alec Kalinin wrote:

> Hello SfePy devels,

>

> Could you, please, explain me some SfePy internals?

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

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?

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?

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?

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.

Sep 7, 2012, 2:07:25 AM9/7/12

to sfepy...@googlegroups.com

Hello,

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

> 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

Sep 7, 2012, 4:02:41 AM9/7/12

to sfepy...@googlegroups.com

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.

Sep 10, 2012, 8:09:23 AM9/10/12

to sfepy...@googlegroups.com

Hello!

I tried to evaluate flux using data_from_qp, but I got a dimensional error:

...

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.

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.

Sep 10, 2012, 6:33:03 PM9/10/12

to sfepy...@googlegroups.com

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.

Sep 11, 2012, 4:11:37 AM9/11/12

to sfepy...@googlegroups.com

Hello!

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

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

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

Sep 11, 2012, 6:16:42 AM9/11/12

to sfepy...@googlegroups.com

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.

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.

Sep 11, 2012, 6:25:55 PM9/11/12

to sfepy...@googlegroups.com

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.

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

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

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.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu