From: Roy Stogner <roystgnr@ic...>  20070521 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 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? 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 nonbroken 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 