Re: [Getfem-users] Interpolate RT0 elements with P0

[Gerry Rizzo]

I've tried to use the interpolator function without success; here what I've done:

getfem::mesh_fem mf_p0_3;

getfem::pfem pf_p0;

pf_p0 = getfem::fem_descriptor("FEM_PK(3,0)");

mf_p0_3.init_with_mesh(mesh);

mf_p0_3.set_qdim(N);

mf_p0_3.set_finite_element(mesh.convex_index(), pf_p0);

std::vector<scalar_type> U_v_int;

U_v_int.resize(mf_p0_3.nb_basic_dof());

getfem::interpolation(mf_v,mf_p0_3,U_v,U_v_int);

mf_v is the RT0 fem and the computed velocity in RT0 is stored in U_v. I have the following error:

terminate called after throwing an instance of 'gmm::gmm_error'

  what():  Error in /usr/local/include/getfem/getfem_fem.h, line 728 : 

Wrong size for coeff vector

Thank you for the availability,

Calogero B. Rizzo

2013/12/2 Gerry Rizzo <address@hidden>

Hi,

I'm trying to make a function to interpolate three dimensional RT0 elements with P0 in order to export the function in vtk; the usual base functions are:

phy_i(x) = (x - a_i) / (d * |K|)

where a_i is a vertex of the tetrahedron, d is the dimension of the space (e.g. 3) and |K| is the volume of  the tetrahedron.

Then I can find the velocity in one point inside the thedraedron:

v(x) = sum ( v_i*phy_i(x) )

This is the function declaration (below there is the full function):

template<typename VEC,typename MAT>

void InterpolatorRT0( const getfem::mesh &mesh,                //mesh

const getfem::mesh_fem &mf_v,         //RT0 fem

const getfem::mesh_fem &mf_p,         //P0 fem

scalar_type &Vx,                                // variable for the computed Vx 

scalar_type &Vy,                                // variable for the computed Vy                                 

scalar_type &Vz,                                // variable for the computed Vz 

const scalar_type &x,                          // x, y, z are the point's coordinates where I want to compute the velocity

const scalar_type &y,

const scalar_type &z,

const VEC &gdl,                                // gdl is the velocity vector with RT0 elements

const size_type &elem,                       // elem is the thetraedron where the point lies

const MAT &A,                                  // A is the mass matrix

const size_type &shift = 0);

This function interpolate the RT0 velocity in one point given by x, y, z.

In order to guess the sign I'm using the mass matrix computed with:

getfem::asm_mass_matrix(A, im_v, mf_v, mf_p);

Unluckly this function doesn't work, I've tested it on a 3D box with constant velocity, I can't figure out the problem;

The problem may be the sign, is there a better/easier way to compute it?

Thank you,

Calogero B. Rizzo

Code snippet:

template<typename VEC,typename MAT>

void InterpolatorRT0( const getfem::mesh &mesh,

const getfem::mesh_fem &mf_v,

const getfem::mesh_fem &mf_p, 

scalar_type &Vx,

scalar_type &Vy,

scalar_type &Vz,

const scalar_type &x,

const scalar_type &y,

const scalar_type &z,

const VEC &gdl,

const size_type &elem,

const MAT &A,

const size_type &shift = 0) {

Vx = 0.;

Vy = 0.;

Vz = 0.;

scalar_type volume = mesh.convex_area_estimate(elem);

int N = mesh.dim();

std::vector<base_node> vertex;

for(int i=0; i<N+1; i++) {

vertex.push_back(mesh.points_of_convex(elem)[i]);

}

// loop over the four faces

for(int i=0; i<mf_v.nb_basic_dof_of_element(elem); i++) {

size_type idof = mf_v.ind_basic_dof_of_face_of_element(elem,i)[0];

scalar_type sign = A(idof,mf_p.ind_basic_dof_of_element(elem)[0]+shift);

sign /= (gmm::abs(sign)>0) ? gmm::abs(sign) : 1;

                        scalar_type area;

if(N==2) {

base_node x1 = mesh.points_of_face_of_convex(elem,i)[0];

base_node x2 = mesh.points_of_face_of_convex(elem,i)[1];

area = pow(pow(x2[0]-x1[0],2)+pow(x2[1]-x1[1],2),0.5);

} else {

assert(N==3);

base_node x1 = mesh.points_of_face_of_convex(elem,i)[0];

base_node x2 = mesh.points_of_face_of_convex(elem,i)[1];

base_node x3 = mesh.points_of_face_of_convex(elem,i)[2];

Geometry::Point3D xx1(x1[0],x1[1],x1[2]);

Geometry::Point3D xx2(x2[0],x2[1],x2[2]);

Geometry::Point3D xx3(x3[0],x3[1],x3[2]);

area = (Geometry::Triangle(xx1,xx2,xx3)).area();

}

// for a tetrahedron the vertex i is in front of the face i

Vx += gdl[idof]*(x-vertex[i][0])/(3*volume)*sign*area;

Vy += gdl[idof]*(y-vertex[i][1])/(3*volume)*sign*area;

Vz += gdl[idof]*(z-vertex[i][2])/(3*volume)*sign*area;

}

}