Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
35template <int dim, int spacedim>
37class DoFHandler;
38#endif
39
40
48{
49public:
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 = ComponentMask());
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,
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
228private:
232 std::vector<IndexSet> boundary_indices;
233
238 std::vector<IndexSet> refinement_edge_indices;
239
244 std::vector<AffineConstraints<double>> level_constraints;
245
249 std::vector<AffineConstraints<double>> user_constraints;
250};
251
252
253template <int dim, int spacedim>
254inline void
256 const DoFHandler<dim, spacedim> &dof,
257 const MGLevelObject<IndexSet> & level_relevant_dofs)
258{
259 boundary_indices.clear();
261 level_constraints.clear();
262 user_constraints.clear();
263
264 const unsigned int nlevels = dof.get_triangulation().n_global_levels();
265 const unsigned int min_level = level_relevant_dofs.min_level();
266 const unsigned int max_level = (level_relevant_dofs.max_level() == 0) ?
267 nlevels - 1 :
268 level_relevant_dofs.max_level();
269 const bool user_level_dofs =
270 (level_relevant_dofs.max_level() == 0) ? false : true;
271
272 // At this point level_constraint and refinement_edge_indices are empty.
273 refinement_edge_indices.resize(nlevels);
274 level_constraints.resize(nlevels);
275 user_constraints.resize(nlevels);
276 for (unsigned int l = min_level; l <= max_level; ++l)
277 {
278 if (user_level_dofs)
279 {
280 level_constraints[l].reinit(level_relevant_dofs[l]);
281 }
282 else
283 {
284 const IndexSet relevant_dofs =
286 level_constraints[l].reinit(relevant_dofs);
287 }
288
289 // Loop through relevant cells and faces finding those which are periodic
290 // neighbors.
291 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(l),
292 endc = dof.end(l);
293 for (; cell != endc; ++cell)
294 if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
295 {
296 for (auto f : cell->face_indices())
297 if (cell->has_periodic_neighbor(f) &&
298 cell->periodic_neighbor(f)->level() == cell->level())
299 {
300 if (cell->is_locally_owned_on_level())
301 {
302 Assert(
303 cell->periodic_neighbor(f)->level_subdomain_id() !=
306 "Periodic neighbor of a locally owned cell must either be owned or ghost."));
307 }
308 // Cell is a level-ghost and its neighbor is a
309 // level-artificial cell nothing to do here
310 else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
312 {
313 Assert(cell->is_locally_owned_on_level() == false,
315 continue;
316 }
317
318 const unsigned int dofs_per_face =
319 dof.get_fe(0).n_dofs_per_face(f);
320 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
321 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
322
323 cell->periodic_neighbor(f)
324 ->face(cell->periodic_neighbor_face_no(f))
325 ->get_mg_dof_indices(l, dofs_1, 0);
326 cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
327 // Store periodicity information in the level
328 // AffineConstraints object. Skip DoFs for which we've
329 // previously entered periodicity constraints already; this
330 // can happen, for example, for a vertex dof at a periodic
331 // boundary that we visit from more than one cell
332 for (unsigned int i = 0; i < dofs_per_face; ++i)
333 if (level_constraints[l].can_store_line(dofs_2[i]) &&
334 level_constraints[l].can_store_line(dofs_1[i]) &&
335 !level_constraints[l].is_constrained(dofs_2[i]) &&
336 !level_constraints[l].is_constrained(dofs_1[i]))
337 {
338 level_constraints[l].add_line(dofs_2[i]);
339 level_constraints[l].add_entry(dofs_2[i],
340 dofs_1[i],
341 1.);
342 }
343 }
344 }
345 level_constraints[l].close();
346
347 // Initialize with empty IndexSet of correct size
349 }
350
352}
353
354
355template <int dim, int spacedim>
356inline void
358 const DoFHandler<dim, spacedim> & dof,
359 const std::set<types::boundary_id> &boundary_ids,
360 const ComponentMask & component_mask)
361{
362 // allocate an IndexSet for each global level. Contents will be
363 // overwritten inside make_boundary_list.
364 const unsigned int n_levels = dof.get_triangulation().n_global_levels();
365 Assert(boundary_indices.size() == 0 || boundary_indices.size() == n_levels,
367 boundary_indices.resize(n_levels);
368
370 boundary_ids,
372 component_mask);
373}
374
375
376
377template <int dim, int spacedim>
378inline void
380 const unsigned int level,
381 const IndexSet &level_boundary_indices)
382{
383 const unsigned int n_levels = dof.get_triangulation().n_global_levels();
384 if (boundary_indices.size() == 0)
385 {
386 boundary_indices.resize(n_levels);
387 for (unsigned int i = 0; i < n_levels; ++i)
388 boundary_indices[i] = IndexSet(dof.n_dofs(i));
389 }
390 AssertDimension(boundary_indices.size(), n_levels);
391 boundary_indices[level].add_indices(level_boundary_indices);
392}
393
394
395
396template <int dim, int spacedim>
397inline void
399 const DoFHandler<dim, spacedim> &dof,
400 const types::boundary_id bid,
401 const unsigned int first_vector_component)
402{
403 // For a given boundary id, find which vector component is on the boundary
404 // and set a zero boundary constraint for those degrees of freedom.
405 const unsigned int n_components = dof.get_fe_collection().n_components();
406 AssertIndexRange(first_vector_component + dim - 1, n_components);
407
408 ComponentMask comp_mask(n_components, false);
409
410
412 face = dof.get_triangulation().begin_face(),
413 endf = dof.get_triangulation().end_face();
414 for (; face != endf; ++face)
415 if (face->at_boundary() && face->boundary_id() == bid)
416 for (unsigned int d = 0; d < dim; ++d)
417 {
418 Tensor<1, dim, double> unit_vec;
419 unit_vec[d] = 1.0;
420
421 const Tensor<1, dim> normal_vec =
422 face->get_manifold().normal_vector(face, face->center());
423
424 if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
425 comp_mask.set(d + first_vector_component, true);
426 else
427 Assert(
428 std::abs(unit_vec * normal_vec) < 1e-10,
430 "We can currently only support no normal flux conditions "
431 "for a specific boundary id if all faces are normal to the "
432 "x, y, or z axis."));
433 }
434
435 Assert(comp_mask.n_selected_components() == 1,
437 "We can currently only support no normal flux conditions "
438 "for a specific boundary id if all faces are facing in the "
439 "same direction, i.e., a boundary normal to the x-axis must "
440 "have a different boundary id than a boundary normal to the "
441 "y- or z-axis and so on. If the mesh here was produced using "
442 "GridGenerator::..., setting colorize=true during mesh generation "
443 "and calling make_no_normal_flux_constraints() for each no normal "
444 "flux boundary will fulfill the condition."));
445
446 this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
447}
448
449
450inline void
452 const unsigned int level,
453 const AffineConstraints<double> &constraints_on_level)
454{
456
457 // Get the relevant DoFs from level_constraints if
458 // the user constraint matrix has not been initialized
459 if (user_constraints[level].get_local_lines().size() == 0)
460 user_constraints[level].reinit(level_constraints[level].get_local_lines());
461
462 user_constraints[level].merge(
463 constraints_on_level,
465 user_constraints[level].close();
466}
467
468
469
470inline void
472{
473 for (auto &constraint : user_constraints)
474 constraint.clear();
475}
476
477
478
479inline void
481{
482 boundary_indices.clear();
484 user_constraints.clear();
485}
486
487
488inline bool
490 const types::global_dof_index index) const
491{
492 if (boundary_indices.size() == 0)
493 return false;
494
496 return boundary_indices[level].is_element(index);
497}
498
499inline bool
501 const types::global_dof_index index) const
502{
504
505 return refinement_edge_indices[level].is_element(index);
506}
507
508inline bool
510 const unsigned int level,
512 const types::global_dof_index j) const
513{
514 const IndexSet &interface_dofs_on_level =
515 this->get_refinement_edge_indices(level);
516
517 return interface_dofs_on_level.is_element(i) // at_refinement_edge(i)
518 && !interface_dofs_on_level.is_element(j) // !at_refinement_edge(j)
519 && !this->is_boundary_index(level, i) // !on_boundary(i)
520 && !this->is_boundary_index(level, j); // !on_boundary(j)
521}
522
523
524
525inline const IndexSet &
527{
529 return boundary_indices[level];
530}
531
532
533
534inline const IndexSet &
536{
539}
540
541
542
543inline bool
545{
546 return boundary_indices.size() != 0;
547}
548
549
550
551inline const AffineConstraints<double> &
553{
555 return level_constraints[level];
556}
557
558
559
560inline const AffineConstraints<double> &
562{
564 return user_constraints[level];
565}
566
567
568
570
571#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 hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
cell_iterator begin(const unsigned int level=0) 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:1776
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:472
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:160
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int level
Definition grid_out.cc:4618
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
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
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition mg_tools.cc:1446
const types::subdomain_id artificial_subdomain_id
Definition types.h:315
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)