From 5e95a73749f0df854ab54074fceaaaf4236a5fa0 Mon Sep 17 00:00:00 2001 From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Tue, 20 Aug 2013 14:46:56 +0000 Subject: [PATCH] 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 --- deal.II/doc/news/changes.h | 8 ++++++++ deal.II/source/fe/mapping_q.cc | 32 +++++++++++++++++++++----------- 2 files changed, 29 insertions(+), 11 deletions(-) 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); -- 2.39.5