]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Interface to MatShell matrix-type in PETSc
authorMartin Steigemann <msteigemann@yahoo.de>
Mon, 24 Sep 2012 21:44:07 +0000 (21:44 +0000)
committerMartin Steigemann <msteigemann@yahoo.de>
Mon, 24 Sep 2012 21:44:07 +0000 (21:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@26689 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/base/quadrature_lib.h
deal.II/include/deal.II/base/table_handler.h
deal.II/include/deal.II/lac/petsc_precondition.h
deal.II/include/deal.II/lac/petsc_solver.h
deal.II/source/base/quadrature_lib.cc
deal.II/source/base/table_handler.cc
deal.II/source/distributed/tria.cc
deal.II/source/lac/petsc_precondition.cc
deal.II/source/lac/petsc_solver.cc

index 678bb73933995b4617c1c4fe94ef960309ceb728..a73a6d2f26d85677df21295b7df9d8cbd7dd1c9e 100644 (file)
@@ -172,7 +172,7 @@ class QTrapez : public Quadrature<dim>
 
 /**
  * Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5.
- * See Stoer: Einfhrung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einführung in die Numerische Mathematik I, p. 102
  */
 template <int dim>
 class QMilne : public Quadrature<dim>
@@ -184,7 +184,7 @@ class QMilne : public Quadrature<dim>
 
 /**
  * Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7.
- * See Stoer: Einfhrung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einführung in die Numerische Mathematik I, p. 102
  */
 template <int dim>
 class QWeddle : public Quadrature<dim>
@@ -430,75 +430,6 @@ class QGaussOneOverR : public Quadrature<dim>
 
 
 
-/**
- * Gauss Quadrature Formula with $1/R^{3/2}$ weighting function. This formula
- * can be used to to integrate $1/R^{3/2} \ f(x)$ on the reference
- * element $[0,1]^2$, where $f$ is a smooth function without
- * singularities, and $R$ is the distance from the point $x$ to the vertex
- * $\xi$, given at construction time by specifying its index. Notice that
- * this distance is evaluated in the reference element.
- *
- * This quadrature formula is a specialization of QGaussOneOverR.
- * We apply a second transformation $R = t^2$ to cancel a singularity
- * of order $R^{1/2}$:
- * \f[
- * \int_0^1 \int_0^1 \frac{1}{R^{1/2}} f(x,y) \biggr|_{x < y} dxdy =
- * \int_0^1 \int_0^{r(\pi/4 v)} \frac{1}{R^{1/2}} f(R \cos(\pi/4 v), R \sin(\pi/4 v)) \frac{\pi}{4} R dr dv =
- * 2 \int_0^1 \int_0^{\sqrt{r(\pi/4 v)}} f(t^2 \cos(\pi/4 v), t^2 \sin(\pi/4 v)) \frac{\pi}{4} t^2 dt dv =
- * 2 \int_0^1 \int_0^1 f(t^2, t^2 \tan(\pi/4 v)) \frac{\pi}{4} \frac{t^2}{(\cos(\pi/4 v))^{3/2}} dt dv
- * \f]
- *
- * Upon construction it is possible to specify wether we want the
- * singularity removed, or not. In other words, this quadrature can be
- * used to integrate $g(x) = 1/R^{3/2}\ f(x)$, or simply $f(x)$, with the $1/R^{3/2}$
- * factor already included in the quadrature weights.
- */
-template<int dim>
-class QGaussOneOverRThreeHalfs : public Quadrature<dim>
-{
-  public:
-    /**
-     * The constructor takes three arguments: the order of the Gauss
-     * formula, the index of the vertex where the singularity is
-     * located, and whether we include the weighting singular function
-     * inside the quadrature, or we leave it in the user function to
-     * be integrated. Notice that this constructor only works for the
-     * vertices of the quadrilateral.
-     *
-     * Traditionally, quadrature formulas include their weighting
-     * function, and the last argument is set to false by
-     * default. There are cases, however, where this is undesirable
-     * (for example when you only know that your singularity has the
-     * same order of 1/R, but cannot be written exactly in this
-     * way).
-     *
-     * In other words, you can use this function in either of
-     * the following way, obtaining the same result:
-     *
-     * @code
-     * QGaussOneOverRThreeHalfs singular_quad(order, vertex_id, false);
-     * // This will produce the integral of f(x)/R^{3/2}
-     * for(unsigned int i=0; i<singular_quad.size(); ++i)
-     *          integral += f(singular_quad.point(i))*singular_quad.weight(i);
-     *
-     * // And the same here
-     * QGaussOneOverRThreeHalfs singular_quad_noR(order, vertex_id, true);
-     *
-     * // This also will produce the integral of f(x)/R^{3/2}, but 1/R^{3/2} has to
-     * // be specified.
-     * for(unsigned int i=0; i<singular_quad.size(); ++i) {
-     *   double R = (singular_quad_noR.point(i)-cell->vertex(vertex_id)).norm();
-     *   integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R^{3/2};
-     * }
-     * @endcode
-     */
-    QGaussOneOverRThreeHalfs(const unsigned int n,
-                             const unsigned int vertex_index,
-                             const bool factor_out_singular_weight=false);
-};
-
-
-
 /*@}*/
 
 /* -------------- declaration of explicit specializations ------------- */
@@ -529,7 +460,6 @@ template <> QWeddle<1>::QWeddle ();
 template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
 template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag);
 template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag);
