100 would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$
for all
103 If we instead consider the
parallelogram given by the convex hull of
104 @f$(0,0)@f$, @f$(1,1)@f$, @f$(1,2)@f$, @f$(0,1)@f$ we can achieve the constraints
105 @f$u(0,y)=u(1,y+1)@f$ by specifying an @p offset:
124 The resulting @p matched_pairs can be used in
126 with periodicity constraints:
131 Apart from
this high level
interface there are also variants of
133 steps (see the variants of DofTools::make_periodicity_constraints).
135 There is also a low level interface to
137 low level variant allows to directly specify two faces that shall be
144 component_mask = <
default value>;
145 face_orientation = <
default value>,
146 face_flip = <
default value>,
147 face_rotation = <
default value>,
148 matrix = <
default value>);
150 Here, we need to specify the orientation of the two faces
using 151 @p face_orientation, @p face_flip and @p face_orientation. For a closer description
153 The remaining parameters are the same as
for the high level
interface apart
154 from the self-explaining @p component_mask and @p constraint_matrix.
156 <a name=
"problem"></a>
157 <a name=
"Apracticalexample"></a><h1>A practical example</h1>
160 In the following, we show how to use the above functions in a more involved
161 example. The task is to enforce rotated periodicity constraints
for the
162 velocity component of a Stokes flow.
164 On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
165 going to solve the Stokes problem
167 -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100*\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
168 -\textrm{div}\; \textbf{u}&=&0,\\
169 \textbf{u}|_{\Gamma_1}&=&{\bf 0},
171 where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1:=\{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
172 For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
174 u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
175 u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
177 * <a name=
"CommProg"></a>
178 * <h1> The commented program</h1>
180 * This example program is a slight modification of @ref step_22
"step-22" running in
parallel 181 *
using Trilinos to demonstrate the usage of periodic boundary conditions in deal.II.
182 * We thus omit to discuss the majority of the source code and only comment on the
183 * parts that deal with periodicity constraints. For the rest have a look at @ref step_22
"step-22" 184 * and the full source code at the bottom.
188 * In order to implement periodic boundary conditions only two functions
189 * have to be modified:
190 * - <code>StokesProblem<dim>::setup_dofs()</code>: To populate a
ConstraintMatrix 191 *
object with periodicity constraints
192 * - <code>StokesProblem<dim>::run()</code>: To supply a distributed triangulation with
193 * periodicity information.
203 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
204 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
208 *
void StokesProblem<dim>::create_mesh()
211 *
const double inner_radius = .5;
212 *
const double outer_radius = 1.;
223 * Before we can prescribe periodicity constraints, we need to ensure that cells
224 * on opposite sides of the domain but connected by periodic faces are part of
225 * the ghost layer
if one of them is stored on the local processor.
226 * At
this point we need to think about how we want to prescribe periodicity.
227 * The vertices @f$\text{vertices}_2@f$ of a face on the left boundary should be
228 * matched to the vertices @f$\text{vertices}_1@f$ of a face on the lower boundary
229 * given by @f$\text{vertices}_2=R\cdot \text{vertices}_1+
b@f$ where the rotation
230 *
matrix @f$R@f$ and the offset @f$b@f$ are given by
236 *
b=\begin{pmatrix}0&0\end{pmatrix}.
238 * The data structure we are saving the reuslitng information into is here based
242 * std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
243 * periodicity_vector;
246 * rotation_matrix[0][1]=1.;
247 * rotation_matrix[1][0]=-1.;
255 * Now telling the triangulation about the desired periodicity is
256 * particularly easy by just calling
260 * triangulation.add_periodicity(periodicity_vector);
262 * triangulation.refine_global (4-dim);
269 * <a name=
"Settingupperiodicityconstraintsondistributedtriangulations"></a>
270 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
274 *
void StokesProblem<dim>::setup_dofs ()
276 * dof_handler.distribute_dofs (fe);
278 * std::vector<unsigned int> block_component (dim+1,0);
279 * block_component[dim] = 1;
282 * std::vector<types::global_dof_index> dofs_per_block (2);
284 *
const unsigned int n_u = dofs_per_block[0],
285 * n_p = dofs_per_block[1];
288 * owned_partitioning.clear();
289 *
IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
290 * owned_partitioning.push_back(locally_owned_dofs.
get_view(0, n_u));
291 * owned_partitioning.push_back(locally_owned_dofs.
get_view(n_u, n_u+n_p));
293 * relevant_partitioning.
clear();
296 * relevant_partitioning.push_back(locally_relevant_dofs.
get_view(0, n_u));
297 * relevant_partitioning.push_back(locally_relevant_dofs.
get_view(n_u, n_u+n_p));
299 * constraints.clear ();
300 * constraints.reinit(locally_relevant_dofs);
309 * BoundaryValues<dim>(),
311 * fe.component_mask(velocities));
315 * BoundaryValues<dim>(),
317 * fe.component_mask(velocities));
321 * After we provided the mesh with the necessary information
for the periodicity
322 * constraints, we are now able to actual create them. For describing the
323 * matching we are
using the same approach as before, i.e., the @f$\text{vertices}_2@f$
324 * of a face on the left boundary should be matched to the vertices
325 * @f$\text{vertices}_1@f$ of a face on the lower boundary given by
326 * @f$\text{vertices}_2=R\cdot \text{vertices}_1+
b@f$ where the rotation
matrix @f$R@f$
327 * and the offset @f$b@f$ are given by
333 *
b=\begin{pmatrix}0&0\end{pmatrix}.
335 * These two objects not only describe how faces should be matched but also
336 * in which sense the solution should be transformed from @f$\text{face}_2@f$ to
337 * @f$\text{face}_1@f$.
341 * rotation_matrix[0][1]=1.;
342 * rotation_matrix[1][0]=-1.;
348 * For setting up the constraints, we first store the periodicity
349 * information in an auxiliary
object of type
351 *
DoFHandler@<dim@>::cell_iterator@> </code>. The periodic boundaries have the
352 * boundary indicators 2 (x=0) and 3 (y=0). All the other parameters we
353 * have
set up before. In
this case the direction does not matter. Due to
354 * @f$\text{vertices}_2=R\cdot \text{vertices}_1+
b@f$
this is exactly what we want.
357 * std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
358 * periodicity_vector;
360 *
const unsigned int direction = 1;
363 * periodicity_vector, offset,
368 * Next we need to provide information on which vector valued components of
369 * the solution should be rotated. Since we choose here to just constraint the
370 * velocity and
this starts at the first component of the solution vector we
374 * std::vector<unsigned int> first_vector_components;
375 * first_vector_components.push_back(0);
379 * After setting up all the information in periodicity_vector all we have
384 * DoFTools::make_periodicity_constraints<DoFHandler<dim> >
385 * (periodicity_vector, constraints, fe.component_mask(velocities),
386 * first_vector_components);
391 * BoundaryValues<dim>(),
393 * fe.component_mask(velocities));
397 * BoundaryValues<dim>(),
399 * fe.component_mask(velocities));
403 * constraints.close ();
407 * (owned_partitioning, owned_partitioning,
408 * relevant_partitioning, mpi_communicator);
411 * (dof_handler, bsp, constraints,
false,
416 * system_matrix.reinit (bsp);
419 * system_rhs.reinit (owned_partitioning,
421 * solution.reinit (owned_partitioning, relevant_partitioning,
428 <a name=
"Results"></a><h1>Results</h1>
431 The created output is not very surprising. We simply see that the solution is
432 periodic with respect to the left and lower boundary:
434 <img src=
"/documentation/images/steps/developer/step-45.periodic.png" alt=
"">
436 Without the periodicity constraints we would have ended up with the following solution:
438 <img src=
"/documentation/images/steps/developer/step-45.non_periodic.png" alt=
"">
439 <a name=
"PlainProg"></a>
440 <h1> The plain program</h1>
441 @include
"step-45.cc" Contents is actually a matrix.
Point< spacedim > point(const gp_Pnt &p, const double &tolerance=1e-10)
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)
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)