From: John Peterson <peterson@cf...>  20070521 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 pinpointed 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 