98 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$
for all
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:
122 The resulting @p matched_pairs can be used in
124 with periodicity constraints:
129 Apart from
this high level
interface there are also variants of
131 steps (see the variants of DofTools::make_periodicity_constraints).
133 There is also a low level interface to
135 low level variant allows to directly specify two faces that shall be
142 component_mask = <
default value>;
143 face_orientation = <
default value>,
144 face_flip = <
default value>,
145 face_rotation = <
default value>,
146 matrix = <
default value>);
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
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.
154 <a name=
"problem"></a>
155 <a name=
"Apracticalexample"></a><h1>A practical example</h1>
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.
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
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},
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.
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].
175 * <a name=
"CommProg"></a>
176 * <h1> The commented program</h1>
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.
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.
201 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
202 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
206 *
void StokesProblem<dim>::create_mesh()
209 *
const double inner_radius = .5;
210 *
const double outer_radius = 1.;
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
234 *
b=
\begin{pmatrix}0&0\end{pmatrix}
.
236 * The data structure we are saving the reuslitng information into is here based
240 * std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
241 * periodicity_vector;
244 * rotation_matrix[0][1]=1.;
245 * rotation_matrix[1][0]=-1.;
253 * Now telling the triangulation about the desired periodicity is
254 * particularly easy by just calling
258 * triangulation.add_periodicity(periodicity_vector);
260 * triangulation.set_boundary(0, boundary);
261 * triangulation.set_boundary(1, boundary);
263 * triangulation.refine_global (4-dim);
270 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
271 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
275 *
void StokesProblem<dim>::setup_dofs ()
277 * dof_handler.distribute_dofs (fe);
279 * std::vector<unsigned int> block_component (dim+1,0);
280 * block_component[dim] = 1;
283 * std::vector<types::global_dof_index> dofs_per_block (2);
285 *
const unsigned int n_u = dofs_per_block[0],
286 * n_p = dofs_per_block[1];
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));
294 * relevant_partitioning.
clear();
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));
300 * constraints.clear ();
301 * constraints.reinit(locally_relevant_dofs);
310 * BoundaryValues<dim>(),
312 * fe.component_mask(velocities));
316 * BoundaryValues<dim>(),
318 * fe.component_mask(velocities));
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
334 *
b=
\begin{pmatrix}0&0\end{pmatrix}
.
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$.
342 * rotation_matrix[0][1]=1.;
343 * rotation_matrix[1][0]=-1.;
349 * For setting up the constraints, we first store the periodicity
350 * information in an auxiliary
object of type
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.
358 * std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
359 * periodicity_vector;
361 *
const unsigned int direction = 1;
364 * periodicity_vector, offset,
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
375 * std::vector<unsigned int> first_vector_components;
376 * first_vector_components.push_back(0);
380 * After setting up all the information in periodicity_vector all we have
385 * DoFTools::make_periodicity_constraints<DoFHandler<dim> >
386 * (periodicity_vector, constraints, fe.component_mask(velocities),
387 * first_vector_components);
392 * BoundaryValues<dim>(),
394 * fe.component_mask(velocities));
398 * BoundaryValues<dim>(),
400 * fe.component_mask(velocities));
404 * constraints.close ();
408 * (owned_partitioning, owned_partitioning,
409 * relevant_partitioning, mpi_communicator);
412 * (dof_handler, bsp, constraints,
false,
417 * system_matrix.reinit (bsp);
420 * system_rhs.reinit (owned_partitioning,
422 * solution.reinit (owned_partitioning, relevant_partitioning,
429 <a name=
"Results"></a><h1>Results</h1>
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:
435 <img src=
"/documentation/images/steps/developer/step-45.periodic.png" alt=
"">
437 Without the periodicity constraints we would have ended up with the following solution:
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"
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(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 > > &)
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, ConstraintMatrix &constraints)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
IndexSet get_view(const size_type begin, const size_type end) const
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)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Point< 3 > point(const gp_Pnt &p)