From: Denis Davydov <davydden@gmail.com>
Date: Sat, 27 Feb 2016 20:07:30 +0000 (+0100)
Subject: extend do_project_boundary_values() to complex numbers with ExcNotImplemented for... 
X-Git-Tag: v8.5.0-rc1~1277^2~6
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b031aa7478606d09bcee5b8f490ef75314e30ef4;p=dealii.git

extend do_project_boundary_values() to complex numbers with ExcNotImplemented for solver
---

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 68d9cc1983..af58c55d22 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -2165,6 +2165,77 @@ namespace VectorTools
 
   namespace
   {
+    // keep first argument is non-reference since we use it
+    // with 1e-8 * number
+    template <typename number1, typename number2>
+    bool real_part_bigger_than(const number1 a,
+                               const number2 &b)
+    {
+      return a > b;
+    }
+
+    template <typename number1, typename number2>
+    bool real_part_bigger_than(const std::complex<number1> a,
+                               const std::complex<number2> &b)
+    {
+      Assert(std::abs(a.imag()) < 1e-15 , ExcInternalError());
+      Assert(std::abs(b.imag()) < 1e-15 , ExcInternalError());
+      return a.real() > b.real();
+    }
+
+    template <typename number>
+    double min_number()
+    {
+      return std::numeric_limits<number>::min();
+    }
+
+    template <typename number>
+    double min_number(const number &dummy)
+    {
+      return std::numeric_limits<number>::min();
+    }
+
+    template <typename number>
+    double min_number(const std::complex<number> &dummy)
+    {
+      return std::numeric_limits<number>::min();
+    }
+
+
+    template <typename number>
+    void invert_mass_matrix(const SparseMatrix<number> &mass_matrix,
+                            const FilteredMatrix<Vector<number> > &filtered_mass_matrix,
+                            FilteredMatrix<Vector<number> > &filtered_precondition,
+                            const Vector<number> &rhs,
+                            Vector<number> &boundary_projection)
+    {
+      // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+      // steps may not be sufficient, since roundoff errors may accumulate for
+      // badly conditioned matrices
+      ReductionControl        control(5*rhs.size(), 0., 1.e-12, false, false);
+      GrowingVectorMemory<Vector<number> > memory;
+      SolverCG<Vector<number> >            cg(control,memory);
+
+      PreconditionSSOR<SparseMatrix<number> > prec;
+      prec.initialize(mass_matrix, 1.2);
+      filtered_precondition.initialize(prec, true);
+      // solve
+      cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
+      filtered_precondition.apply_constraints(boundary_projection, true);
+      filtered_precondition.clear();
+    }
+
+    template <typename number>
+    void invert_mass_matrix(const SparseMatrix<std::complex<number> > &mass_matrix,
+                            const FilteredMatrix<Vector<std::complex<number> > > &filtered_mass_matrix,
+                            FilteredMatrix<Vector<std::complex<number> > > &filtered_precondition,
+                            const Vector<std::complex<number> > &rhs,
+                            Vector<std::complex<number> > &boundary_projection)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+
     template <int dim, int spacedim, template <int, int> class DoFHandlerType,
               template <int,int> class M_or_MC, template <int> class Q_or_QC, typename number>
     void
@@ -2292,13 +2363,18 @@ namespace VectorTools
       FilteredMatrix<Vector<number> > filtered_precondition;
       std::vector<bool> excluded_dofs(mass_matrix.m(), false);
 
+      // we assemble mass matrix with unit weight,
+      // thus it will be real-valued irrespectively of the underlying algebra
+      // with positive elements on diagonal.
+      // Thus in order to extend this filtering to complex-algebra simply take
+      // the real-part of element.
       number max_element = 0.;
       for (unsigned int i=0; i<mass_matrix.m(); ++i)
-        if (mass_matrix.diag_element(i) > max_element)
+        if (real_part_bigger_than(mass_matrix.diag_element(i),max_element))
           max_element = mass_matrix.diag_element(i);
 
       for (unsigned int i=0; i<mass_matrix.m(); ++i)
-        if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
+        if (real_part_bigger_than(1.e-8 * max_element,mass_matrix.diag_element(i)))
           {
             filtered_mass_matrix.add_constraint(i, 0.);
             filtered_precondition.add_constraint(i, 0.);
@@ -2310,24 +2386,15 @@ namespace VectorTools
 
       // cannot reduce residual in a useful way if we are close to the square
       // root of the minimal double value
-      if (rhs.norm_sqr() < 1e28 * std::numeric_limits<number>::min())
+      if (rhs.norm_sqr() < 1e28 * min_number(number()))
         boundary_projection = 0;
       else
         {
-          // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
-          // steps may not be sufficient, since roundoff errors may accumulate for
-          // badly conditioned matrices
-          ReductionControl        control(5*rhs.size(), 0., 1.e-12, false, false);
-          GrowingVectorMemory<Vector<number> > memory;
-          SolverCG<Vector<number> >              cg(control,memory);
-
-          PreconditionSSOR<SparseMatrix<number> > prec;
-          prec.initialize(mass_matrix, 1.2);
-          filtered_precondition.initialize(prec, true);
-          // solve
-          cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
-          filtered_precondition.apply_constraints(boundary_projection, true);
-          filtered_precondition.clear();
+          invert_mass_matrix(mass_matrix,
+                             filtered_mass_matrix,
+                             filtered_precondition,
+                             rhs,
+                             boundary_projection);
         }
       // fill in boundary values
       for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)