Reference documentation for deal.II version GIT a189bc2bdf 2022-12-07 02:45: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 - 2022 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 class DoFHandler;
37 #endif
38 
39 
47 {
48 public:
49  using size_dof = std::vector<std::set<types::global_dof_index>>::size_type;
75  template <int dim, int spacedim>
76  void
78  const MGLevelObject<IndexSet> & level_relevant_dofs =
80 
91  template <int dim, int spacedim>
92  void
94  const DoFHandler<dim, spacedim> & dof,
95  const std::set<types::boundary_id> &boundary_ids,
96  const ComponentMask & component_mask = ComponentMask());
97 
104  template <int dim, int spacedim>
105  void
107  const unsigned int level,
108  const IndexSet & boundary_indices);
109 
124  void
125  add_user_constraints(const unsigned int level,
126  const AffineConstraints<double> &constraints_on_level);
127 
142  template <int dim, int spacedim>
143  void
145  const types::boundary_id bid,
146  const unsigned int first_vector_component);
147 
151  void
153 
157  void
158  clear();
159 
163  bool
164  is_boundary_index(const unsigned int level,
165  const types::global_dof_index index) const;
166 
170  bool
171  at_refinement_edge(const unsigned int level,
172  const types::global_dof_index index) const;
173 
174 
181  bool
182  is_interface_matrix_entry(const unsigned int level,
183  const types::global_dof_index i,
184  const types::global_dof_index j) const;
185 
192  const IndexSet &
193  get_boundary_indices(const unsigned int level) const;
194 
195 
200  const IndexSet &
201  get_refinement_edge_indices(unsigned int level) const;
202 
203 
207  bool
208  have_boundary_indices() const;
209 
215  get_level_constraints(const unsigned int level) const;
216 
225  get_user_constraint_matrix(const unsigned int level) const;
226 
227 private:
231  std::vector<IndexSet> boundary_indices;
232 
237  std::vector<IndexSet> refinement_edge_indices;
238 
243  std::vector<AffineConstraints<double>> level_constraints;
244 
248  std::vector<AffineConstraints<double>> user_constraints;
249 };
250 
251 
252 template <int dim, int spacedim>
253 inline void
255  const DoFHandler<dim, spacedim> &dof,
256  const MGLevelObject<IndexSet> & level_relevant_dofs)
257 {
258  boundary_indices.clear();
259  refinement_edge_indices.clear();
260  level_constraints.clear();
261  user_constraints.clear();
262 
263  const unsigned int nlevels = dof.get_triangulation().n_global_levels();
264  const unsigned int min_level = level_relevant_dofs.min_level();
265  const unsigned int max_level = (level_relevant_dofs.max_level() == 0) ?
266  nlevels - 1 :
267  level_relevant_dofs.max_level();
268  const bool user_level_dofs =
269  (level_relevant_dofs.max_level() == 0) ? false : true;
270 
271  // At this point level_constraint and refinement_edge_indices are empty.
272  refinement_edge_indices.resize(nlevels);
273  level_constraints.resize(nlevels);
274  user_constraints.resize(nlevels);
275  for (unsigned int l = min_level; l <= max_level; ++l)
276  {
277  if (user_level_dofs)
278  {
279  level_constraints[l].reinit(level_relevant_dofs[l]);
280  }
281  else
282  {
283  const IndexSet relevant_dofs =
285  level_constraints[l].reinit(relevant_dofs);
286  }
287 
288  // Loop through relevant cells and faces finding those which are periodic
289  // neighbors.
290  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
291  endc = dof.end(l);
292  for (; cell != endc; ++cell)
293  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
294  {
295  for (auto f : cell->face_indices())
296  if (cell->has_periodic_neighbor(f) &&
297  cell->periodic_neighbor(f)->level() == cell->level())
298  {
299  if (cell->is_locally_owned_on_level())
300  {
301  Assert(
302  cell->periodic_neighbor(f)->level_subdomain_id() !=
304  ExcMessage(
305  "Periodic neighbor of a locally owned cell must either be owned or ghost."));
306  }
307  // Cell is a level-ghost and its neighbor is a
308  // level-artificial cell nothing to do here
309  else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
311  {
312  Assert(cell->is_locally_owned_on_level() == false,
313  ExcInternalError());
314  continue;
315  }
316 
317  const unsigned int dofs_per_face =
318  dof.get_fe(0).n_dofs_per_face(f);
319  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
320  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
321 
322  cell->periodic_neighbor(f)
323  ->face(cell->periodic_neighbor_face_no(f))
324  ->get_mg_dof_indices(l, dofs_1, 0);
325  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
326  // Store periodicity information in the level
327  // AffineConstraints object. Skip DoFs for which we've
328  // previously entered periodicity constraints already; this
329  // can happen, for example, for a vertex dof at a periodic
330  // boundary that we visit from more than one cell
331  for (unsigned int i = 0; i < dofs_per_face; ++i)
332  if (level_constraints[l].can_store_line(dofs_2[i]) &&
333  level_constraints[l].can_store_line(dofs_1[i]) &&
334  !level_constraints[l].is_constrained(dofs_2[i]) &&
335  !level_constraints[l].is_constrained(dofs_1[i]))
336  {
337  level_constraints[l].add_line(dofs_2[i]);
338  level_constraints[l].add_entry(dofs_2[i],
339  dofs_1[i],
340  1.);
341  }
342  }
343  }
344  level_constraints[l].close();
345 
346  // Initialize with empty IndexSet of correct size
348  }
349 
351 }
352 
353 
354 template <int dim, int spacedim>
355 inline void
357  const DoFHandler<dim, spacedim> & dof,
358  const std::set<types::boundary_id> &boundary_ids,
359  const ComponentMask & component_mask)
360 {
361  // allocate an IndexSet for each global level. Contents will be
362  // overwritten inside make_boundary_list.
363  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
364  Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
365  ExcInternalError());
366  boundary_indices.resize(n_levels);
367 
369  boundary_ids,
371  component_mask);
372 }
373 
374 
375 
376 template <int dim, int spacedim>
377 inline void
379  const unsigned int level,
380  const IndexSet &level_boundary_indices)
381 {
382  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
383  if (boundary_indices.size() == 0)
384  {
385  boundary_indices.resize(n_levels);
386  for (unsigned int i = 0; i < n_levels; ++i)
387  boundary_indices[i] = IndexSet(dof.n_dofs(i));
388  }
389  AssertDimension(boundary_indices.size(), n_levels);
390  boundary_indices[level].add_indices(level_boundary_indices);
391 }
392 
393 
394 
395 template <int dim, int spacedim>
396 inline void
398  const DoFHandler<dim, spacedim> &dof,
399  const types::boundary_id bid,
400  const unsigned int first_vector_component)
401 {
402  // For a given boundary id, find which vector component is on the boundary
403  // and set a zero boundary constraint for those degrees of freedom.
404  const unsigned int n_components = dof.get_fe_collection().n_components();
405  AssertIndexRange(first_vector_component + dim - 1, n_components);
406 
407  ComponentMask comp_mask(n_components, false);
408 
409 
411  face = dof.get_triangulation().begin_face(),
412  endf = dof.get_triangulation().end_face();
413  for (; face != endf; ++face)
414  if (face->at_boundary() && face->boundary_id() == bid)
415  for (unsigned int d = 0; d < dim; ++d)
416  {
417  Tensor<1, dim, double> unit_vec;
418  unit_vec[d] = 1.0;
419 
420  const Tensor<1, dim> normal_vec =
421  face->get_manifold().normal_vector(face, face->center());
422 
423  if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
424  comp_mask.set(d + first_vector_component, true);
425  else
426  Assert(
427  std::abs(unit_vec * normal_vec) < 1e-10,
428  ExcMessage(
429  "We can currently only support no normal flux conditions "
430  "for a specific boundary id if all faces are normal to the "
431  "x, y, or z axis."));
432  }
433 
434  Assert(comp_mask.n_selected_components() == 1,
435  ExcMessage(
436  "We can currently only support no normal flux conditions "
437  "for a specific boundary id if all faces are facing in the "
438  "same direction, i.e., a boundary normal to the x-axis must "
439  "have a different boundary id than a boundary normal to the "
440  "y- or z-axis and so on. If the mesh here was produced using "
441  "GridGenerator::..., setting colorize=true during mesh generation "
442  "and calling make_no_normal_flux_constraints() for each no normal "
443  "flux boundary will fulfill the condition."));
444 
445  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
446 }
447 
448 
449 inline void
451  const unsigned int level,
452  const AffineConstraints<double> &constraints_on_level)
453 {
455 
456  // Get the relevant DoFs from level_constraints if
457  // the user constraint matrix has not been initialized
458  if (user_constraints[level].get_local_lines().size() == 0)
459  user_constraints[level].reinit(level_constraints[level].get_local_lines());
460 
461  user_constraints[level].merge(
462  constraints_on_level,
464  user_constraints[level].close();
465 }
466 
467 
468 
469 inline void
471 {
472  for (auto &constraint : user_constraints)
473  constraint.clear();
474 }
475 
476 
477 
478 inline void
480 {
481  boundary_indices.clear();
482  refinement_edge_indices.clear();
483  user_constraints.clear();
484 }
485 
486 
487 inline bool
489  const types::global_dof_index index) const
490 {
491  if (boundary_indices.size() == 0)
492  return false;
493 
495  return boundary_indices[level].is_element(index);
496 }
497 
498 inline bool
500  const types::global_dof_index index) const
501 {
503 
504  return refinement_edge_indices[level].is_element(index);
505 }
506 
507 inline bool
509  const unsigned int level,
510  const types::global_dof_index i,
511  const types::global_dof_index j) const
512 {
513  const IndexSet &interface_dofs_on_level =
514  this->get_refinement_edge_indices(level);
515 
516  return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
517  && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
518  && !this->is_boundary_index(level, i) // !on_boundary(i)
519  && !this->is_boundary_index(level, j); // !on_boundary(j)
520 }
521 
522 
523 
524 inline const IndexSet &
526 {
528  return boundary_indices[level];
529 }
530 
531 
532 
533 inline const IndexSet &
535 {
538 }
539 
540 
541 
542 inline bool
544 {
545  return boundary_indices.size() != 0;
546 }
547 
548 
549 
550 inline const AffineConstraints<double> &
552 {
554  return level_constraints[level];
555 }
556 
557 
558 
559 inline const AffineConstraints<double> &
561 {
563  return user_constraints[level];
564 }
565 
566 
567 
569 
570 #endif
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 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:1734
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
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
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
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const std::set< types::boundary_id > &boundary_ids, const ComponentMask &component_mask=ComponentMask())
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
unsigned int level
Definition: grid_out.cc:4608
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
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:1190
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:1446
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=ComponentMask())
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:308
unsigned int global_dof_index
Definition: types.h:81
unsigned int boundary_id
Definition: types.h:134