Reference documentation for deal.II version 9.0.0
step-45.h
1  0,
96  /*b_id2*/ 1,
97  /*direction*/ 0,
98  matched_pairs);
99 @endcode
100 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$ for all
101 @f$y\in[0,1]@f$.
102 
103 If we instead consider the parallelogram given by the convex hull of
104 @f$(0,0)@f$, @f$(1,1)@f$, @f$(1,2)@f$, @f$(0,1)@f$ we can achieve the constraints
105 @f$u(0,y)=u(1,y+1)@f$ by specifying an @p offset:
106 @code
108  /*b_id1*/ 0,
109  /*b_id2*/ 1,
110  /*direction*/ 0,
111  matched_pairs,
112  Tensor<1, 2>(0.,1.));
113 @endcode
114 or
115 @code
117  /*b_id1*/ 0,
118  /*b_id2*/ 1,
119  /*arbitrary direction*/ 0,
120  matched_pairs,
121  Tensor<1, 2>(1.,1.));
122 @endcode
123 
124 The resulting @p matched_pairs can be used in
126 with periodicity constraints:
127 @code
128 DoFTools::make_periodicity_constraints(matched_pairs, constraints);
129 @endcode
130 
131 Apart from this high level interface there are also variants of
132 DoFTools::make_periodicity_constraints available that combine those two
133 steps (see the variants of DofTools::make_periodicity_constraints).
134 
135 There is also a low level interface to
136 DoFTools::make_periodicity_constraints if more flexibility is needed. The
137 low level variant allows to directly specify two faces that shall be
138 constrained:
139 @code
140 using namespace DoFTools;
142  face_2,
143  constraint_matrix,
144  component_mask = <default value>;
145  face_orientation = <default value>,
146  face_flip = <default value>,
147  face_rotation = <default value>,
148  matrix = <default value>);
149 @endcode
150 Here, we need to specify the orientation of the two faces using
151 @p face_orientation, @p face_flip and @p face_orientation. For a closer description
152 have a look at the documentation of DoFTools::make_periodicity_constraints.
153 The remaining parameters are the same as for the high level interface apart
154 from the self-explaining @p component_mask and @p constraint_matrix.
155 
156 <a name="problem"></a>
157 <a name="Apracticalexample"></a><h1>A practical example</h1>
158 
159 
160 In the following, we show how to use the above functions in a more involved
161 example. The task is to enforce rotated periodicity constraints for the
162 velocity component of a Stokes flow.
163 
164 On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
165 going to solve the Stokes problem
166 @f{eqnarray*}
167  -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100*\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
168  -\textrm{div}\; \textbf{u}&=&0,\\
169  \textbf{u}|_{\Gamma_1}&=&{\bf 0},
170 @f}
171 where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1:=\{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
172 For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
173 @f{align*}
174  u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
175  u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
176 @f}
177  * <a name="CommProg"></a>
178  * <h1> The commented program</h1>
179  *
180  * This example program is a slight modification of @ref step_22 "step-22" running in parallel
181  * using Trilinos to demonstrate the usage of periodic boundary conditions in deal.II.
182  * We thus omit to discuss the majority of the source code and only comment on the
183  * parts that deal with periodicity constraints. For the rest have a look at @ref step_22 "step-22"
184  * and the full source code at the bottom.
185  *
186 
187  *
188  * In order to implement periodic boundary conditions only two functions
189  * have to be modified:
190  * - <code>StokesProblem<dim>::setup_dofs()</code>: To populate a ConstraintMatrix
191  * object with periodicity constraints
192  * - <code>StokesProblem<dim>::run()</code>: To supply a distributed triangulation with
193  * periodicity information.
194  *
195 
196  *
197  *
198 
199  *
200  *
201  *
202  *
203  * <a name="Settingupperiodicityconstraintsondistributedtriangulations"></a>
204  * <h3>Setting up periodicity constraints on distributed triangulations</h3>
205  *
206  * @code
207  * template <int dim>
208  * void StokesProblem<dim>::create_mesh()
209  * {
210  * Point<dim> center;
211  * const double inner_radius = .5;
212  * const double outer_radius = 1.;
213  *
214  * GridGenerator::quarter_hyper_shell (triangulation,
215  * center,
216  * inner_radius,
217  * outer_radius,
218  * 0,
219  * true);
220  *
221  * @endcode
222  *
223  * Before we can prescribe periodicity constraints, we need to ensure that cells
224  * on opposite sides of the domain but connected by periodic faces are part of
225  * the ghost layer if one of them is stored on the local processor.
226  * At this point we need to think about how we want to prescribe periodicity.
227  * The vertices @f$\text{vertices}_2@f$ of a face on the left boundary should be
228  * matched to the vertices @f$\text{vertices}_1@f$ of a face on the lower boundary
229  * given by @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ where the rotation
230  * matrix @f$R@f$ and the offset @f$b@f$ are given by
231  * @f{align*}
232  * R=\begin{pmatrix}
233  * 0&1\\-1&0
234  * \end{pmatrix},
235  * \quad
236  * b=\begin{pmatrix}0&0\end{pmatrix}.
237  * @f}
238  * The data structure we are saving the reuslitng information into is here based
239  * on the Triangulation.
240  *
241  * @code
242  * std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
243  * periodicity_vector;
244  *
245  * FullMatrix<double> rotation_matrix(dim);
246  * rotation_matrix[0][1]=1.;
247  * rotation_matrix[1][0]=-1.;
248  *
249  * GridTools::collect_periodic_faces(triangulation, 2, 3, 1,
250  * periodicity_vector, Tensor<1, dim>(),
251  * rotation_matrix);
252  *
253  * @endcode
254  *
255  * Now telling the triangulation about the desired periodicity is
256  * particularly easy by just calling
258  *
259  * @code
260  * triangulation.add_periodicity(periodicity_vector);
261  *
262  * triangulation.refine_global (4-dim);
263  * }
264  *
265  *
266  * @endcode
267  *
268  *
269  * <a name="Settingupperiodicityconstraintsondistributedtriangulations"></a>
270  * <h3>Setting up periodicity constraints on distributed triangulations</h3>
271  *
272  * @code
273  * template <int dim>
274  * void StokesProblem<dim>::setup_dofs ()
275  * {
276  * dof_handler.distribute_dofs (fe);
277  *
278  * std::vector<unsigned int> block_component (dim+1,0);
279  * block_component[dim] = 1;
280  * DoFRenumbering::component_wise (dof_handler, block_component);
281  *
282  * std::vector<types::global_dof_index> dofs_per_block (2);
283  * DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
284  * const unsigned int n_u = dofs_per_block[0],
285  * n_p = dofs_per_block[1];
286  *
287  * {
288  * owned_partitioning.clear();
289  * IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
290  * owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
291  * owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u+n_p));
292  *
293  * relevant_partitioning.clear();
294  * IndexSet locally_relevant_dofs;
295  * DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
296  * relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
297  * relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u, n_u+n_p));
298  *
299  * constraints.clear ();
300  * constraints.reinit(locally_relevant_dofs);
301  *
302  * FEValuesExtractors::Vector velocities(0);
303  *
305  * constraints);
307  * dof_handler,
308  * 0,
309  * BoundaryValues<dim>(),
310  * constraints,
311  * fe.component_mask(velocities));
313  * dof_handler,
314  * 1,
315  * BoundaryValues<dim>(),
316  * constraints,
317  * fe.component_mask(velocities));
318  *
319  * @endcode
320  *
321  * After we provided the mesh with the necessary information for the periodicity
322  * constraints, we are now able to actual create them. For describing the
323  * matching we are using the same approach as before, i.e., the @f$\text{vertices}_2@f$
324  * of a face on the left boundary should be matched to the vertices
325  * @f$\text{vertices}_1@f$ of a face on the lower boundary given by
326  * @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ where the rotation matrix @f$R@f$
327  * and the offset @f$b@f$ are given by
328  * @f{align*}
329  * R=\begin{pmatrix}
330  * 0&1\\-1&0
331  * \end{pmatrix},
332  * \quad
333  * b=\begin{pmatrix}0&0\end{pmatrix}.
334  * @f}
335  * These two objects not only describe how faces should be matched but also
336  * in which sense the solution should be transformed from @f$\text{face}_2@f$ to
337  * @f$\text{face}_1@f$.
338  *
339  * @code
340  * FullMatrix<double> rotation_matrix(dim);
341  * rotation_matrix[0][1]=1.;
342  * rotation_matrix[1][0]=-1.;
343  *
344  * Tensor<1,dim> offset;
345  *
346  * @endcode
347  *
348  * For setting up the constraints, we first store the periodicity
349  * information in an auxiliary object of type
350  * <code>std::vector@<GridTools::PeriodicFacePair<typename
351  * DoFHandler@<dim@>::cell_iterator@> </code>. The periodic boundaries have the
352  * boundary indicators 2 (x=0) and 3 (y=0). All the other parameters we
353  * have set up before. In this case the direction does not matter. Due to
354  * @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ this is exactly what we want.
355  *
356  * @code
357  * std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
358  * periodicity_vector;
359  *
360  * const unsigned int direction = 1;
361  *
362  * GridTools::collect_periodic_faces(dof_handler, 2, 3, direction,
363  * periodicity_vector, offset,
364  * rotation_matrix);
365  *
366  * @endcode
367  *
368  * Next we need to provide information on which vector valued components of
369  * the solution should be rotated. Since we choose here to just constraint the
370  * velocity and this starts at the first component of the solution vector we
371  * simply insert a 0:
372  *
373  * @code
374  * std::vector<unsigned int> first_vector_components;
375  * first_vector_components.push_back(0);
376  *
377  * @endcode
378  *
379  * After setting up all the information in periodicity_vector all we have
380  * to do is to tell make_periodicity_constraints to create the desired
381  * constraints.
382  *
383  * @code
384  * DoFTools::make_periodicity_constraints<DoFHandler<dim> >
385  * (periodicity_vector, constraints, fe.component_mask(velocities),
386  * first_vector_components);
387  *
389  * dof_handler,
390  * 0,
391  * BoundaryValues<dim>(),
392  * constraints,
393  * fe.component_mask(velocities));
395  * dof_handler,
396  * 1,
397  * BoundaryValues<dim>(),
398  * constraints,
399  * fe.component_mask(velocities));
400  *
401  * }
402  *
403  * constraints.close ();
404  *
405  * {
407  * (owned_partitioning, owned_partitioning,
408  * relevant_partitioning, mpi_communicator);
409  *
411  * (dof_handler, bsp, constraints, false,
412  * Utilities::MPI::this_mpi_process(mpi_communicator));
413  *
414  * bsp.compress();
415  *
416  * system_matrix.reinit (bsp);
417  * }
418  *
419  * system_rhs.reinit (owned_partitioning,
420  * mpi_communicator);
421  * solution.reinit (owned_partitioning, relevant_partitioning,
422  * mpi_communicator);
423  * }
424  *
425  *
426  * @endcode
427  *
428 <a name="Results"></a><h1>Results</h1>
429 
430 
431 The created output is not very surprising. We simply see that the solution is
432 periodic with respect to the left and lower boundary:
433 
434 <img src="/documentation/images/steps/developer/step-45.periodic.png" alt="">
435 
436 Without the periodicity constraints we would have ended up with the following solution:
437 
438 <img src="/documentation/images/steps/developer/step-45.non_periodic.png" alt="">
439 <a name="PlainProg"></a>
440 <h1> The plain program</h1>
441 @include "step-45.cc"
442  */
Contents is actually a matrix.
Point< spacedim > point(const gp_Pnt &p, const double &tolerance=1e-10)
Definition: utilities.cc:183
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)
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator > > &)
Definition: tria.cc:3427
void clear()
Definition: index_set.h:1549
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, ConstraintMatrix &constraints)
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 >())
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1087
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:254
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:1836
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const ConstraintMatrix &constraints=ConstraintMatrix(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: mpi.h:53
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:75
void make_periodicity_constraints(const FaceIterator &face_1, const typename identity< FaceIterator >::type &face_2, ::ConstraintMatrix &constraint_matrix, 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 >())