From 11893c3b76883a1a902eed8afba43296b00753eb Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 4 Aug 2015 17:39:13 -0500 Subject: [PATCH] Redesign step-45 to support existing functions for periodic boundary conditions --- doc/doxygen/headers/glossary.h | 4 +- examples/step-45/doc/intro.dox | 203 ++++-- examples/step-45/doc/results.dox | 155 +---- examples/step-45/step-45.cc | 954 ++++++++++++++++++--------- include/deal.II/distributed/tria.h | 4 + include/deal.II/dofs/dof_tools.h | 2 +- include/deal.II/grid/grid_tools.h | 2 +- source/dofs/dof_tools_constraints.cc | 4 +- source/grid/grid_tools.cc | 2 +- 9 files changed, 807 insertions(+), 523 deletions(-) diff --git a/doc/doxygen/headers/glossary.h b/doc/doxygen/headers/glossary.h index 219d825a21..6ec0bba515 100644 --- a/doc/doxygen/headers/glossary.h +++ b/doc/doxygen/headers/glossary.h @@ -1233,7 +1233,9 @@ * 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. + * * *
@anchor GlossPrimitive Primitive finite * elements
diff --git a/examples/step-45/doc/intro.dox b/examples/step-45/doc/intro.dox index f878982f21..475c33de26 100644 --- a/examples/step-45/doc/intro.dox +++ b/examples/step-45/doc/intro.dox @@ -1,72 +1,149 @@
-This program was contributed by Markus Bürg. -
- +This program was contributed by Daniel Arndt and Matthias Maier.

Introduction

-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: - - - -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. + + +

Procedure

+ +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 std::vector 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 = , + matrix = , + first_vector_components = ); +@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 = ; + face_orientation = , + face_flip = , + face_rotation = , + matrix = ); +@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 practical example

+ +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. diff --git a/examples/step-45/doc/results.dox b/examples/step-45/doc/results.dox index 4affd89176..900ae319cf 100644 --- a/examples/step-45/doc/results.dox +++ b/examples/step-45/doc/results.dox @@ -1,155 +1,10 @@

Results

-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 + -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: - - -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. - - - -

Possibilities for extensions

