Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_constrained_dofs.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 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 template <int dim, int spacedim>
33 class DoFHandler;
34 
35 
43 {
44 public:
45  using size_dof = std::vector<std::set<types::global_dof_index>>::size_type;
65  template <int dim, int spacedim>
66  void
68 
79  template <int dim, int spacedim>
80  DEAL_II_DEPRECATED void
82  const std::map<types::boundary_id, const Function<spacedim> *>
83  & function_map,
84  const ComponentMask &component_mask = ComponentMask());
85 
96  template <int dim, int spacedim>
97  void
99  const DoFHandler<dim, spacedim> & dof,
100  const std::set<types::boundary_id> &boundary_ids,
101  const ComponentMask & component_mask = ComponentMask());
102 
117  template <int dim, int spacedim>
118  void
120  const types::boundary_id bid,
121  const unsigned int first_vector_component);
122 
126  void
127  clear();
128 
132  bool
133  is_boundary_index(const unsigned int level,
134  const types::global_dof_index index) const;
135 
139  bool
140  at_refinement_edge(const unsigned int level,
141  const types::global_dof_index index) const;
142 
143 
150  bool
151  is_interface_matrix_entry(const unsigned int level,
152  const types::global_dof_index i,
153  const types::global_dof_index j) const;
154 
161  const IndexSet &
162  get_boundary_indices(const unsigned int level) const;
163 
164 
169  const IndexSet &
170  get_refinement_edge_indices(unsigned int level) const;
171 
172 
176  bool
177  have_boundary_indices() const;
178 
184  get_level_constraints(const unsigned int level) const;
185 
192  DEAL_II_DEPRECATED
194  get_level_constraint_matrix(const unsigned int level) const;
195 
196 private:
200  std::vector<IndexSet> boundary_indices;
201 
206  std::vector<IndexSet> refinement_edge_indices;
207 
212  std::vector<AffineConstraints<double>> level_constraints;
213 };
214 
215 
216 template <int dim, int spacedim>
217 inline void
219 {
220  boundary_indices.clear();
221  refinement_edge_indices.clear();
222  level_constraints.clear();
223 
224  const unsigned int nlevels = dof.get_triangulation().n_global_levels();
225 
226  // At this point level_constraint and refinement_edge_indices are empty.
227  level_constraints.resize(nlevels);
228  refinement_edge_indices.resize(nlevels);
229  for (unsigned int l = 0; l < nlevels; ++l)
230  {
231  IndexSet relevant_dofs;
232  DoFTools::extract_locally_relevant_level_dofs(dof, l, relevant_dofs);
233  level_constraints[l].reinit(relevant_dofs);
234 
235  // Loop through relevant cells and faces finding those which are periodic
236  // neighbors.
237  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
238  endc = dof.end(l);
239  for (; cell != endc; ++cell)
240  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
241  {
242  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
243  if (cell->has_periodic_neighbor(f) &&
244  cell->periodic_neighbor(f)->level() == cell->level())
245  {
246  if (cell->is_locally_owned_on_level())
247  {
248  Assert(
249  cell->periodic_neighbor(f)->level_subdomain_id() !=
251  ExcMessage(
252  "Periodic neighbor of a locally owned cell must either be owned or ghost."));
253  }
254  // Cell is a level-ghost and its neighbor is a
255  // level-artificial cell nothing to do here
256  else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
258  {
259  Assert(cell->is_locally_owned_on_level() == false,
260  ExcInternalError());
261  continue;
262  }
263 
264  const unsigned int dofs_per_face =
265  cell->face(f)->get_fe(0).dofs_per_face;
266  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
267  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
268 
269  cell->periodic_neighbor(f)
270  ->face(cell->periodic_neighbor_face_no(f))
271  ->get_mg_dof_indices(l, dofs_1, 0);
272  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
273  // Store periodicity information in the level
274  // AffineConstraints object. Skip DoFs for which we've
275  // previously entered periodicity constraints already; this
276  // can happen, for example, for a vertex dof at a periodic
277  // boundary that we visit from more than one cell
278  for (unsigned int i = 0; i < dofs_per_face; ++i)
279  if (level_constraints[l].can_store_line(dofs_2[i]) &&
280  level_constraints[l].can_store_line(dofs_1[i]) &&
281  !level_constraints[l].is_constrained(dofs_2[i]) &&
282  !level_constraints[l].is_constrained(dofs_1[i]))
283  {
284  level_constraints[l].add_line(dofs_2[i]);
285  level_constraints[l].add_entry(dofs_2[i],
286  dofs_1[i],
287  1.);
288  }
289  }
290  }
291  level_constraints[l].close();
292 
293  // Initialize with empty IndexSet of correct size
295  }
296 
298 }
299 
300 
301 template <int dim, int spacedim>
302 inline void
304  const DoFHandler<dim, spacedim> & dof,
305  const std::map<types::boundary_id, const Function<spacedim> *> &function_map,
306  const ComponentMask &component_mask)
307 {
308  initialize(dof);
309 
310  // allocate an IndexSet for each global level. Contents will be
311  // overwritten inside make_boundary_list.
312  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
313  // At this point boundary_indices is empty.
314  boundary_indices.resize(n_levels);
315 
317  function_map,
319  component_mask);
320 }
321 
322 
323 template <int dim, int spacedim>
324 inline void
326  const DoFHandler<dim, spacedim> & dof,
327  const std::set<types::boundary_id> &boundary_ids,
328  const ComponentMask & component_mask)
329 {
330  // allocate an IndexSet for each global level. Contents will be
331  // overwritten inside make_boundary_list.
332  const unsigned int n_levels = dof.get_triangulation().n_global_levels();
333  Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
334  ExcInternalError());
335  boundary_indices.resize(n_levels);
336 
338  boundary_ids,
340  component_mask);
341 }
342 
343 
344 template <int dim, int spacedim>
345 inline void
347  const DoFHandler<dim, spacedim> &dof,
348  const types::boundary_id bid,
349  const unsigned int first_vector_component)
350 {
351  // For a given boundary id, find which vector component is on the boundary
352  // and set a zero boundary constraint for those degrees of freedom.
353  const unsigned int n_components = DoFTools::n_components(dof);
354  Assert(first_vector_component + dim <= n_components,
355  ExcIndexRange(first_vector_component, 0, n_components - dim + 1));
356 
357  ComponentMask comp_mask(n_components, false);
358 
359 
361  face = dof.get_triangulation().begin_face(),
362  endf = dof.get_triangulation().end_face();
363  for (; face != endf; ++face)
364  if (face->at_boundary() && face->boundary_id() == bid)
365  for (unsigned int d = 0; d < dim; ++d)
366  {
367  Tensor<1, dim, double> unit_vec;
368  unit_vec[d] = 1.0;
369 
370  const Tensor<1, dim> normal_vec =
371  face->get_manifold().normal_vector(face, face->center());
372 
373  if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
374  comp_mask.set(d + first_vector_component, true);
375  else
376  Assert(
377  std::abs(unit_vec * normal_vec) < 1e-10,
378  ExcMessage(
379  "We can currently only support no normal flux conditions "
380  "for a specific boundary id if all faces are normal to the "
381  "x, y, or z axis."));
382  }
383 
384  Assert(comp_mask.n_selected_components() == 1,
385  ExcMessage(
386  "We can currently only support no normal flux conditions "
387  "for a specific boundary id if all faces are facing in the "
388  "same direction, i.e., a boundary normal to the x-axis must "
389  "have a different boundary id than a boundary normal to the "
390  "y- or z-axis and so on. If the mesh here was produced using "
391  "GridGenerator::..., setting colorize=true during mesh generation "
392  "and calling make_no_normal_flux_constraints() for each no normal "
393  "flux boundary will fulfill the condition."));
394 
395  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
396 }
397 
398 
399 inline void
401 {
402  boundary_indices.clear();
403  refinement_edge_indices.clear();
404 }
405 
406 
407 inline bool
408 MGConstrainedDoFs::is_boundary_index(const unsigned int level,
409  const types::global_dof_index index) const
410 {
411  if (boundary_indices.size() == 0)
412  return false;
413 
414  AssertIndexRange(level, boundary_indices.size());
415  return boundary_indices[level].is_element(index);
416 }
417 
418 inline bool
419 MGConstrainedDoFs::at_refinement_edge(const unsigned int level,
420  const types::global_dof_index index) const
421 {
423 
424  return refinement_edge_indices[level].is_element(index);
425 }
426 
427 inline bool
429  const unsigned int level,
430  const types::global_dof_index i,
431  const types::global_dof_index j) const
432 {
433  const IndexSet &interface_dofs_on_level =
434  this->get_refinement_edge_indices(level);
435 
436  return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
437  && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
438  && !this->is_boundary_index(level, i) // !on_boundary(i)
439  && !this->is_boundary_index(level, j); // !on_boundary(j)
440 }
441 
442 
443 
444 inline const IndexSet &
445 MGConstrainedDoFs::get_boundary_indices(const unsigned int level) const
446 {
447  AssertIndexRange(level, boundary_indices.size());
448  return boundary_indices[level];
449 }
450 
451 
452 
453 inline const IndexSet &
455 {
457  return refinement_edge_indices[level];
458 }
459 
460 
461 
462 inline bool
464 {
465  return boundary_indices.size() != 0;
466 }
467 
468 
469 
470 inline const AffineConstraints<double> &
471 MGConstrainedDoFs::get_level_constraints(const unsigned int level) const
472 {
473  AssertIndexRange(level, level_constraints.size());
474  return level_constraints[level];
475 }
476 
477 
478 
479 inline const AffineConstraints<double> &
480 MGConstrainedDoFs::get_level_constraint_matrix(const unsigned int level) const
481 {
482  return get_level_constraints(level);
483 }
484 
485 
486 
487 DEAL_II_NAMESPACE_CLOSE
488 
489 #endif
const Triangulation< dim, spacedim > & get_triangulation() const
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:930
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:959
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
std::vector< IndexSet > boundary_indices
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
void set(const unsigned int index, const bool value)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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:1179
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
#define Assert(cond, exc)
Definition: exceptions.h:1407
types::global_dof_index n_dofs() const
const IndexSet & get_refinement_edge_indices(unsigned int level) const
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
unsigned int global_dof_index
Definition: types.h:89
const types::subdomain_id artificial_subdomain_id
Definition: types.h:275
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:1255
Definition: mpi.h:90
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:1536
bool is_element(const size_type index) const
Definition: index_set.h:1665
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:231
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:352
bool have_boundary_indices() const
std::vector< AffineConstraints< double > > level_constraints
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int boundary_id
Definition: types.h:111
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:10342
static ::ExceptionBase & ExcInternalError()
void initialize(const DoFHandler< dim, spacedim > &dof)