1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_mg_constrained_dofs_h
17 #define dealii_mg_constrained_dofs_h
18
19 #include <deal.II/base/config.h>
20
21 #include <deal.II/base/subscriptor.h>
22
23 #include <deal.II/lac/affine_constraints.h>
24
25 #include <deal.II/multigrid/mg_tools.h>
26
27 #include <set>
28 #include <vector>
29
30 DEAL_II_NAMESPACE_OPEN
31
32 // Forward declaration
33 #ifndef DOXYGEN
34 template <int dim, int spacedim>
35 class DoFHandler;
36 #endif
37
38
39 /**
40 * Collection of boundary constraints and refinement edge constraints for
41 * level vectors.
42 *
43 * @ingroup mg
44 */
45 class MGConstrainedDoFs : public Subscriptor
46 {
47 public:
48 using size_dof = std::vector<std::set<types::global_dof_index>>::size_type;
49 /**
50 * Fill the internal data structures with hanging node constraints extracted
51 * from the dof handler object. Works with natural boundary conditions only.
52 * There exists a sister function setting up boundary constraints as well.
53 *
54 * This function ensures that on every level, degrees of freedom at interior
55 * edges of a refinement level are treated corrected but leaves degrees of
56 * freedom at the boundary of the domain untouched assuming that no
57 * Dirichlet boundary conditions for them exist.
58 *
59 * Furthermore, this call sets up an AffineConstraints object on each
60 * level that contains possible periodicity constraints in case those
61 * have been added to the underlying triangulation. The AffineConstraints
62 * object can be queried by get_level_constraints(level). Note that the
63 * current implementation of periodicity constraints in this class does
64 * not support rotation matrices in the periodicity definition, i.e., the
65 * respective argument in the GridTools::collect_periodic_faces() may not
66 * be different from the identity matrix.
67 */
68 template <int dim, int spacedim>
69 void
70 initialize(const DoFHandler<dim, spacedim> &dof);
71
72 /**
73 * Fill the internal data structures with information
74 * about Dirichlet boundary dofs.
75 *
76 * The initialize() function has to be called before
77 * to set hanging node constraints.
78 *
79 * This function can be called multiple times to allow considering
80 * different sets of boundary_ids for different components.
81 */
82 template <int dim, int spacedim>
83 void
84 make_zero_boundary_constraints(
85 const DoFHandler<dim, spacedim> & dof,
86 const std::set<types::boundary_id> &boundary_ids,
87 const ComponentMask & component_mask = ComponentMask());
88
89 /**
90 * Add user defined constraints to be used on level @p level.
91 *
92 * The user can call this function multiple times and any new,
93 * conflicting constraints will overwrite the previous constraints
94 * for that DoF.
95 *
96 * Before the transfer, the user defined constraints will be distributed
97 * to the source vector, and then any DoF index set using
98 * make_zero_boundary_constraints() will be overwritten with
99 * value zero.
100 *
101 * @note This is currently only implemented for MGTransferMatrixFree.
102 */
103 void
104 add_user_constraints(const unsigned int level,
105 const AffineConstraints<double> &constraints_on_level);
106
107 /**
108 * Fill the internal data structures with information
109 * about no normal flux boundary dofs.
110 *
111 * This function is limited to meshes whose no normal flux boundaries
112 * have faces which are normal to the x-, y-, or z-axis. Also, for a
113 * specific boundary id, all faces must be facing in the same direction,
114 * i.e., a boundary normal to the x-axis must have a different boundary
115 * id than a boundary normal to the y- or z-axis and so on. If the mesh
116 * was produced, for example, using the <tt>GridGenerator::hyper_cube()</tt>
117 * function, setting <tt>colorize=true</tt> during mesh generation and calling
118 * make_no_normal_flux_constraints() for each no normal flux boundary is
119 * sufficient.
120 */
121 template <int dim, int spacedim>
122 void
123 make_no_normal_flux_constraints(const DoFHandler<dim, spacedim> &dof,
124 const types::boundary_id bid,
125 const unsigned int first_vector_component);
126
127 /**
128 * Reset the data structures.
129 */
130 void
131 clear();
132
133 /**
134 * Determine whether a dof index is subject to a boundary constraint.
135 */
136 bool
137 is_boundary_index(const unsigned int level,
138 const types::global_dof_index index) const;
139
140 /**
141 * Determine whether a dof index is at the refinement edge.
142 */
143 bool
144 at_refinement_edge(const unsigned int level,
145 const types::global_dof_index index) const;
146
147
148 /**
149 * Determine whether the (i,j) entry of the interface matrix
150 * on a given level should be set. This is taken in terms of
151 * dof i, that is, return true if i is at a refinement edge,
152 * j is not, and both are not on the external boundary.
153 */
154 bool
155 is_interface_matrix_entry(const unsigned int level,
156 const types::global_dof_index i,
157 const types::global_dof_index j) const;
158
159 /**
160 * Return the indices of level dofs on the given level that are subject to
161 * Dirichlet boundary conditions (as set by the @p function_map parameter in
162 * initialize()). The indices are restricted to the set of locally relevant
163 * level dofs.
164 */
165 const IndexSet &
166 get_boundary_indices(const unsigned int level) const;
167
168
169 /**
170 * Return the indices of dofs on the given level that lie on an refinement
171 * edge (dofs on faces to neighbors that are coarser).
172 */
173 const IndexSet &
174 get_refinement_edge_indices(unsigned int level) const;
175
176
177 /**
178 * Return if Dirichlet boundary indices are set in initialize().
179 */
180 bool
181 have_boundary_indices() const;
182
183 /**
184 * Return the AffineConstraints object for a given level, containing
185 * periodicity constraints (if enabled on the triangulation).
186 */
187 const AffineConstraints<double> &
188 get_level_constraints(const unsigned int level) const;
189
190 /**
191 * Return the user defined constraint matrix for a given level. These
192 * constraints are set using the function add_user_constraints() and
193 * should not contain constraints for DoF indices set in
194 * make_zero_boundary_constraints() as they will be overwritten during
195 * the transfer.
196 */
197 const AffineConstraints<double> &
198 get_user_constraint_matrix(const unsigned int level) const;
199
200 private:
201 /**
202 * The indices of boundary dofs for each level.
203 */
204 std::vector<IndexSet> boundary_indices;
205
206 /**
207 * The degrees of freedom on a given level that live on the refinement edge
208 * between the level and cells on a coarser level.
209 */
210 std::vector<IndexSet> refinement_edge_indices;
211
212 /**
213 * Constraint matrices containing information regarding potential
214 * periodic boundary conditions for each level .
215 */
216 std::vector<AffineConstraints<double>> level_constraints;
217
218 /**
219 * Constraint matrices defined by user.
220 */
221 std::vector<AffineConstraints<double>> user_constraints;
222 };
223
224
225 template <int dim, int spacedim>
226 inline void
initialize(const DoFHandler<dim,spacedim> & dof)227 MGConstrainedDoFs::initialize(const DoFHandler<dim, spacedim> &dof)
228 {
229 boundary_indices.clear();
230 refinement_edge_indices.clear();
231 level_constraints.clear();
232 user_constraints.clear();
233
234 const unsigned int nlevels = dof.get_triangulation().n_global_levels();
235
236 // At this point level_constraint and refinement_edge_indices are empty.
237 refinement_edge_indices.resize(nlevels);
238 level_constraints.resize(nlevels);
239 user_constraints.resize(nlevels);
240 for (unsigned int l = 0; l < nlevels; ++l)
241 {
242 IndexSet relevant_dofs;
243 DoFTools::extract_locally_relevant_level_dofs(dof, l, relevant_dofs);
244 level_constraints[l].reinit(relevant_dofs);
245
246 // Loop through relevant cells and faces finding those which are periodic
247 // neighbors.
248 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
249 endc = dof.end(l);
250 for (; cell != endc; ++cell)
251 if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
252 {
253 for (auto f : cell->face_indices())
254 if (cell->has_periodic_neighbor(f) &&
255 cell->periodic_neighbor(f)->level() == cell->level())
256 {
257 if (cell->is_locally_owned_on_level())
258 {
259 Assert(
260 cell->periodic_neighbor(f)->level_subdomain_id() !=
261 numbers::artificial_subdomain_id,
262 ExcMessage(
263 "Periodic neighbor of a locally owned cell must either be owned or ghost."));
264 }
265 // Cell is a level-ghost and its neighbor is a
266 // level-artificial cell nothing to do here
267 else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
268 numbers::artificial_subdomain_id)
269 {
270 Assert(cell->is_locally_owned_on_level() == false,
271 ExcInternalError());
272 continue;
273 }
274
275 const unsigned int dofs_per_face =
276 cell->face(f)->get_fe(0).n_dofs_per_face();
277 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
278 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
279
280 cell->periodic_neighbor(f)
281 ->face(cell->periodic_neighbor_face_no(f))
282 ->get_mg_dof_indices(l, dofs_1, 0);
283 cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
284 // Store periodicity information in the level
285 // AffineConstraints object. Skip DoFs for which we've
286 // previously entered periodicity constraints already; this
287 // can happen, for example, for a vertex dof at a periodic
288 // boundary that we visit from more than one cell
289 for (unsigned int i = 0; i < dofs_per_face; ++i)
290 if (level_constraints[l].can_store_line(dofs_2[i]) &&
291 level_constraints[l].can_store_line(dofs_1[i]) &&
292 !level_constraints[l].is_constrained(dofs_2[i]) &&
293 !level_constraints[l].is_constrained(dofs_1[i]))
294 {
295 level_constraints[l].add_line(dofs_2[i]);
296 level_constraints[l].add_entry(dofs_2[i],
297 dofs_1[i],
298 1.);
299 }
300 }
301 }
302 level_constraints[l].close();
303
304 // Initialize with empty IndexSet of correct size
305 refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
306 }
307
308 MGTools::extract_inner_interface_dofs(dof, refinement_edge_indices);
309 }
310
311
312
313 template <int dim, int spacedim>
314 inline void
make_zero_boundary_constraints(const DoFHandler<dim,spacedim> & dof,const std::set<types::boundary_id> & boundary_ids,const ComponentMask & component_mask)315 MGConstrainedDoFs::make_zero_boundary_constraints(
316 const DoFHandler<dim, spacedim> & dof,
317 const std::set<types::boundary_id> &boundary_ids,
318 const ComponentMask & component_mask)
319 {
320 // allocate an IndexSet for each global level. Contents will be
321 // overwritten inside make_boundary_list.
322 const unsigned int n_levels = dof.get_triangulation().n_global_levels();
323 Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
324 ExcInternalError());
325 boundary_indices.resize(n_levels);
326
327 MGTools::make_boundary_list(dof,
328 boundary_ids,
329 boundary_indices,
330 component_mask);
331 }
332
333
334 template <int dim, int spacedim>
335 inline void
make_no_normal_flux_constraints(const DoFHandler<dim,spacedim> & dof,const types::boundary_id bid,const unsigned int first_vector_component)336 MGConstrainedDoFs::make_no_normal_flux_constraints(
337 const DoFHandler<dim, spacedim> &dof,
338 const types::boundary_id bid,
339 const unsigned int first_vector_component)
340 {
341 // For a given boundary id, find which vector component is on the boundary
342 // and set a zero boundary constraint for those degrees of freedom.
343 const unsigned int n_components = dof.get_fe_collection().n_components();
344 AssertIndexRange(first_vector_component + dim - 1, n_components);
345
346 ComponentMask comp_mask(n_components, false);
347
348
349 typename Triangulation<dim>::face_iterator
350 face = dof.get_triangulation().begin_face(),
351 endf = dof.get_triangulation().end_face();
352 for (; face != endf; ++face)
353 if (face->at_boundary() && face->boundary_id() == bid)
354 for (unsigned int d = 0; d < dim; ++d)
355 {
356 Tensor<1, dim, double> unit_vec;
357 unit_vec[d] = 1.0;
358
359 const Tensor<1, dim> normal_vec =
360 face->get_manifold().normal_vector(face, face->center());
361
362 if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
363 comp_mask.set(d + first_vector_component, true);
364 else
365 Assert(
366 std::abs(unit_vec * normal_vec) < 1e-10,
367 ExcMessage(
368 "We can currently only support no normal flux conditions "
369 "for a specific boundary id if all faces are normal to the "
370 "x, y, or z axis."));
371 }
372
373 Assert(comp_mask.n_selected_components() == 1,
374 ExcMessage(
375 "We can currently only support no normal flux conditions "
376 "for a specific boundary id if all faces are facing in the "
377 "same direction, i.e., a boundary normal to the x-axis must "
378 "have a different boundary id than a boundary normal to the "
379 "y- or z-axis and so on. If the mesh here was produced using "
380 "GridGenerator::..., setting colorize=true during mesh generation "
381 "and calling make_no_normal_flux_constraints() for each no normal "
382 "flux boundary will fulfill the condition."));
383
384 this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
385 }
386
387
388 inline void
add_user_constraints(const unsigned int level,const AffineConstraints<double> & constraints_on_level)389 MGConstrainedDoFs::add_user_constraints(
390 const unsigned int level,
391 const AffineConstraints<double> &constraints_on_level)
392 {
393 AssertIndexRange(level, user_constraints.size());
394
395 // Get the relevant DoFs from level_constraints if
396 // the user constraint matrix has not been initialized
397 if (user_constraints[level].get_local_lines().size() == 0)
398 user_constraints[level].reinit(level_constraints[level].get_local_lines());
399
400 user_constraints[level].merge(
401 constraints_on_level,
402 AffineConstraints<double>::MergeConflictBehavior::right_object_wins);
403 user_constraints[level].close();
404 }
405
406
407 inline void
clear()408 MGConstrainedDoFs::clear()
409 {
410 boundary_indices.clear();
411 refinement_edge_indices.clear();
412 user_constraints.clear();
413 }
414
415
416 inline bool
is_boundary_index(const unsigned int level,const types::global_dof_index index)417 MGConstrainedDoFs::is_boundary_index(const unsigned int level,
418 const types::global_dof_index index) const
419 {
420 if (boundary_indices.size() == 0)
421 return false;
422
423 AssertIndexRange(level, boundary_indices.size());
424 return boundary_indices[level].is_element(index);
425 }
426
427 inline bool
at_refinement_edge(const unsigned int level,const types::global_dof_index index)428 MGConstrainedDoFs::at_refinement_edge(const unsigned int level,
429 const types::global_dof_index index) const
430 {
431 AssertIndexRange(level, refinement_edge_indices.size());
432
433 return refinement_edge_indices[level].is_element(index);
434 }
435
436 inline bool
is_interface_matrix_entry(const unsigned int level,const types::global_dof_index i,const types::global_dof_index j)437 MGConstrainedDoFs::is_interface_matrix_entry(
438 const unsigned int level,
439 const types::global_dof_index i,
440 const types::global_dof_index j) const
441 {
442 const IndexSet &interface_dofs_on_level =
443 this->get_refinement_edge_indices(level);
444
445 return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
446 && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
447 && !this->is_boundary_index(level, i) // !on_boundary(i)
448 && !this->is_boundary_index(level, j); // !on_boundary(j)
449 }
450
451
452
453 inline const IndexSet &
get_boundary_indices(const unsigned int level)454 MGConstrainedDoFs::get_boundary_indices(const unsigned int level) const
455 {
456 AssertIndexRange(level, boundary_indices.size());
457 return boundary_indices[level];
458 }
459
460
461
462 inline const IndexSet &
get_refinement_edge_indices(unsigned int level)463 MGConstrainedDoFs::get_refinement_edge_indices(unsigned int level) const
464 {
465 AssertIndexRange(level, refinement_edge_indices.size());
466 return refinement_edge_indices[level];
467 }
468
469
470
471 inline bool
have_boundary_indices()472 MGConstrainedDoFs::have_boundary_indices() const
473 {
474 return boundary_indices.size() != 0;
475 }
476
477
478
479 inline const AffineConstraints<double> &
get_level_constraints(const unsigned int level)480 MGConstrainedDoFs::get_level_constraints(const unsigned int level) const
481 {
482 AssertIndexRange(level, level_constraints.size());
483 return level_constraints[level];
484 }
485
486
487
488 inline const AffineConstraints<double> &
get_user_constraint_matrix(const unsigned int level)489 MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const
490 {
491 AssertIndexRange(level, user_constraints.size());
492 return user_constraints[level];
493 }
494
495
496
497 DEAL_II_NAMESPACE_CLOSE
498
499 #endif
500