]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Redesign step-45 to support existing functions for periodic boundary conditions 1421/head
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Tue, 4 Aug 2015 22:39:13 +0000 (17:39 -0500)
committerMatthias Maier <tamiko@43-1.org>
Thu, 15 Oct 2015 16:45:32 +0000 (11:45 -0500)
doc/doxygen/headers/glossary.h
examples/step-45/doc/intro.dox
examples/step-45/doc/results.dox
examples/step-45/step-45.cc
include/deal.II/distributed/tria.h
include/deal.II/dofs/dof_tools.h
include/deal.II/grid/grid_tools.h
source/dofs/dof_tools_constraints.cc
source/grid/grid_tools.cc

index 219d825a21d68928c8e9ac49e113c4e03e9b17f0..6ec0bba5158ef97be16406a50363fa0e0d2dc1a1 100644 (file)
  * using parallel::distributed::Triangulation::add_periodicity()
  * -# Gather the periodic faces using GridTools::collect_periodic_faces() (DoFHandler)
  * -# Add periodicity constraints using DoFTools::make_periodicity_constraints()
- *
+ * 
+ * An example for this can be found in step-45. 
+ * 
  *
  * <dt class="glossary">@anchor GlossPrimitive <b>Primitive finite
  * elements</b></dt>
index f878982f21c448b229d2d5b8a85953131fbedb06..475c33de267788fa51b575fff05cf6471c1a2443 100644 (file)
 <br>
 
-<i>This program was contributed by Markus B&uuml;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.
index 4affd891768da72dd3ca3db1e3940cc711ee20d1..900ae319cf7c4e394b4e7444de27b0f6413030f9 100644 (file)
 <h1>Results</h1>
 
-The textual output the program generates is not very surprising. It just
-prints out the usual information on degrees of freedom and active cells,
-in much the same way as step-3 did:
+The created output is not very surprising. We simply see that the solution is
+periodic with respect to the left and lower boundary:
 
-@code
-Number of active cells: 1024
-Degrees of freedom: 1089
-@endcode
+<img src="http://www.dealii.org/images/steps/developer/step-45.periodic.png" alt="">
 
-The plot of the solution can be found in the figure below. We can see that the
-solution is constant zero on the upper and the lower part of the boundary as
-required by the homogeneous Dirichlet boundary conditions. On the left and
-right parts the values coincide with each other, just as we wanted:
+Without the periodicity constraints we would have ended up with the following solution:
 
