Reference documentation for deal.II version 9.3.0
\(\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 - 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 
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 
112  void
113  add_user_constraints(const unsigned int level,
114  const AffineConstraints<double> &constraints_on_level);
115 
130  template <int dim, int spacedim>
131  void
133  const types::boundary_id bid,
134  const unsigned int first_vector_component);
135 
139  void
140  clear();
141 
145  bool
146  is_boundary_index(const unsigned int level,
147  const types::global_dof_index index) const;
148 
152  bool
153  at_refinement_edge(const unsigned int level,
154  const types::global_dof_index index) const;
155 
156 
163  bool
164  is_interface_matrix_entry(const unsigned int level,
165  const types::global_dof_index i,
166  const types::global_dof_index j) const;
167 
174  const IndexSet &
175  get_boundary_indices(const unsigned int level) const;
176 
177 
182  const IndexSet &
183  get_refinement_edge_indices(unsigned int level) const;
184 
185 
189  bool
190  have_boundary_indices() const;
191 
197  get_level_constraints(const unsigned int level) const;
198 
207  get_user_constraint_matrix(const unsigned int level) const;
208 
209 private:
213  std::vector<IndexSet> boundary_indices;
214 
219  std::vector<IndexSet> refinement_edge_indices;
220 
225  std::vector<AffineConstraints<double>> level_constraints;
226 
230  std::vector<AffineConstraints<double>> user_constraints;
231 };
232 
233 
234 template <int dim, int spacedim>
235 inline void
237  const DoFHandler<dim, spacedim> &dof,
238  const MGLevelObject<IndexSet> & level_relevant_dofs)
239 {
240  boundary_indices.clear();
241  refinement_edge_indices.clear();
242  level_constraints.clear();
243  user_constraints.clear();
244 
245  const unsigned int nlevels = dof.get_triangulation().n_global_levels();
246  const unsigned int min_level = level_relevant_dofs.min_level();
247  const unsigned int max_level = (level_relevant_dofs.max_level() == 0) ?
248  nlevels - 1 :
249  level_relevant_dofs.max_level();
250  const bool user_level_dofs =
251  (level_relevant_dofs.max_level() == 0) ? false : true;
252 
253  // At this point level_constraint and refinement_edge_indices are empty.
254  refinement_edge_indices.resize(nlevels);
255  level_constraints.resize(nlevels);
256  user_constraints.resize(nlevels);
257  for (unsigned int l = min_level; l <= max_level; ++l)
258  {
259  if (user_level_dofs)
260  {
261  level_constraints[l].reinit(level_relevant_dofs[l]);
262  }
263  else
264  {
265  IndexSet relevant_dofs;
267  level_constraints[l].reinit(relevant_dofs);
268  }
269 
270  // Loop through relevant cells and faces finding those which are periodic
271  // neighbors.
272  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
273  endc = dof.end(l);
274  for (; cell != endc; ++cell)
275  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
276  {
277  for (auto f : cell->face_indices())
278  if (cell->has_periodic_neighbor(f) &&
279  cell->periodic_neighbor(f)->level() == cell->level())
280  {
281  if (cell->is_locally_owned_on_level())
282  {
283  Assert(
284  cell->periodic_neighbor(f)->level_subdomain_id() !=
286  ExcMessage(
287  "Periodic neighbor of a locally owned cell must either be owned or ghost."));
288  }
289  // Cell is a level-ghost and its neighbor is a
290  // level-artificial cell nothing to do here
291  else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
293  {
294  Assert(cell->is_locally_owned_on_level() == false,
295  ExcInternalError());
296  continue;
297  }
298 
299  const unsigned int dofs_per_face =
300  dof.get_fe(0).n_dofs_per_face(f);
301  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
302  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
303 
304  cell->periodic_neighbor(f)
305  ->face(cell->periodic_neighbor_face_no(f))
306  ->get_mg_dof_indices(l, dofs_1, 0);
307  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
308  // Store periodicity information in the level
309  // AffineConstraints object. Skip DoFs for which we've
310  // previously entered periodicity constraints already; this
311  // can happen, for example, for a vertex dof at a periodic
312  // boundary that we visit from more than one cell
313  for (unsigned int i = 0; i < dofs_per_face; ++i)
314  if (level_constraints[l].can_store_line(dofs_2[i]) &&
315  level_constraints[l].can_store_line(dofs_1[i]) &&
316  !level_constraints[l].is_constrained(dofs_2[i]) &&
317  !level_constraints[l].is_constrained(dofs_1[i]))
318  {
319  level_constraints[l].add_line(dofs_2[i]);
320  level_constraints[l].add_entry(dofs_2[i],
321  dofs_1[i],
322  1.);
323  }
324  }
325  }
326  level_constraints[l].close();
327 
328  // Initialize with empty IndexSet of correct size
330  }
331 
333 }
334 
335 
336 template <int dim, int spacedim>
337 inline void
339  const DoFHandler<dim, spacedim> & dof,
340  const std::set<types::boundary_id> &boundary_ids,
341  const ComponentMask & component_mask)
342 {
343  // allocate an IndexSet for each global level. Contents will be
344  // overwritten inside make_boundary_list.
345  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
346  Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
347  ExcInternalError());
348  boundary_indices.resize(n_levels);
349 
351  boundary_ids,
353  component_mask);
354 }
355 
356 
357 
358 template <int dim, int spacedim>
359 inline void
361  const DoFHandler<dim, spacedim> &dof,
362  const types::boundary_id bid,
363  const unsigned int first_vector_component)
364 {
365  // For a given boundary id, find which vector component is on the boundary
366  // and set a zero boundary constraint for those degrees of freedom.
367  const unsigned int n_components = dof.get_fe_collection().n_components();
368  AssertIndexRange(first_vector_component + dim - 1, n_components);
369 
370  ComponentMask comp_mask(n_components, false);
371 
372 
374  face = dof.get_triangulation().begin_face(),
375  endf = dof.get_triangulation().end_face();
376  for (; face != endf; ++face)
377  if (face->at_boundary() && face->boundary_id() == bid)
378  for (unsigned int d = 0; d < dim; ++d)
379  {
380  Tensor<1, dim, double> unit_vec;
381  unit_vec[d] = 1.0;
382 
383  const Tensor<1, dim> normal_vec =
384  face->get_manifold().normal_vector(face, face->center());
385 
386  if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
387  comp_mask.set(d + first_vector_component, true);
388  else
389  Assert(
390  std::abs(unit_vec * normal_vec) < 1e-10,
391  ExcMessage(
392  "We can currently only support no normal flux conditions "
393  "for a specific boundary id if all faces are normal to the "
394  "x, y, or z axis."));
395  }
396 
397  Assert(comp_mask.n_selected_components() == 1,
398  ExcMessage(
399  "We can currently only support no normal flux conditions "
400  "for a specific boundary id if all faces are facing in the "
401  "same direction, i.e., a boundary normal to the x-axis must "
402  "have a different boundary id than a boundary normal to the "
403  "y- or z-axis and so on. If the mesh here was produced using "
404  "GridGenerator::..., setting colorize=true during mesh generation "
405  "and calling make_no_normal_flux_constraints() for each no normal "
406  "flux boundary will fulfill the condition."));
407 
408  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
409 }
410 
411 
412 inline void
414  const unsigned int level,
415  const AffineConstraints<double> &constraints_on_level)
416 {
417  AssertIndexRange(level, user_constraints.size());
418 
419  // Get the relevant DoFs from level_constraints if
420  // the user constraint matrix has not been initialized
421  if (user_constraints[level].get_local_lines().size() == 0)
422  user_constraints[level].reinit(level_constraints[level].get_local_lines());
423 
424  user_constraints[level].merge(
425  constraints_on_level,
427  user_constraints[level].close();
428 }
429 
430 
431 inline void
433 {
434  boundary_indices.clear();
435  refinement_edge_indices.clear();
436  user_constraints.clear();
437 }
438 
439 
440 inline bool
442  const types::global_dof_index index) const
443 {
444  if (boundary_indices.size() == 0)
445  return false;
446 
447  AssertIndexRange(level, boundary_indices.size());
448  return boundary_indices[level].is_element(index);
449 }
450 
451 inline bool
453  const types::global_dof_index index) const
454 {
456 
457  return refinement_edge_indices[level].is_element(index);
458 }
459 
460 inline bool
462  const unsigned int level,
463  const types::global_dof_index i,
464  const types::global_dof_index j) const
465 {
466  const IndexSet &interface_dofs_on_level =
467  this->get_refinement_edge_indices(level);
468 
469  return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
470  && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
471  && !this->is_boundary_index(level, i) // !on_boundary(i)
472  && !this->is_boundary_index(level, j); // !on_boundary(j)
473 }
474 
475 
476 
477 inline const IndexSet &
479 {
480  AssertIndexRange(level, boundary_indices.size());
481  return boundary_indices[level];
482 }
483 
484 
485 
486 inline const IndexSet &
488 {
491 }
492 
493 
494 
495 inline bool
497 {
498  return boundary_indices.size() != 0;
499 }
500 
501 
502 
503 inline const AffineConstraints<double> &
505 {
506  AssertIndexRange(level, level_constraints.size());
507  return level_constraints[level];
508 }
509 
510 
511 
512 inline const AffineConstraints<double> &
514 {
515  AssertIndexRange(level, user_constraints.size());
516  return user_constraints[level];
517 }
518 
519 
520 
522 
523 #endif
unsigned int max_level() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
cell_iterator begin(const unsigned int level=0) const
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
void make_no_normal_flux_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id bid, const unsigned int first_vector_component)
cell_iterator end() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
std::vector< IndexSet > boundary_indices
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
std::vector< std::set< types::global_dof_index > >::size_type size_dof
void set(const unsigned int index, const bool value)
unsigned int min_level() const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1252
bool is_interface_matrix_entry(const unsigned int level, const types::global_dof_index i, const types::global_dof_index j) const
static ::ExceptionBase & ExcMessage(std::string arg1)
const IndexSet & get_boundary_indices(const unsigned int level) const
void add_user_constraints(const unsigned int level, const AffineConstraints< double > &constraints_on_level)
#define Assert(cond, exc)
Definition: exceptions.h:1465
types::global_dof_index n_dofs() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
unsigned int level
Definition: grid_out.cc:4590
const IndexSet & get_refinement_edge_indices(unsigned int level) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const std::set< types::boundary_id > &boundary_ids, const ComponentMask &component_mask=ComponentMask())
void initialize(const DoFHandler< dim, spacedim > &dof, const MGLevelObject< IndexSet > &level_relevant_dofs=MGLevelObject< IndexSet >())
std::vector< IndexSet > refinement_edge_indices
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
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:1270
Definition: tensor.h:462
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1480
bool is_element(const size_type index) const
Definition: index_set.h:1765
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:239
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
bool have_boundary_indices() const
std::vector< AffineConstraints< double > > level_constraints
const AffineConstraints< double > & get_user_constraint_matrix(const unsigned int level) const
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:10259
std::vector< AffineConstraints< double > > user_constraints
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()