Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/reduced_dof_map.h"
19 :
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/int_range.h"
24 : #include "libmesh/mesh_base.h"
25 : #include "libmesh/system.h"
26 : #include "timpi/parallel_sync.h"
27 :
28 : #include <algorithm>
29 : #include <set>
30 :
31 : namespace libMesh
32 : {
33 96 : ReducedDofMap::ReducedDofMap(MeshBase & mesh,
34 : System & system,
35 : const DofMap & dof_map,
36 96 : std::string reduced_system_name)
37 : : DofMapBase(dof_map.comm()),
38 32 : _mesh(mesh),
39 32 : _system(system),
40 32 : _dof_map(dof_map),
41 64 : _reduced_system_name(std::move(reduced_system_name)),
42 32 : _parent(nullptr),
43 32 : _initialized(false),
44 96 : _reduced_system(nullptr)
45 : {
46 96 : }
47 :
48 96 : ReducedDofMap::~ReducedDofMap() = default;
49 :
50 4338 : void ReducedDofMap::retained_dof_added(const Elem &, dof_id_type) {}
51 :
52 780 : void ReducedDofMap::unretained_dof_found(const Elem &, dof_id_type) {}
53 :
54 2483979 : bool ReducedDofMap::retains_field_dof(const Elem & elem,
55 : const unsigned int node_num,
56 : const unsigned int var_num,
57 : const dof_id_type field_dof) const
58 : {
59 4967628 : return (!_parent || _parent->retains_field_dof(elem, node_num, var_num, field_dof)) &&
60 3311642 : this->retain_field_dof(elem, node_num, var_num, field_dof);
61 : }
62 :
63 108 : bool ReducedDofMap::retains_scalar_dof(const dof_id_type scalar_dof) const
64 : {
65 216 : return (!_parent || _parent->retains_scalar_dof(scalar_dof)) &&
66 144 : this->retain_scalar_dof(scalar_dof);
67 : }
68 :
69 594801 : void ReducedDofMap::retain_dof(const dof_id_type full_dof)
70 : {
71 594801 : if (_dof_map.local_index(full_dof))
72 193075 : _local_retained_set.insert(full_dof);
73 : else
74 25960 : _nonlocal_retained[_dof_map.dof_owner(full_dof)].insert(full_dof);
75 594801 : }
76 :
77 670080 : void ReducedDofMap::retain_dof_and_constraining_dofs(const Elem * const elem,
78 : const dof_id_type full_dof,
79 : std::unordered_set<dof_id_type> & elem_seen,
80 : std::vector<dof_id_type> & elem_var_dofs)
81 : {
82 670080 : if (!elem_seen.insert(full_dof).second)
83 25097 : return;
84 :
85 : // A constraining dof belongs to this space only if it belongs to the parent space; one the
86 : // parent eliminates still passes its own constraining dofs on
87 594789 : if (_parent && _parent->find_reduced_dof(full_dof) == DofObject::invalid_id)
88 : {
89 0 : this->retain_constraining_dofs(elem, full_dof, elem_seen, elem_var_dofs);
90 0 : return;
91 : }
92 :
93 594789 : this->retain_dof(full_dof);
94 594789 : elem_var_dofs.push_back(full_dof);
95 594789 : this->retained_dof_added(*elem, full_dof);
96 :
97 594789 : this->retain_constraining_dofs(elem, full_dof, elem_seen, elem_var_dofs);
98 : }
99 :
100 727599 : void ReducedDofMap::retain_constraining_dofs(const Elem * const elem,
101 : const dof_id_type full_dof,
102 : std::unordered_set<dof_id_type> & elem_seen,
103 : std::vector<dof_id_type> & elem_var_dofs)
104 : {
105 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
106 727599 : const auto & constraints = _dof_map.get_dof_constraints();
107 727599 : if (const auto it = constraints.find(full_dof); it != constraints.end())
108 236835 : for (const auto & [constraining_dof, weight] : it->second)
109 : {
110 53127 : libmesh_ignore(weight);
111 : // Constraining dofs may themselves be constrained
112 159381 : this->retain_dof_and_constraining_dofs(elem, constraining_dof, elem_seen, elem_var_dofs);
113 : }
114 : #else
115 : libmesh_ignore(elem, full_dof, elem_seen, elem_var_dofs);
116 : #endif
117 727599 : }
118 :
119 402 : void ReducedDofMap::reinit()
120 : {
121 402 : if (this->initialized())
122 312 : this->clear();
123 :
124 : // SCALAR dofs are always retained. They are owned by the last rank
125 894 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
126 492 : if (const auto & vg_description = _dof_map.variable_group(vg);
127 492 : vg_description.type().family == SCALAR)
128 : {
129 8 : std::vector<dof_id_type> scalar_dof_indices;
130 24 : for (const auto vg_vn : make_range(vg_description.n_variables()))
131 : {
132 16 : _dof_map.SCALAR_dof_indices(scalar_dof_indices, vg_description.number(vg_vn));
133 24 : for (const auto dof : scalar_dof_indices)
134 12 : if (this->retains_scalar_dof(dof))
135 12 : this->retain_dof(dof);
136 : }
137 : }
138 :
139 268 : std::vector<dof_id_type> elem_dofs; // only used to satisfy the dof_indices API
140 402 : std::vector<dof_id_type> * elem_var_dofs = nullptr;
141 268 : std::unordered_set<dof_id_type> elem_seen;
142 268 : std::set<unsigned int> full_vars_present;
143 :
144 16 : auto scalar_dofs_functor = [this, &elem_seen, &elem_var_dofs](
145 : const Elem & elem,
146 : std::vector<dof_id_type> & dof_indices,
147 160 : const std::vector<dof_id_type> & scalar_dof_indices) {
148 64 : dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
149 96 : for (const auto dof : scalar_dof_indices)
150 48 : if (this->retains_scalar_dof(dof))
151 48 : this->retain_dof_and_constraining_dofs(&elem, dof, elem_seen, *elem_var_dofs);
152 : else
153 : {
154 0 : if (!_dof_map.local_index(dof))
155 0 : _nonlocal_unretained.insert(dof);
156 0 : this->retain_constraining_dofs(&elem, dof, elem_seen, *elem_var_dofs);
157 : }
158 316 : };
159 :
160 214487 : auto field_dofs_functor = [this, &elem_seen, &elem_var_dofs](const Elem & elem,
161 : const unsigned int node_num,
162 : const unsigned int var_num,
163 : std::vector<dof_id_type> & dof_indices,
164 2587487 : const dof_id_type field_dof) {
165 643461 : dof_indices.push_back(field_dof);
166 643461 : if (this->retains_field_dof(elem, node_num, var_num, field_dof))
167 510651 : this->retain_dof_and_constraining_dofs(&elem, field_dof, elem_seen, *elem_var_dofs);
168 : else
169 : {
170 : // Only dofs of the parent space belong to this space's problem
171 132810 : if (!_parent || _parent->retains_field_dof(elem, node_num, var_num, field_dof))
172 132645 : this->unretained_dof_found(elem, field_dof);
173 132810 : if (!_dof_map.local_index(field_dof))
174 6 : _nonlocal_unretained.insert(field_dof);
175 132810 : this->retain_constraining_dofs(&elem, field_dof, elem_seen, *elem_var_dofs);
176 : }
177 643863 : };
178 :
179 80040 : for (const auto * const elem : _mesh.active_local_element_ptr_range())
180 : {
181 47622 : auto & elem_dofs_by_var = _elem_reduced_dofs[elem->id()];
182 47622 : elem_dofs_by_var.resize(_dof_map.n_variables());
183 15874 : elem_seen.clear();
184 :
185 31748 : const auto sub_id = elem->subdomain_id();
186 95925 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
187 : {
188 48303 : const auto & var_group = _dof_map.variable_group(vg);
189 32202 : if (!var_group.active_on_subdomain(sub_id))
190 0 : continue;
191 97182 : for (const auto v : make_range(var_group.n_variables()))
192 : {
193 48879 : const auto var_num = var_group.number(v);
194 48879 : elem_var_dofs = &elem_dofs_by_var[var_num];
195 81465 : _dof_map.dof_indices(
196 32586 : elem, elem_dofs, var_num, scalar_dofs_functor, field_dofs_functor, elem->p_level());
197 65172 : if (!elem_var_dofs->empty())
198 32074 : full_vars_present.insert(var_num);
199 : }
200 : }
201 134 : }
202 :
203 : //
204 : // A rank may own a retained dof that none of its local elements use, for instance a mid-edge
205 : // node of a coarse element that is a vertex of a finer element owned elsewhere. Every rank that
206 : // retains a nonlocal dof therefore tells the dof's owner
207 : //
208 268 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> nonlocal_retained;
209 672 : for (const auto & [pid, set] : _nonlocal_retained)
210 180 : nonlocal_retained[pid].assign(set.begin(), set.end());
211 134 : _nonlocal_retained.clear();
212 :
213 670 : TIMPI::push_parallel_vector_data(
214 268 : this->comm(),
215 : nonlocal_retained,
216 314 : [this](processor_id_type, const std::vector<dof_id_type> & dofs) {
217 270 : _local_retained_set.insert(dofs.begin(), dofs.end());
218 180 : });
219 :
220 402 : _local_retained_dofs.assign(_local_retained_set.begin(), _local_retained_set.end());
221 134 : _local_retained_set.clear();
222 402 : std::sort(_local_retained_dofs.begin(), _local_retained_dofs.end());
223 :
224 : //
225 : // Number the reduced space
226 : //
227 536 : this->compute_dof_info(_local_retained_dofs.size());
228 :
229 402 : const auto full_first = _dof_map.first_dof();
230 536 : _local_full_to_reduced.assign(_dof_map.n_local_dofs(), DofObject::invalid_id);
231 258483 : for (const auto i : index_range(_local_retained_dofs))
232 430135 : _local_full_to_reduced[_local_retained_dofs[i] - full_first] = this->first_dof() + i;
233 :
234 : // Every nonlocal dof on a local element gets its reduced index, which is invalid_id for an
235 : // eliminated dof, so that find_reduced_dof() can tell eliminated dofs from unknown ones
236 420 : for (const auto dof : _nonlocal_unretained)
237 18 : nonlocal_retained[_dof_map.dof_owner(dof)].push_back(dof);
238 134 : _nonlocal_unretained.clear();
239 :
240 90 : auto gather_reduced = [this](processor_id_type,
241 : const std::vector<dof_id_type> & full_dofs,
242 9879 : std::vector<dof_id_type> & reduced_dofs) {
243 360 : reduced_dofs.resize(full_dofs.size());
244 9969 : for (const auto i : index_range(full_dofs))
245 12932 : reduced_dofs[i] = this->find_reduced_dof(full_dofs[i]);
246 538 : };
247 90 : auto store_reduced = [this](processor_id_type,
248 : const std::vector<dof_id_type> & full_dofs,
249 6646 : const std::vector<dof_id_type> & reduced_dofs) {
250 9969 : for (const auto i : index_range(full_dofs))
251 12932 : _nonlocal_full_to_reduced[full_dofs[i]] = reduced_dofs[i];
252 448 : };
253 402 : TIMPI::pull_parallel_vector_data(this->comm(),
254 : nonlocal_retained,
255 : gather_reduced,
256 : store_reduced,
257 : &DofObject::invalid_id);
258 :
259 : // The parent indices of the local retained dofs, and the ghosts of a vector on this space
260 536 : _local_retained_parent_dofs.resize(_local_retained_dofs.size());
261 258483 : for (const auto i : index_range(_local_retained_dofs))
262 344108 : _local_retained_parent_dofs[i] = this->parent_dof(_local_retained_dofs[i]);
263 10101 : for (const auto & [full_dof, reduced] : _nonlocal_full_to_reduced)
264 : {
265 3233 : libmesh_ignore(full_dof);
266 9699 : if (reduced != DofObject::invalid_id)
267 9681 : _send_list.push_back(reduced);
268 : }
269 402 : std::sort(_send_list.begin(), _send_list.end());
270 :
271 : //
272 : // Determine the variables with any dofs in the reduced space
273 : //
274 402 : this->comm().set_union(full_vars_present);
275 536 : _reduced_vars.reserve(full_vars_present.size());
276 134 : unsigned int first_local_number = 0;
277 870 : for (const auto full_var_num : full_vars_present)
278 : {
279 468 : const auto & full_var = _dof_map.variable(full_var_num);
280 936 : _reduced_vars.push_back(Variable{nullptr,
281 156 : full_var.name(),
282 : cast_int<unsigned int>(_reduced_vars.size()),
283 : first_local_number,
284 : full_var.type()});
285 468 : first_local_number += _reduced_vars.back().n_components(_mesh);
286 : }
287 :
288 : // Renumber the element dofs, keeping only the variables in the reduced space
289 48024 : for (auto & [elem_id, elem_dofs_by_var] : _elem_reduced_dofs)
290 : {
291 15874 : libmesh_ignore(elem_id);
292 47622 : std::vector<std::vector<dof_id_type>> reduced_dofs_by_var;
293 47622 : reduced_dofs_by_var.reserve(full_vars_present.size());
294 95733 : for (const auto full_var_num : full_vars_present)
295 : {
296 16037 : auto & var_dofs = reduced_dofs_by_var.emplace_back(
297 64148 : std::move(elem_dofs_by_var[full_var_num]));
298 642900 : for (auto & dof : var_dofs)
299 594789 : dof = this->reduced_dof(dof);
300 : }
301 15874 : elem_dofs_by_var = std::move(reduced_dofs_by_var);
302 15874 : }
303 :
304 : // The sparsity pattern of the reduced space, with rows indexed by the local full dofs
305 670 : _reduced_sp = _dof_map.build_sparsity(_mesh, /*calculate_constrained=*/false, this);
306 134 : const auto & nnz = _reduced_sp->get_n_nz();
307 134 : const auto & noz = _reduced_sp->get_n_oz();
308 134 : libmesh_assert_equal_to(nnz.size(), noz.size());
309 :
310 : // Exact per-row preallocation counts. These are critical when SCALAR dofs introduce dense rows,
311 : // since the maximum row count would then allocate a dense matrix
312 536 : _reduced_nnz.resize(_local_retained_dofs.size());
313 536 : _reduced_noz.resize(_local_retained_dofs.size());
314 258483 : for (const auto i : index_range(_local_retained_dofs))
315 : {
316 258081 : const dof_id_type local_full_i = _local_retained_dofs[i] - full_first;
317 86027 : libmesh_assert_less(local_full_i, nnz.size());
318 344108 : _reduced_nnz[i] = nnz[local_full_i];
319 430135 : _reduced_noz[i] = noz[local_full_i];
320 : }
321 :
322 : //
323 : // The stand-in system gives the reduced variables a system number. Nodes hold no reduced dofs,
324 : // since reduced dofs are queried by element only
325 : //
326 536 : const std::vector<unsigned int> nvpg(_reduced_vars.size(), 1);
327 :
328 : // add_system returns the existing system if there is one, so this is safe on every reinit
329 402 : _reduced_system = &_system.get_equation_systems().add_system("Basic", _reduced_system_name);
330 402 : if (!_reduced_system->is_initialized())
331 : {
332 84 : _reduced_system->init();
333 34076 : for (auto * const nd : _mesh.active_node_ptr_range())
334 : {
335 25452 : nd->set_n_vars_per_group(_reduced_system->number(), nvpg);
336 53946 : for (const auto g : index_range(nvpg))
337 28494 : nd->set_n_comp_group(_reduced_system->number(), g, 0);
338 28 : }
339 : }
340 :
341 : // The stand-in system holds no solution data worth writing
342 402 : _reduced_system->hide_output() = true;
343 :
344 402 : _initialized = true;
345 402 : }
346 :
347 934674 : dof_id_type ReducedDofMap::find_reduced_dof(const dof_id_type full_dof) const
348 : {
349 934674 : if (_dof_map.local_index(full_dof))
350 1187824 : return _local_full_to_reduced[full_dof - _dof_map.first_dof()];
351 43806 : return libmesh_map_find(_nonlocal_full_to_reduced, full_dof);
352 : }
353 :
354 980397 : dof_id_type ReducedDofMap::parent_dof(const dof_id_type full_dof) const
355 : {
356 980397 : return _parent ? _parent->find_reduced_dof(full_dof) : full_dof;
357 : }
358 :
359 720 : dof_id_type ReducedDofMap::parent_n_dofs() const
360 : {
361 720 : return _parent ? _parent->n_dofs() : _dof_map.n_dofs();
362 : }
363 :
364 720 : dof_id_type ReducedDofMap::parent_n_local_dofs() const
365 : {
366 720 : return _parent ? _parent->n_local_dofs() : _dof_map.n_local_dofs();
367 : }
368 :
369 0 : dof_id_type ReducedDofMap::parent_first_dof() const
370 : {
371 0 : return _parent ? _parent->first_dof() : _dof_map.first_dof();
372 : }
373 :
374 0 : dof_id_type ReducedDofMap::parent_end_dof() const
375 : {
376 0 : return _parent ? _parent->end_dof() : _dof_map.end_dof();
377 : }
378 :
379 360 : const std::vector<dof_id_type> & ReducedDofMap::parent_send_list() const
380 : {
381 360 : return _parent ? _parent->send_list() : _dof_map.get_send_list();
382 : }
383 :
384 595344 : dof_id_type ReducedDofMap::reduced_dof(const dof_id_type full_dof) const
385 : {
386 595344 : const auto reduced = this->find_reduced_dof(full_dof);
387 198448 : libmesh_assert_msg(reduced != DofObject::invalid_id,
388 : "Full dof " << full_dof << " is not retained in the reduced space");
389 595344 : return reduced;
390 : }
391 :
392 : const std::vector<std::vector<dof_id_type>> &
393 178113 : ReducedDofMap::elem_dof_indices(const dof_id_type elem_id) const
394 : {
395 178113 : return libmesh_map_find(_elem_reduced_dofs, elem_id);
396 : }
397 :
398 0 : unsigned int ReducedDofMap::n_variables() const { return _reduced_vars.size(); }
399 :
400 0 : const Variable & ReducedDofMap::variable(const unsigned int c) const
401 : {
402 0 : libmesh_assert_less(c, _reduced_vars.size());
403 0 : return _reduced_vars[c];
404 : }
405 :
406 0 : void ReducedDofMap::dof_indices(const Elem * const elem,
407 : std::vector<dof_id_type> & di,
408 : const unsigned int vn,
409 : int /*p_level*/) const
410 : {
411 0 : di = libmesh_map_find(_elem_reduced_dofs, elem->id())[vn];
412 0 : }
413 :
414 0 : void ReducedDofMap::dof_indices(const Node * const, std::vector<dof_id_type> &, const unsigned int) const
415 : {
416 0 : libmesh_error_msg("ReducedDofMap dof indices are only meant to be queried with elements, not "
417 : "nodes");
418 : }
419 :
420 408 : void ReducedDofMap::clear()
421 : {
422 408 : DofMapBase::clear();
423 136 : _local_retained_dofs.clear();
424 136 : _local_retained_parent_dofs.clear();
425 136 : _send_list.clear();
426 136 : _local_full_to_reduced.clear();
427 136 : _nonlocal_full_to_reduced.clear();
428 136 : _local_retained_set.clear();
429 136 : _nonlocal_retained.clear();
430 136 : _nonlocal_unretained.clear();
431 136 : _elem_reduced_dofs.clear();
432 272 : _reduced_vars.clear();
433 136 : _reduced_sp = nullptr;
434 136 : _reduced_nnz.clear();
435 136 : _reduced_noz.clear();
436 408 : _initialized = false;
437 408 : }
438 : }
|
