Dear Calogero B. Rizzo,
If you can give me a simple compilable test program with the call of
getfem::interpolation I can try to check what is wrong.
Yves.
Le 03/12/2013 15:21, Gerry Rizzo a écrit :
> 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 <[email protected]
> <mailto:[email protected]>>
>
> 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
