Re: [Libmesh-users] Quadrature points From: John Peterson - 2009-09-30 14:44 ```On Wed, Sep 30, 2009 at 9:23 AM, Petry Stefan 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 .0100975 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 ```

[Libmesh-users] Quadrature points Petry Stefan <petryste@hs...>
 Re: [Libmesh-users] Quadrature points From: John Peterson - 2009-09-30 14:44 ```On Wed, Sep 30, 2009 at 9:23 AM, Petry Stefan 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 .0100975 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 ```

 Re: [Libmesh-users] Quadrature points From: Kirk, Benjamin (JSC-EG311) - 2009-09-30 14:57 ```>> 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> .0100975> 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 ```

 Re: [Libmesh-users] Quadrature points From: Roy Stogner - 2009-09-30 15:10 ```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 .0100975 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 ```

 Re: [Libmesh-users] Quadrature points From: John Peterson - 2009-09-30 15:35 ```On Wed, Sep 30, 2009 at 9:23 AM, Petry Stefan 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 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 ```

 Re: [Libmesh-users] Quadrature points From: Roy Stogner - 2009-10-01 14:15 ```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 ```

 Re: [Libmesh-users] Quadrature points From: John Peterson - 2009-10-01 14:24 ```On Thu, Oct 1, 2009 at 1:46 AM, Petry Stefan 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

 Re: [Libmesh-users] Quadrature points From: Roy Stogner - 2009-10-02 13:28 ```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 ```

 Re: [Libmesh-users] Quadrature points From: Roy Stogner - 2009-10-05 13:45 ```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 ```