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 <[email protected]> > > 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]> > >> >> >> 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]> >> >>> 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]> >>> >>>> 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 >> [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 >> >> --------- >> >> > > > -- > > 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
