/**
* Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5.
- * See Stoer: Einf�hrung 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>
/**
* Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7.
- * See Stoer: Einf�hrung 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>
-/**
- * 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 ------------- */
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);
*/
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
}
}
-
-
template <typename T>
void TableHandler::add_value (const std::string &key,
const T value)
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;
- };
-
/**
* initialize() before using this
* object.
*/
- PreconditionParaSails ();
+ PreconditionBoomerAMG ();
/**
* Constructor. Take the matrix which
* and additional flags if there are
* any.
*/
- PreconditionParaSails (const MatrixBase &matrix,
+ PreconditionBoomerAMG (const MatrixBase &matrix,
const AdditionalData &additional_data = AdditionalData());
/**
*/
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;
- };
}
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
<< "An error with error number " << arg1
<< " occurred while calling a PETSc function");
- DeclException1 (ExcPETScSolverError,
- char*,
- << "PETSc solver failed: " << arg1);
-
protected:
/**
}
-
-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
}
-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)
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());
}
}
-/* ----------------- 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::
}
- 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)
{
// 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);
}
// 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);
}
projects / dealii-svn.git / commitdiff