From dec0650b26beb1d7373113908867a496691ef1d9 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
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 <deal.II/base/numbers.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/quadrature_lib.h>
@@ -476,7 +477,7 @@ namespace internal
         for (unsigned int component=0; component<n_components; ++component)
           if (parallel_data.component_mask[component] == true)
             for (unsigned int p=0; p<n_q_points; ++p)
-              face_integral[n] += Utilities::fixed_power<2>(parallel_data.phi[n][p][component]) *
+              face_integral[n] += numbers::NumberTraits<number>::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<const parallel::distributed::Triangulation<1,spacedim>*>
       (&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<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > >
+  std::vector<std::vector<std::vector<Tensor<1,spacedim,number> > > >
   gradients_here (n_solution_vectors,
-                  std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > >(2, std::vector<Tensor<1,spacedim,typename InputVector::value_type> >(n_components)));
-  std::vector<std::vector<std::vector<Tensor<1,spacedim,typename InputVector::value_type> > > >
+                  std::vector<std::vector<Tensor<1,spacedim,number> > >(2, std::vector<Tensor<1,spacedim,number> >(n_components)));
+  std::vector<std::vector<std::vector<Tensor<1,spacedim,number> > > >
   gradients_neighbor (gradients_here);
-  std::vector<Vector<typename InputVector::value_type> >
-  grad_neighbor (n_solution_vectors, Vector<typename InputVector::value_type>(n_components));
+  std::vector<Vector<number> >
+  grad_neighbor (n_solution_vectors, Vector<number>(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<number>::abs_square(jump) * cell->diameter();
                   }
           }
 
-- 
2.39.5