step-45.h

Go to the documentation of this file.

 0,
                                  /*b_id2*/ 1,
                                  /*direction*/ 0,
                                  matched_pairs);
@endcode
would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$ for all
@f$y\in[0,1]@f$.
 
If we instead consider the parallelogram given by the convex hull of
@f$(0,0)@f$, @f$(1,1)@f$, @f$(1,2)@f$, @f$(0,1)@f$ we can achieve the constraints
@f$u(0,y)=u(1,y+1)@f$ 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
Here, again, the assignment of boundary indicators 0 and 1 stems from
what GridGenerator::parallelogram() documents.
 
The resulting @p matched_pairs can be used in
DoFTools::make_periodicity_constraints for populating an AffineConstraints
object 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,
                             affine_constraints,
                             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 affine_constraints.
 
 
<a name="problem"></a>
<a name="Apracticalexample"></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 @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ 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}
where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
@f{align*}
  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}
 
The mesh will be generated by GridGenerator::quarter_hyper_shell(),
which also documents how it assigns boundary indicators to its various
boundaries if its `colorize` argument is set to `true`.
 *
 *
 * <a name="CommProg"></a>
 * <h1> The commented program</h1>
 * 
 * This example program is a slight modification of @ref step_22 "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 @ref step_22 "step-22" and the full source code at the bottom.
 * 
 
 * 
 * In order to implement periodic boundary conditions only two functions
 * have to be modified:
 * - <code>StokesProblem<dim>::setup_dofs()</code>:
 * To populate an AffineConstraints object with periodicity constraints
 * - <code>StokesProblem<dim>::create_mesh()</code>:
 * To supply a distributed triangulation with periodicity information.
 * 
 
 * 
 * The rest of the program is identical to @ref step_22 "step-22", so let us skip this part
 * and only show these two functions in the following. (The full program can be
 * found in the "Plain program" section below, though.)
 * 
 
 * 
 * 
 
 * 
 *   
 
 * 
 * 
 * <a name="Settingupperiodicityconstraintsondistributedtriangulations"></a> 
 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
 * 
 * @code
 *   template <int dim>
 *   void StokesProblem<dim>::create_mesh()
 *   {
 *     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);
 * 
 * @endcode
 * 
 * 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 @f$\text{vertices}_2@f$ of a face on the left
 * boundary should be matched to the vertices @f$\text{vertices}_1@f$ of a face
 * on the lower boundary given by @f$\text{vertices}_2=R\cdot
 * \text{vertices}_1+b@f$ where the rotation matrix @f$R@f$ and the offset @f$b@f$ 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 resulting information into is here
 * based on the Triangulation.
 * 
 * @code
 *     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);
 * 
 * @endcode
 * 
 * Now telling the triangulation about the desired periodicity is
 * particularly easy by just calling
 * parallel::distributed::Triangulation::add_periodicity.
 * 
 * @code
 *     triangulation.add_periodicity(periodicity_vector);
 * 
 *     triangulation.refine_global(4 - dim);
 *   }
 * 
 * 
 *   template <int dim>
 *   void StokesProblem<dim>::setup_dofs()
 *   {
 *     dof_handler.distribute_dofs(fe);
 * 
 *     std::vector<unsigned int> block_component(dim + 1, 0);
 *     block_component[dim] = 1;
 *     DoFRenumbering::component_wise(dof_handler, block_component);
 * 
 *     const std::vector<types::global_dof_index> dofs_per_block =
 *       DoFTools::count_dofs_per_fe_block(dof_handler, 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();
 *       const IndexSet locally_relevant_dofs =
 *         DoFTools::extract_locally_relevant_dofs(dof_handler);
 *       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);
 * 
 *       const 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));
 * 
 * @endcode
 * 
 * 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 @f$\text{vertices}_2@f$ of a face on the left boundary should be
 * matched to the vertices
 * @f$\text{vertices}_1@f$ of a face on the lower boundary given by
 * @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ where the rotation
 * matrix @f$R@f$ and the offset @f$b@f$ 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
 * @f$\text{face}_2@f$ to
 * @f$\text{face}_1@f$.
 * 
 * @code
 *       FullMatrix<double> rotation_matrix(dim);
 *       rotation_matrix[0][1] = 1.;
 *       rotation_matrix[1][0] = -1.;
 * 
 *       Tensor<1, dim> offset;
 * 
 * @endcode
 * 
 * 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 @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ this is
 * exactly what we want.
 * 
 * @code
 *       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);
 * 
 * @endcode
 * 
 * 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:
 * 
 * @code
 *       std::vector<unsigned int> first_vector_components;
 *       first_vector_components.push_back(0);
 * 
 * @endcode
 * 
 * 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.
 * 
 * @code
 *       DoFTools::make_periodicity_constraints<dim, dim>(periodicity_vector,
 *                                                        constraints,
 *                                                        fe.component_mask(
 *                                                          velocities),
 *                                                        first_vector_components);
 *     }
 * 
 *     constraints.close();
 * 
 *     {
 *       TrilinosWrappers::BlockSparsityPattern bsp(owned_partitioning,
 *                                                  owned_partitioning,
 *                                                  relevant_partitioning,
 *                                                  mpi_communicator);
 * 
 *       Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
 *       for (unsigned int c = 0; c < dim + 1; ++c)
 *         for (unsigned int d = 0; d < dim + 1; ++d)
 *           if (!((c == dim) && (d == dim)))
 *             coupling[c][d] = DoFTools::always;
 *           else
 *             coupling[c][d] = DoFTools::none;
 * 
 *       DoFTools::make_sparsity_pattern(dof_handler,
 *                                       coupling,
 *                                       bsp,
 *                                       constraints,
 *                                       false,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         mpi_communicator));
 * 
 *       bsp.compress();
 * 
 *       system_matrix.reinit(bsp);
 *     }
 * 
 *     {
 *       TrilinosWrappers::BlockSparsityPattern preconditioner_bsp(
 *         owned_partitioning,
 *         owned_partitioning,
 *         relevant_partitioning,
 *         mpi_communicator);
 * 
 *       Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
 *       for (unsigned int c = 0; c < dim + 1; ++c)
 *         for (unsigned int d = 0; d < dim + 1; ++d)
 *           if ((c == dim) && (d == dim))
 *             preconditioner_coupling[c][d] = DoFTools::always;
 *           else
 *             preconditioner_coupling[c][d] = DoFTools::none;
 * 
 *       DoFTools::make_sparsity_pattern(dof_handler,
 *                                       preconditioner_coupling,
 *                                       preconditioner_bsp,
 *                                       constraints,
 *                                       false,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         mpi_communicator));
 * 
 *       preconditioner_bsp.compress();
 * 
 *       preconditioner_matrix.reinit(preconditioner_bsp);
 *     }
 * 
 *     system_rhs.reinit(owned_partitioning, mpi_communicator);
 *     solution.reinit(owned_partitioning,
 *                     relevant_partitioning,
 *                     mpi_communicator);
 *   }
 * 
 * @endcode
 * 
 * The rest of the program is then again identical to @ref step_22 "step-22". We will omit
 * it here now, but as before, you can find these parts in the "Plain program"
 * section below.
 * 
 
 * 
<a name="Results"></a><h1>Results</h1>
 
 
The created output is not very surprising. We simply see that the solution is
periodic with respect to the left and lower boundary:
 
<img src="https://www.dealii.org/images/steps/developer/step-45.periodic.png" alt="">
 
Without the periodicity constraints we would have ended up with the following solution:
 
<img src="https://www.dealii.org/images/steps/developer/step-45.non_periodic.png" alt="">
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-45.cc"
*/