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
