From c609fd8c74420f366e0a8cc5a60453c510d2d3d5 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 6 Apr 2001 16:30:46 +0000 Subject: [PATCH] Disable 3x3 eigenvalues again, since it doesn't work with Maple... :-( git-svn-id: https://svn.dealii.org/trunk@4396 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/derivative_approximation.cc | 141 +++++------------- 1 file changed, 41 insertions(+), 100 deletions(-) diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc index 56aa426b90..1493891c9d 100644 --- a/deal.II/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc @@ -159,106 +159,39 @@ double DerivativeApproximation::SecondDerivative<3>:: derivative_norm (const Derivative &d) { - // compute the three eigenvalues of - // the tensor @p{d} and take the - // largest: - const double t1 = d[1][2]*d[1][2]; - const double t2 = d[0][0]*t1; - const double t3 = d[0][1]*d[0][1]; - const double t4 = t3*d[2][2]; - const double t5 = d[0][2]*d[0][2]; - const double t6 = t5*d[1][1]; - const double t7 = d[0][0]*d[1][1]; - const double t8 = t7*d[2][2]; - const double t9 = d[0][1]*d[0][2]; - const double t10 = t9*d[1][2]; - const double t11 = t3*d[0][0]; - const double t12 = t3*d[1][1]; - const double t13 = d[0][0]*d[0][0]; - const double t14 = t13*d[1][1]; - const double t15 = d[1][1]*d[1][1]; - const double t16 = d[0][0]*t15; - const double t17 = t13*d[2][2]; - const double t19 = d[2][2]*d[2][2]; - const double t24 = t5*d[2][2]; - const double t25 = t1*d[2][2]; - const double t27 = t13*d[0][0]; - const double t28 = t19*d[2][2]; - const double t29 = t15*d[1][1]; - const double t30 = t1*t1; - const double t32 = d[1][1]*t28; - const double t35 = t15*t19; - const double t37 = t29*d[2][2]; - const double t39 = t13*t15; - const double t41 = d[0][0]*d[2][2]; - const double t43 = t19*t19; - const double t45 = t3*t1; - const double t47 = t3*t13; - const double t49 = t27*d[1][1]; - const double t51 = t27*d[2][2]; - const double t53 = t15*t15; - const double t57 = -3.0*t30*t15-6.0*t32*t5+6.0*t11*t29-6.0*t35*t5+24.0*t37*t5-6.0*t39* -t5-24.0*t41*t30+6.0*t7*t43-6.0*t45*t15-24.0*t47*t15+24.0*t49*t1-6.0*t51*t15 --12.0*t5*t53-36.0*t5*t30-3.0*t30*t19; - const double t59 = t13*t13; - const double t61 = t5*t5; - const double t64 = t3*t5; - const double t71 = d[0][0]*t29; - const double t78 = 24.0*t30*t13-12.0*t1*t59-3.0*t61*t19+6.0*t27*t28-6.0*t64*t13-60.0* -t5*t13*t1-60.0*t64*t19-60.0*t45*t19-6.0*t39*t1+24.0*t71*t5-6.0*t49*t19-3.0*t13* -t43+6.0*t29*t28-36.0*t61*t1-3.0*t61*t13; - const double t81 = t3*t3; - const double t82 = t81*d[0][0]; - const double t95 = d[0][0]*t28; - const double t99 = 18.0*t39*t19-24.0*t82*d[2][2]-3.0*t59*t15+6.0*t27*t29+6.0*t59*d[1][1]*d[2][2] --6.0*t14*t28-6.0*t71*t19-6.0*t16*t28-3.0*t13*t53+24.0*t81*t19-36.0*t81*t5-3.0* -t81*t13+6.0*t95*t5-6.0*t95*t1+30.0*t41*t61; - const double t105 = t13*t19; - const double t111 = d[1][1]*d[2][2]; - const double t114 = t5*t1; - const double t119 = -3.0*t81*t15-12.0*t3*t43-6.0*t71*t1-24.0*t7*t61-24.0*t7*t30-6.0* -t105*t1+24.0*t51*t1-24.0*t105*t5-6.0*t17*t29+6.0*t51*t5-24.0*t111*t61+30.0*t111 -*t30-6.0*t114*t19-60.0*t114*t15+6.0*t37*t1+6.0*t41*t114; - const double t140 = 30.0*t7*t1*t19+114.0*t7*t114+30.0*t7*t5*t19+30.0*t16*t25+30.0*t14* -t24-60.0*t7*t3*t19+30.0*t16*t4-60.0*t14*t25+30.0*t14*t4-60.0*t16*t24+6.0*t12* -t25-12.0*t81*t3+6.0*t11*t24+6.0*t12*t2+114.0*t4*t2; - const double t141 = t1*d[1][2]; - const double t163 = -216.0*d[0][0]*t141*t9+6.0*t11*t6+114.0*t4*t6+24.0*t61*t15+6.0*t32*t1 --3.0*t53*t19+6.0*d[0][0]*t53*d[2][2]-6.0*t49*t5-3.0*t15*t43-3.0*t59*t19-24.0*t81*d[1][1]* -d[2][2]-36.0*t3*t30-36.0*t3*t61-36.0*t81*t1+24.0*t9*d[1][2]*t27+24.0*t9*d[1][2]*t28; - const double t169 = t5*d[0][2]; - const double t170 = d[0][1]*t169; - const double t173 = d[1][2]*d[0][0]; - const double t182 = d[1][2]*d[1][1]; - const double t187 = d[1][2]*t13; - const double t191 = t3*d[0][1]; - const double t192 = t191*d[0][2]; - const double t198 = 24.0*t9*d[1][2]*t29+108.0*t9*t141*d[1][1]+108.0*t170*d[1][2]*d[2][2]+108.0*t170* -t173+108.0*t9*t141*d[2][2]-36.0*t9*d[1][2]*t15*d[2][2]-36.0*t9*t173*t19-36.0*t9*t182*t19 --36.0*t9*t173*t15-36.0*t9*t187*d[2][2]-60.0*t47*t1+108.0*t192*t182-36.0*t9*t187*d[1][1] -+108.0*t192*t173+144.0*t8*t10; - const double t209 = t3*t27; - const double t222 = -216.0*t169*d[1][1]*d[0][1]*d[1][2]-216.0*t191*d[2][2]*d[0][2]*d[1][2]-12.0*t30*t1-12.0* -t61*t5+6.0*t111*t114-6.0*t47*t19+6.0*t209*d[1][1]-60.0*t64*t15+30.0*t82*d[1][1]-6.0* -t209*d[2][2]-6.0*t3*t15*t19-24.0*t35*t1+24.0*t11*t28+24.0*t12*t28-6.0*t3*t29*d[2][2]+ -252.0*t64*t1; - const double t226 = sqrt(t57+t78+t99+t119+t140+t163+t198+t222); - const double t227 = -12.0*d[0][0]*t19-12.0*d[1][1]*t19-12.0*t15*d[2][2]+36.0*t5*d[0][0]+36.0*t24+36.0* -t25+36.0*t1*d[1][1]+8.0*t27+8.0*t28+8.0*t29+12.0*t226; - const double t229 = pow(-72.0*t2-72.0*t4-72.0*t6+48.0*t8+216.0*t10+36.0*t11+36.0*t12 --12.0*t14-12.0*t16-12.0*t17+t227,1.0/3.0); - const double t232 = (-t3/3+t7/9+t41/9+t111/9-t5/3-t1/3-t13/9-t19/9-t15/9)/t229; - const double t234 = sqrt(3.0); - const double t236 = t234*(t229/6+6.0*t232); - - const double eigenvalues[3] - = { t229/6-6.0*t232+d[0][0]/3+d[2][2]/3+d[1][1]/3, - -t229/12+3.0*t232+d[0][0]/3+d[2][2]/3+d[1][1]/3+sqrt(-1.0)*t236/2, - -t229/12+3.0*t232+d[0][0]/3+d[2][2]/3+d[1][1]/3-sqrt(-1.0)*t236/2 }; + cout << "**************** d" << endl + << d + << endl; + +/* + compute the three eigenvalues of the tensor @p{d} and take the + largest. one could use the following maple script to generate C + code: - return std::max (std::fabs (eigenvalues[0]), - std::max (std::fabs (eigenvalues[1]), - std::fabs (eigenvalues[2]))); + with(linalg); + readlib(C); + A:=matrix(3,3,[[a00,a01,a02],[a01,a11,a12],[a02,a12,a22]]); + E:=eigenvals(A); + EE:=vector(3,[E[1],E[2],E[3]]); + C(EE); + + Unfortunately, with both optimized and non-optimized output, at some + places cthe code `sqrt(-1.0)' is emitted, and I don't know what + Maple intends to do with it. This happens both with Maple4 and + Maple5. + + So, if someone has a handy way to compute the three eigenvalues of a + 3x3 matrix, send it to us. The trick is probably to tell Maple or + some other code generator that the matrix is symmetric and the + eigenvalues thus real, but how to do that? +*/ + + Assert (false, ExcNotImplemented()); + + const double EE[3] = { 0, 0, 0 }; + return std::max (std::fabs (EE[0]), + std::max (std::fabs (EE[1]), + std::fabs (EE[2]))); }; #endif @@ -295,12 +228,18 @@ void DerivativeApproximation::SecondDerivative::symmetrize (Derivative &d) { // symmetrize non-diagonal entries + cout << "----------1 d" << endl + << d + << endl; for (unsigned int i=0; i &mapping, // gradient AssertThrow (determinant(Y) != 0, ExcInsufficientDirections()); - + cout << "--------------- Y" << endl + << Y << endl; + // first symmetrize g DerivativeDescription::symmetrize (projected_derivative); -- 2.39.5