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;
}
}