deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
step-45.h
Go to the documentation of this file.
1 0,
106 /*b_id2*/ 1,
107 /*direction*/ 0,
108 matched_pairs);
109@endcode
110would yield periodicity constraints such that @f$u(0,y)=u(1,y)@f$ for all
111@f$y\in[0,1]@f$.
112
113If we instead consider the parallelogram given by the convex hull of
114@f$(0,0)@f$, @f$(1,1)@f$, @f$(1,2)@f$, @f$(0,1)@f$ we can achieve the constraints
115@f$u(0,y)=u(1,y+1)@f$ by specifying an @p offset:
116@code
118 /*b_id1*/ 0,
119 /*b_id2*/ 1,
120 /*direction*/ 0,
121 matched_pairs,
122 Tensor<1, 2>(0.,1.));
123@endcode
124or
125@code
127 /*b_id1*/ 0,
128 /*b_id2*/ 1,
129 /*arbitrary direction*/ 0,
130 matched_pairs,
131 Tensor<1, 2>(1.,1.));
132@endcode
133Here, again, the assignment of boundary indicators 0 and 1 stems from
134what GridGenerator::parallelogram() documents.
135
136The resulting @p matched_pairs can be used in
137DoFTools::make_periodicity_constraints for populating an AffineConstraints
138object with periodicity constraints:
139@code
140DoFTools::make_periodicity_constraints(matched_pairs, constraints);
141@endcode
142
143Apart from this high level interface there are also variants of
144DoFTools::make_periodicity_constraints available that combine those two
145steps (see the variants of DofTools::make_periodicity_constraints).
146
147There is also a low level interface to
148DoFTools::make_periodicity_constraints if more flexibility is needed. The
149low level variant allows to directly specify two faces that shall be
150constrained:
151@code
152using namespace DoFTools;
153make_periodicity_constraints(face_1,
154 face_2,
155 affine_constraints,
156 component_mask = <default value>;
157 face_orientation = <default value>,
158 face_flip = <default value>,
159 face_rotation = <default value>,
160 matrix = <default value>);
161@endcode
162Here, we need to specify the orientation of the two faces using
163@p face_orientation, @p face_flip and @p face_orientation. For a closer description
164have a look at the documentation of DoFTools::make_periodicity_constraints.
165The remaining parameters are the same as for the high level interface apart
166from the self-explaining @p component_mask and @p affine_constraints.
167
168
169<a name="step-45-problem"></a>
170<a name="step_45-Apracticalexample"></a><h1>A practical example</h1>
171
172
173In the following, we show how to use the above functions in a more involved
174example. The task is to enforce rotated periodicity constraints for the
175velocity component of a Stokes flow.
176
177On a quarter-circle defined by @f$\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}@f$ we are
178going to solve the Stokes problem
179@f{eqnarray*}{
180 -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
181 -\textrm{div}\; \textbf{u}&=&0,\\
182 \textbf{u}|_{\Gamma_1}&=&{\bf 0},
183@f}
184where the boundary @f$\Gamma_1@f$ is defined as @f$\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}@f$.
185For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
186@f{align*}{
187 u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
188 u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
189@f}
190
191The mesh will be generated by GridGenerator::quarter_hyper_shell(),
192which also documents how it assigns boundary indicators to its various
193boundaries if its `colorize` argument is set to `true`.
194 *
195 *
196 * <a name="step_45-CommProg"></a>
197 * <h1> The commented program</h1>
198 *
199 * This example program is a slight modification of @ref step_22 "step-22" running in parallel
200 * using Trilinos to demonstrate the usage of periodic boundary conditions in
201 * deal.II. We thus omit to discuss the majority of the source code and only
202 * comment on the parts that deal with periodicity constraints. For the rest
203 * have a look at @ref step_22 "step-22" and the full source code at the bottom.
204 *
205
206 *
207 * In order to implement periodic boundary conditions only two functions
208 * have to be modified:
209 * - <code>StokesProblem<dim>::setup_dofs()</code>:
210 * To populate an AffineConstraints object with periodicity constraints
211 * - <code>StokesProblem<dim>::create_mesh()</code>:
212 * To supply a distributed triangulation with periodicity information.
213 *
214
215 *
216 * The rest of the program is identical to @ref step_22 "step-22", so let us skip this part
217 * and only show these two functions in the following. (The full program can be
218 * found in the "Plain program" section below, though.)
219 *
220
221 *
222 *
223
224 *
225 *
226
227 *
228 *
229 * <a name="step_45-Settingupperiodicityconstraintsondistributedtriangulations"></a>
230 * <h3>Setting up periodicity constraints on distributed triangulations</h3>
231 *
232 * @code
233 *   template <int dim>
234 *   void StokesProblem<dim>::create_mesh()
235 *   {
236 *   Point<dim> center;
237 *   const double inner_radius = .5;
238 *   const double outer_radius = 1.;
239 *  
241 *   triangulation, center, inner_radius, outer_radius, 0, true);
242 *  
243 * @endcode
244 *
245 * Before we can prescribe periodicity constraints, we need to ensure that
246 * cells on opposite sides of the domain but connected by periodic faces are
247 * part of the ghost layer if one of them is stored on the local processor.
248 * At this point we need to think about how we want to prescribe
249 * periodicity. The vertices @f$\text{vertices}_2@f$ of a face on the left
250 * boundary should be matched to the vertices @f$\text{vertices}_1@f$ of a face
251 * on the lower boundary given by @f$\text{vertices}_2=R\cdot
252 * \text{vertices}_1+b@f$ where the rotation matrix @f$R@f$ and the offset @f$b@f$ are
253 * given by
254 * @f{align*}{
255 * R=\begin{pmatrix}
256 * 0&1\\-1&0
257 * \end{pmatrix},
258 * \quad
259 * b=\begin{pmatrix}0&0\end{pmatrix}.
260 * @f}
261 * The data structure we are saving the resulting information into is here
262 * based on the Triangulation.
263 *
264 * @code
265 *   std::vector<GridTools::PeriodicFacePair<
267 *   periodicity_vector;
268 *  
269 *   FullMatrix<double> rotation_matrix(dim);
270 *   rotation_matrix[0][1] = 1.;
271 *   rotation_matrix[1][0] = -1.;
272 *  
274 *   2,
275 *   3,
276 *   1,
277 *   periodicity_vector,
278 *   Tensor<1, dim>(),
279 *   rotation_matrix);
280 *  
281 * @endcode
282 *
283 * Now telling the triangulation about the desired periodicity is
284 * particularly easy by just calling
286 *
287 * @code
288 *   triangulation.add_periodicity(periodicity_vector);
289 *  
290 *   triangulation.refine_global(4 - dim);
291 *   }
292 *  
293 *  
294 *   template <int dim>
295 *   void StokesProblem<dim>::setup_dofs()
296 *   {
297 *   dof_handler.distribute_dofs(fe);
298 *  
299 *   std::vector<unsigned int> block_component(dim + 1, 0);
300 *   block_component[dim] = 1;
301 *   DoFRenumbering::component_wise(dof_handler, block_component);
302 *  
303 *   const std::vector<types::global_dof_index> dofs_per_block =
304 *   DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
305 *   const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
306 *  
307 *   {
308 *   owned_partitioning.clear();
309 *   const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
310 *   owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
311 *   owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u + n_p));
312 *  
313 *   relevant_partitioning.clear();
314 *   const IndexSet locally_relevant_dofs =
316 *   relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
317 *   relevant_partitioning.push_back(
318 *   locally_relevant_dofs.get_view(n_u, n_u + n_p));
319 *  
320 *   constraints.clear();
321 *   constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
322 *  
323 *   const FEValuesExtractors::Vector velocities(0);
324 *  
325 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
327 *   dof_handler,
328 *   0,
329 *   BoundaryValues<dim>(),
330 *   constraints,
331 *   fe.component_mask(velocities));
333 *   dof_handler,
334 *   1,
335 *   BoundaryValues<dim>(),
336 *   constraints,
337 *   fe.component_mask(velocities));
338 *  
339 * @endcode
340 *
341 * After we provided the mesh with the necessary information for the
342 * periodicity constraints, we are now able to actual create them. For
343 * describing the matching we are using the same approach as before, i.e.,
344 * the @f$\text{vertices}_2@f$ of a face on the left boundary should be
345 * matched to the vertices
346 * @f$\text{vertices}_1@f$ of a face on the lower boundary given by
347 * @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ where the rotation
348 * matrix @f$R@f$ and the offset @f$b@f$ are given by
349 * @f{align*}{
350 * R=\begin{pmatrix}
351 * 0&1\\-1&0
352 * \end{pmatrix},
353 * \quad
354 * b=\begin{pmatrix}0&0\end{pmatrix}.
355 * @f}
356 * These two objects not only describe how faces should be matched but
357 * also in which sense the solution should be transformed from
358 * @f$\text{face}_2@f$ to
359 * @f$\text{face}_1@f$.
360 *
361 * @code
362 *   FullMatrix<double> rotation_matrix(dim);
363 *   rotation_matrix[0][1] = 1.;
364 *   rotation_matrix[1][0] = -1.;
365 *  
366 *   Tensor<1, dim> offset;
367 *  
368 * @endcode
369 *
370 * For setting up the constraints, we first store the periodicity
371 * information in an auxiliary object of type
372 * <code>std::vector@<GridTools::PeriodicFacePair<typename
373 * DoFHandler@<dim@>::%cell_iterator@> </code>. The periodic boundaries
374 * have the boundary indicators 2 (x=0) and 3 (y=0). All the other
375 * parameters we have set up before. In this case the direction does not
376 * matter. Due to @f$\text{vertices}_2=R\cdot \text{vertices}_1+b@f$ this is
377 * exactly what we want.
378 *
379 * @code
380 *   std::vector<
382 *   periodicity_vector;
383 *  
384 *   const unsigned int direction = 1;
385 *  
386 *   GridTools::collect_periodic_faces(dof_handler,
387 *   2,
388 *   3,
389 *   direction,
390 *   periodicity_vector,
391 *   offset,
392 *   rotation_matrix);
393 *  
394 * @endcode
395 *
396 * Next, we need to provide information on which vector valued components
397 * of the solution should be rotated. Since we choose here to just
398 * constraint the velocity and this starts at the first component of the
399 * solution vector, we simply insert a 0:
400 *
401 * @code
402 *   std::vector<unsigned int> first_vector_components;
403 *   first_vector_components.push_back(0);
404 *  
405 * @endcode
406 *
407 * After setting up all the information in periodicity_vector all we have
408 * to do is to tell make_periodicity_constraints to create the desired
409 * constraints.
410 *
411 * @code
412 *   DoFTools::make_periodicity_constraints<dim, dim>(periodicity_vector,
413 *   constraints,
414 *   fe.component_mask(
415 *   velocities),
416 *   first_vector_components);
417 *   }
418 *  
419 *   constraints.close();
420 *  
421 *   {
422 *   TrilinosWrappers::BlockSparsityPattern bsp(owned_partitioning,
423 *   owned_partitioning,
424 *   relevant_partitioning,
425 *   mpi_communicator);
426 *  
427 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
428 *   for (unsigned int c = 0; c < dim + 1; ++c)
429 *   for (unsigned int d = 0; d < dim + 1; ++d)
430 *   if (!((c == dim) && (d == dim)))
431 *   coupling[c][d] = DoFTools::always;
432 *   else
433 *   coupling[c][d] = DoFTools::none;
434 *  
435 *   DoFTools::make_sparsity_pattern(dof_handler,
436 *   coupling,
437 *   bsp,
438 *   constraints,
439 *   false,
441 *   mpi_communicator));
442 *  
443 *   bsp.compress();
444 *  
445 *   system_matrix.reinit(bsp);
446 *   }
447 *  
448 *   {
449 *   TrilinosWrappers::BlockSparsityPattern preconditioner_bsp(
450 *   owned_partitioning,
451 *   owned_partitioning,
452 *   relevant_partitioning,
453 *   mpi_communicator);
454 *  
455 *   Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
456 *   for (unsigned int c = 0; c < dim + 1; ++c)
457 *   for (unsigned int d = 0; d < dim + 1; ++d)
458 *   if ((c == dim) && (d == dim))
459 *   preconditioner_coupling[c][d] = DoFTools::always;
460 *   else
461 *   preconditioner_coupling[c][d] = DoFTools::none;
462 *  
463 *   DoFTools::make_sparsity_pattern(dof_handler,
464 *   preconditioner_coupling,
465 *   preconditioner_bsp,
466 *   constraints,
467 *   false,
469 *   mpi_communicator));
470 *  
471 *   preconditioner_bsp.compress();
472 *  
473 *   preconditioner_matrix.reinit(preconditioner_bsp);
474 *   }
475 *  
476 *   system_rhs.reinit(owned_partitioning, mpi_communicator);
477 *   solution.reinit(owned_partitioning,
478 *   relevant_partitioning,
479 *   mpi_communicator);
480 *   }
481 *  
482 * @endcode
483 *
484 * The rest of the program is then again identical to @ref step_22 "step-22". We will omit
485 * it here now, but as before, you can find these parts in the "Plain program"
486 * section below.
487 *
488
489 *
490<a name="step_45-Results"></a><h1>Results</h1>
491
492
493The created output is not very surprising. We simply see that the solution is
494periodic with respect to the left and lower boundary:
495
496<img src="https://www.dealii.org/images/steps/developer/step-45.periodic.png" alt="">
497
498Without the periodicity constraints we would have ended up with the following solution:
499
500<img src="https://www.dealii.org/images/steps/developer/step-45.non_periodic.png" alt="">
501 *
502 *
503<a name="step_45-PlainProg"></a>
504<h1> The plain program</h1>
505@include "step-45.cc"
506*/
void clear()
Definition index_set.h:1752
Definition point.h:111
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator > > &) override
Definition tria.cc:3768
Point< 3 > center
Point< 3 > vertices[4]
bool colorize
Definition grid_out.cc:4625
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition tria.h:288
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, 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 >())
void make_periodicity_constraints(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, AffineConstraints< number > &constraints, const ComponentMask &component_mask={}, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation(), const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >(), const number periodicity_factor=1.)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
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 >())
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 unsigned 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.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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:107
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={})
int(&) functions(const void *v1, const void *v2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation