Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30:01+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
mg_constrained_dofs.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
19
22
23#include <set>
24
25
27
28
29template <int dim, int spacedim>
30void
33 const MGLevelObject<IndexSet> &level_relevant_dofs)
34{
35 boundary_indices.clear();
37 level_constraints.clear();
38 user_constraints.clear();
39
40 const unsigned int nlevels = dof.get_triangulation().n_global_levels();
41 const unsigned int min_level = level_relevant_dofs.min_level();
42 const unsigned int max_level = (level_relevant_dofs.max_level() == 0) ?
43 nlevels - 1 :
44 level_relevant_dofs.max_level();
45 const bool use_provided_level_relevant_dofs =
46 (level_relevant_dofs.max_level() > 0);
47
48 // At this point level_constraint and refinement_edge_indices are empty.
49 refinement_edge_indices.resize(nlevels);
50 level_constraints.resize(nlevels);
51 user_constraints.resize(nlevels);
52 for (unsigned int l = min_level; l <= max_level; ++l)
53 {
54 if (use_provided_level_relevant_dofs)
55 {
57 level_relevant_dofs[l]);
58 }
59 else
60 {
61 const IndexSet relevant_dofs =
64 relevant_dofs);
65 }
66 }
67
68 // TODO: currently we only consider very basic periodic constraints
69 const IdentityMatrix transformation(dof.get_fe().n_dofs_per_face());
70 const ComponentMask component_mask;
71 const double periodicity_factor = 1.0;
72
73 for (const auto &[first_cell, second_cell] :
75 {
76 // only consider non-artificial cells
77 if (first_cell.first->is_artificial_on_level())
78 continue;
79 if (second_cell.first.first->is_artificial_on_level())
80 continue;
81
82 // consider cell pairs with the same level
83 if (first_cell.first->level() != second_cell.first.first->level())
84 continue;
85
87 first_cell.first->as_dof_handler_level_iterator(dof)->face(
88 first_cell.second),
89 second_cell.first.first->as_dof_handler_level_iterator(dof)->face(
90 second_cell.first.second),
91 transformation,
92 level_constraints[first_cell.first->level()],
93 component_mask,
94 second_cell.second,
95 periodicity_factor,
96 first_cell.first->level());
97 }
98
99 for (unsigned int l = min_level; l <= max_level; ++l)
100 {
101 level_constraints[l].close();
102
103 // Initialize with empty IndexSet of correct size
105 }
106
108}
109
110
111
112template <int dim, int spacedim>
113void
115 const DoFHandler<dim, spacedim> &dof,
116 const std::set<types::boundary_id> &boundary_ids,
117 const ComponentMask &component_mask)
118{
119 // allocate an IndexSet for each global level. Contents will be
120 // overwritten inside make_boundary_list.
121 const unsigned int n_levels = dof.get_triangulation().n_global_levels();
122 Assert(boundary_indices.empty() || boundary_indices.size() == n_levels,
124 boundary_indices.resize(n_levels);
125
127 boundary_ids,
129 component_mask);
130}
131
132
133
134template <int dim, int spacedim>
135void
137 const unsigned int level,
138 const IndexSet &level_boundary_indices)
139{
140 const unsigned int n_levels = dof.get_triangulation().n_global_levels();
141 if (boundary_indices.empty())
142 {
143 boundary_indices.resize(n_levels);
144 for (unsigned int i = 0; i < n_levels; ++i)
145 boundary_indices[i] = IndexSet(dof.n_dofs(i));
146 }
147 AssertDimension(boundary_indices.size(), n_levels);
148 boundary_indices[level].add_indices(level_boundary_indices);
149}
150
151
152
153template <int dim, int spacedim>
154void
156 const DoFHandler<dim, spacedim> &dof,
157 const types::boundary_id bid,
158 const unsigned int first_vector_component)
159{
160 // For a given boundary id, find which vector component is on the boundary
161 // and set a zero boundary constraint for those degrees of freedom.
162 const unsigned int n_components = dof.get_fe_collection().n_components();
163 AssertIndexRange(first_vector_component + dim - 1, n_components);
164
165 ComponentMask comp_mask(n_components, false);
166
167
169 face = dof.get_triangulation().begin_face(),
170 endf = dof.get_triangulation().end_face();
171 for (; face != endf; ++face)
172 if (face->at_boundary() && face->boundary_id() == bid)
173 for (unsigned int d = 0; d < dim; ++d)
174 {
175 Tensor<1, dim, double> unit_vec;
176 unit_vec[d] = 1.0;
177
178 const Tensor<1, dim> normal_vec =
179 face->get_manifold().normal_vector(face, face->center());
180
181 if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
182 comp_mask.set(d + first_vector_component, true);
183 else
184 Assert(
185 std::abs(unit_vec * normal_vec) < 1e-10,
187 "We can currently only support no normal flux conditions "
188 "for a specific boundary id if all faces are normal to the "
189 "x, y, or z axis."));
190 }
191
192 Assert(comp_mask.n_selected_components() == 1,
194 "We can currently only support no normal flux conditions "
195 "for a specific boundary id if all faces are facing in the "
196 "same direction, i.e., a boundary normal to the x-axis must "
197 "have a different boundary id than a boundary normal to the "
198 "y- or z-axis and so on. If the mesh here was produced using "
199 "GridGenerator::..., setting colorize=true during mesh generation "
200 "and calling make_no_normal_flux_constraints() for each no normal "
201 "flux boundary will fulfill the condition."));
202
203 this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
204}
205
206
207
208void
210 const unsigned int level,
211 const AffineConstraints<double> &constraints_on_level)
212{
214
215 // Get the relevant DoFs from level_constraints if
216 // the user constraint matrix has not been initialized
217 if (user_constraints[level].get_local_lines().size() == 0)
218 user_constraints[level].reinit(
219 level_constraints[level].get_locally_owned_indices(),
220 level_constraints[level].get_local_lines());
221
222 user_constraints[level].merge(
223 constraints_on_level,
225 user_constraints[level].close();
226}
227
228
229
230void
232{
233 for (auto &constraint : user_constraints)
234 constraint.clear();
235}
236
237
238
239void
246
247
248
249template <typename Number>
250void
252 const unsigned int level,
253 const bool add_boundary_indices,
254 const bool add_refinement_edge_indices,
255 const bool add_level_constraints,
256 const bool add_user_constraints) const
257{
258 constraints.clear();
259
260 // determine local lines
261 IndexSet index_set(this->get_refinement_edge_indices(level).size());
262
264 index_set.add_indices(this->get_boundary_indices(level));
265
266 if (add_refinement_edge_indices)
267 index_set.add_indices(this->get_refinement_edge_indices(level));
268
269 if (add_level_constraints)
270 index_set.add_indices(this->get_level_constraints(level).get_local_lines());
271
273 index_set.add_indices(
274 this->get_user_constraint_matrix(level).get_local_lines());
275
276 constraints.reinit(level_constraints[level].get_locally_owned_indices(),
277 index_set);
278
279 // merge constraints
281 for (const auto i : this->get_boundary_indices(level))
282 constraints.constrain_dof_to_zero(i);
283
284 if (add_refinement_edge_indices)
285 for (const auto i : this->get_refinement_edge_indices(level))
286 constraints.constrain_dof_to_zero(i);
287
288 if (add_level_constraints)
289 constraints.merge(this->get_level_constraints(level),
291 true);
292
294 constraints.merge(this->get_user_constraint_matrix(level),
296 true);
297
298 // finalize setup
299 constraints.close();
300}
301
302#include "multigrid/mg_constrained_dofs.inst"
303
304
void merge(const AffineConstraints< other_number > &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
void constrain_dof_to_zero(const size_type constrained_dof)
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
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
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)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const std::set< types::boundary_id > &boundary_ids, const ComponentMask &component_mask={})
std::vector< IndexSet > boundary_indices
void add_boundary_indices(const DoFHandler< dim, spacedim > &dof, const unsigned int level, const IndexSet &boundary_indices)
void initialize(const DoFHandler< dim, spacedim > &dof, const MGLevelObject< IndexSet > &level_relevant_dofs=MGLevelObject< IndexSet >())
void merge_constraints(AffineConstraints< Number > &constraints, const unsigned int level, const bool add_boundary_indices, const bool add_refinement_edge_indices, const bool add_level_constraints, const bool add_user_constraints) 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
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
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, types::geometric_orientation > > & get_periodic_face_map() const
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
unsigned int level
Definition grid_out.cc:4635
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::size_t size
Definition mpi.cc:745
void set_periodicity_constraints(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, const FullMatrix< double > &transformation, AffineConstraints< number > &affine_constraints, const ComponentMask &component_mask, const types::geometric_orientation combined_orientation, const number periodicity_factor, const unsigned int level=numbers::invalid_unsigned_int)
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={})
Definition mg_tools.cc:1241
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition mg_tools.cc:1445
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)