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