<br>
-<i>This program was contributed by Markus Bürg.</i>
-<br>
-
+<i>This program was contributed by Daniel Arndt and Matthias Maier.</i>
<a name="Intro"></a>
<h1>Introduction</h1>
-In this example we consider how to use periodic boundary conditions in
-deal.II. Periodic boundary conditions are a typical approach when one wants to
-solve some equation on a representative piece of a larger domain that repeats
-in one or more direction; an example is the simulation of the electronic
-structure of photonic
+In this example we present how to use periodic boundary conditions in
+deal.II. Periodic boundary conditions are algebraic constraints that
+typically occur in computations on representative regions of a larger
+domain that repeat in one or more directions.
+
+An example is the simulation of the electronic structure of photonic
crystals, because they have a lattice-like structure and, thus, it often
-suffices to do the actual computation on only one cell of the lattice. To be
-able to proceed this way one has to assume that the computation can be periodically
-extended to the other cells. This requires the solution to be periodic with respect
-to the cells. Hence the solution has to obtain the same nodal values on opposite
-parts of the boundary. In the figure below we show this
-concept in two space-dimensions. There, all dashed faces with the same color should
-have the same boundary values:
-
-<img src="http://www.dealii.org/images/steps/developer/step-45.periodic_cells.png" alt="">
-
-To keep things simple, in this tutorial we will consider an academic,
-simplified problem that allows us to focus on only that part that we are
-interested in here, namely how to set up periodic boundary
-conditions. Specifically, we solve the Poisson problem on a domain, where the
-left and right parts of the boundary are identified. Let $\Omega=(0,1)^2$ and
-consider the problem
-@f{align*}
- -\Delta u &=
- \cos(2\pi x)e^{-2x}\cos(2\pi y)e^{-2y} \qquad &&\text{in }\Omega
- \\
- u(x,0) &= 0 \qquad &&\text{for }x\in(0,1)\qquad &&\text{(bottom boundary)}
- \\
- u(x,1) &= 0 \qquad &&\text{for }x\in(0,1)\qquad &&\text{(top boundary)}
- \\
- u(0,y) &= u(1,y) \qquad &&\text{for }y\in(0,1)
- \qquad && \text{(left and right boundaries)}
+suffices to do the actual computation on only one box of the lattice. To
+be able to proceed this way one has to assume that the computation can be
+periodically extended to the other boxes; this requires the solution to
+have a periodic structure.
+
+<a name="Procedure"></a>
+<h1>Procedure</h1>
+
+deal.II provides a number of high level entry points to impose periodic
+boundary conditions.
+The general approach to apply periodic boundary conditions consists of
+three steps (see also the
+@ref GlossPeriodicConstraints "Glossary entry on periodic boundary conditions"):
+-# Create a mesh
+-# Identify those pairs of faces on different parts of the boundary across which
+ the solution should be symmetric, using GridTools::collect_periodic_faces()
+-# Add the periodicity information to the mesh
+ using parallel::distributed::Triangulation::add_periodicity()
+-# Add periodicity constraints using DoFTools::make_periodicity_constraints()
+
+The second and third step are necessary for distributed meshes
+to ensure that cells on opposite sides of the domain but connected by periodic
+faces are part of the ghost layer if one of them is stored on the local processor.
+If the Triangulation is not a parallel::distributed::Triangulation,
+these steps have to be omitted.
+
+The first step consists of collecting matching periodic faces and storing them in
+a <code>std::vector</code> of GridTools::PeriodicFacePair. This is done with the
+function GridTools::collect_periodic_faces() that can be invoked for example
+like this:
+@code
+GridTools::collect_periodic_faces(dof_handler,
+ b_id1,
+ b_id2,
+ direction,
+ matched_pairs,
+ offset = <default value>,
+ matrix = <default value>,
+ first_vector_components = <default value>);
+@endcode
+
+This call loops over all faces of the container dof_handler on the opposing
+boundaries with boundary indicator @p b_id1 and @p b_id2,
+respecitvely. If $\text{vertices}_{1/2}$ are the vertices of $\text{face}_{1/2}$,
+it matches pairs of faces (and dofs) such that the difference between $\text{vertices}_2$
+and $matrix\cdot \text{vertices}_1+\text{offset}$ vanishes in every component apart from direction
+and stores the resulting pairs with associated data in @p matched_pairs. (See
+GridTools::orthogonal_equality() for detailed information about the
+matching process.)
+
+Consider, for example, the colored unit square $\Omega=[0,1]^2$ with boundary
+indicator 0 on the left, 1 on the right, 2 on the bottom and 3 on the top
+faces. Then,
+@code
+GridTools::collect_periodic_faces(dof_handler,
+ /*b_id1*/ 0,
+ /*b_id2*/ 1,
+ /*direction*/ 0,
+ matched_pairs);
+@endcode
+would yield periodicity constraints such that $u(0,y)=u(1,y)$ for all
+$y\in[0,1]$.
+
+If we instead consider the parallelogram given by the convex hull of
+$(0,0)$, $(1,1)$, $(1,2)$, $(0,1)$ we can achieve the constraints
+$u(0,y)=u(1,y+1)$ by specifying an @p offset:
+@code
+GridTools::collect_periodic_faces(dof_handler,
+ /*b_id1*/ 0,
+ /*b_id2*/ 1,
+ /*direction*/ 0,
+ matched_pairs,
+ Tensor<1, 2>(0.,1.));
+@endcode
+or
+@code
+GridTools::collect_periodic_faces(dof_handler,
+ /*b_id1*/ 0,
+ /*b_id2*/ 1,
+ /*arbitrary direction*/ 0,
+ matched_pairs,
+ Tensor<1, 2>(1.,1.));
+@endcode
+
+The resulting @p matched_pairs can be used in
+DoFTools::make_periodicity_constraints for populating a ConstraintMatrix
+with periodicity constraints:
+@code
+DoFTools::make_periodicity_constraints(matched_pairs, constraints);
+@endcode
+
+Apart from this high level interface there are also variants of
+DoFTools::make_periodicity_constraints available that combine those two
+steps (see the variants of DofTools::make_periodicity_constraints).
+
+There is also a low level interface to
+DoFTools::make_periodicity_constraints if more flexibility is needed. The
+low level variant allows to directly specify two faces that shall be
+constrained:
+@code
+using namespace DoFTools;
+make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask = <default value>;
+ face_orientation = <default value>,
+ face_flip = <default value>,
+ face_rotation = <default value>,
+ matrix = <default value>);
+@endcode
+Here, we need to specify the orientation of the two faces using
+@p face_orientation, @p face_flip and @p face_orientation. For a closer description
+have a look at the documentation of DoFTools::make_periodicity_constraints.
+The remaining parameters are the same as for the high level interface apart
+from the self-explaining @p component_mask and @p constraint_matrix.
+
+<a name="problem"></a>
+<h1>A practical example</h1>
+
+In the following, we show how to use the above functions in a more involved
+example. The task is to enforce rotated periodicity constraints for the
+velocity component of a Stokes flow.
+
+On a quarter-circle defined by $\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}$ we are
+going to solve the Stokes problem
+@f{eqnarray*}
+ -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100*\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
+ -\textrm{div}\; \textbf{u}&=&0,\\
+ \textbf{u}|_{\Gamma_1}&=&{\bf 0},
@f}
-Note that the source term is not symmetric and so the solution would not be
-periodic unless this is explicitly enforced.
-
-The way one has to see these periodic boundary conditions $u(x,0) = u(x,1)$ is
-as follows: Assume for a moment (as we do in this program) that we have a
-uniformly refined mesh. Then, after discretization there are a number of nodes
-(degrees of freedom) with indices $i \in {\cal I}_l$ on the left boundary of
-the domain, and a second set of nodes at the right boundary $j \in {\cal
-I}_r$. Since we have assumed that the mesh is uniformly refined, there is
-exactly one node $j \in {\cal I}_r$ for each $i \in {\cal I}_l$ so that
-${\mathrm x}_j = {\mathrm x}_i + (1,0)^T$, i.e. the two of them match with
-respect to the periodicity. We will then write that $j=\text{periodic}(i)$
-(and, if you want, $i=\text{periodic}(j)$).
-If now $U_k, k=0,\ldots,N-1,$ are the unknowns of our discretized problem, then
-the periodic boundary condition boils down to the following set of
-constraints:
+where the boundary $\Gamma_1$ is defined as $\Gamma_1:=\{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}$.
+For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
@f{align*}
- U_{\text{periodic}(i)} = U_i, \qquad \forall i \in {\cal I}_l.
+ u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
+ u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
@f}
-Now, this is exactly the sort of constraint that the ConstraintMatrix class,
-first introduced in step-6,
-handles and can enforce in a linear system. Consequently, the main point of this
-program is how we fill the ConstraintMatrix object that stores these
-constraints, and how this is applied to the resulting linear system.
-
-The code for solving this problem is simple and based on step-3 since we want
-to focus on the implementation of the periodic boundary conditions. The code
-could be much more sophisticated, of course. For example, we could want to
-enforce periodic boundary conditions for adaptively refined meshes in which
-there is no longer a one-to-one relationship between degrees of freedom. We
-will discuss this at the end of the results section of this program.
/* ---------------------------------------------------------------------
*
- * Copyright (C) 2010 - 2015 by the deal.II authors
+ * Copyright (C) 2008 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* ---------------------------------------------------------------------
*
- * Author: Markus Buerg, University of Karlsruhe, 2010
+ * Author: Daniel Arndt, Matthias Maier, 2015
+ *
+ * Based on step-22 by Wolfgang Bangerth and Martin Kronbichler
*/
+// This example program is a slight modification of step-22 running in parallel
+// using Trilinos to demonstrate the usage of periodic boundary conditions in deal.II.
+// We thus omit to discuss the majority of the source code and only comment on the
+// parts that deal with periodicity constraints. For the rest have a look at step-22
+// and the full source code at the bottom.
-// @sect3{Include files}
+// In order to implement periodic boundary conditions only two functions
+// have to be modified:
+// - <code>StokesProblem<dim>::setup_dofs()</code>: To populate a ConstraintMatrix
+// object with periodicity constraints
+// - <code>StokesProblem<dim>::run()</code>: To supply a distributed triangulation with
+// periodicity information.
-// The include files are already known. The one critical for the current
-// program is the one that contains the ConstraintMatrix in the
-// <code>lac/</code> directory:
-#include <deal.II/base/function.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/precondition.h>
+// @cond SKIP
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+
#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/sparsity_pattern.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_tools.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_system.h>
-#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-
-#include <iostream>
-#include <fstream>
-
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
namespace Step45
{
using namespace dealii;
- // @sect3{The <code>LaplaceProblem</code> class}
-
- // The class <code>LaplaceProblem</code> is the main class of this
- // problem. As mentioned in the introduction, it is fashioned after the
- // corresponding class in step-3. Correspondingly, the documentation from
- // that tutorial program applies here as well. The only new member variable
- // is the <code>constraints</code> variables that will hold the constraints
- // from the periodic boundary condition. We will initialize it in the
- // <code>make_periodicity_constraints()</code> function which we call from
- // <code>make_grid_and_dofs()</code>.
- class LaplaceProblem
+ template <int dim>
+ class StokesProblem
{
public:
- LaplaceProblem ();
+ StokesProblem (const unsigned int degree);
void run ();
private:
- Triangulation<2> triangulation;
+ void create_mesh();
+ void setup_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results (const unsigned int refinement_cycle) const;
+ void refine_mesh ();
- FE_Q<2> fe;
- DoFHandler<2> dof_handler;
+ const unsigned int degree;
- ConstraintMatrix constraints;
+ MPI_Comm mpi_communicator;
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
- Vector<double> system_rhs;
- Vector<double> solution;
+ HyperShellBoundary<dim> boundary;
+ parallel::distributed::Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
- void assemble_system ();
- void output_results ();
- void make_grid_and_dofs ();
- void make_periodicity_constraints ();
- void solve ();
+ ConstraintMatrix constraints;
+ std::vector<IndexSet> owned_partitioning;
+ std::vector<IndexSet> relevant_partitioning;
+
+ TrilinosWrappers::BlockSparsityPattern sparsity_pattern;
+ TrilinosWrappers::BlockSparseMatrix system_matrix;
+
+ TrilinosWrappers::MPI::BlockVector solution;
+ TrilinosWrappers::MPI::BlockVector system_rhs;
+
+ ConditionalOStream pcout;
+
+ MappingQ<dim> mapping;
};
- // @sect3{The <code>RightHandSide</code> class}
- // The following implements the right hand side function discussed in the
- // introduction. Its implementation is obvious given what has been shown in
- // step-4 before:
- class RightHandSide: public Function<2>
+ template <int dim>
+ class BoundaryValues : public Function<dim>
{
public:
- RightHandSide ();
+ BoundaryValues () : Function<dim>(dim+1) {}
- virtual double value (const Point<2> &p,
- const unsigned int component = 0) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &value) const;
};
- RightHandSide::RightHandSide ()
- :
- Function<2> ()
- {}
+ template <int dim>
+ double
+ BoundaryValues<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int component) const
+ {
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+
+ return 0;
+ }
+
+
+ template <int dim>
+ void
+ BoundaryValues<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ {
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = BoundaryValues<dim>::value (p, c);
+ }
+
+
+ template <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>(dim+1) {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &value) const;
+
+ };
+ template <int dim>
double
- RightHandSide::value (const Point<2> &p,
- const unsigned int) const
+ RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
{
- return (std::cos (2 * numbers::PI * p(0)) *
- std::exp (- 2 * p(0)) *
- std::cos (2 * numbers::PI * p(1)) *
- std::exp (- 2 * p(1)));
+ const Point<dim> center(0.75, 0.1);
+ const double r = (p-center).norm();
+
+ if (component==0)
+ return std::exp(-100.*r*r);
+ return 0;
}
- // @sect3{Implementation of the <code>LaplaceProblem</code> class}
- // The first part of implementing the main class is the constructor. It is
- // unchanged from step-3 and step-4:
- LaplaceProblem::LaplaceProblem ()
+ template <int dim>
+ void
+ RightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ {
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = RightHandSide<dim>::value (p, c);
+ }
+
+
+
+ template <class Matrix, class Preconditioner>
+ class InverseMatrix : public Subscriptor
+ {
+ public:
+ InverseMatrix (const Matrix &m,
+ const Preconditioner &preconditioner,
+ const IndexSet &locally_owned,
+ const MPI_Comm &mpi_communicator);
+
+ void vmult (TrilinosWrappers::MPI::Vector &dst,
+ const TrilinosWrappers::MPI::Vector &src) const;
+
+ private:
+ const SmartPointer<const Matrix> matrix;
+ const SmartPointer<const Preconditioner> preconditioner;
+
+ const MPI_Comm *mpi_communicator;
+ mutable TrilinosWrappers::MPI::Vector tmp;
+ };
+
+
+
+ template <class Matrix, class Preconditioner>
+ InverseMatrix<Matrix,Preconditioner>::InverseMatrix
+ (const Matrix &m,
+ const Preconditioner &preconditioner,
+ const IndexSet &locally_owned,
+ const MPI_Comm &mpi_communicator)
:
- fe (1),
- dof_handler (triangulation)
+ matrix (&m),
+ preconditioner (&preconditioner),
+ mpi_communicator (&mpi_communicator),
+ tmp(locally_owned, mpi_communicator)
{}
- // @sect4{LaplaceProblem::make_grid_and_dofs}
-
- // The following is the first function to be called in
- // <code>run()</code>. It sets up the mesh and degrees of freedom.
- //
- // We start by creating the usual square mesh and changing the boundary
- // indicator on the parts of the boundary where we have Dirichlet boundary
- // conditions (top and bottom, i.e. faces two and three of the reference
- // cell as defined by GeometryInfo), so that we can distinguish between the
- // parts of the boundary where periodic and where Dirichlet boundary
- // conditions hold. We then refine the mesh a fixed number of times, with
- // child faces inheriting the boundary indicators previously set on the
- // coarse mesh from their parents.
- void LaplaceProblem::make_grid_and_dofs ()
+
+ template <class Matrix, class Preconditioner>
+ void InverseMatrix<Matrix,Preconditioner>::vmult
+ (TrilinosWrappers::MPI::Vector &dst,
+ const TrilinosWrappers::MPI::Vector &src) const
{
- GridGenerator::hyper_cube (triangulation);
- triangulation.begin_active ()->face (2)->set_boundary_id (1);
- triangulation.begin_active ()->face (3)->set_boundary_id (1);
- triangulation.refine_global (5);
+ SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
+ TrilinosWrappers::SolverCG cg (solver_control,
+ TrilinosWrappers::SolverCG::AdditionalData());
+
+ tmp = 0.;
+ cg.solve (*matrix, tmp, src, *preconditioner);
+ dst = tmp;
+ }
+
+
+
+ template <class Preconditioner>
+ class SchurComplement : public TrilinosWrappers::SparseMatrix
+ {
+ public:
+ SchurComplement ( const TrilinosWrappers::BlockSparseMatrix &system_matrix,
+ const InverseMatrix<TrilinosWrappers::SparseMatrix,
+ Preconditioner> &A_inverse,
+ const IndexSet &owned_pres,
+ const MPI_Comm &mpi_communicator);
+
+ void vmult (TrilinosWrappers::MPI::Vector &dst,
+ const TrilinosWrappers::MPI::Vector &src) const;
+
+ private:
+ const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> system_matrix;
+ const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
+ Preconditioner> > A_inverse;
+ mutable TrilinosWrappers::MPI::Vector tmp1, tmp2;
+ };
+
+
+
+ template <class Preconditioner>
+ SchurComplement<Preconditioner>::
+ SchurComplement (const TrilinosWrappers::BlockSparseMatrix &system_matrix,
+ const InverseMatrix<TrilinosWrappers::SparseMatrix,
+ Preconditioner> &A_inverse,
+ const IndexSet &owned_vel,
+ const MPI_Comm &mpi_communicator)
+ :
+ system_matrix (&system_matrix),
+ A_inverse (&A_inverse),
+ tmp1 (owned_vel, mpi_communicator),
+ tmp2 (tmp1)
+ {}
+
- // The next step is to distribute the degrees of freedom and produce a
- // little bit of graphical output:
- dof_handler.distribute_dofs (fe);
- std::cout << "Number of active cells: "
- << triangulation.n_active_cells ()
- << std::endl
- << "Degrees of freedom: " << dof_handler.n_dofs ()
- << std::endl;
-
- // Now it is the time for the constraints that come from the periodicity
- // constraints. We do this in the following, separate function, after
- // clearing any possible prior content from the constraints object:
- constraints.clear ();
- make_periodicity_constraints ();
-
- // We also incorporate the homogeneous Dirichlet boundary conditions on
- // the upper and lower parts of the boundary (i.e. the ones with boundary
- // indicator 1) and close the ConstraintMatrix object:
- VectorTools::interpolate_boundary_values (dof_handler, 1,
- ZeroFunction<2> (),
- constraints);
- constraints.close ();
- // Then we create the sparsity pattern and the system matrix and
- // initialize the solution and right-hand side vectors. This is again as
- // in step-3 or step-6, for example:
- DynamicSparsityPattern dsp (dof_handler.n_dofs(),
- dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler,
- dsp,
- constraints,
- false);
- dsp.compress ();
- sparsity_pattern.copy_from (dsp);
-
- system_matrix.reinit (sparsity_pattern);
- system_rhs.reinit (dof_handler.n_dofs());
- solution.reinit (dof_handler.n_dofs());
+ template <class Preconditioner>
+ void SchurComplement<Preconditioner>::vmult
+ (TrilinosWrappers::MPI::Vector &dst,
+ const TrilinosWrappers::MPI::Vector &src) const
+ {
+ system_matrix->block(0,1).vmult (tmp1, src);
+ A_inverse->vmult (tmp2, tmp1);
+ system_matrix->block(1,0).vmult (dst, tmp2);
}
- // @sect4{LaplaceProblem::make_periodicity_constraints}
-
- // This is the function that provides the new material of this tutorial
- // program. The general outline of the algorithm is as follows: we first
- // loop over all the degrees of freedom on the right boundary and record
- // their $y$-locations in a map together with their global indices. Then we
- // go along the left boundary, find matching $y$-locations for each degree
- // of freedom, and then add constraints that identify these matched degrees
- // of freedom.
- //
- // In this function, we make use of the fact that we have a scalar element
- // (i.e. the only valid vector component that can be passed to
- // DoFAccessor::vertex_dof_index is zero) and that we have a $Q_1$ element
- // for which all degrees of freedom live in the vertices of the
- // cell. Furthermore, we have assumed that we are in 2d and that meshes were
- // not refined adaptively — the latter assumption would imply that
- // there may be vertices that aren't matched one-to-one and for which we
- // won't be able to compute constraints this easily. We will discuss in the
- // "outlook" part of the results section below other strategies to write the
- // current function that can work in cases like this as well.
- void LaplaceProblem::make_periodicity_constraints ()
+ template <int dim>
+ StokesProblem<dim>::StokesProblem (const unsigned int degree)
+ :
+ degree (degree),
+ mpi_communicator (MPI_COMM_WORLD),
+ triangulation (mpi_communicator),
+ fe (FE_Q<dim>(degree+1), dim,
+ FE_Q<dim>(degree) , 1),
+ dof_handler (triangulation),
+ pcout (std::cout,
+ Utilities::MPI::this_mpi_process(mpi_communicator) == 0),
+ mapping(degree+1)
+ {}
+// @endcond
+//
+// @sect3{Setting up periodicity constraints on distributed triangulations}
+ template <int dim>
+ void StokesProblem<dim>::create_mesh()
{
- // To start with the actual implementation, we loop over all active cells
- // and check whether the cell is located at the right boundary (i.e. face
- // 1 — the one at the right end of the cell — is at the
- // boundary). If that is so, then we use that for the currently used
- // finite element, each degree of freedom of the face is located on one
- // vertex, and store their $y$-coordinate along with the global number of
- // this degree of freedom in the following map:
- std::map<unsigned int, double> dof_locations;
-
- for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
- cell != dof_handler.end (); ++cell)
- if (cell->at_boundary ()
- &&
- cell->face(1)->at_boundary ())
- {
- dof_locations[cell->face(1)->vertex_dof_index(0, 0)]
- = cell->face(1)->vertex(0)[1];
- dof_locations[cell->face(1)->vertex_dof_index(1, 0)]
- = cell->face(1)->vertex(1)[1];
- }
- // Note that in the above block, we add vertices zero and one of the
- // affected face to the map. This means that we will add each vertex
- // twice, once from each of the two adjacent cells (unless the vertex is a
- // corner of the domain). Since the coordinates of the vertex are the same
- // both times of course, there is no harm: we replace one value in the map
- // with itself the second time we visit an entry.
- //
- // The same will be true below where we add the same constraint twice to
- // the ConstraintMatrix — again, we will overwrite the constraint
- // with itself, and no harm is done.
-
- // Now we have to find the corresponding degrees of freedom on the left
- // part of the boundary. Therefore we loop over all cells again and choose
- // the ones where face 0 is at the boundary:
- for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active ();
- cell != dof_handler.end (); ++cell)
- if (cell->at_boundary ()
- &&
- cell->face (0)->at_boundary ())
- {
- // Every degree of freedom on this face needs to have a
- // corresponding one on the right side of the face, and our goal is
- // to add a constraint for the one on the left in terms of the one
- // on the right. To this end we first add a new line to the
- // constraint matrix for this one degree of freedom. Then we
- // identify it with the corresponding degree of freedom on the right
- // part of the boundary by constraining the degree of freedom on the
- // left with the one on the right times a weight of 1.0.
- //
- // Consequently, we loop over the two vertices of each face we find
- // and then loop over all the $y$-locations we've previously
- // recorded to find which degree of freedom on the right boundary
- // corresponds to the one we currently look at. Note that we have
- // entered these into a map, and when looping over the iterators
- // <code>p</code> of this map, <code>p-@>first</code> corresponds to
- // the "key" of an entry (the global number of the degree of
- // freedom), whereas <code>p-@>second</code> is the "value" (the
- // $y$-location we have entered above).
- //
- // We are quite sure here that we should be finding such a
- // corresponding degree of freedom. However, sometimes stuff happens
- // and so the bottom of the block contains an assertion that our
- // assumption was indeed correct and that a vertex was found.
- for (unsigned int face_vertex = 0; face_vertex<2; ++face_vertex)
- {
- constraints.add_line (cell->face(0)->vertex_dof_index (face_vertex, 0));
-
- std::map<unsigned int, double>::const_iterator p = dof_locations.begin();
- for (; p != dof_locations.end(); ++p)
- if (std::fabs(p->second - cell->face(0)->vertex(face_vertex)[1]) < 1e-8)
- {
- constraints.add_entry (cell->face(0)->vertex_dof_index (face_vertex, 0),
- p->first, 1.0);
- break;
- }
- Assert (p != dof_locations.end(),
- ExcMessage ("No corresponding degree of freedom was found!"));
- }
- }
+ Point<dim> center;
+ const double inner_radius = .5;
+ const double outer_radius = 1.;
+
+ GridGenerator::quarter_hyper_shell (triangulation,
+ center,
+ inner_radius,
+ outer_radius,
+ 0,
+ true);
+
+// Before we can prescribe periodicity constraints, we need to ensure that cells
+// on opposite sides of the domain but connected by periodic faces are part of
+// the ghost layer if one of them is stored on the local processor.
+// At this point we need to think about how we want to prescribe periodicity.
+// The vertices $\text{vertices}_2}$ of a face on the left boundary should be
+// matched to the vertices $\text{vertices}_1$ of a face on the lower boundary
+// given by $\text{vertices}_2=R\cdot \text{vertices}_1+b$ where the rotation
+// matrix $R$ and the offset $b$ are given by
+// @f{align*}
+// R=\begin{pmatrix}
+// 0&1\\-1&0
+// \end{pmatrix},
+// \quad
+// b=\begin{pmatrix}0&0\end{pmatrix}.
+// @f}
+// The data structure we are saving the reuslitng information into is here based
+// on the Triangulation.
+ std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
+ periodicity_vector;
+
+ FullMatrix<double> rotation_matrix(dim);
+ rotation_matrix[0][1]=1.;
+ rotation_matrix[1][0]=-1.;
+
+ GridTools::collect_periodic_faces(triangulation, 2, 3, 1,
+ periodicity_vector, Tensor<1, dim>(),
+ rotation_matrix);
+
+// Now telling the triangulation about the desired periodicity is
+// particularly easy by just calling
+// parallel::distributed::Triangulation::add_periodicity.
+ triangulation.add_periodicity(periodicity_vector);
+
+ triangulation.set_boundary(0, boundary);
+ triangulation.set_boundary(1, boundary);
+
+ triangulation.refine_global (4-dim);
}
+// @sect3{Setting up periodicity constraints on distributed triangulations}
+ template <int dim>
+ void StokesProblem<dim>::setup_dofs ()
+ {
+ dof_handler.distribute_dofs (fe);
- // @sect4{LaplaceProblem::assemble_system}
+ std::vector<unsigned int> block_component (dim+1,0);
+ block_component[dim] = 1;
+ DoFRenumbering::component_wise (dof_handler, block_component);
- // Assembling the system matrix and the right-hand side vector is done as in
- // other tutorials before.
- //
- // The only difference here is that we don't copy elements from local
- // contributions into the global matrix and later fix up constrained degrees
- // of freedom, but that we let the ConstraintMatrix do this job in one swoop
- // for us using the ConstraintMatrix::distribute_local_to_global
- // function(). This was previously already demonstrated in step-16, step-22,
- // for example, along with a discussion in the introduction of step-27.
- void LaplaceProblem::assemble_system ()
+ std::vector<types::global_dof_index> dofs_per_block (2);
+ DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
+ const unsigned int n_u = dofs_per_block[0],
+ n_p = dofs_per_block[1];
+
+ {
+ owned_partitioning.clear();
+ IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+ owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
+ owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u+n_p));
+
+ relevant_partitioning.clear();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
+ relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
+ relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u, n_u+n_p));
+
+ constraints.clear ();
+ constraints.reinit(locally_relevant_dofs);
+
+ FEValuesExtractors::Vector velocities(0);
+
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ VectorTools::interpolate_boundary_values (mapping,
+ dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ constraints,
+ fe.component_mask(velocities));
+ VectorTools::interpolate_boundary_values (mapping,
+ dof_handler,
+ 1,
+ BoundaryValues<dim>(),
+ constraints,
+ fe.component_mask(velocities));
+
+// After we provided the mesh with the necessary information for the periodicity
+// constraints, we are now able to actual create them. For describing the
+// matching we are using the same approach as before, i.e the $\text{vertices}_2}$
+// of a face on the left boundary should be matched to the vertices
+// $\text{vertices}_1$ of a face on the lower boundary given by
+// $\text{vertices}_2=R\cdot \text{vertices}_1+b$ where the rotation matrix $R$
+// and the offset $b$ are given by
+// @f{align*}
+// R=\begin{pmatrix}
+// 0&1\\-1&0
+// \end{pmatrix},
+// \quad
+// b=\begin{pmatrix}0&0\end{pmatrix}.
+// @f}
+// These two objects not only describe how faces should be matched but also
+// in which sense the solution should be transformed from $\text{face}_2$ to
+// $\text{face}_1$.
+ FullMatrix<double> rotation_matrix(dim);
+ rotation_matrix[0][1]=1.;
+ rotation_matrix[1][0]=-1.;
+
+ Tensor<1,dim> offset;
+
+// For setting up the constraints, we first store the periodicity
+// information in an auxiliary object of type
+// <code>std::vector@<GridTools::PeriodicFacePair<typename
+// DoFHandler@<dim@>::cell_iterator@> </code>. The periodic boundaries have the
+// boundary indicators 2 (x=0) and 3 (y=0). All the other parameters we
+// have set up before. In this case the direction does not matter. Due to
+// $\text{vertices}_2=R\cdot \text{vertices}_1+b$ this is exactly what we want.
+ std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
+ periodicity_vector;
+
+ const unsigned int direction = 1;
+
+ GridTools::collect_periodic_faces(dof_handler, 2, 3, direction,
+ periodicity_vector, offset,
+ rotation_matrix);
+
+// Next we need to provide information on which vector valued components of
+// the solution should be rotated. Since we choose here to just constraint the
+// velocity and this starts at the first component of the solution vector we
+// simply insert a 0:
+ std::vector<unsigned int> first_vector_components;
+ first_vector_components.push_back(0);
+
+// After setting up all the information in periodicity_vector all we have
+// to do is to tell make_periodicity_constraints to create the desired
+// constraints.
+ DoFTools::make_periodicity_constraints<DoFHandler<dim> >
+ (periodicity_vector, constraints, fe.component_mask(velocities),
+ first_vector_components);
+
+ VectorTools::interpolate_boundary_values (mapping,
+ dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ constraints,
+ fe.component_mask(velocities));
+ VectorTools::interpolate_boundary_values (mapping,
+ dof_handler,
+ 1,
+ BoundaryValues<dim>(),
+ constraints,
+ fe.component_mask(velocities));
+
+ }
+
+ constraints.close ();
+
+ {
+ TrilinosWrappers::BlockSparsityPattern bsp
+ (owned_partitioning, owned_partitioning,
+ relevant_partitioning, mpi_communicator);
+
+ DoFTools::make_sparsity_pattern
+ (dof_handler, bsp, constraints, false,
+ Utilities::MPI::this_mpi_process(mpi_communicator));
+
+ bsp.compress();
+
+ system_matrix.reinit (bsp);
+ }
+
+ system_rhs.reinit (owned_partitioning,
+ mpi_communicator);
+ solution.reinit (owned_partitioning, relevant_partitioning,
+ mpi_communicator);
+ }
+
+
+// @cond SKIP
+ template <int dim>
+ void StokesProblem<dim>::assemble_system ()
{
- QGauss<2> quadrature_formula(2);
- FEValues<2> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
+ system_matrix=0.;
+ system_rhs=0.;
+
+ QGauss<dim> quadrature_formula(degree+2);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.size();
+ FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values |
+ update_gradients);
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_rhs (dofs_per_cell);
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
- const RightHandSide right_hand_side;
+ const RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim+1));
+
+ const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+
+ std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
+ std::vector<double> div_phi_u (dofs_per_cell);
+ std::vector<double> phi_p (dofs_per_cell);
- DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
- {
- fe_values.reinit (cell);
- cell_matrix = 0;
- cell_rhs = 0;
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+ local_matrix = 0;
+ local_rhs = 0;
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+ rhs_values);
+
+ for (unsigned int q=0; q<n_q_points; ++q)
{
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
- fe_values.shape_grad (j, q_point) *
- fe_values.JxW (q_point));
-
- cell_rhs(i) += (fe_values.shape_value (i, q_point) *
- right_hand_side.value (fe_values.quadrature_point (q_point)) *
- fe_values.JxW (q_point));
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ {
+ symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+ div_phi_u[k] = fe_values[velocities].divergence (k, q);
+ phi_p[k] = fe_values[pressure].value (k, q);
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<=i; ++j)
+ {
+ local_matrix(i,j) += (symgrad_phi_u[i] * symgrad_phi_u[j]
+ - div_phi_u[i] * phi_p[j]
+ - phi_p[i] * div_phi_u[j]
+ + phi_p[i] * phi_p[j])
+ * fe_values.JxW(q);
+ }
+
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ local_rhs(i) += fe_values.shape_value(i,q) *
+ rhs_values[q](component_i) *
+ fe_values.JxW(q);
+ }
}
- cell->get_dof_indices (local_dof_indices);
- constraints.distribute_local_to_global (cell_matrix, cell_rhs,
- local_dof_indices,
- system_matrix, system_rhs);
- }
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+ local_matrix(i,j) = local_matrix(j,i);
+
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global (local_matrix, local_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+
+ system_matrix.compress (VectorOperation::add);
+ system_rhs.compress (VectorOperation::add);
+
+ pcout << " Computing preconditioner..." << std::endl << std::flush;
}
- // @sect4{LaplaceProblem::solve}
- // To solve the linear system of equations $Au=b$ we use the CG solver with
- // an SSOR-preconditioner. This is, again, copied almost verbatim from
- // step-6. As in step-6, we need to make sure that constrained degrees of
- // freedom get their correct values after solving by calling the
- // ConstraintMatrix::distribute function:
- void LaplaceProblem::solve ()
+ template <int dim>
+ void StokesProblem<dim>::solve ()
{
- SolverControl solver_control (dof_handler.n_dofs (), 1e-12);
- PreconditionSSOR<SparseMatrix<double> > precondition;
+ TrilinosWrappers::PreconditionJacobi A_preconditioner;
+ A_preconditioner.initialize(system_matrix.block(0,0));
- precondition.initialize (system_matrix);
+ const InverseMatrix<TrilinosWrappers::SparseMatrix,
+ TrilinosWrappers::PreconditionJacobi>
+ A_inverse (system_matrix.block(0,0),
+ A_preconditioner,
+ owned_partitioning[0],
+ mpi_communicator);
- SolverCG<> cg (solver_control);
+ TrilinosWrappers::MPI::BlockVector tmp (owned_partitioning,
+ mpi_communicator);
- cg.solve (system_matrix, solution, system_rhs, precondition);
- constraints.distribute (solution);
+ {
+ TrilinosWrappers::MPI::Vector schur_rhs (owned_partitioning[1],
+ mpi_communicator);
+ A_inverse.vmult (tmp.block(0), system_rhs.block(0));
+ system_matrix.block(1,0).vmult (schur_rhs, tmp.block(0));
+ schur_rhs -= system_rhs.block(1);
+
+ SchurComplement<TrilinosWrappers::PreconditionJacobi>
+ schur_complement (system_matrix, A_inverse,
+ owned_partitioning[0],
+ mpi_communicator);
+
+ SolverControl solver_control (solution.block(1).size(),
+ 1e-6*schur_rhs.l2_norm());
+ SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
+
+ TrilinosWrappers::PreconditionAMG preconditioner;
+ preconditioner.initialize (system_matrix.block(1,1));
+
+ InverseMatrix<TrilinosWrappers::SparseMatrix,
+ TrilinosWrappers::PreconditionAMG>
+ m_inverse (system_matrix.block(1,1), preconditioner,
+ owned_partitioning[1], mpi_communicator);
+
+ cg.solve (schur_complement,
+ tmp.block(1),
+ schur_rhs,
+ preconditioner);
+
+ constraints.distribute (tmp);
+ solution.block(1)=tmp.block(1);
+ }
+
+ {
+ system_matrix.block(0,1).vmult (tmp.block(0), tmp.block(1));
+ tmp.block(0) *= -1;
+ tmp.block(0) += system_rhs.block(0);
+
+ A_inverse.vmult (tmp.block(0), tmp.block(0));
+
+ constraints.distribute (tmp);
+ solution.block(0)=tmp.block(0);
+ }
}
- // @sect4{LaplaceProblem::output_results}
- // This is another function copied from previous tutorial programs. It
- // generates graphical output in VTK format:
- void LaplaceProblem::output_results ()
+ template <int dim>
+ void
+ StokesProblem<dim>::output_results (const unsigned int refinement_cycle) const
{
- DataOut<2> data_out;
+ std::vector<std::string> solution_names (dim, "velocity");
+ solution_names.push_back ("pressure");
+
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation
+ (dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_component_interpretation
+ .push_back (DataComponentInterpretation::component_is_scalar);
+ DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (solution, "u");
- data_out.build_patches ();
+ data_out.add_data_vector (solution, solution_names,
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ Vector<float> subdomain (triangulation.n_active_cells());
+ for (unsigned int i=0; i<subdomain.size(); ++i)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+ data_out.add_data_vector (subdomain, "subdomain");
+ data_out.build_patches (mapping, degree+1);
+
+ std::ostringstream filename;
+ filename << "solution-"
+ << Utilities::int_to_string (refinement_cycle, 2)
+ << "."
+ << Utilities::int_to_string (triangulation.locally_owned_subdomain(),2)
+ << ".vtu";
+
+ std::ofstream output (filename.str().c_str());
+ data_out.write_vtu (output);
+
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ {
+ std::vector<std::string> filenames;
+ for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+ filenames.push_back (std::string("solution-") +
+ Utilities::int_to_string (refinement_cycle, 2) +
+ "." +
+ Utilities::int_to_string(i, 2) +
+ ".vtu");
+ const std::string
+ pvtu_master_filename = ("solution-" +
+ Utilities::int_to_string (refinement_cycle, 2) +
+ ".pvtu");
+ std::ofstream pvtu_master (pvtu_master_filename.c_str());
+ data_out.write_pvtu_record (pvtu_master, filenames);
+ }
+ }
- std::ofstream output ("solution.vtk");
- data_out.write_vtk (output);
- }
+ template <int dim>
+ void
+ StokesProblem<dim>::refine_mesh ()
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ FEValuesExtractors::Scalar pressure(dim);
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1>(degree+1),
+ typename FunctionMap<dim>::type(),
+ solution,
+ estimated_error_per_cell,
+ fe.component_mask(pressure));
+
+ parallel::distributed::GridRefinement::
+ refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.0);
+ triangulation.execute_coarsening_and_refinement ();
+ }
- // @sect4{LaplaceProblem::run}
- // And another function copied from previous programs:
- void LaplaceProblem::run ()
+ template <int dim>
+ void StokesProblem<dim>::run ()
{
- make_grid_and_dofs();
- assemble_system ();
- solve ();
- output_results ();
+ create_mesh();
+
+ for (unsigned int refinement_cycle = 0; refinement_cycle<9;
+ ++refinement_cycle)
+ {
+ pcout << "Refinement cycle " << refinement_cycle << std::endl;
+
+ if (refinement_cycle > 0)
+ refine_mesh ();
+
+ setup_dofs ();
+
+ pcout << " Assembling..." << std::endl << std::flush;
+ assemble_system ();
+
+ pcout << " Solving..." << std::flush;
+ solve ();
+
+ output_results (refinement_cycle);
+
+ pcout << std::endl;
+ }
}
}
-// @sect3{The <code>main</code> function}
-// And at the end we have the main function as usual, this time copied from
-// step-6:
-int main ()
+int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step45;
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
deallog.depth_console (0);
- LaplaceProblem laplace_problem;
- laplace_problem.run ();
+ {
+ StokesProblem<2> flow_problem(1);
+ flow_problem.run ();
+ }
}
catch (std::exception &exc)
{
return 0;
}
+// @endcond