Hello,
I encountered a strange problem with the use of quadrature points. In a function I interpolate values of a given variable onto the quadrature points. The following code worked fine for elements of type TRI:
unsigned int quadrature_points = 0;
double **Mx_on_qpoints, **My_on_qpoints;
double rVal = 0., phiVal = 0., MrTmp = 0., MphiTmp = 0.;
MeshBase::const_element_iterator el = mesh>elements_begin(); // prepare iterator for elements
const MeshBase::const_element_iterator end_el = mesh>elements_end();
Mx_on_qpoints = new double *[mesh>n_elem()];
My_on_qpoints = new double *[mesh>n_elem()];
FEType fe_type; // create and initialize finite element with appropriate approximation order and quadrature rule
fe_type.order = (libMeshEnums::Order)VAL_ORDER_ES;
fe_type.family = (libMeshEnums::FEFamily)VAL_FAMILY;
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, (libMeshEnums::Order)VAL_ORDER_GAUSS);
fe>attach_quadrature_rule(&qrule);
/*
* we are computing all quadrature points first before we use it,
* because then we have the right number for e.g. using triangles and quads
*/
Compute_Quadrature_Points(quadrature_points);
const vector<Point>& q_point = fe>get_xyz(); // physical coordinates of quadrature points
ofstream Magnetization_On_Quadrature_Points_in_file("./output/Magnetization_On_Quadrature_Points.txt");
for ( ; el != end_el ; ++el)
{
const Elem* elem = *el;
fe>reinit(elem); // update values for current element
unsigned int elemId = elem>id();
Mx_on_qpoints[elemId] = new double[quadrature_points];
My_on_qpoints[elemId] = new double[quadrature_points];
if(Material_Number[elem>subdomain_id()] == MAGNET_LINEAR) // in magnets
{
///////////////////// inserted only for debugging purposes
Point Zentrum=elem>centroid();
double xz=Zentrum(0);
double yz=Zentrum(1);
double rz=sqrt(xz*xz+yz*yz);
double xmin,xmax,ymin,ymax;
Node *nd;
xmin=1.0e300;
ymin=1.0e300;
xmax=1.0e300;
ymax=1.0e300;
for(unsigned int n=0;n<elem>n_nodes();n++)
{
nd=elem>get_node(n);
xz=(*nd)(0);
if(xz<xmin)
xmin=xz;
if(xz>xmax)
xmax=xz;
yz=(*nd)(1);
if(yz<ymin)
ymin=yz;
if(yz>ymax)
ymax=yz;
}
cout << xmin << "<x<" << xmax << endl;
cout << ymin << "<y<" << ymax << endl;
///////////////////////////
for(unsigned int i=0; i<quadrature_points; ++i)
{
const Real x = q_point[i](0), y = q_point[i](1);
///////////////////// inserted only for debugging purposes
cout << "Quadrature point " << i << ": ";
cout << "x=" << x << " ; y=" << y << endl;
cout << "Quadrature point " << i;
if(elem>contains_point(q_point[i]))
cout << " is in element!" << endl;
else
cout << " is not in element!" << endl;
///////////////////////////
phiVal = atan2(y, x); // convert to polar coordinates
if(phiVal < 0)
phiVal += 2*M_PI; // angular coordinates are from the interval [0,2*pi] instead [pi,pi]
rVal = sqrt(pow(x, 2) + pow(y, 2));
if(rVal<R_Mag[0]  rVal>R_Mag[N_MagR1]) // outside the magnet
{
Mx_on_qpoints[elemId][i] = 0.;
My_on_qpoints[elemId][i] = 0.;
}
else
{
Interpolation(rVal, phiVal, Mr, Mphi, R_Mag, Phi_Mag, N_MagR, N_MagPhi, MrTmp, MphiTmp); // bilinear interpolation for MrTmp and MphiTmp
Mx_on_qpoints[elemId][i] = MrTmp*cos(phiVal)  MphiTmp*sin(phiVal); // conversion to cartesian vector components
My_on_qpoints[elemId][i] = MrTmp*sin(phiVal) + MphiTmp*cos(phiVal);
}
Magnetization_On_Quadrature_Points_in_file << x << "\t\t" << y << "\t\t" <<
Mx_on_qpoints[elemId][i] << "\t\t" << My_on_qpoints[elemId][i] << endl;
}
}
else
{
for(unsigned int i=0; i<quadrature_points; ++i) // no magnet
{
Mx_on_qpoints[elemId][i] = 0.;
My_on_qpoints[elemId][i] = 0.;
}
}
}
Transfer_Magnetization(Mx_on_qpoints, My_on_qpoints, mesh>n_elem(), quadrature_points); // transfer values to Magnetostatic_Simualtion_Private
Clean_Up(My_on_qpoints, mesh>n_elem());
Clean_Up(Mx_on_qpoints, mesh>n_elem());
} //end Magnetization_On_Quadrature_Points()
Part of the code has only been added for debugging purposes. For elements of type QUAD4 the quadrature points are no longer inside the element, when I check this by determining the minimum x and y coordinate directly from the nodes of the element. However, the libmesh function elem>contains_point(q_point[i]) states that the point is inside the element. Here is a sample output
0.0101701<x<0.010009
.0100975<y<0.0102585
Quadrature point 0: x=0.00839908 ; y=0.00843638
Quadrature point 0 is in element!
Quadrature point 1: x=0.00958875 ; y=0.00972795
Quadrature point 1 is in element!
Quadrature point 2: x=0.00381006 ; y=0.00384751
Quadrature point 2 is in element!
Quadrature point 3: x=0.00838955 ; y=0.00852933
Quadrature point 3 is in element!
As you can see, although the xvalues for the element vary between 0.0101701 and 0.010009, the xvalue of the first quadrature point is 0.00839908. I guess that I do something wrong with the initialization of the finite element or the quadrature rule, but I can not figure out what. The values of the initialization variable are
dim=2
VAL_ORDER_ES=1
VAL_FAMILY=0
VAL_ORDER_GAUSS=3
Thanks for your help in advance,
Werner
