|
From: | Gerry Rizzo |
Subject: | Re: [Getfem-users] Interpolate RT0 elements with P0 |
Date: | Thu, 5 Dec 2013 12:46:05 +0100 |
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 vectorThank you for the availability,Calogero B. Rizzo2013/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, //meshconst getfem::mesh_fem &mf_v, //RT0 femconst getfem::mesh_fem &mf_p, //P0 femscalar_type &Vx, // variable for the computed Vxscalar_type &Vy, // variable for the computed Vyscalar_type &Vz, // variable for the computed Vzconst scalar_type &x, // x, y, z are the point's coordinates where I want to compute the velocityconst scalar_type &y,const scalar_type &z,const VEC &gdl, // gdl is the velocity vector with RT0 elementsconst size_type &elem, // elem is the thetraedron where the point liesconst MAT &A, // A is the mass matrixconst 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. RizzoCode 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 facesfor(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 iVx += 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;}}
[Prev in Thread] | Current Thread | [Next in Thread] |