Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 
22 
24 
26 
27 #include <set>
28 #include <vector>
29 
31 
32 // Forward declaration
33 #ifndef DOXYGEN
34 template <int dim, int spacedim>
35 class DoFHandler;
36 #endif
37 
38 
46 {
47 public:
48  using size_dof = std::vector<std::set<types::global_dof_index>>::size_type;
68  template <int dim, int spacedim>
69  void
71 
82  template <int dim, int spacedim>
85  const std::map<types::boundary_id, const Function<spacedim> *>
86  & function_map,
87  const ComponentMask &component_mask = ComponentMask());
88 
99  template <int dim, int spacedim>
100  void
102  const DoFHandler<dim, spacedim> & dof,
103  const std::set<types::boundary_id> &boundary_ids,
104  const ComponentMask & component_mask = ComponentMask());
105 
120  void
121  add_user_constraints(const unsigned int level,
122  const AffineConstraints<double> &constraints_on_level);
123 
138  template <int dim, int spacedim>
139  void
141  const types::boundary_id bid,
142  const unsigned int first_vector_component);
143 
147  void
148  clear();
149 
153  bool
154  is_boundary_index(const unsigned int level,
155  const types::global_dof_index index) const;
156 
160  bool
161  at_refinement_edge(const unsigned int level,
162  const types::global_dof_index index) const;
163 
164 
171  bool
172  is_interface_matrix_entry(const unsigned int level,
173  const types::global_dof_index i,
174  const types::global_dof_index j) const;
175 
182  const IndexSet &
183  get_boundary_indices(const unsigned int level) const;
184 
185 
190  const IndexSet &
191  get_refinement_edge_indices(unsigned int level) const;
192 
193 
197  bool
198  have_boundary_indices() const;
199 
205  get_level_constraints(const unsigned int level) const;
206 
215  get_level_constraint_matrix(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 {
256  boundary_indices.clear();
257  refinement_edge_indices.clear();
258  level_constraints.clear();
259  user_constraints.clear();
260 
261  const unsigned int nlevels = dof.get_triangulation().n_global_levels();
262 
263  // At this point level_constraint and refinement_edge_indices are empty.
264  refinement_edge_indices.resize(nlevels);
265  level_constraints.resize(nlevels);
266  user_constraints.resize(nlevels);
267  for (unsigned int l = 0; l < nlevels; ++l)
268  {
269  IndexSet relevant_dofs;
271  level_constraints[l].reinit(relevant_dofs);
272 
273  // Loop through relevant cells and faces finding those which are periodic
274  // neighbors.
275  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
276  endc = dof.end(l);
277  for (; cell != endc; ++cell)
278  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
279  {
280  for (auto f : GeometryInfo<dim>::face_indices())
281  if (cell->has_periodic_neighbor(f) &&
282  cell->periodic_neighbor(f)->level() == cell->level())
283  {
284  if (cell->is_locally_owned_on_level())
285  {
286  Assert(
287  cell->periodic_neighbor(f)->level_subdomain_id() !=
289  ExcMessage(
290  "Periodic neighbor of a locally owned cell must either be owned or ghost."));
291  }
292  // Cell is a level-ghost and its neighbor is a
293  // level-artificial cell nothing to do here
294  else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
296  {
297  Assert(cell->is_locally_owned_on_level() == false,
298  ExcInternalError());
299  continue;
300  }
301 
302  const unsigned int dofs_per_face =
303  cell->face(f)->get_fe(0).dofs_per_face;
304  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
305  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
306 
307  cell->periodic_neighbor(f)
308  ->face(cell->periodic_neighbor_face_no(f))
309  ->get_mg_dof_indices(l, dofs_1, 0);
310  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
311  // Store periodicity information in the level
312  // AffineConstraints object. Skip DoFs for which we've
313  // previously entered periodicity constraints already; this
314  // can happen, for example, for a vertex dof at a periodic
315  // boundary that we visit from more than one cell
316  for (unsigned int i = 0; i < dofs_per_face; ++i)
317  if (level_constraints[l].can_store_line(dofs_2[i]) &&
318  level_constraints[l].can_store_line(dofs_1[i]) &&
319  !level_constraints[l].is_constrained(dofs_2[i]) &&
320  !level_constraints[l].is_constrained(dofs_1[i]))
321  {
322  level_constraints[l].add_line(dofs_2[i]);
323  level_constraints[l].add_entry(dofs_2[i],
324  dofs_1[i],
325  1.);
326  }
327  }
328  }
329  level_constraints[l].close();
330 
331  // Initialize with empty IndexSet of correct size
333  }
334 
336 }
337 
338 
339 template <int dim, int spacedim>
340 inline void
342  const DoFHandler<dim, spacedim> & dof,
343  const std::map<types::boundary_id, const Function<spacedim> *> &function_map,
344  const ComponentMask &component_mask)
345 {
346  initialize(dof);
347 
348  // allocate an IndexSet for each global level. Contents will be
349  // overwritten inside make_boundary_list.
350  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
351  // At this point boundary_indices is empty.
352  boundary_indices.resize(n_levels);
353 
355  function_map,
357  component_mask);
358 }
359 
360 
361 template <int dim, int spacedim>
362 inline void
364  const DoFHandler<dim, spacedim> & dof,
365  const std::set<types::boundary_id> &boundary_ids,
366  const ComponentMask & component_mask)
367 {
368  // allocate an IndexSet for each global level. Contents will be
369  // overwritten inside make_boundary_list.
370  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
371  Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
372  ExcInternalError());
373  boundary_indices.resize(n_levels);
374 
376  boundary_ids,
378  component_mask);
379 }
380 
381 
382 template <int dim, int spacedim>
383 inline void
385  const DoFHandler<dim, spacedim> &dof,
386  const types::boundary_id bid,
387  const unsigned int first_vector_component)
388 {
389  // For a given boundary id, find which vector component is on the boundary
390  // and set a zero boundary constraint for those degrees of freedom.
391  const unsigned int n_components = DoFTools::n_components(dof);
392  AssertIndexRange(first_vector_component + dim - 1, n_components);
393 
394  ComponentMask comp_mask(n_components, false);
395 
396 
398  face = dof.get_triangulation().begin_face(),
399  endf = dof.get_triangulation().end_face();
400  for (; face != endf; ++face)
401  if (face->at_boundary() && face->boundary_id() == bid)
402  for (unsigned int d = 0; d < dim; ++d)
403  {
404  Tensor<1, dim, double> unit_vec;
405  unit_vec[d] = 1.0;
406 
407  const Tensor<1, dim> normal_vec =
408  face->get_manifold().normal_vector(face, face->center());
409 
410  if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
411  comp_mask.set(d + first_vector_component, true);
412  else
413  Assert(
414  std::abs(unit_vec * normal_vec) < 1e-10,
415  ExcMessage(
416  "We can currently only support no normal flux conditions "
417  "for a specific boundary id if all faces are normal to the "
418  "x, y, or z axis."));
419  }
420 
421  Assert(comp_mask.n_selected_components() == 1,
422  ExcMessage(
423  "We can currently only support no normal flux conditions "
424  "for a specific boundary id if all faces are facing in the "
425  "same direction, i.e., a boundary normal to the x-axis must "
426  "have a different boundary id than a boundary normal to the "
427  "y- or z-axis and so on. If the mesh here was produced using "
428  "GridGenerator::..., setting colorize=true during mesh generation "
429  "and calling make_no_normal_flux_constraints() for each no normal "
430  "flux boundary will fulfill the condition."));
431 
432  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
433 }
434 
435 
436 inline void
438  const unsigned int level,
439  const AffineConstraints<double> &constraints_on_level)
440 {
441  AssertIndexRange(level, user_constraints.size());
442 
443  // Get the relevant DoFs from level_constraints if
444  // the user constraint matrix has not been initialized
445  if (user_constraints[level].get_local_lines().size() == 0)
446  user_constraints[level].reinit(level_constraints[level].get_local_lines());
447 
448  user_constraints[level].merge(
449  constraints_on_level,
451  user_constraints[level].close();
452 }
453 
454 
455 inline void
457 {
458  boundary_indices.clear();
459  refinement_edge_indices.clear();
460  user_constraints.clear();
461 }
462 
463 
464 inline bool
466  const types::global_dof_index index) const
467 {
468  if (boundary_indices.size() == 0)
469  return false;
470 
471  AssertIndexRange(level, boundary_indices.size());
472  return boundary_indices[level].is_element(index);
473 }
474 
475 inline bool
477  const types::global_dof_index index) const
478 {
480 
481  return refinement_edge_indices[level].is_element(index);
482 }
483 
484 inline bool
486  const unsigned int level,
487  const types::global_dof_index i,
488  const types::global_dof_index j) const
489 {
490  const IndexSet &interface_dofs_on_level =
491  this->get_refinement_edge_indices(level);
492 
493  return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
494  && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
495  && !this->is_boundary_index(level, i) // !on_boundary(i)
496  && !this->is_boundary_index(level, j); // !on_boundary(j)
497 }
498 
499 
500 
501 inline const IndexSet &
503 {
504  AssertIndexRange(level, boundary_indices.size());
505  return boundary_indices[level];
506 }
507 
508 
509 
510 inline const IndexSet &
512 {
515 }
516 
517 
518 
519 inline bool
521 {
522  return boundary_indices.size() != 0;
523 }
524 
525 
526 
527 inline const AffineConstraints<double> &
529 {
530  AssertIndexRange(level, level_constraints.size());
531  return level_constraints[level];
532 }
533 
534 
535 
536 inline const AffineConstraints<double> &
538 {
539  return get_level_constraints(level);
540 }
541 
542 
543 
544 inline const AffineConstraints<double> &
546 {
547  AssertIndexRange(level, user_constraints.size());
548  return user_constraints[level];
549 }
550 
551 
552 
554 
555 #endif
const Triangulation< dim, spacedim > & get_triangulation() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:918
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
Definition: dof_handler.cc:947
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
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)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
const AffineConstraints< double > & get_level_constraint_matrix(const unsigned int level) const
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:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
unsigned int level
Definition: grid_out.cc:4355
types::global_dof_index n_dofs() const
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())
std::vector< IndexSet > refinement_edge_indices
const types::subdomain_id artificial_subdomain_id
Definition: types.h:302
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:1243
Definition: tensor.h:450
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
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:1528
bool is_element(const size_type index) const
Definition: index_set.h:1766
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:237
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
bool have_boundary_indices() const
std::vector< AffineConstraints< double > > level_constraints
const AffineConstraints< double > & get_user_constraint_matrix(const unsigned int level) const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
#define DEAL_II_DEPRECATED
Definition: config.h:151
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:10321
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()
void initialize(const DoFHandler< dim, spacedim > &dof)