]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Disable 3x3 eigenvalues again, since it doesn't work with Maple... :-(
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 6 Apr 2001 16:30:46 +0000 (16:30 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 6 Apr 2001 16:30:46 +0000 (16:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@4396 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/numerics/derivative_approximation.cc

index 56aa426b90b2e196ceba359948050c42a2c7570d..1493891c9daa949c497ecdab039a3c38a040e19b 100644 (file)
@@ -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<dim>::symmetrize (Derivative &d)
 {
                                   // symmetrize non-diagonal entries
+  cout << "----------1 d" << endl
+       << d
+       << endl;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=i+1; j<dim; ++j)
       {
        const double s = (d[i][j] + d[j][i]) / 2;
        d[i][j] = d[j][i] = s;
       };
+  cout << "----------2 d" << endl
+       << d
+       << endl;
 };
 
 
@@ -607,7 +546,9 @@ DerivativeApproximation::approximate (const Mapping<dim>    &mapping,
                                       // gradient
       AssertThrow (determinant(Y) != 0,
                   ExcInsufficientDirections());
-
+      cout << "--------------- Y" << endl
+          << Y << endl;
+      
                                       // first symmetrize g
       DerivativeDescription::symmetrize (projected_derivative);
       

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.