]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Juan Carlos Araujo Cabarcas: use higher precision arithmetic in computing...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Aug 2013 14:46:56 +0000 (14:46 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Aug 2013 14:46:56 +0000 (14:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@30357 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/fe/mapping_q.cc

index 760fb61e2676bea90ef5111649daf8124196459c..c3781bca612388b00f8ca51c19f33eef26216588 100644 (file)
@@ -56,6 +56,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  Fixed: Some operations in the MappingQ class are now done in higher
+  precision arithmetic to mitigate the ill-conditioning that appears
+  when using mappings of high order (say, order 6 or 8 or 10).
+  <br>
+  (Juan Carlos Araujo Cabarcas, 2013/08/20)
+  </li>
+
   <li>
   Fixed: The SLEPcWrappers classes could not be compiled for 64-bit
   indices. This is now fixed.
index 18cb2b111045fd59d847b1ffcd54f81ffe349603..8b41ecab316713a77ec02d2b7026be969c88736c 100644 (file)
@@ -742,7 +742,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
   // one. check this
   for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
     Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
-                                     lohvs[unit_point].end(),0.) - 1)<1e-13,
+                                     lohvs[unit_point].end(),0.) - 1)<1e-12*this->degree*this->degree,
            ExcInternalError());
 }
 
@@ -790,29 +790,39 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
 
   // Compute the stiffness matrix of
   // the inner dofs
-  FullMatrix<double> S(n_inner);
+  FullMatrix<long double> S(n_inner);
   for (unsigned int point=0; point<n_q_points; ++point)
     for (unsigned int i=0; i<n_inner; ++i)
       for (unsigned int j=0; j<n_inner; ++j)
-        S(i,j) += contract(quadrature_data.derivative(point, n_outer+i),
-                           quadrature_data.derivative(point, n_outer+j))
-                  * quadrature.weight(point);
+       {
+         long double res = 0.;
+         for (unsigned int l=0; l<dim; ++l)
+           res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                  (long double)quadrature_data.derivative(point, n_outer+j)[l];
+
+         S(i,j) += res * (long double)quadrature.weight(point);
+       }
 
   // Compute the components of T to be the
   // product of gradients of inner and
   // outer shape functions.
-  FullMatrix<double> T(n_inner, n_outer);
+  FullMatrix<long double> T(n_inner, n_outer);
   for (unsigned int point=0; point<n_q_points; ++point)
     for (unsigned int i=0; i<n_inner; ++i)
       for (unsigned int k=0; k<n_outer; ++k)
-        T(i,k) += contract(quadrature_data.derivative(point, n_outer+i),
-                           quadrature_data.derivative(point, k))
-                  *quadrature.weight(point);
+       {
+         long double res = 0.;
+         for (unsigned int l=0; l<dim; ++l)
+           res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                  (long double)quadrature_data.derivative(point, k)[l];
+
+         T(i,k) += res *(long double)quadrature.weight(point);
+       }
 
-  FullMatrix<double> S_1(n_inner);
+  FullMatrix<long double> S_1(n_inner);
   S_1.invert(S);
 
-  FullMatrix<double> S_1_T(n_inner, n_outer);
+  FullMatrix<long double> S_1_T(n_inner, n_outer);
 
   // S:=S_1*T
   S_1.mmult(S_1_T,T);

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.