extend error estimator to complex-valued PETSc

diff --git a/include/deal.II/numerics/error_estimator.templates.h b/include/deal.II/numerics/error_estimator.templates.h
index e4f8c4b047929e54e5a50747abcbffa0c6f35391..908752922006ea97812fe3f73fe34a4b6e5b9a46 100644 (file)
--- 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 3420168aa76f54e34884f0199429d373f80d1466..788488570aa75facaea73f3bd4e31b93e9459373 100644 (file)
--- 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();
                   }
           }