112 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$
for all
115 If 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:
135 Here, again, the assignment of boundary indicators 0 and 1 stems from
138 The resulting @p matched_pairs can be used in
140 object with periodicity constraints:
145 Apart from this high
level interface there are also variants of
149 There is also a low
level interface to
151 low
level variant allows to directly specify two faces that shall be
158 component_mask = <default
value>;
159 face_orientation = <default
value>,
160 face_flip = <default
value>,
161 face_rotation = <default
value>,
164 Here, 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
167 The remaining parameters are the same as for the high
level interface apart
168 from the self-explaining @p component_mask and @p affine_constraints.
171 <a name="problem"></a>
172 <a name="Apracticalexample"></a><h1>
A practical example</h1>
175 In the following, we show how to use the above
functions in a more involved
176 example. The task is to enforce rotated periodicity constraints for the
177 velocity component of a Stokes flow.
179 On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
180 going to solve the Stokes problem
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},
186 where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
187 For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
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].
194 which also documents how it assigns boundary indicators to its various
195 boundaries
if its `
colorize` argument is
set to `
true`.
198 * <a name=
"CommProg"></a>
199 * <h1> The commented program</h1>
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.
209 * In order to implement periodic boundary conditions only two
functions
210 * have to be modified:
211 * - <code>StokesProblem<dim>::setup_dofs()</code>:
213 * - <code>StokesProblem<dim>::create_mesh()</code>:
214 * To supply a distributed
triangulation with periodicity information.
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.)
231 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
232 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
236 *
void StokesProblem<dim>::create_mesh()
239 *
const double inner_radius = .5;
240 *
const double outer_radius = 1.;
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
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
263 * The data structure we are saving the resulting information into is here
269 * periodicity_vector;
272 * rotation_matrix[0][1] = 1.;
273 * rotation_matrix[1][0] = -1.;
279 * periodicity_vector,
285 * Now telling the
triangulation about the desired periodicity is
286 * particularly easy by just calling
297 *
void StokesProblem<dim>::setup_dofs()
299 * dof_handler.distribute_dofs(fe);
301 * std::vector<unsigned int> block_component(dim + 1, 0);
302 * block_component[dim] = 1;
305 *
const std::vector<types::global_dof_index> dofs_per_block =
307 *
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
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));
315 * relevant_partitioning.
clear();
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));
323 * constraints.clear();
324 * constraints.reinit(locally_relevant_dofs);
332 * BoundaryValues<dim>(),
334 * fe.component_mask(velocities));
338 * BoundaryValues<dim>(),
340 * fe.component_mask(velocities));
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
349 * @f$\text{
vertices}_1@f$ of a face on the lower boundary given by
351 *
matrix @f$R@f$ and the offset @f$b@f$ are given by
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$.
366 * rotation_matrix[0][1] = 1.;
367 * rotation_matrix[1][0] = -1.;
373 * For setting up the constraints, we
first store the periodicity
374 * information in an auxiliary
object of type
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
380 * exactly what we want.
385 * periodicity_vector;
387 *
const unsigned int direction = 1;
393 * periodicity_vector,
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:
405 * std::vector<unsigned int> first_vector_components;
406 * first_vector_components.push_back(0);
410 * After setting up all the information in periodicity_vector all we have
415 * DoFTools::make_periodicity_constraints<DoFHandler<dim>>(
416 * periodicity_vector,
418 * fe.component_mask(velocities),
419 * first_vector_components);
424 * BoundaryValues<dim>(),
426 * fe.component_mask(velocities));
430 * BoundaryValues<dim>(),
432 * fe.component_mask(velocities));
435 * constraints.close();
439 * owned_partitioning,
440 * relevant_partitioning,
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)))
457 * mpi_communicator));
461 * system_matrix.reinit(bsp);
466 * owned_partitioning,
467 * owned_partitioning,
468 * relevant_partitioning,
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))
480 * preconditioner_coupling,
481 * preconditioner_bsp,
485 * mpi_communicator));
487 * preconditioner_bsp.compress();
489 * preconditioner_matrix.reinit(preconditioner_bsp);
492 * system_rhs.reinit(owned_partitioning, mpi_communicator);
493 * solution.reinit(owned_partitioning,
494 * relevant_partitioning,
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"
506 <a name=
"Results"></a><h1>Results</h1>
509 The created output is not very surprising. We simply see that the solution is
510 periodic with respect to the left and lower boundary:
512 <img src=
"https://www.dealii.org/images/steps/developer/step-45.periodic.png" alt=
"">
514 Without the periodicity constraints we would have ended up with the following solution:
516 <img src=
"https://www.dealii.org/images/steps/developer/step-45.non_periodic.png" alt=
"">
519 <a name=
"PlainProg"></a>
520 <h1> The plain program</h1>
521 @include
"step-45.cc"