Reference documentation for deal.II version GIT 3779fa9eb4 2023-09-28 13:00:02+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\}}\)
mg_constrained_dofs.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2023 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 
23 
25 
27 
28 #include <set>
29 #include <vector>
30 
32 
33 // Forward declaration
34 #ifndef DOXYGEN
35 template <int dim, int spacedim>
36 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
37 class DoFHandler;
38 #endif
39 
40 
48 {
49 public:
50  using size_dof = std::vector<std::set<types::global_dof_index>>::size_type;
76  template <int dim, int spacedim>
77  void
79  const MGLevelObject<IndexSet> &level_relevant_dofs =
81 
92  template <int dim, int spacedim>
93  void
95  const DoFHandler<dim, spacedim> &dof,
96  const std::set<types::boundary_id> &boundary_ids,
97  const ComponentMask &component_mask = {});
98 
105  template <int dim, int spacedim>
106  void
108  const unsigned int level,
109  const IndexSet &boundary_indices);
110 
125  void
126  add_user_constraints(const unsigned int level,
127  const AffineConstraints<double> &constraints_on_level);
128 
143  template <int dim, int spacedim>
144  void
146  const types::boundary_id bid,
147  const unsigned int first_vector_component);
148 
152  void
154 
158  void
159  clear();
160 
164  bool
165  is_boundary_index(const unsigned int level,
166  const types::global_dof_index index) const;
167 
171  bool
172  at_refinement_edge(const unsigned int level,
173  const types::global_dof_index index) const;
174 
175 
182  bool
183  is_interface_matrix_entry(const unsigned int level,
184  const types::global_dof_index i,
185  const types::global_dof_index j) const;
186 
193  const IndexSet &
194  get_boundary_indices(const unsigned int level) const;
195 
196 
201  const IndexSet &
202  get_refinement_edge_indices(unsigned int level) const;
203 
204 
208  bool
209  have_boundary_indices() const;
210 
216  get_level_constraints(const unsigned int level) const;
217 
226  get_user_constraint_matrix(const unsigned int level) const;
227 
240  template <typename Number>
241  void
243  const unsigned int level,
244  const bool add_boundary_indices,
245  const bool add_refinement_edge_indices,
246  const bool add_level_constraints,
247  const bool add_user_constraints) const;
248 
249 private:
253  std::vector<IndexSet> boundary_indices;
254 
259  std::vector<IndexSet> refinement_edge_indices;
260 
265  std::vector<AffineConstraints<double>> level_constraints;
266 
270  std::vector<AffineConstraints<double>> user_constraints;
271 };
272 
273 
274 template <int dim, int spacedim>
275 inline void
277  const DoFHandler<dim, spacedim> &dof,
278  const MGLevelObject<IndexSet> &level_relevant_dofs)
279 {
280  boundary_indices.clear();
281  refinement_edge_indices.clear();
282  level_constraints.clear();
283  user_constraints.clear();
284 
285  const unsigned int nlevels = dof.get_triangulation().n_global_levels();
286  const unsigned int min_level = level_relevant_dofs.min_level();
287  const unsigned int max_level = (level_relevant_dofs.max_level() == 0) ?
288  nlevels - 1 :
289  level_relevant_dofs.max_level();
290  const bool user_level_dofs =
291  (level_relevant_dofs.max_level() == 0) ? false : true;
292 
293  // At this point level_constraint and refinement_edge_indices are empty.
294  refinement_edge_indices.resize(nlevels);
295  level_constraints.resize(nlevels);
296  user_constraints.resize(nlevels);
297  for (unsigned int l = min_level; l <= max_level; ++l)
298  {
299  if (user_level_dofs)
300  {
302  level_relevant_dofs[l]);
303  }
304  else
305  {
306  const IndexSet relevant_dofs =
309  relevant_dofs);
310  }
311 
312  // Loop through relevant cells and faces finding those which are periodic
313  // neighbors.
314  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
315  endc = dof.end(l);
316  for (; cell != endc; ++cell)
317  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
318  {
319  for (auto f : cell->face_indices())
320  if (cell->has_periodic_neighbor(f) &&
321  cell->periodic_neighbor(f)->level() == cell->level())
322  {
323  if (cell->is_locally_owned_on_level())
324  {
325  Assert(
326  cell->periodic_neighbor(f)->level_subdomain_id() !=
328  ExcMessage(
329  "Periodic neighbor of a locally owned cell must either be owned or ghost."));
330  }
331  // Cell is a level-ghost and its neighbor is a
332  // level-artificial cell nothing to do here
333  else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
335  {
336  Assert(cell->is_locally_owned_on_level() == false,
337  ExcInternalError());
338  continue;
339  }
340 
341  const unsigned int dofs_per_face =
342  dof.get_fe(0).n_dofs_per_face(f);
343  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
344  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
345 
346  cell->periodic_neighbor(f)
347  ->face(cell->periodic_neighbor_face_no(f))
348  ->get_mg_dof_indices(l, dofs_1, 0);
349  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
350  // Store periodicity information in the level
351  // AffineConstraints object. Skip DoFs for which we've
352  // previously entered periodicity constraints already; this
353  // can happen, for example, for a vertex dof at a periodic
354  // boundary that we visit from more than one cell
355  for (unsigned int i = 0; i < dofs_per_face; ++i)
356  if (level_constraints[l].can_store_line(dofs_2[i]) &&
357  level_constraints[l].can_store_line(dofs_1[i]) &&
358  !level_constraints[l].is_constrained(dofs_2[i]) &&
359  !level_constraints[l].is_constrained(dofs_1[i]))
360  {
361  level_constraints[l].add_line(dofs_2[i]);
362  level_constraints[l].add_entry(dofs_2[i],
363  dofs_1[i],
364  1.);
365  }
366  }
367  }
368  level_constraints[l].close();
369 
370  // Initialize with empty IndexSet of correct size
372  }
373 
375 }
376 
377 
378 template <int dim, int spacedim>
379 inline void
381  const DoFHandler<dim, spacedim> &dof,
382  const std::set<types::boundary_id> &boundary_ids,
383  const ComponentMask &component_mask)
384 {
385  // allocate an IndexSet for each global level. Contents will be
386  // overwritten inside make_boundary_list.
387  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
388  Assert(boundary_indices.empty() || boundary_indices.size() == n_levels,
389  ExcInternalError());
390  boundary_indices.resize(n_levels);
391 
393  boundary_ids,
395  component_mask);
396 }
397 
398 
399 
400 template <int dim, int spacedim>
401 inline void
403  const unsigned int level,
404  const IndexSet &level_boundary_indices)
405 {
406  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
407  if (boundary_indices.empty())
408  {
409  boundary_indices.resize(n_levels);
410  for (unsigned int i = 0; i < n_levels; ++i)
411  boundary_indices[i] = IndexSet(dof.n_dofs(i));
412  }
413  AssertDimension(boundary_indices.size(), n_levels);
414  boundary_indices[level].add_indices(level_boundary_indices);
415 }
416 
417 
418 
419 template <int dim, int spacedim>
420 inline void
422  const DoFHandler<dim, spacedim> &dof,
423  const types::boundary_id bid,
424  const unsigned int first_vector_component)
425 {
426  // For a given boundary id, find which vector component is on the boundary
427  // and set a zero boundary constraint for those degrees of freedom.
428  const unsigned int n_components = dof.get_fe_collection().n_components();
429  AssertIndexRange(first_vector_component + dim - 1, n_components);
430 
431  ComponentMask comp_mask(n_components, false);
432 
433 
435  face = dof.get_triangulation().begin_face(),
436  endf = dof.get_triangulation().end_face();
437  for (; face != endf; ++face)
438  if (face->at_boundary() && face->boundary_id() == bid)
439  for (unsigned int d = 0; d < dim; ++d)
440  {
441  Tensor<1, dim, double> unit_vec;
442  unit_vec[d] = 1.0;
443 
444  const Tensor<1, dim> normal_vec =
445  face->get_manifold().normal_vector(face, face->center());
446 
447  if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
448  comp_mask.set(d + first_vector_component, true);
449  else
450  Assert(
451  std::abs(unit_vec * normal_vec) < 1e-10,
452  ExcMessage(
453  "We can currently only support no normal flux conditions "
454  "for a specific boundary id if all faces are normal to the "
455  "x, y, or z axis."));
456  }
457 
458  Assert(comp_mask.n_selected_components() == 1,
459  ExcMessage(
460  "We can currently only support no normal flux conditions "
461  "for a specific boundary id if all faces are facing in the "
462  "same direction, i.e., a boundary normal to the x-axis must "
463  "have a different boundary id than a boundary normal to the "
464  "y- or z-axis and so on. If the mesh here was produced using "
465  "GridGenerator::..., setting colorize=true during mesh generation "
466  "and calling make_no_normal_flux_constraints() for each no normal "
467  "flux boundary will fulfill the condition."));
468 
469  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
470 }
471 
472 
473 inline void
475  const unsigned int level,
476  const AffineConstraints<double> &constraints_on_level)
477 {
479 
480  // Get the relevant DoFs from level_constraints if
481  // the user constraint matrix has not been initialized
482  if (user_constraints[level].get_local_lines().size() == 0)
483  user_constraints[level].reinit(
484  level_constraints[level].get_locally_owned_indices(),
485  level_constraints[level].get_local_lines());
486 
487  user_constraints[level].merge(
488  constraints_on_level,
490  user_constraints[level].close();
491 }
492 
493 
494 
495 inline void
497 {
498  for (auto &constraint : user_constraints)
499  constraint.clear();
500 }
501 
502 
503 
504 inline void
506 {
507  boundary_indices.clear();
508  refinement_edge_indices.clear();
509  user_constraints.clear();
510 }
511 
512 
513 inline bool
515  const types::global_dof_index index) const
516 {
517  if (boundary_indices.empty())
518  return false;
519 
521  return boundary_indices[level].is_element(index);
522 }
523 
524 inline bool
526  const types::global_dof_index index) const
527 {
529 
530  return refinement_edge_indices[level].is_element(index);
531 }
532 
533 inline bool
535  const unsigned int level,
536  const types::global_dof_index i,
537  const types::global_dof_index j) const
538 {
539  const IndexSet &interface_dofs_on_level =
540  this->get_refinement_edge_indices(level);
541 
542  return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
543  && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
544  && !this->is_boundary_index(level, i) // !on_boundary(i)
545  && !this->is_boundary_index(level, j); // !on_boundary(j)
546 }
547 
548 
549 
550 inline const IndexSet &
552 {
554  return boundary_indices[level];
555 }
556 
557 
558 
559 inline const IndexSet &
561 {
564 }
565 
566 
567 
568 inline bool
570 {
571  return boundary_indices.size() != 0;
572 }
573 
574 
575 
576 inline const AffineConstraints<double> &
578 {
580  return level_constraints[level];
581 }
582 
583 
584 
585 inline const AffineConstraints<double> &
587 {
589  return user_constraints[level];
590 }
591 
592 
593 template <typename Number>
594 inline void
596  const unsigned int level,
597  const bool add_boundary_indices,
598  const bool add_refinement_edge_indices,
599  const bool add_level_constraints,
600  const bool add_user_constraints) const
601 {
602  constraints.clear();
603 
604  // determine local lines
605  IndexSet index_set(this->get_refinement_edge_indices(level).size());
606 
608  index_set.add_indices(this->get_boundary_indices(level));
609 
610  if (add_refinement_edge_indices)
611  index_set.add_indices(this->get_refinement_edge_indices(level));
612 
613  if (add_level_constraints)
614  index_set.add_indices(this->get_level_constraints(level).get_local_lines());
615 
617  index_set.add_indices(
618  this->get_user_constraint_matrix(level).get_local_lines());
619 
620  constraints.reinit(level_constraints[level].get_locally_owned_indices(),
621  index_set);
622 
623  // merge constraints
625  constraints.add_lines(this->get_boundary_indices(level));
626 
627  if (add_refinement_edge_indices)
628  constraints.add_lines(this->get_refinement_edge_indices(level));
629 
630  if (add_level_constraints)
631  constraints.merge(this->get_level_constraints(level),
633  true);
634 
636  constraints.merge(this->get_user_constraint_matrix(level),
638  true);
639 
640  // finalize setup
641  constraints.close();
642 }
643 
644 
645 
647 
648 #endif
void merge(const AffineConstraints< other_number > &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
void add_lines(const std::vector< bool > &lines)
void set(const unsigned int index, const bool value)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
cell_iterator end() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
cell_iterator begin(const unsigned int level=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
bool is_element(const size_type index) const
Definition: index_set.h:1814
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1751
std::vector< IndexSet > refinement_edge_indices
std::vector< AffineConstraints< double > > user_constraints
void make_no_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id bid, const unsigned int first_vector_component)
bool is_interface_matrix_entry(const unsigned int level, const types::global_dof_index i, const types::global_dof_index j) const
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const std::set< types::boundary_id > &boundary_ids, const ComponentMask &component_mask={})
std::vector< IndexSet > boundary_indices
void add_boundary_indices(const DoFHandler< dim, spacedim > &dof, const unsigned int level, const IndexSet &boundary_indices)
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
void initialize(const DoFHandler< dim, spacedim > &dof, const MGLevelObject< IndexSet > &level_relevant_dofs=MGLevelObject< IndexSet >())
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
void merge_constraints(AffineConstraints< Number > &constraints, const unsigned int level, const bool add_boundary_indices, const bool add_refinement_edge_indices, const bool add_level_constraints, const bool add_user_constraints) const
std::vector< AffineConstraints< double > > level_constraints
bool have_boundary_indices() const
void add_user_constraints(const unsigned int level, const AffineConstraints< double > &constraints_on_level)
const IndexSet & get_refinement_edge_indices(unsigned int level) const
const AffineConstraints< double > & get_user_constraint_matrix(const unsigned int level) const
std::vector< std::set< types::global_dof_index > >::size_type size_dof
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
const IndexSet & get_boundary_indices(const unsigned int level) const
unsigned int max_level() const
unsigned int min_level() const
face_iterator end_face() const
virtual unsigned int n_global_levels() const
face_iterator begin_face() const
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int level
Definition: grid_out.cc:4617
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcMessage(std::string arg1)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition: dof_tools.cc:1184
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1444
void make_boundary_list(const DoFHandler< dim, spacedim > &mg_dof, const std::map< types::boundary_id, const Function< spacedim > * > &function_map, std::vector< std::set< types::global_dof_index >> &boundary_indices, const ComponentMask &component_mask={})
Definition: mg_tools.cc:1240
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const types::subdomain_id artificial_subdomain_id
Definition: types.h:315
unsigned int global_dof_index
Definition: types.h:82
unsigned int boundary_id
Definition: types.h:141