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