-<img src="http://www.dealii.org/images/steps/developer/step-45.solution.png" alt="">
-
-Note also that the solution is clearly not left-right symmetric and so would
-not likely have been periodic had we prescribed, for example, homogeneous
-Neumann boundary condition. However, it is periodic thanks to the constraints
-imposed.
-
-
-<a name="extensions"></a>
-<h3>Possibilities for extensions</h3>
-
-The function <code>LaplaceProblem::make_periodicity_constraints</code> is
-relatively simple in that it just matches the location of degrees of
-freedom. This makes it flexible when the periodicity boundary conditions are
-posed not just on opposite faces of the unit rectangle but on separate parts
-of a possibly more complicated domain. Or, if one wanted to "twist" the
-boundary condition by prescribing, for example,
-@f{align*}
-   u(0,y) &= u(1,1-y) \qquad &\text{for }y\in(0,1).
-@f}
-
-On the other hand, the function is somewhat limited by the assumption that the
-domain is two-dimensional and that we only use $Q_1$ elements. The former
-assumption is easily lifted by looping over all four vertices of a face in 3d,
-but the latter is somewhat more complicated to lift because we have assumed
-that degrees of freedom are only located in vertices. In the following,
-therefore, let us describe a function that computes the same constraints but
-in a dimension-independent way and for any finite element one may want to
-consider.
-
-@note The discussion below is meant as an explanation of how one might
-approach this kind of problem. If you need this functionality, you may also
-want to take a look at DoFTools::make_periodicity_constraints() that already
-provides some of it.
-
-The idea is to work recursively on pairs of faces. For example, let us start
-with the left and right face of the (single) coarse mesh cell. They need to
-match, but they are not active (i.e. they are further refined) and so there
-are no degrees of freedom associated with these faces. However, if the two
-current faces are periodic, then so are the zeroth children of the two as well
-as the respective first children, etc. We can then in turn work on each of
-these pairs of faces. If they are not active, we may recurse further into this
-refinement tree until we find a pair of active faces. In that case, we enter
-equivalences between matching degrees of freedom on the two active faces.
-
-An implementation of this idea would look like follows (with the
-<code>make_periodicity_constraint_recursively()</code> function &mdash; an
-implementation detail, not an external interface &mdash; put into an anonymous
-namespace):
-@code
-namespace
-{
-  template <int dim>
-  void
-  make_periodicity_constraints_recursively
-  (const typename DoFHandler<dim>::face_iterator &face_1,
-   const typename DoFHandler<dim>::face_iterator &face_2,
-   ConstraintMatrix &constraints)
-  {
-    Assert (face_1->n_children() == face_2->n_children(),
-           ExcNotImplemented());
-    if (face_1->has_children())
-      {
-       for (unsigned int c=0; c<face_1->n_children(); ++c)
-         make_periodicity_constraints_recursively<dim> (face_1->child(c),
-                                                        face_2->child(c),
-                                                        constraints);
-      }
-    else
-      {
-       const unsigned int dofs_per_face
-         = face_1->get_dof_handler().get_fe().dofs_per_face;
-
-       std::vector<unsigned int> local_dof_indices_1 (dofs_per_face);
-       face_1->get_dof_indices (local_dof_indices_1);
-
-       std::vector<unsigned int> local_dof_indices_2 (dofs_per_face);
-       face_2->get_dof_indices (local_dof_indices_2);
-
-       for (unsigned int i=0; i<dofs_per_face; ++i)
-         {
-           constraints.add_line (local_dof_indices_1[i]);
-           constraints.add_entry (local_dof_indices_1[i],
-                                  local_dof_indices_2[i],
-                                  1.0);
-         }
-      }
-  }
-}
-
-
-void LaplaceProblem::make_periodicity_constraints ()
-{
-  make_periodicity_constraints_recursively<2> (dof_handler.begin(0)->face(0),
-                                               dof_handler.begin(0)->face(1),
-                                               constraints);
-}
-@endcode
-
-The implementation of the recursive function should be mostly self explanatory
-given the discussion above. The
-<code>LaplaceProblem::make_periodicity_constraints()</code> function simply
-calls the former with matching faces of the first (and only) coarse mesh cell
-on refinement level 0. Note that when calling the recursive function we have
-to explicitly specify the template argument since the compiler can not deduce
-it (the template argument is only used in a "non-deducible context").
-
-This function is now dimension and finite element independent, but it still
-has the restriction that it assumes that the mesh is uniformly refined (or, in
-fact, only that matching periodic faces are refined equally). We check this at
-the beginning by asserting that both faces have the same number of children
-(that includes that neither have any children, i.e. that both are active).
-On the other hand, the function above can be extended to also allow this sort
-of thing. In that case, if we encounter the situation that only one cell is
-refined, we would have to recurse into its children and interpolate their
-degrees of freedom with respect to the degrees of freedom to the coarser
-matching face. This can use the same facilities the finite element classes
-already provide for computing constraints based on hanging nodes. We leave
-implementing this as an exercise, however.
-
-@note The functions above make one assumption, namely that the degrees of
-freedom on one face match one-to-one to the corresponding other face. This is
-sometimes difficult to establish, especially in 3d. For example, consider a
-long string of cube cells where we want to match the far left face with the
-far right face for periodicity. If the cells are all undistorted cubes, then
-everything will work as expected. But imagine we have twisted our cells so
-that the string as a whole now has a 90 degree twist; in that case, the
-coordinate systems of the far left and far right face are also rotated
-relative to each other, and the first DoF on the far left face will no longer
-be at the same location as the first DoF on the far right face. To make things
-a bit worse, the 3d case also allows for pathological cases where mesh cells
-are no longer orientable in the standard order (see the
-@ref GlossFaceOrientation "Face orientation" glossary entry), making matching
-coordinate systems difficult. Ultimately, whether you will encounter these
-cases depends on the kind of mesh you have: if your coarse mesh is just the
-unit cube or a subdivided hyper rectangle, you are definitely on the safe
-side. In other cases, if in doubt, verify the locations of degrees of freedom
-by printing the kind of information we have used in the tutorial's own
-implementation.
+<img src="http://www.dealii.org/images/steps/developer/step-45.non_periodic.png" alt="">
index 4499e22af62952e76bf26693d95c70363dc7d661..4fe52954ab71ce343cbd89df7ab3130f2ff1094c 100644 (file)
@@ -1,6 +1,6 @@
 /* ---------------------------------------------------------------------
  *
- * 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 &mdash; 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 &mdash; the one at the right end of the cell &mdash; 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 &mdash; 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)
     {
@@ -453,3 +798,4 @@ int main ()
 
   return 0;
 }
+// @endcond
index 7638ab9adeec1fcb94bc74bc0db20d1f25b1e8a1..291101d33349e8a9cc3a03402abfdc60ae3e343f 100644 (file)
@@ -809,6 +809,10 @@ namespace parallel
        * The vector can be filled by the function
        * GridTools::collect_periodic_faces.
        *
+       * For more information on periodic boundary conditions see
+       * GridTools::collect_periodic_faces, DoFTools::make_periodicity_constraints
+       * and step-45.
+       *
        * @note Before this function can be used the Triangulation has to be
        * initialized and must not be refined. Calling this function more than
        * once is possible, but not recommended: The function destroys and
index f87c9a47a2ce93d79f77d13cc89564e6ef7a7f28..e27154340e3c1ffbdb5d1835ec29e669f5e6020c 100644 (file)
@@ -947,7 +947,7 @@ namespace DoFTools
    *
    * @see
    * @ref GlossPeriodicConstraints "Glossary entry on periodic boundary conditions"
-   * for further information.
+   * and step-45 for further information.
    *
    * @author Daniel Arndt, Matthias Maier, 2013 - 2015
    */
