Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-45.h
Go to the documentation of this file.
1  0,
108  /*b_id2*/ 1,
109  /*direction*/ 0,
110  matched_pairs);
111 @endcode
112 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$ for all
113 @f$y\in[0,1]@f$.
114 
115 If we instead consider the parallelogram given by the convex hull of
116 @f$(0,0)@f$, @f$(1,1)@f$, @f$(1,2)@f$, @f$(0,1)@f$ we can achieve the constraints
117 @f$u(0,y)=u(1,y+1)@f$ by specifying an @p offset:
118 @code
120  /*b_id1*/ 0,
121  /*b_id2*/ 1,
122  /*direction*/ 0,
123  matched_pairs,
124  Tensor<1, 2>(0.,1.));
125 @endcode
126 or
127 @code
129  /*b_id1*/ 0,
130  /*b_id2*/ 1,
131  /*arbitrary direction*/ 0,
132  matched_pairs,
133  Tensor<1, 2>(1.,1.));
134 @endcode
135 Here, again, the assignment of boundary indicators 0 and 1 stems from
136 what GridGenerator::parallelogram() documents.
137 
138 The resulting @p matched_pairs can be used in
140 object with periodicity constraints:
141 @code
142 DoFTools::make_periodicity_constraints(matched_pairs, constraints);
143 @endcode
144 
145 Apart from this high level interface there are also variants of
146 DoFTools::make_periodicity_constraints available that combine those two
147 steps (see the variants of DofTools::make_periodicity_constraints).
148 
149 There is also a low level interface to
150 DoFTools::make_periodicity_constraints if more flexibility is needed. The
151 low level variant allows to directly specify two faces that shall be
152 constrained:
153 @code
154 using namespace DoFTools;
156  face_2,
157  affine_constraints,
158  component_mask = <default value>;
159  face_orientation = <default value>,
160  face_flip = <default value>,
161  face_rotation = <default value>,
162  matrix = <default value>);
163 @endcode
164 Here, we need to specify the orientation of the two faces using
165 @p face_orientation, @p face_flip and @p face_orientation. For a closer description
166 have a look at the documentation of DoFTools::make_periodicity_constraints.
167 The remaining parameters are the same as for the high level interface apart
168 from the self-explaining @p component_mask and @p affine_constraints.
169 
170 
171 <a name="problem"></a>
172 <a name="Apracticalexample"></a><h1>A practical example</h1>
173 
174 
175 In the following, we show how to use the above functions in a more involved
176 example. The task is to enforce rotated periodicity constraints for the
177 velocity component of a Stokes flow.
178 
179 On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
180 going to solve the Stokes problem
181 @f{eqnarray*}
182  -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
183  -\textrm{div}\; \textbf{u}&=&0,\\
184  \textbf{u}|_{\Gamma_1}&=&{\bf 0},
185 @f}
186 where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
187 For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
188 @f{align*}
189  u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
190  u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
191 @f}
192 
193 The mesh will be generated by GridGenerator::quarter_hyper_shell(),
194 which also documents how it assigns boundary indicators to its various
195 boundaries if its `colorize` argument is set to `true`.
196  *
197  *
198  * <a name="CommProg"></a>
199  * <h1> The commented program</h1>
200  *
201  * This example program is a slight modification of @ref step_22 "step-22" running in parallel
202  * using Trilinos to demonstrate the usage of periodic boundary conditions in
203  * deal.II. We thus omit to discuss the majority of the source code and only
204  * comment on the parts that deal with periodicity constraints. For the rest
205  * have a look at @ref step_22 "step-22" and the full source code at the bottom.
206  *
207 
208  *
209  * In order to implement periodic boundary conditions only two functions
210  * have to be modified:
211  * - <code>StokesProblem<dim>::setup_dofs()</code>:
212  * To populate an AffineConstraints object with periodicity constraints
213  * - <code>StokesProblem<dim>::create_mesh()</code>:
214  * To supply a distributed triangulation with periodicity information.
215  *
216 
217  *
218  * The rest of the program is identical to @ref step_22 "step-22", so let us skip this part
219  * and only show these two functions in the following. (The full program can be
220  * found in the "Plain program" section below, though.)
221  *
222 
223  *
224  *
225 
226  *
227  *
228 
229  *
230  *
231  * <a name="Settingupperiodicityconstraintsondistributedtriangulations"></a>
232  * <h3>Setting up periodicity constraints on distributed triangulations</h3>
233  *
234  * @code
235  * template <int dim>
236  * void StokesProblem<dim>::create_mesh()
237  * {
238  * Point<dim> center;
239  * const double inner_radius = .5;
240  * const double outer_radius = 1.;
241  *
243  * triangulation, center, inner_radius, outer_radius, 0, true);
244  *
245  * @endcode
246  *
247  * Before we can prescribe periodicity constraints, we need to ensure that
248  * cells on opposite sides of the domain but connected by periodic faces are
249  * part of the ghost layer if one of them is stored on the local processor.
250  * At this point we need to think about how we want to prescribe
251  * periodicity. The vertices @f$\text{vertices}_2@f$ of a face on the left
252  * boundary should be matched to the vertices @f$\text{vertices}_1@f$ of a face
253  * on the lower boundary given by @f$\text{vertices}_2=R\cdot
254  * \text{vertices}_1+b@f$ where the rotation matrix @f$R@f$ and the offset @f$b@f$ are
255  * given by
256  * @f{align*}
257  * R=\begin{pmatrix}
258  * 0&1\\-1&0
259  * \end{pmatrix},
260  * \quad
261  * b=\begin{pmatrix}0&0\end{pmatrix}.
262  * @f}
263  * The data structure we are saving the resulting information into is here
264  * based on the Triangulation.
265  *
266  * @code
267  * std::vector<GridTools::PeriodicFacePair<
269  * periodicity_vector;
270  *
271  * FullMatrix<double> rotation_matrix(dim);
272  * rotation_matrix[0][1] = 1.;
273  * rotation_matrix[1][0] = -1.;
274  *
276  * 2,
277  * 3,
278  * 1,
279  * periodicity_vector,
280  * Tensor<1, dim>(),
281  * rotation_matrix);
282  *
283  * @endcode
284  *
285  * Now telling the triangulation about the desired periodicity is
286  * particularly easy by just calling
288  *
289  * @code
290  * triangulation.add_periodicity(periodicity_vector);
291  *
292  * triangulation.refine_global(4 - dim);
293  * }
294  *
295  *
296  * template <int dim>
297  * void StokesProblem<dim>::setup_dofs()
298  * {
299  * dof_handler.distribute_dofs(fe);
300  *
301  * std::vector<unsigned int> block_component(dim + 1, 0);
302  * block_component[dim] = 1;
303  * DoFRenumbering::component_wise(dof_handler, block_component);
304  *
305  * const std::vector<types::global_dof_index> dofs_per_block =
306  * DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
307  * const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
308  *
309  * {
310  * owned_partitioning.clear();
311  * IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
312  * owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
313  * owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u + n_p));
314  *
315  * relevant_partitioning.clear();
316  * IndexSet locally_relevant_dofs;
318  * locally_relevant_dofs);
319  * relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
320  * relevant_partitioning.push_back(
321  * locally_relevant_dofs.get_view(n_u, n_u + n_p));
322  *
323  * constraints.clear();
324  * constraints.reinit(locally_relevant_dofs);
325  *
326  * FEValuesExtractors::Vector velocities(0);
327  *
328  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
330  * dof_handler,
331  * 0,
332  * BoundaryValues<dim>(),
333  * constraints,
334  * fe.component_mask(velocities));
336  * dof_handler,
337  * 1,
338  * BoundaryValues<dim>(),
339  * constraints,
340  * fe.component_mask(velocities));
341  *
342  * @endcode
343  *
344  * After we provided the mesh with the necessary information for the
345  * periodicity constraints, we are now able to actual create them. For
346  * describing the matching we are using the same approach as before, i.e.,
347  * the @f$\text{vertices}_2@f$ of a face on the left boundary should be
348  * matched to the vertices
349  * @f$\text{vertices}_1@f$ of a face on the lower boundary given by
350  * @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ where the rotation
351  * matrix @f$R@f$ and the offset @f$b@f$ are given by
352  * @f{align*}
353  * R=\begin{pmatrix}
354  * 0&1\\-1&0
355  * \end{pmatrix},
356  * \quad
357  * b=\begin{pmatrix}0&0\end{pmatrix}.
358  * @f}
359  * These two objects not only describe how faces should be matched but
360  * also in which sense the solution should be transformed from
361  * @f$\text{face}_2@f$ to
362  * @f$\text{face}_1@f$.
363  *
364  * @code
365  * FullMatrix<double> rotation_matrix(dim);
366  * rotation_matrix[0][1] = 1.;
367  * rotation_matrix[1][0] = -1.;
368  *
369  * Tensor<1, dim> offset;
370  *
371  * @endcode
372  *
373  * For setting up the constraints, we first store the periodicity
374  * information in an auxiliary object of type
375  * <code>std::vector@<GridTools::PeriodicFacePair<typename
376  * DoFHandler@<dim@>::%cell_iterator@> </code>. The periodic boundaries
377  * have the boundary indicators 2 (x=0) and 3 (y=0). All the other
378  * parameters we have set up before. In this case the direction does not
379  * matter. Due to @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ this is
380  * exactly what we want.
381  *
382  * @code
383  * std::vector<
385  * periodicity_vector;
386  *
387  * const unsigned int direction = 1;
388  *
389  * GridTools::collect_periodic_faces(dof_handler,
390  * 2,
391  * 3,
392  * direction,
393  * periodicity_vector,
394  * offset,
395  * rotation_matrix);
396  *
397  * @endcode
398  *
399  * Next, we need to provide information on which vector valued components
400  * of the solution should be rotated. Since we choose here to just
401  * constraint the velocity and this starts at the first component of the
402  * solution vector, we simply insert a 0:
403  *
404  * @code
405  * std::vector<unsigned int> first_vector_components;
406  * first_vector_components.push_back(0);
407  *
408  * @endcode
409  *
410  * After setting up all the information in periodicity_vector all we have
411  * to do is to tell make_periodicity_constraints to create the desired
412  * constraints.
413  *
414  * @code
415  * DoFTools::make_periodicity_constraints<DoFHandler<dim>>(
416  * periodicity_vector,
417  * constraints,
418  * fe.component_mask(velocities),
419  * first_vector_components);
420  *
422  * dof_handler,
423  * 0,
424  * BoundaryValues<dim>(),
425  * constraints,
426  * fe.component_mask(velocities));
428  * dof_handler,
429  * 1,
430  * BoundaryValues<dim>(),
431  * constraints,
432  * fe.component_mask(velocities));
433  * }
434  *
435  * constraints.close();
436  *
437  * {
438  * TrilinosWrappers::BlockSparsityPattern bsp(owned_partitioning,
439  * owned_partitioning,
440  * relevant_partitioning,
441  * mpi_communicator);
442  *
443  * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
444  * for (unsigned int c = 0; c < dim + 1; ++c)
445  * for (unsigned int d = 0; d < dim + 1; ++d)
446  * if (!((c == dim) && (d == dim)))
447  * coupling[c][d] = DoFTools::always;
448  * else
449  * coupling[c][d] = DoFTools::none;
450  *
451  * DoFTools::make_sparsity_pattern(dof_handler,
452  * coupling,
453  * bsp,
454  * constraints,
455  * false,
457  * mpi_communicator));
458  *
459  * bsp.compress();
460  *
461  * system_matrix.reinit(bsp);
462  * }
463  *
464  * {
465  * TrilinosWrappers::BlockSparsityPattern preconditioner_bsp(
466  * owned_partitioning,
467  * owned_partitioning,
468  * relevant_partitioning,
469  * mpi_communicator);
470  *
471  * Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
472  * for (unsigned int c = 0; c < dim + 1; ++c)
473  * for (unsigned int d = 0; d < dim + 1; ++d)
474  * if ((c == dim) && (d == dim))
475  * preconditioner_coupling[c][d] = DoFTools::always;
476  * else
477  * preconditioner_coupling[c][d] = DoFTools::none;
478  *
479  * DoFTools::make_sparsity_pattern(dof_handler,
480  * preconditioner_coupling,
481  * preconditioner_bsp,
482  * constraints,
483  * false,
485  * mpi_communicator));
486  *
487  * preconditioner_bsp.compress();
488  *
489  * preconditioner_matrix.reinit(preconditioner_bsp);
490  * }
491  *
492  * system_rhs.reinit(owned_partitioning, mpi_communicator);
493  * solution.reinit(owned_partitioning,
494  * relevant_partitioning,
495  * mpi_communicator);
496  * }
497  *
498  * @endcode
499  *
500  * The rest of the program is then again identical to @ref step_22 "step-22". We will omit
501  * it here now, but as before, you can find these parts in the "Plain program"
502  * section below.
503  *
504 
505  *
506 <a name="Results"></a><h1>Results</h1>
507 
508 
509 The created output is not very surprising. We simply see that the solution is
510 periodic with respect to the left and lower boundary:
511 
512 <img src="https://www.dealii.org/images/steps/developer/step-45.periodic.png" alt="">
513 
514 Without the periodicity constraints we would have ended up with the following solution:
515 
516 <img src="https://www.dealii.org/images/steps/developer/step-45.non_periodic.png" alt="">
517  *
518  *
519 <a name="PlainProg"></a>
520 <h1> The plain program</h1>
521 @include "step-45.cc"
522 */
GridTools::collect_periodic_faces
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
Definition: grid_tools_dof_handlers.cc:2196
IndexSet
Definition: index_set.h:74
parallel::distributed::Triangulation::add_periodicity
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator >> &) override
Definition: tria.cc:4575
Triangulation
Definition: tria.h:1109
DoFTools::make_periodicity_constraints
void make_periodicity_constraints(const FaceIterator &face_1, const typename identity< FaceIterator >::type &face_2, AffineConstraints< number > &constraints, const ComponentMask &component_mask=ComponentMask(), const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >(), const number periodicity_factor=1.)
Definition: dof_tools_constraints.cc:2204
IndexSet::get_view
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:211
GridGenerator::quarter_hyper_shell
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
DoFTools::always
@ always
Definition: dof_tools.h:236
GridGenerator::parallelogram
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler
Definition: dof_handler.h:205
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Table
Definition: table.h:699
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
level
unsigned int level
Definition: grid_out.cc:4355
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
LAPACKSupport::T
static const char T
Definition: lapack_support.h:163
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
VectorizedArray::exp
VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5368
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
DoFTools::none
@ none
Definition: dof_tools.h:232
Tensor
Definition: tensor.h:450
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
TrilinosWrappers::BlockSparsityPattern
Definition: block_sparsity_pattern.h:638
IndexSet::clear
void clear()
Definition: index_set.h:1611
DoFTools::count_dofs_per_fe_block
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandlerType &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:2012
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
AffineConstraints
Definition: affine_constraints.h:180
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
GridTools::PeriodicFacePair
Definition: grid_tools.h:2463
Point< dim >
colorize
bool colorize
Definition: grid_out.cc:4354
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
DoFTools
Definition: dof_tools.h:214
first
Point< 2 > first
Definition: grid_out.cc:4352
DoFRenumbering::component_wise
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition: dof_renumbering.cc:633
TriaIterator
Definition: tria_iterator.h:578
center
Point< 3 > center
Definition: data_out_base.cc:169
parallel
Definition: distributed.h:416