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_MagR-1])
// 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 x-values for the element vary between -0.0101701
and -0.010009, the x-value 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
------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
