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

[Yves Renard]

Dear Gerry,

I order to interpolate the RTO on given points, you can still directly use the 
interpolation functions using the mesh_trans_inv objects where you store the 
points on which you need to
interpolate (see the file getfem/getfem_interpolation.h).

Yves.

----- Mail original -----
De: "Gerry Rizzo" <address@hidden>
À: "Yves Renard" <address@hidden>
Cc: address@hidden
Envoyé: Jeudi 12 Décembre 2013 11:18:56
Objet: Re: [Getfem-users] Interpolate RT0 elements with P0


Dear Yves, 


I've tested the fixed code and now it works like a charm! 
I have one more question concerning RT0 interpolation, I couldn't find a method 
to compute the velocity in one point given an RT0 vector; I thought in building 
a new fem with only one point, then use the interpolation function. 


Gerry 



2013/12/9 Yves Renard < address@hidden > 





Dear Gerry, 

It seems to come from an error on line 320 of getfem_interpolation.h 
coeff[qq].resize(nb_s*qdim); 
have to be replaced by 
coeff[qq].resize(cvnbdof); 
the loop which is just after is from 0 to cvnbdof. I don't understand why this 
error has not been detected yet ! 

I commited the changes. 


Yves. 



Le 09/12/2013 17:07, Gerry Rizzo a écrit : 





Hi, 
here I've attached the example, please let me know if you receive it correctly. 


Thank you again, 
Gerry 



2013/12/9 Yves Renard < address@hidden > 






Dear Gerry, 

Sorry, but it seems that your attached file has not been transmitted to the 
list. Could you send me it directly ? 

Yves. 




Le 05/12/2013 12:46, Gerry Rizzo a écrit : 





Hi, 
I've written a very simple test case for the RT0 interpolation, you can find it 
attached to this email. 
It loads a vector (example.txt) containing the velocity of the solution of 
Darcy equations over a 10x10x10 mesh with constant permeability; the velocity 
is a positive constant along x-axis and zero along y and z axes. 


I still have the same error I've showed you in the previous email. 
I've tested this code with getfem 4.2, ubuntu 12.10 32bit and gcc version 
4.7.2. 


Thank you, 
Calogero B. Rizzo 



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



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; 
} 
} 





_______________________________________________
Getfem-users mailing list address@hidden 
https://mail.gna.org/listinfo/getfem-users 

-- 

  Yves Renard ( address@hidden )       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE http://math.univ-lyon1.fr/~renard --------- 


-- 

  Yves Renard ( address@hidden )       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE http://math.univ-lyon1.fr/~renard ---------