Dear Gerry Rizzo,
Why don't you use Getfem interpolation functions in getfem_interpolation.h ?
They should work in that case also.
You can interpolate with them on a Lagrange fem (for instance P0) or on
a given cloud of point.
Yves.
Le 02/12/2013 16:49, Gerry Rizzo a écrit :
> 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
> [email protected]
> https://mail.gna.org/listinfo/getfem-users
--
Yves Renard ([email protected]) 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
---------
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
