From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 20 Aug 2013 14:46:56 +0000 (+0000)
Subject: Patch by Juan Carlos Araujo Cabarcas: use higher precision arithmetic in computing... 
X-Git-Tag: v8.1.0~1031
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=55a9b189eb16bcdf1f9b1de231fbb384aab71aa8;p=dealii.git

Patch by Juan Carlos Araujo Cabarcas: use higher precision arithmetic in computing some data for high order mappings.

git-svn-id: https://svn.dealii.org/trunk@30357 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 760fb61e26..c3781bca61 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -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.
diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc
index 18cb2b1110..8b41ecab31 100644
--- a/deal.II/source/fe/mapping_q.cc
+++ b/deal.II/source/fe/mapping_q.cc
@@ -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);