## [Libmesh-devel] A possible bug in 3D?

 [Libmesh-devel] A possible bug in 3D? From: - 2007-05-21 19:06:02 ```Hi, I'd like to report a possible bug in the 3D element matrix calculation. The Hamiltonian in my test problem is: H = -1/2*laplacian + 1/2*x^2+2*x^4+1/2*x^6 + 1/2*y^2+2*y^4+1/2*y^6 + 1/2*z^2+2*z^4+1/2*z^6 + x*y + x*z + y*z So, the code for the element stiffness matrix is const Real x = q_point[qp](0); const Real y = q_point[qp](1); const Real z = q_point[qp](2); Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + ( .5*pow(x, 2.)+2.*pow(x, 4.)+.5*pow(x, 6.) +.5*pow(y, 2.)+2.*pow(y, 4.)+.5*pow(y, 6.) +.5*pow(z, 2.)+2.*pow(z, 4.)+.5*pow(z, 6.) +x*y+x*z+y*z) *phi[i][qp]*phi[j][qp]); The problem is that the matrix assembly stalled when I had 6th power functions in my code, I pin-pointed this problem by simply replacing pow(x, 6.) with x*x*x*x*x*x and the code ran just fine. The odd thing is that pow(x, 6) works for a 2D test case. I'd appreciate it if anyone will look into this and let me know. Thanks a lot in advance, -DX ```

 [Libmesh-devel] A possible bug in 3D? From: - 2007-05-21 19:06:02 ```Hi, I'd like to report a possible bug in the 3D element matrix calculation. The Hamiltonian in my test problem is: H = -1/2*laplacian + 1/2*x^2+2*x^4+1/2*x^6 + 1/2*y^2+2*y^4+1/2*y^6 + 1/2*z^2+2*z^4+1/2*z^6 + x*y + x*z + y*z So, the code for the element stiffness matrix is const Real x = q_point[qp](0); const Real y = q_point[qp](1); const Real z = q_point[qp](2); Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] + ( .5*pow(x, 2.)+2.*pow(x, 4.)+.5*pow(x, 6.) +.5*pow(y, 2.)+2.*pow(y, 4.)+.5*pow(y, 6.) +.5*pow(z, 2.)+2.*pow(z, 4.)+.5*pow(z, 6.) +x*y+x*z+y*z) *phi[i][qp]*phi[j][qp]); The problem is that the matrix assembly stalled when I had 6th power functions in my code, I pin-pointed this problem by simply replacing pow(x, 6.) with x*x*x*x*x*x and the code ran just fine. The odd thing is that pow(x, 6) works for a 2D test case. I'd appreciate it if anyone will look into this and let me know. Thanks a lot in advance, -DX ```
 [Libmesh-devel] [Libmesh-users] A possible bug in 3D? From: John Peterson - 2007-05-21 19:11:48 ```dxu@... writes: > Hi, I'd like to report a possible bug in the 3D element matrix calculation. > > The Hamiltonian in my test problem is: > > H = -1/2*laplacian > + 1/2*x^2+2*x^4+1/2*x^6 > + 1/2*y^2+2*y^4+1/2*y^6 > + 1/2*z^2+2*z^4+1/2*z^6 > + x*y + x*z + y*z > > So, the code for the element stiffness matrix is > > const Real x = q_point[qp](0); > const Real y = q_point[qp](1); > const Real z = q_point[qp](2); > > Ke(i,j) += JxW[qp]*(.5*dphi[i][qp]*dphi[j][qp] > + ( .5*pow(x, 2.)+2.*pow(x, 4.)+.5*pow(x, 6.) > +.5*pow(y, 2.)+2.*pow(y, 4.)+.5*pow(y, 6.) > +.5*pow(z, 2.)+2.*pow(z, 4.)+.5*pow(z, 6.) > +x*y+x*z+y*z) > *phi[i][qp]*phi[j][qp]); > > The problem is that the matrix assembly stalled when I had 6th power > functions in my code, I pin-pointed this problem by simply replacing > pow(x, 6.) with x*x*x*x*x*x and the code ran just fine. The odd thing > is that pow(x, 6) works for a 2D test case. "Stalled"? What does that mean, exactly? And why is a problem with "pow" a problem with libmesh? -John > I'd appreciate it if anyone will look into this and let me know. > > Thanks a lot in advance, > > -DX ```
 Re: [Libmesh-devel] [Libmesh-users] A possible bug in 3D? From: Roy Stogner - 2007-05-21 19:48:55 ```On Mon, 21 May 2007, John Peterson wrote: > dxu@... writes: > > > The problem is that the matrix assembly stalled when I had 6th power > > functions in my code, I pin-pointed this problem by simply replacing > > pow(x, 6.) with x*x*x*x*x*x and the code ran just fine. The odd thing > > is that pow(x, 6) works for a 2D test case. > > "Stalled"? What does that mean, exactly? And why is a problem with "pow" > a problem with libmesh? I think some old math libraries implement pow(x,y) by evaluating exp(y*log(x)), which can fail if x is zero or negative even when y is a positive integer. I can't see why it would work in 2D, though, unless xi can be negative in the 3D elements being used (e.g. hexes) but can't in the 2D elements (e.g. triangles). I also don't see why it would stall, unless the broken std::pow goes into an infinite loop instead of at least returning after setting EDOM. Anyway, the solution is to get a non-broken math library, use the templated Utility::pow<6>(x) code in libMesh, or write out x*x*... by hand. I'm not replying to repeat the obvious solution, though, but to ask a curious question: what sort of problem requires integrating 6th powers of master element coordinates? Perhaps I've just been working in Eulerian frames of reference too long, but I've always thought of the master element coordinates as just an arbitrary system we use to make numerical integration easier. If we changed the underlying tensor product elements from [-1,1]^d to [0,1]^d, would the application code here have to change to match? --- Roy ```
 Re: [Libmesh-devel] [Libmesh-users] A possible bug in 3D? From: Roy Stogner - 2007-05-21 19:52:32 ```On Mon, 21 May 2007, Roy Stogner wrote: > curious question: what sort of problem requires integrating 6th powers > of master element coordinates? Perhaps I've just been working in > Eulerian frames of reference too long, but I've always thought of the > master element coordinates as just an arbitrary system we use to make > numerical integration easier. If we changed the underlying tensor > product elements from [-1,1]^d to [0,1]^d, would the application code > here have to change to match? And of course, just as I was hitting send I realized that you're integrating physical coordinates, not local coordinates. Never mind. --- Roy ```