From: Petry S. <pet...@hs...> - 2009-09-30 14:37:19
|
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 |
From: John P. <pet...@cf...> - 2009-09-30 14:44:57
|
On Wed, Sep 30, 2009 at 9:23 AM, Petry Stefan <pet...@hs...> wrote: > 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: [code snipped] 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 Hi there, Are you confusing points in reference space (which is where the quadrature points are defined) and points in physical space (which is where the element is defined)? If you want to know the physical space locations of your quadrature points, just catch the reference returned by fe->get_xyz(); I only gave your code a cursory examination, if this is not the real issue I'm sorry. -- John |
From: Kirk, B. (JSC-EG311) <ben...@na...> - 2009-09-30 14:57:49
|
>> 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: > > [code snipped] > > 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 > > > Hi there, > > Are you confusing points in reference space (which is where the > quadrature points are defined) and points in physical space (which is > where the element is defined)? If you want to know the physical space > locations of your quadrature points, just catch the reference returned > by > > fe->get_xyz(); > > I only gave your code a cursory examination, if this is not the real > issue I'm sorry. Also, by default we are using a quadrature rule with nonnegative weights, and that ends up with a 4-point rule. the distribution thus "squishes" two points toward a single corner and are not generally symmetric - is that what is bothering you? The Gaussian quadrature points points are located *within* the element, and thus will be inside the bounds of the element... -Ben |
From: Roy S. <roy...@ic...> - 2009-09-30 15:10:09
|
I don't think John or Ben has hit on the problem yet; I likewise haven't had a chance to do more than skim your code, but it did look like you're correctly setting up a non-exterior quadrature rule and pulling xyz coordinates. On Wed, 30 Sep 2009, Petry Stefan wrote: > -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! Can you output all the nodal coordinates for this example, not just the min/max range? I'm not sure what's wrong, but more information might help us replicate or diagnose the problem. --- Roy |
From: John P. <pet...@cf...> - 2009-09-30 15:35:56
|
On Wed, Sep 30, 2009 at 9:23 AM, Petry Stefan <pet...@hs...> wrote: > 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); Is the value of "quadrature_points" still zero when you call this function? Is it passed by reference? What does this Compute_Quadrature_Points function do? -- John |
From: Roy S. <roy...@ic...> - 2009-10-01 14:15:55
|
On Thu, 1 Oct 2009, Petry Stefan wrote: > node 0: x=-0.0100975 ; y=0.0100975 > node 1: x=-0.010009 ; y=0.0101852 > node 2: x=-0.010081 ; y=0.0102585 > node 3: x=-0.0101701 ; y=0.0101701 This appears to be an inverted element by libMesh ordering. Have you run in debug or devel mode? Try recompiling and rerunning with METHOD=dbg and see if you trip any assertions. An inverted element by itself shouldn't cause the garbage quadrature points you're seeing, but it should have thrown an assertion once you tried to reinitialize the FE object on that element. Perhaps there's some other problem which the debugging assertions might catch. --- Roy |
From: John P. <pet...@cf...> - 2009-10-01 14:24:01
|
On Thu, Oct 1, 2009 at 1:46 AM, Petry Stefan <pet...@hs...> wrote: > Thank you for your help. > > The function Compute_Quadrature_Points evaluates the maximum number of quadrature points over all elements. This is done because there might be different types of elements within the mesh. In this example the variable quadrature_points has the value 4 on return. If it's the maximum over all elements, then I would be wary of using the following loop for(unsigned int i=0; i<quadrature_points; ++i) for *all* elements, since it will likely be going past the end of an array for elements which have fewer quadrature points... -- John |
From: Roy S. <roy...@ic...> - 2009-10-02 13:28:42
|
On Fri, 2 Oct 2009, Petry Stefan wrote: > Unfortunately, I am not able to use configure, because I am working > with libmesh on Windows. A student created an eclipse workspace for > the compilation of libmesh, Interesting. What are you guys doing about the POSIX dependencies that libMesh has? sys/time.h, unistd.h, stdlib.h... we've gotten libMesh to compile with Cygwin on Windows, but even that took a little doing. > but I am no experienced user of gcc. Do you know which compiler > switch of gcc has to be set to show debugging assertions? We do some different things depending on which version of gcc you're using, but the important parts for assertions are: Don't use -DNDEBUG, which turns off all our assertions. Use -DDEBUG, which turns on some of our debugging-only code. Use "-D_GLIBCXX_DEBUG -D_GLIBCXX_DEBUG_PEDANTIC", which turns on the libstdc++ container tests. (All this will make your code much too slow for production use, but is worth it when something's going wrong) Finally, recompile everything. We use extensions like dbg.o and opt.o to let debug-mode and optimized-mode object files live side-by-side, but your Eclipse setup may not be doing the same. And if it's not, then depending on how smart it is about recognizing compiler flags as a dependency, you may need to tell it delete your non-debugging object files before recompiling. --- Roy |
From: Roy S. <roy...@ic...> - 2009-10-05 13:45:42
|
On Mon, 5 Oct 2009, Petry Stefan wrote: > I recompiled libmesh with the compiler flags "-DDEBUG > -D_GLIBCXX_DEBUG_PEDANTIC" and linked it to our program. > Unfortunately, I had to skip the flag "-D_GLIBCXX_DEBUG", because > linker errors occured, when I tried to link a libmesh dll compiled > with this flag to our program. But If i understood your remarks > correctly the most important switch for debug assertions is "DDEBUG > ". However, rerunning our program with this debug version of libmesh > produced no assertions. Turning on DEBUG will activate the libMesh assertions. It won't catch array bounds errors, though. It does change the libstdc++ ABI, so to use it you'd have to recompile your own program with the same flag on. I think that's worth a try. There's not much else we can do for you unless you can create a test case that reproduces the problem and is small enough to post. > Concerning our use of libmesh on Windows I have to explain that we > only employ a very limited subset of libmesh's capabilities. > Especially we are only using a single processor. We were not able to > port a parallel-processing version of libmesh to Windows. For our 2d > application we can live with such a restriction, but I think in > general this will be a knock-out criterion. The tools that we use > are mingw and gcc. How far does the port attempt get in parallel? If PETSc works but libMesh doesn't, it shouldn't take too much effort to go the rest of the way. --- Roy |