-template <> QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs (const unsigned int n, const unsigned int index, const bool flag);
 
 
 
index 8bcceffebff8f41d414299c3bf2374769f6f6f52..0c6366717e335171cd51ed0528d91b712a711c46 100644 (file)
@@ -505,23 +505,6 @@ class TableHandler
                                       */
     void write_tex (std::ostream &file, const bool with_header=true) const;
 
-                                     /**
-                                      * Write table as a tex file
-                                      * and use the booktabs package.
-                                      * Values are written in
-                                      * math mode: $value$
-                                      * If with_header is set to false
-                                      * (it is true by default), then
-                                      * no "\documentclass{...}",
-                                      * "\begin{document}" and
-                                      * "\end{document}" are used. In
-                                      * this way the file can be
-                                      * included into an existing tex
-                                      * file using a command like
-                                      * "\input{table_file}".
-                                      */
-    void write_results (std::ostream &file, const bool with_header=true) const;
-
                                      /**
                                       * Read or write the data of this
                                       * object to or from a stream for
@@ -888,8 +871,6 @@ namespace internal
   }
 }
 
-
-
 template <typename T>
 void TableHandler::add_value (const std::string &key,
                               const T            value)
index 43730abcf156bea60c54f49dcbb3e7f4e7b52f3a..af3e05acab0c42863a1fc44aca4c8d9a76ab6ce9 100644 (file)
@@ -804,107 +804,6 @@ namespace PETScWrappers
           bool output_details;
       };
 
-                                       /**
-                                        * Empty Constructor. You need to call
-                                        * initialize() before using this
-                                        * object.
-                                        */
-      PreconditionBoomerAMG ();
-
-                                       /**
-                                        * Constructor. Take the matrix which
-                                        * is used to form the preconditioner,
-                                        * and additional flags if there are
-                                        * any.
-                                        */
-      PreconditionBoomerAMG (const MatrixBase     &matrix,
-                             const AdditionalData &additional_data = AdditionalData());
-
-                                       /**
-                                        * Initializes the preconditioner
-                                        * object and calculate all data that
-                                        * is necessary for applying it in a
-                                        * solver. This function is
-                                        * automatically called when calling
-                                        * the constructor with the same
-                                        * arguments and is only used if you
-                                        * create the preconditioner without
-                                        * arguments.
-                                        */
-      void initialize (const MatrixBase     &matrix,
-                       const AdditionalData &additional_data = AdditionalData());
-
-    protected:
-                                       /**
-                                        * Store a copy of the flags for this
-                                        * particular preconditioner.
-                                        */
-      AdditionalData additional_data;
-  };
-
-
-
-/**
- * A class that implements the interface to use the ParaSails sparse
- * approximate inverse preconditioner from the HYPRE suite. Note that
- * PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
- *
- * @ingroup PETScWrappers
- * @author Martin Steigemann, 2011
- */
-  class PreconditionParaSails : public PreconditionerBase
-  {
-    public:
-                                       /**
-                                        * Standardized data struct to
-                                        * pipe additional flags to the
-                                        * preconditioner.
-                                        */
-      struct AdditionalData
-      {
-                                           /**
-                                            * Constructor.
-                                            */
-          AdditionalData (
-            const unsigned int symmetric = 0,
-            const unsigned int n_levels = 1,
-            const double threshold = 0.1,
-            const double filter = 0.05,
-            const double load_bal = 0.,
-            const bool output_details = false
-          );
-
-                                           /**
-                                            * This parameter has the following meanings,
-                                            * to indicate the symmetry and definiteness
-                                            * of the problem, and to specify the type
-                                            * of the preconditioner to construct:
-                                            * <ul>
-                                            * <li> @p 0: nonsymmetric and/or indefinite problem, and nonsymmetric preconditioner
-                                            * <li> @p 1: SPD problem, and SPD (factored) preconditioner
-                                            * <li> @p 2: nonsymmetric, definite problem, and SPD (factored) preconditioner
-                                            * </ul>
-                                            *
-                                            */
-          unsigned int symmetric;
-
-          unsigned int n_levels;
-
-          double threshold;
-
-          double filter;
-
-          double load_bal;
-
-                                           /**
-                                            * Setting this flag to true
-                                            * produces output from HYPRE,
-                                            * when the preconditioner
-                                            * is constructed.
-                                            */
-         bool output_details;
-      };
-
 
 
                                        /**
@@ -912,7 +811,7 @@ namespace PETScWrappers
                                         * initialize() before using this
                                         * object.
                                         */
-      PreconditionParaSails ();
+      PreconditionBoomerAMG ();
 
                                        /**
                                         * Constructor. Take the matrix which
@@ -920,7 +819,7 @@ namespace PETScWrappers
                                         * and additional flags if there are
                                         * any.
                                         */
-      PreconditionParaSails (const MatrixBase     &matrix,
+      PreconditionBoomerAMG (const MatrixBase     &matrix,
                              const AdditionalData &additional_data = AdditionalData());
 
                                        /**
@@ -944,160 +843,6 @@ namespace PETScWrappers
                                         */
       AdditionalData additional_data;
   };
-
-
-
-/**
- * A class that implements the interface to use the scalable implementation
- * of the Parallel ILU algorithm in the Euclid library from the HYPRE suite.
- * Note that PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
- *
- * @ingroup PETScWrappers
- * @author Martin Steigemann, 2011
- */
-  class PreconditionEuclid : public PreconditionerBase
-  {
-    public:
-                                       /**
-                                        * Standardized data struct to
-                                        * pipe additional flags to the
-                                        * preconditioner.
-                                        */
-      struct AdditionalData
-      {
-                                           /**
-                                            * Constructor.
-                                            */
-          AdditionalData (
-            const unsigned int level = 1,
-            const bool use_block_jacobi = false,
-            const bool output_details = false
-          );
-
-                                           /**
-                                            * Factorization level for ILU(k). Default is 1.
-                                            * For 2D convection-diffusion and similar problems,
-                                            * fastest solution time is typically obtained with
-                                            * levels 4 through 8. For 3D problems, fastest solution
-                                            * time is typically obtained with level 1.
-                                            */
-          unsigned int level;
-
-                                           /**
-                                            * Use Block Jacobi ILU preconditioning
-                                            * instead of PILU. Default is false.
-                                            * If subdomains contain relatively few nodes
-                                            * (less than 1000), or the problem is not
-                                            * well partitioned, Block Jacobi ILU
-                                            * may give faster solution time then PILU.
-                                            */
-          unsigned int use_block_jacobi;
-
-                                           /**
-                                            * Setting this flag to true
-                                            * produces debug output from
-                                            * HYPRE, when the preconditioner
-                                            * is constructed.
-                                            */
-          bool output_details;
-      };
-
-
-
-                                       /**
-                                        * Empty Constructor. You need to call
-                                        * initialize() before using this
-                                        * object.
-                                        */
-      PreconditionEuclid ();
-
-                                       /**
-                                        * Constructor. Take the matrix which
-                                        * is used to form the preconditioner,
-                                        * and additional flags if there are
-                                        * any.
-                                        */
-      PreconditionEuclid (const MatrixBase     &matrix,
-                         const AdditionalData &additional_data = AdditionalData());
-
-                                       /**
-                                        * Initializes the preconditioner
-                                        * object and calculate all data that
-                                        * is necessary for applying it in a
-                                        * solver. This function is
-                                        * automatically called when calling
-                                        * the constructor with the same
-                                        * arguments and is only used if you
-                                        * create the preconditioner without
-                                        * arguments.
-                                        */
-      void initialize (const MatrixBase     &matrix,
-                       const AdditionalData &additional_data = AdditionalData());
-
-    protected:
-                                       /**
-                                        * Store a copy of the flags for this
-                                        * particular preconditioner.
-                                        */
-      AdditionalData additional_data;
-  };
-
-
-
-/**
- * A class that implements a non-preconditioned Krylov method.
- *
- * @ingroup PETScWrappers
- * @author Martin Steigemann, 2011
- */
-  class PreconditionNone : public PreconditionerBase
-  {
-    public:
-                                       /**
-                                        * Standardized data struct to
-                                        * pipe additional flags to the
-                                        * preconditioner.
-                                        */
-      struct AdditionalData
-      {};
-
-                                       /**
-                                        * Empty Constructor. You need to call
-                                        * initialize() before using this
-                                        * object.
-                                        */
-      PreconditionNone ();
-
-                                       /**
-                                        * Constructor. Take the matrix which
-                                        * is used to form the preconditioner,
-                                        * and additional flags if there are
-                                        * any.
-                                        */
-      PreconditionNone (const MatrixBase     &matrix,
-                       const AdditionalData &additional_data = AdditionalData());
-
-                                       /**
-                                        * Initializes the preconditioner
-                                        * object and calculate all data that
-                                        * is necessary for applying it in a
-                                        * solver. This function is
-                                        * automatically called when calling
-                                        * the constructor with the same
-                                        * arguments and is only used if you
-                                        * create the preconditioner without
-                                        * arguments.
-                                        */
-      void initialize (const MatrixBase     &matrix,
-                       const AdditionalData &additional_data = AdditionalData());
-
-    protected:
-                                       /**
-                                        * Store a copy of the flags for this
-                                        * particular preconditioner.
-                                        */
-      AdditionalData additional_data;
-  };
 }
 
 
