Processing math: 100%
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
step-45.h
Go to the documentation of this file.
1 0,
108 /*b_id2*/ 1,
109 /*direction*/ 0,
110 matched_pairs);
111@endcode
112would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$ for all
113@f$y\in[0,1]@f$.
114
115If 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
126or
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
135Here, again, the assignment of boundary indicators 0 and 1 stems from
136what GridGenerator::parallelogram() documents.
137
138The resulting @p matched_pairs can be used in
140object with periodicity constraints:
141@code
142DoFTools::make_periodicity_constraints(matched_pairs, constraints);
143@endcode
144
145Apart from this high level interface there are also variants of
146DoFTools::make_periodicity_constraints available that combine those two
147steps (see the variants of DofTools::make_periodicity_constraints).
148
149There is also a low level interface to
150DoFTools::make_periodicity_constraints if more flexibility is needed. The
151low level variant allows to directly specify two faces that shall be
152constrained:
153@code
154using 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
164Here, 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
166have a look at the documentation of DoFTools::make_periodicity_constraints.
167The remaining parameters are the same as for the high level interface apart
168from 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
175In the following, we show how to use the above functions in a more involved
176example. The task is to enforce rotated periodicity constraints for the
177velocity component of a Stokes flow.
178
179On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
180going 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}
186where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
187For 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
193The mesh will be generated by GridGenerator::quarter_hyper_shell(),
194which also documents how it assigns boundary indicators to its various
195boundaries 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 * {
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 *
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<dim, dim>(periodicity_vector,
416 * constraints,
417 * fe.component_mask(
418 * 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 *
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 *
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
509The created output is not very surprising. We simply see that the solution is
510periodic 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
514Without 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*/
void clear()
Definition: index_set.h:1610
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:211
Definition: point.h:111
Definition: tensor.h:472
VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator > > &) override
Definition: tria.cc:3048
Point< 3 > center
Point< 3 > vertices[4]
bool colorize
Definition: grid_out.cc:4589
Point< 2 > first
Definition: grid_out.cc:4587
unsigned int level
Definition: grid_out.cc:4590
__global__ void set(Number *val, const Number s, const size_type N)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &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)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:2047
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.)
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1210
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
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 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 >())
@ matrix
Contents is actually a matrix.
static const char A
static const char T
static const types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< 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())
int(&) functions(const void *v1, const void *v2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation