113 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$
for all
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:
136 Here, again, the assignment of boundary indicators 0 and 1 stems from
139 The resulting @p matched_pairs can be used in
141 object with periodicity constraints:
143 DoFTools::make_periodicity_constraints(matched_pairs, constraints);
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).
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
156 make_periodicity_constraints(face_1,
159 component_mask = <default value>;
160 face_orientation = <default value>,
161 face_flip = <default value>,
162 face_rotation = <default value>,
163 matrix = <default value>);
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.
172 <a name="problem"></a>
173 <a name="Apracticalexample"></a><h1>A practical example</h1>
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.
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
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},
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.
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].
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>
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.
208 * In order to implement periodic boundary conditions only two functions
209 * have to be modified:
210 * - <code>StokesProblem<dim>::setup_dofs()</code>:
212 * - <code>StokesProblem<dim>::run()</code>:
213 * To supply a distributed triangulation with periodicity information.
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.)
230 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
231 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
235 *
void StokesProblem<dim>::create_mesh()
238 *
const double inner_radius = .5;
239 *
const double outer_radius = 1.;
242 * triangulation, center, inner_radius, outer_radius, 0,
true);
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
260 *
b=
\begin{pmatrix}0&0\end{pmatrix}.
262 * The data structure we are saving the resulting information into is here
268 * periodicity_vector;
271 * rotation_matrix[0][1] = 1.;
272 * rotation_matrix[1][0] = -1.;
278 * periodicity_vector,
284 * Now telling the triangulation about the desired periodicity is
285 * particularly easy by just calling
289 * triangulation.add_periodicity(periodicity_vector);
291 * triangulation.refine_global(4 - dim);
298 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
299 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
303 *
void StokesProblem<dim>::setup_dofs()
305 * dof_handler.distribute_dofs(fe);
307 * std::vector<unsigned int> block_component(dim + 1, 0);
308 * block_component[dim] = 1;
311 * std::vector<types::global_dof_index> dofs_per_block(2);
315 *
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
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));
323 * relevant_partitioning.
clear();
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));
331 * constraints.clear();
332 * constraints.reinit(locally_relevant_dofs);
340 * BoundaryValues<dim>(),
342 * fe.component_mask(velocities));
346 * BoundaryValues<dim>(),
348 * fe.component_mask(velocities));
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
365 *
b=
\begin{pmatrix}0&0\end{pmatrix}.
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$.
374 * rotation_matrix[0][1] = 1.;
375 * rotation_matrix[1][0] = -1.;
381 * For setting up the constraints, we first store the periodicity
382 * information in an auxiliary
object of type
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.
393 * periodicity_vector;
395 *
const unsigned int direction = 1;
401 * periodicity_vector,
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:
413 * std::vector<unsigned int> first_vector_components;
414 * first_vector_components.push_back(0);
418 * After setting up all the information in periodicity_vector all we have
423 * DoFTools::make_periodicity_constraints<DoFHandler<dim>>(
424 * periodicity_vector,
426 * fe.component_mask(velocities),
427 * first_vector_components);
432 * BoundaryValues<dim>(),
434 * fe.component_mask(velocities));
438 * BoundaryValues<dim>(),
440 * fe.component_mask(velocities));
443 * constraints.close();
447 * owned_partitioning,
448 * relevant_partitioning,
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)))
465 * mpi_communicator));
469 * system_matrix.reinit(bsp);
474 * owned_partitioning,
475 * owned_partitioning,
476 * relevant_partitioning,
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))
488 * preconditioner_coupling,
489 * preconditioner_bsp,
493 * mpi_communicator));
495 * preconditioner_bsp.compress();
497 * preconditioner_matrix.reinit(preconditioner_bsp);
500 * system_rhs.reinit(owned_partitioning, mpi_communicator);
501 * solution.reinit(owned_partitioning,
502 * relevant_partitioning,
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" 514 <a name=
"Results"></a><h1>Results</h1>
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:
520 <img src=
"https://www.dealii.org/images/steps/developer/step-45.periodic.png" alt=
"">
522 Without the periodicity constraints we would have ended up with the following solution:
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" Contents is actually a matrix.
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, 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)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
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
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator >> &) override