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
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?