- -The function LaplaceProblem::make_periodicity_constraints 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 -make_periodicity_constraint_recursively() function — an -implementation detail, not an external interface — put into an anonymous -namespace): -@code -namespace -{ - template - void - make_periodicity_constraints_recursively - (const typename DoFHandler::face_iterator &face_1, - const typename DoFHandler::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; cn_children(); ++c) - make_periodicity_constraints_recursively (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 local_dof_indices_1 (dofs_per_face); - face_1->get_dof_indices (local_dof_indices_1); - - std::vector local_dof_indices_2 (dofs_per_face); - face_2->get_dof_indices (local_dof_indices_2); - - for (unsigned int i=0; i (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 -LaplaceProblem::make_periodicity_constraints() 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. + diff --git a/examples/step-45/step-45.cc b/examples/step-45/step-45.cc index 4499e22af6..4fe52954ab 100644 --- a/examples/step-45/step-45.cc +++ b/examples/step-45/step-45.cc @@ -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. * @@ -14,417 +14,762 @@ * --------------------------------------------------------------------- * - * 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: +// - StokesProblem::setup_dofs(): To populate a ConstraintMatrix +// object with periodicity constraints +// - StokesProblem::run(): 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 -// lac/ directory: -#include -#include -#include -#include +// @cond SKIP +#include + +#include + #include -#include -#include -#include -#include +#include + +#include +#include +#include +#include +#include #include -#include +#include +#include -#include -#include +#include #include #include -#include +#include -#include #include - -#include -#include - +#include +#include namespace Step45 { using namespace dealii; - // @sect3{The LaplaceProblem class} - - // The class LaplaceProblem 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 constraints variables that will hold the constraints - // from the periodic boundary condition. We will initialize it in the - // make_periodicity_constraints() function which we call from - // make_grid_and_dofs(). - class LaplaceProblem + template + 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 system_matrix; - Vector system_rhs; - Vector solution; + HyperShellBoundary boundary; + parallel::distributed::Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; - void assemble_system (); - void output_results (); - void make_grid_and_dofs (); - void make_periodicity_constraints (); - void solve (); + ConstraintMatrix constraints; + std::vector owned_partitioning; + std::vector relevant_partitioning; + + TrilinosWrappers::BlockSparsityPattern sparsity_pattern; + TrilinosWrappers::BlockSparseMatrix system_matrix; + + TrilinosWrappers::MPI::BlockVector solution; + TrilinosWrappers::MPI::BlockVector system_rhs; + + ConditionalOStream pcout; + + MappingQ mapping; }; - // @sect3{The RightHandSide 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 + class BoundaryValues : public Function { public: - RightHandSide (); + BoundaryValues () : Function(dim+1) {} - virtual double value (const Point<2> &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; }; - RightHandSide::RightHandSide () - : - Function<2> () - {} + template + double + BoundaryValues::value (const Point &/*p*/, + const unsigned int component) const + { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + + return 0; + } + + + template + void + BoundaryValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = BoundaryValues::value (p, c); + } + + + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function(dim+1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + + }; + template double - RightHandSide::value (const Point<2> &p, - const unsigned int) const + RightHandSide::value (const Point &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 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 LaplaceProblem class} - // The first part of implementing the main class is the constructor. It is - // unchanged from step-3 and step-4: - LaplaceProblem::LaplaceProblem () + template + void + RightHandSide::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = RightHandSide::value (p, c); + } + + + + template + 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 matrix; + const SmartPointer preconditioner; + + const MPI_Comm *mpi_communicator; + mutable TrilinosWrappers::MPI::Vector tmp; + }; + + + + template + InverseMatrix::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 - // run(). 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 + void InverseMatrix::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 SchurComplement : public TrilinosWrappers::SparseMatrix + { + public: + SchurComplement ( const TrilinosWrappers::BlockSparseMatrix &system_matrix, + const InverseMatrix &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 system_matrix; + const SmartPointer > A_inverse; + mutable TrilinosWrappers::MPI::Vector tmp1, tmp2; + }; + + + + template + SchurComplement:: + SchurComplement (const TrilinosWrappers::BlockSparseMatrix &system_matrix, + const InverseMatrix &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 + void SchurComplement::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 + StokesProblem::StokesProblem (const unsigned int degree) + : + degree (degree), + mpi_communicator (MPI_COMM_WORLD), + triangulation (mpi_communicator), + fe (FE_Q(degree+1), dim, + FE_Q(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 + void StokesProblem::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 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 - // p of this map, p-@>first corresponds to - // the "key" of an entry (the global number of the degree of - // freedom), whereas p-@>second 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::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 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::cell_iterator> > + periodicity_vector; + + FullMatrix 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 + void StokesProblem::setup_dofs () + { + dof_handler.distribute_dofs (fe); - // @sect4{LaplaceProblem::assemble_system} + std::vector 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 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(), + constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values (mapping, + dof_handler, + 1, + BoundaryValues(), + 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 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 +// std::vector@::cell_iterator@> . 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::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 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 > + (periodicity_vector, constraints, fe.component_mask(velocities), + first_vector_components); + + VectorTools::interpolate_boundary_values (mapping, + dof_handler, + 0, + BoundaryValues(), + constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values (mapping, + dof_handler, + 1, + BoundaryValues(), + 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 + void StokesProblem::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 quadrature_formula(degree+2); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + FEValues fe_values (mapping, fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector 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 local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - const RightHandSide right_hand_side; + const RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim+1)); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); - DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::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_pointget_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; iget_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 + void StokesProblem::solve () { - SolverControl solver_control (dof_handler.n_dofs (), 1e-12); - PreconditionSSOR > precondition; + TrilinosWrappers::PreconditionJacobi A_preconditioner; + A_preconditioner.initialize(system_matrix.block(0,0)); - precondition.initialize (system_matrix); + const InverseMatrix + 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 + 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 cg(solver_control); + + TrilinosWrappers::PreconditionAMG preconditioner; + preconditioner.initialize (system_matrix.block(1,1)); + + InverseMatrix + 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 + void + StokesProblem::output_results (const unsigned int refinement_cycle) const { - DataOut<2> data_out; + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + DataOut 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::type_dof_data, + data_component_interpretation); + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i + void + StokesProblem::refine_mesh () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + FEValuesExtractors::Scalar pressure(dim); + KellyErrorEstimator::estimate (dof_handler, + QGauss(degree+1), + typename FunctionMap::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 + void StokesProblem::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 main 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 diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 7638ab9ade..291101d333 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -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 diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index f87c9a47a2..e27154340e 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -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 */ diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index f154b32673..c8909f8d22 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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 diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index e3b26a5dfd..7859ed3ab5 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -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 */ template -- 2.39.5