index f154b32673b1450d370ff0f6fdc44d6b99e75577..c8909f8d22e991333b6c60c518be585bfba7cddb 100644 (file)
@@ -1183,7 +1183,7 @@ namespace GridTools
    * (For more details see DoFTools::make_periodicity_constraints(), the
    * glossary
    * @ref GlossPeriodicConstraints "glossary entry on periodic conditions"
-   * and @ref step_45 "step-45"). Second, @p matrix will be stored in the
+   * and step-45). Second, @p matrix will be stored in the
    * PeriodicFacePair collection @p matched_pairs for further use.
    *
    * @tparam Container A type that satisfies the requirements of a mesh
index e3b26a5dfde7d27fa061f39ed1c7be09714b4e87..7859ed3ab5ee8d50cbfcdebf9453078e1aada6f1 100644 (file)
@@ -1946,8 +1946,8 @@ namespace DoFTools
                   unsigned int inverse_constraint_target = numbers::invalid_unsigned_int;
                   if (is_inverse_constrained)
                     for (unsigned int jj=0; jj<dofs_per_face; ++jj)
-                      if (((std::abs(transformation(i,jj)) < eps)) ||
-                          (std::abs(transformation(i,jj)+1) < eps) == false)
+                      if (((std::abs(transformation(i,jj)) < eps) ||
+                           (std::abs(transformation(i,jj)+1) < eps)) == false)
                         {
                           is_inverse_constrained = false;
                           break;
index 72ecb6e84ccbd09202e086fd75958a8973a5a0e2..af5fecc92be68999feeb4583695d4efa30c3d228 100644 (file)
@@ -2926,7 +2926,7 @@ next_cell:
    * An orthogonal equality test for points:
    *
    * point1 and point2 are considered equal, if
-   *   matrix.(point1 + offset) - point2
+   *   matrix.point1 + offset - point2
    * is parallel to the unit vector in <direction>
    */
   template<int spacedim>

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.