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 <[email protected] > <mailto:[email protected]>> > > > > 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 <[email protected] >> <mailto:[email protected]>> >> >> 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] <mailto:[email protected]> >> https://mail.gna.org/listinfo/getfem-users > > > -- > > Yves Renard ([email protected] > <mailto:[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 <http://math.univ-lyon1.fr/%7Erenard> > > --------- > > -- 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