index 0f3b3b5ab49186b2d29c1d554257cc82836d9071..984ea7b6b5ebd18c7e2818a4bf38671d9144b326 100644 (file)
@@ -138,14 +138,6 @@ namespace PETScWrappers
              const PreconditionerBase &preconditioner);
 
 
-
-      void
-      solve (const MatrixBase               &A,
-             VectorBase                     &x,
-             const VectorBase               &b,
-             const PreconditionerBase       &preconditioner,
-             const std::vector<VectorBase>  &nullspace);
-
                                        /**
                                         * Resets the contained preconditioner
                                         * and solver object. See class
@@ -177,10 +169,6 @@ namespace PETScWrappers
                       << "An error with error number " << arg1
                       << " occurred while calling a PETSc function");
 
-      DeclException1 (ExcPETScSolverError,
-                      char*,
-                      << "PETSc solver failed: " << arg1);
-
     protected:
 
                                        /**
index 498f45f66b3db5c3ab458a9593fb2ddb363fb5b1..a38e566433b5fa6c463fbbd722655322f0911c7c 100644 (file)
@@ -925,93 +925,6 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
 }
 
 
-
-template<>
-QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs(const unsigned int n,
-                                                      const unsigned int vertex_index,
-                                                      const bool factor_out_singularity) :
-    Quadrature<2>(2*n*n)
-{
-                                   // This version of the constructor
-                                   // works only for the 4
-                                   // vertices. If you need a more
-                                   // general one, you should use the
-                                   // one with the Point<2> in the
-                                   // constructor.
-    Assert(vertex_index <4, ExcIndexRange(vertex_index, 0, 4));
-
-                                   // Start with the gauss quadrature
-                                   // formula on the (u,v) reference
-                                   // element.
-    QGauss<2> gauss(n);
-
-    Assert(gauss.size() == n*n, ExcInternalError());
-    Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
-
-    // We create only the first one. All other pieces are rotation of
-    // this one.
-    // In this case the transformation is
-    //
-    // (x,y) = (u*u, u*u tan(pi/4 v))
-    //
-    // with Jacobian
-    //
-    // J = 2 pi/4 R*R sqrt(cos(pi/4 v))
-    //
-    // And we get rid of R to take into account the singularity,
-    // unless specified differently in the constructor.
-    std::vector<Point<2> >  &ps = this->quadrature_points;
-    std::vector<double>     &ws = this->weights;
-    double pi4 = numbers::PI/4;
-
-    for(unsigned int q=0; q<gauss.size(); ++q) {
-        const Point<2> &gp = gauss.point(q);
-        ps[q][0] = gp[0]*gp[0];
-        ps[q][1] = gp[0]*gp[0]*std::tan(pi4 *gp[1]);
-        ws[q]    = 2.*gauss.weight(q)*pi4/std::sqrt(std::cos(pi4 *gp[1]));
-        if(factor_out_singularity) {
-            const double  abs_value = (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm();
-            ws[q] *= abs_value;
-            ws[q] *= std::sqrt(abs_value);
-        }
-        // The other half of the quadrilateral is symmetric with
-        // respect to xy plane.
-        ws[gauss.size()+q]    = ws[q];
-        ps[gauss.size()+q][0] = ps[q][1];
-        ps[gauss.size()+q][1] = ps[q][0];
-    }
-
-    // Now we distribute these vertices in the correct manner
-    double theta = 0;
-    switch(vertex_index) {
-    case 0:
-        theta = 0;
-        break;
-    case 1:
-        theta = numbers::PI/2;
-        break;
-    case 2:
-        theta = -numbers::PI/2;
-        break;
-    case 3:
-        theta = numbers::PI;
-        break;
-    }
-
-    double R00 =  std::cos(theta), R01 = -std::sin(theta);
-    double R10 =  std::sin(theta), R11 =  std::cos(theta);
-
-    if(vertex_index != 0)
-        for(unsigned int q=0; q<size(); ++q) {
-            double x = ps[q][0]-.5,  y = ps[q][1]-.5;
-
-            ps[q][0] = R00*x + R01*y + .5;
-            ps[q][1] = R10*x + R11*y + .5;
-        }
-}
-
-
-
 // construct the quadrature formulae in higher dimensions by
 // tensor product of lower dimensions
 
index 51973470119dcc5b63b62dd2fce587958c723913..674651b356699f6983746bd5878e3e61890af783 100644 (file)
@@ -635,136 +635,6 @@ void TableHandler::write_tex (std::ostream &out, const bool with_header) const
 }
 
 
-void TableHandler::write_results (std::ostream &out, const bool with_header) const
-{
-  bool  math_mode = true;
-
-  AssertThrow (out, ExcIO());
-  if (with_header)
-    out << "\\documentclass[10pt]{report}" << std::endl
-        << "\\usepackage{float}" << std::endl
-        << "\\usepackage{booktabs}" << std::endl
-        << "\\usepackage{amsmath}" << std::endl << std::endl
-        << "\\begin{document}" << std::endl;
-
-  out << "\\begin{table}" << std::endl
-      << "\\begin{center}" << std::endl
-      << "\\begin{tabular}{@{}";
-
-  std::vector<std::string> sel_columns;
-  get_selected_columns(sel_columns);
-
-                                   // write the column formats
-  for (unsigned int j=0; j<column_order.size(); ++j)
-    {
-      std::string key=column_order[j];
-                                       // avoid `supercolumns[key]'
-      const std::map<std::string, std::vector<std::string> >::const_iterator 
-        super_iter=supercolumns.find(key);
-
-      if (super_iter!=supercolumns.end())
-        {
-          const unsigned int n_subcolumns=super_iter->second.size();
-          for (unsigned int k=0; k<n_subcolumns; ++k)
-            {
-                                       // avoid `columns[supercolumns[key]]'
-              const std::map<std::string, Column>::const_iterator
-                col_iter=columns.find(super_iter->second[k]);
-              Assert(col_iter!=columns.end(), ExcInternalError());
-
-              out << col_iter->second.tex_format;
-            }
-        }
-      else
-        {
-                                     // avoid `columns[key]';
-          const std::map<std::string, Column>::const_iterator
-            col_iter=columns.find(key);
-          Assert(col_iter!=columns.end(), ExcInternalError());
-          out << col_iter->second.tex_format;
-        }
-    }
-  out << "@{}} \\toprule" << std::endl;
-
-                                   // write the caption line of the table
-  for (unsigned int j=0; j<column_order.size(); ++j)
-    {
-      std::string key=column_order[j];
-      const std::map<std::string, std::vector<std::string> >::const_iterator 
-        super_iter=supercolumns.find(key);
-
-      if (super_iter!=supercolumns.end())
-        {
-          const unsigned int n_subcolumns=super_iter->second.size();
-                                       // avoid use of `tex_supercaptions[key]'
-          std::map<std::string,std::string>::const_iterator 
-            tex_super_cap_iter=tex_supercaptions.find(key);
-          out << std::endl << "\\multicolumn{" << n_subcolumns << "}{c}{" 
-              << tex_super_cap_iter->second << "}";
-        }
-      else
-        {
-                                       // col_iter->second=columns[col];
-          const std::map<std::string, Column>::const_iterator
-            col_iter=columns.find(key);
-          Assert(col_iter!=columns.end(), ExcInternalError());
-          out << col_iter->second.tex_caption;
-        }
-      if (j<column_order.size()-1)
-        out << " & ";
-    }
-  out << "\\\\ \\midrule" << std::endl;
-
-                                   // write the n rows
-  const unsigned int nrows=n_rows();
-  for (unsigned int i=0; i<nrows; ++i)
-    {
-      const unsigned int n_cols=sel_columns.size();
-
-      for (unsigned int j=0; j<n_cols; ++j)
-        {
-          std::string key=sel_columns[j];
-                                       // avoid `column[key]'
-          const std::map<std::string, Column>::const_iterator
-            col_iter=columns.find(key);
-          Assert(col_iter!=columns.end(), ExcInternalError());
-
-          const Column &column=col_iter->second;
-
-          out << std::setprecision(column.precision);
-
-          if (col_iter->second.scientific)
-            out.setf(std::ios::scientific, std::ios::floatfield);
-          else
-            out.setf(std::ios::fixed, std::ios::floatfield);
-
-          if (math_mode)
-            out << "$";
-
-          out << column.entries[i].value;
-
-          if (math_mode)
-            out << "$";
-
-          if (j<n_cols-1)
-            out << " & ";
-        }
-
-      out << (i == nrows-1 ? "\\\\ \\bottomrule" : "\\\\") << std::endl;
-    }
-
-  out << "\\end{tabular}" << std::endl
-      << "\\end{center}" << std::endl;
-  if(tex_table_caption!="")
-    out << "\\caption{"  << tex_table_caption << "}" << std::endl;
-  if(tex_table_label!="")
-    out << "\\caption{"   << tex_table_label << "}" << std::endl;
-    out << "\\end{table}" << std::endl;
-  if (with_header)
-    out << "\\end{document}" << std::endl;
-}
-
-
 unsigned int TableHandler::n_rows() const
 {
   if (columns.size() == 0)
index 8d794272cc70e171de6b4f2aede8eec0de2c8de3..e025665818c3256f8a525cce5e267ceabc1a07ef 100644 (file)
@@ -3005,79 +3005,9 @@ namespace parallel
     template <int dim, int spacedim>
     void
     Triangulation<dim,spacedim>::
-    copy_triangulation (const dealii::Triangulation<dim, spacedim> &old_tria)
+    copy_triangulation (const dealii::Triangulation<dim, spacedim> &)
     {
-      clear();
-
-      try
-        {
-          dealii::Triangulation<dim,spacedim>::
-            copy_triangulation (old_tria);
-    }
-      catch (const typename dealii::Triangulation<dim,spacedim>::DistortedCellList &)
-        {
-                                       // the underlying
-                                       // triangulation should not
-                                       // be checking for
-                                       // distorted cells
-          AssertThrow (false, ExcInternalError());
-        }
-
-                                       // note that now we have some content in
-                                       // the p4est objects and call the
-                                       // functions that do the actual work
-                                       // (which are dimension dependent, so
-                                       // separate)
-      triangulation_has_content = true;
-
-      Assert (old_tria.n_levels() == 1,
-              ExcMessage ("Parallel distributed triangulations can only be copied, "
-                          "if they are not refined!"));
-
-      if (dynamic_cast<const dealii::parallel::distributed::Triangulation<dim,spacedim> *>(&old_tria) != 0)
-        {
-          Assert (!(dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
-                   (old_tria).refinement_in_progress),
-                  ExcMessage ("Parallel distributed triangulations can only "
-                              "be copied, if no refinement is in progress!"));
-
-          coarse_cell_to_p4est_tree_permutation =
-            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
-            (old_tria).coarse_cell_to_p4est_tree_permutation;
-
-          p4est_tree_to_coarse_cell_permutation =
-            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
-            (old_tria).p4est_tree_to_coarse_cell_permutation;
-
-          attached_data_size =
-            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
-            (old_tria).attached_data_size;
-
-          n_attached_datas   =
-            dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
-            (old_tria).n_attached_datas;
-        }
-      else
-        {
-          setup_coarse_cell_to_p4est_tree_permutation ();
-        };
-
-      copy_new_triangulation_to_p4est (dealii::internal::int2type<dim>());
-
-      try
-        {
-          copy_local_forest_to_triangulation ();
-        }
-      catch (const typename Triangulation<dim>::DistortedCellList &)
-        {
-                                       // the underlying
-                                       // triangulation should not
-                                       // be checking for
-                                       // distorted cells
-          AssertThrow (false, ExcInternalError());
-        }
-
-      update_number_cache ();
+      Assert (false, ExcNotImplemented());
     }
 
 
index 370d5c99b5de874b140a55313bf3e1fba25b8068..2df82399d2060b90fd170068603724aea65bcdc4 100644 (file)
@@ -483,223 +483,6 @@ namespace PETScWrappers
   }
 
 
-/* ----------------- PreconditionParaSails -------------------- */
-
-  PreconditionParaSails::AdditionalData::
-  AdditionalData(const unsigned int symmetric,
-                 const unsigned int n_levels,
-                 const double threshold,
-                 const double filter,
-                 const double load_bal,
-                 const bool output_details)
-                  :
-                  symmetric(symmetric),
-                  n_levels(n_levels),
-                  threshold(threshold),
-                  filter(filter),
-                  load_bal(load_bal),
-                  output_details(output_details)
-  {}
-
-
-  PreconditionParaSails::PreconditionParaSails ()
-  {}
-
-
-  PreconditionParaSails::PreconditionParaSails (const MatrixBase     &matrix,
-                                                const AdditionalData &additional_data)
-  {
-    initialize(matrix, additional_data);
-  }
-
-
-  void
-  PreconditionParaSails::initialize (const MatrixBase     &matrix_,
-                                     const AdditionalData &additional_data_)
-  {
-    matrix = static_cast<Mat>(matrix_);
-    additional_data = additional_data_;
-
-#ifdef PETSC_HAVE_HYPRE
-    create_pc();
-
-    int ierr;
-    ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCHYPRESetType(pc, "parasails");
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    if (additional_data.output_details)
-      PetscOptionsSetValue("-pc_hypre_parasails_logging","1");
-
-    Assert ((additional_data.symmetric == 0 ||
-             additional_data.symmetric == 1 ||
-             additional_data.symmetric == 2),
-            ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
-
-    std::stringstream ssStream;
-
-    switch (additional_data.symmetric)
-      {
-        case 0:
-        {
-          ssStream << "nonsymmetric";
-          break;
-        }
-
-        case 1:
-        {
-          ssStream << "SPD";
-          break;
-        }
-
-        case 2:
-        {
-          ssStream << "nonsymmetric,SPD";
-          break;
-        }
-
-        default:
-          Assert (false,
-                  ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
-      };
-
-    PetscOptionsSetValue("-pc_hypre_parasails_sym",ssStream.str().c_str());
-
-    PetscOptionsSetValue("-pc_hypre_parasails_nlevels",
-                         Utilities::int_to_string(
-                           additional_data.n_levels
-                         ).c_str());
-
-    ssStream.str(""); // empty the stringstream
-    ssStream << additional_data.threshold;
-    PetscOptionsSetValue("-pc_hypre_parasails_thresh", ssStream.str().c_str());
-
-    ssStream.str(""); // empty the stringstream
-    ssStream << additional_data.filter;
-    PetscOptionsSetValue("-pc_hypre_parasails_filter", ssStream.str().c_str());
-
-    ssStream.str(""); // empty the stringstream
-    ssStream << additional_data.load_bal;
-    PetscOptionsSetValue("-pc_hypre_parasails_loadbal", ssStream.str().c_str());
-
-    ierr = PCSetFromOptions (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCSetUp (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-#else // PETSC_HAVE_HYPRE
-    (void)pc;
-    Assert (false,
-            ExcMessage ("Your PETSc installation does not include a copy of "
-                        "the hypre package necessary for this preconditioner."));
-#endif
-  }
-
-
-/* ----------------- PreconditionEuclid ----------------------- */
-
-  PreconditionEuclid::AdditionalData::
-  AdditionalData(const unsigned int level,
-                 const bool use_block_jacobi,
-                 const bool output_details)
-                  :
-                  level(level),
-                  use_block_jacobi(use_block_jacobi),
-                  output_details(output_details)
-  {}
-
-
-  PreconditionEuclid::PreconditionEuclid ()
-  {}
-
-
-  PreconditionEuclid::PreconditionEuclid (const MatrixBase     &matrix,
-                                                const AdditionalData &additional_data)
-  {
-    initialize(matrix, additional_data);
-  }
-
-
-  void
-  PreconditionEuclid::initialize (const MatrixBase     &matrix_,
-                                  const AdditionalData &additional_data_)
-  {
-    matrix = static_cast<Mat>(matrix_);
-    additional_data = additional_data_;
-
-#ifdef PETSC_HAVE_HYPRE
-    create_pc();
-
-    int ierr;
-    ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCHYPRESetType(pc, "euclid");
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    PetscOptionsSetValue("-pc_hypre_euclid_levels",
-                         Utilities::int_to_string(
-                           additional_data.level
-                         ).c_str());
-
-    if (additional_data.use_block_jacobi)
-      PetscOptionsSetValue("-pc_hypre_euclid_bj","1");
-
-    if (additional_data.output_details)
-      PetscOptionsSetValue("-pc_hypre_euclid_print_statistics","1");
-
-    ierr = PCSetFromOptions (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCSetUp (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-#else // PETSC_HAVE_HYPRE
-    (void)pc;
-    Assert (false,
-            ExcMessage ("Your PETSc installation does not include a copy of "
-                        "the hypre package necessary for this preconditioner."));
-#endif
-  }
-
-
-/* ----------------- PreconditionNone ------------------------- */
-
-  PreconditionNone::PreconditionNone ()
-  {}
-
-
-  PreconditionNone::PreconditionNone (const MatrixBase     &matrix,
-                                      const AdditionalData &additional_data)
-  {
-    initialize(matrix, additional_data);
-  }
-
-
-  void
-  PreconditionNone::initialize (const MatrixBase     &matrix_,
-                                const AdditionalData &additional_data_)
-  {
-    matrix = static_cast<Mat>(matrix_);
-    additional_data = additional_data_;
-
-    create_pc();
-
-    int ierr;
-    ierr = PCSetType (pc, const_cast<char *>(PCNONE));
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCSetFromOptions (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    ierr = PCSetUp (pc);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-  }
-
-
 /* ----------------- PreconditionLU -------------------- */
 
   PreconditionLU::AdditionalData::
index 0cbb09d2bc5f0a23dd493173daf5d0ed9a9384ec..d01dae47f093932a1cea97103968cb4552016564 100644 (file)
@@ -129,138 +129,6 @@ namespace PETScWrappers
   }
 
 
-  void
-  SolverBase::solve (const MatrixBase               &A,
-                     VectorBase                     &x,
-                     const VectorBase               &b,
-                     const PreconditionerBase       &preconditioner,
-                     const std::vector<VectorBase>  &nullspace)
-  {
-    int ierr;
-                                     // first create a solver object if this
-                                     // is necessary
-    if (solver_data.get() == 0)
-      {
-        solver_data.reset (new SolverData());
-
-        ierr = KSPCreate (mpi_communicator, &solver_data->ksp);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-                                         // set the matrices involved. the
-                                         // last argument is irrelevant here,
-                                         // since we use the solver only once
-                                         // anyway
-        ierr = KSPSetOperators (solver_data->ksp, A, preconditioner,
-                                SAME_PRECONDITIONER);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-                                         // let derived classes set the solver
-                                         // type, and the preconditioning
-                                         // object set the type of
-                                         // preconditioner
-        set_solver_type (solver_data->ksp);
-
-       ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc());
-       AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        PC  prec;
-        ierr = KSPGetPC (solver_data->ksp, &prec);
-        ierr = PCFactorSetShiftType (prec, MAT_SHIFT_POSITIVE_DEFINITE);
-
-                                         // then a convergence monitor
-                                         // function. that function simply
-                                         // checks with the solver_control
-                                         // object we have in this object for
-                                         // convergence
-#if DEAL_II_PETSC_VERSION_LT(3,0,0)
-        KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
-                               reinterpret_cast<void *>(&solver_control));
-#else
-        KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
-                               reinterpret_cast<void *>(&solver_control),
-                               PETSC_NULL);
-#endif
-
-      }
-
-                                     // then do the real work: set up solver
-                                     // internal data and solve the
-                                     // system. 
-    ierr = KSPSetUp (solver_data->ksp);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    unsigned int dim_nullsp = nullspace.size();
-
-    if (dim_nullsp > 0)
-      {
-        Vec * nullsp_basis;
-
-        PetscMalloc (dim_nullsp*sizeof(Vec), &nullsp_basis);
-
-        for (unsigned int i=0; i<dim_nullsp; ++i)
-          nullsp_basis[i] = &(*nullspace[i]);
-
-        MatNullSpace  nullsp;
-
-        ierr = MatNullSpaceCreate (mpi_communicator, PETSC_FALSE, dim_nullsp, nullsp_basis, &nullsp);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        PetscBool  is_nullspace;
-
-        ierr = MatNullSpaceTest (nullsp, A,  &is_nullspace);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        ierr = KSPSetNullSpace (solver_data->ksp, nullsp);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-      };
-
-    ierr = KSPSolve (solver_data->ksp, b, x);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    KSPConvergedReason  reason;
-    ierr = KSPGetConvergedReason (solver_data->ksp, &reason);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    if (reason < 0)
-      {
-        std::ostringstream error_code;
-
-        if (reason == KSP_DIVERGED_NULL)
-          error_code << "KSP_DIVERGED_NULL";
-        else if (reason == KSP_DIVERGED_ITS)
-          error_code << "KSP_DIVERGED_ITS";
-        else if (reason == KSP_DIVERGED_ITS)
-          error_code << "KSP_DIVERGED_DTOL";
-        else if (reason == KSP_DIVERGED_BREAKDOWN)
-          error_code << "KSP_DIVERGED_BREAKDOWN";
-        else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
-          error_code << "KSP_DIVERGED_BREAKDOWN_BICG";
-        else if (reason == KSP_DIVERGED_NONSYMMETRIC)
-          error_code << "KSP_DIVERGED_NONSYMMETRIC";
-        else if (reason == KSP_DIVERGED_INDEFINITE_PC)
-          error_code << "KSP_DIVERGED_INDEFINITE_PC";
-        else if (reason == KSP_DIVERGED_NAN)
-          error_code << "KSP_DIVERGED_NAN";
-        else if (reason == KSP_DIVERGED_INDEFINITE_MAT)
-          error_code << "KSP_DIVERGED_INDEFINITE_MAT";
-        else
-          error_code << "Unknown Error";
-
-        AssertThrow (false, ExcPETScSolverError(error_code.str().c_str()));
-      };
-
-                                     // do not destroy solver object
-//    solver_data.reset ();
-
-                                     // in case of failure: throw
-                                     // exception
-    if (solver_control.last_check() != SolverControl::success)
-      throw SolverControl::NoConvergence (solver_control.last_step(),
-                                          solver_control.last_value());
-                                     // otherwise exit as normal
-  }
-
-
   void
   SolverBase::set_prefix(const std::string &prefix)
   {
@@ -434,9 +302,6 @@ namespace PETScWrappers
                                      // honor the initial guess in the
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
-
-    KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
-                     PETSC_DEFAULT, this->solver_control.max_steps()+1);
   }
 
 
@@ -539,9 +404,6 @@ namespace PETScWrappers
                                      // honor the initial guess in the
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
-
-    KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
-                     PETSC_DEFAULT, this->solver_control.max_steps()+1);
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.