From dec0650b26beb1d7373113908867a496691ef1d9 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 22 Feb 2016 13:39:15 +0100 Subject: [PATCH] extend error estimator to complex-valued PETSc --- .../numerics/error_estimator.templates.h | 3 ++- source/numerics/error_estimator_1d.cc | 17 +++++++++-------- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h index e4f8c4b047..9087529220 100644 --- a/include/deal.II/numerics/error_estimator.templates.h +++ b/include/deal.II/numerics/error_estimator.templates.h @@ -13,6 +13,7 @@ // // --------------------------------------------------------------------- +#include #include #include #include @@ -476,7 +477,7 @@ namespace internal for (unsigned int component=0; component(parallel_data.phi[n][p][component]) * + face_integral[n] += numbers::NumberTraits::abs_square(parallel_data.phi[n][p][component]) * parallel_data.JxW_values[p]; return face_integral; diff --git a/source/numerics/error_estimator_1d.cc b/source/numerics/error_estimator_1d.cc index 3420168aa7..788488570a 100644 --- a/source/numerics/error_estimator_1d.cc +++ b/source/numerics/error_estimator_1d.cc @@ -219,6 +219,7 @@ estimate (const Mapping<1,spacedim> &mapping, const types::subdomain_id subdomain_id_, const types::material_id material_id) { + typedef typename InputVector::value_type number; #ifdef DEAL_II_WITH_P4EST if (dynamic_cast*> (&dof_handler.get_triangulation()) @@ -301,13 +302,13 @@ estimate (const Mapping<1,spacedim> &mapping, // // for the neighbor gradient, we need several auxiliary fields, depending on // the way we get it (see below) - std::vector > > > + std::vector > > > gradients_here (n_solution_vectors, - std::vector > >(2, std::vector >(n_components))); - std::vector > > > + std::vector > >(2, std::vector >(n_components))); + std::vector > > > gradients_neighbor (gradients_here); - std::vector > - grad_neighbor (n_solution_vectors, Vector(n_components)); + std::vector > + grad_neighbor (n_solution_vectors, Vector(n_components)); // reserve some space for coefficient values at one point. if there is no // coefficient, then we fill it by unity once and for all and don't set it @@ -434,12 +435,12 @@ estimate (const Mapping<1,spacedim> &mapping, if (component_mask[component] == true) { // get gradient here - const double grad_here = gradients_here[s][n][component] + const number grad_here = gradients_here[s][n][component] * normal; - const double jump = ((grad_here - grad_neighbor[s](component)) * + const number jump = ((grad_here - grad_neighbor[s](component)) * coefficient_values(component)); - (*errors[s])(cell->active_cell_index()) += jump*jump * cell->diameter(); + (*errors[s])(cell->active_cell_index()) += numbers::NumberTraits::abs_square(jump) * cell->diameter(); } } -- 2.39.5