Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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/static_condensation_dof_map.h"
19 :
20 : #include "libmesh/mesh_base.h"
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/int_range.h"
24 : #include "libmesh/system.h"
25 : #include "libmesh/equation_systems.h"
26 : #include "timpi/parallel_sync.h"
27 : #include <unordered_set>
28 :
29 : namespace libMesh
30 : {
31 350 : StaticCondensationDofMap::StaticCondensationDofMap(MeshBase & mesh,
32 : System & system,
33 350 : const DofMap & dof_map)
34 : : DofMapBase(dof_map.comm()),
35 330 : _mesh(mesh),
36 330 : _system(system),
37 330 : _dof_map(dof_map),
38 350 : _sc_is_initialized(false)
39 : {
40 : libmesh_experimental();
41 350 : }
42 :
43 1320 : StaticCondensationDofMap::~StaticCondensationDofMap() = default;
44 :
45 30932 : void StaticCondensationDofMap::add_uncondensed_dof(
46 : const dof_id_type full_dof_number,
47 : const bool involved_in_constraints,
48 : std::unordered_map<dof_id_type, dof_id_type> & uncondensed_global_to_local_map,
49 : std::unordered_set<dof_id_type> & local_uncondensed_dofs_set,
50 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> &
51 : nonlocal_uncondensed_dofs,
52 : std::vector<dof_id_type> & elem_uncondensed_dofs,
53 : dof_id_type & uncondensed_local_dof_number,
54 : std::unordered_set<dof_id_type> & constraint_dofs)
55 : {
56 5624 : if (uncondensed_global_to_local_map.count(full_dof_number))
57 : // We've already seen this dof on this element
58 2320 : return;
59 :
60 28380 : if (_dof_map.local_index(full_dof_number))
61 2458 : local_uncondensed_dofs_set.insert(full_dof_number);
62 : else
63 7826 : nonlocal_uncondensed_dofs[_dof_map.dof_owner(full_dof_number)].insert(full_dof_number);
64 :
65 28380 : elem_uncondensed_dofs.push_back(full_dof_number);
66 28380 : (uncondensed_global_to_local_map)[full_dof_number] = uncondensed_local_dof_number++;
67 28380 : if (involved_in_constraints)
68 232 : constraint_dofs.insert(full_dof_number);
69 : }
70 :
71 30932 : void StaticCondensationDofMap::add_uncondensed_dof_plus_constraint_dofs(
72 : const dof_id_type full_dof_number,
73 : bool involved_in_constraints,
74 : std::unordered_map<dof_id_type, dof_id_type> & uncondensed_global_to_local_map,
75 : std::unordered_set<dof_id_type> & local_uncondensed_dofs_set,
76 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> &
77 : nonlocal_uncondensed_dofs,
78 : std::vector<dof_id_type> & elem_uncondensed_dofs,
79 : dof_id_type & uncondensed_local_dof_number,
80 : std::unordered_set<dof_id_type> & constraint_dofs)
81 : {
82 30932 : const auto & full_dof_constraints = _dof_map.get_dof_constraints();
83 2812 : auto it = full_dof_constraints.find(full_dof_number);
84 2812 : const bool is_constrained = it != full_dof_constraints.end();
85 30932 : involved_in_constraints = involved_in_constraints || is_constrained;
86 :
87 30932 : this->add_uncondensed_dof(full_dof_number,
88 : involved_in_constraints,
89 : uncondensed_global_to_local_map,
90 : local_uncondensed_dofs_set,
91 : nonlocal_uncondensed_dofs,
92 : elem_uncondensed_dofs,
93 : uncondensed_local_dof_number,
94 : constraint_dofs);
95 30932 : if (is_constrained)
96 5104 : for (const auto [full_constraining_dof, weight] : it->second)
97 : {
98 348 : libmesh_ignore(weight);
99 : // Our constraining dofs may themselves be constrained
100 3828 : this->add_uncondensed_dof_plus_constraint_dofs(full_constraining_dof,
101 : /*involved_in_constraints=*/true,
102 : uncondensed_global_to_local_map,
103 : local_uncondensed_dofs_set,
104 : nonlocal_uncondensed_dofs,
105 : elem_uncondensed_dofs,
106 : uncondensed_local_dof_number,
107 : constraint_dofs);
108 : }
109 30932 : }
110 :
111 490 : void StaticCondensationDofMap::reinit()
112 : {
113 490 : if (this->initialized())
114 140 : this->clear();
115 :
116 28 : std::vector<dof_id_type> elem_dofs, elem_uncondensed_dofs; // only used to satisfy API
117 490 : dof_id_type condensed_local_dof_number = 0, uncondensed_local_dof_number = 0;
118 490 : std::unordered_map<dof_id_type, dof_id_type> *condensed_global_to_local_map = nullptr,
119 490 : *uncondensed_global_to_local_map = nullptr;
120 28 : std::set<unsigned int> full_vars_present_in_reduced_sys;
121 28 : std::unordered_set<dof_id_type> local_uncondensed_dofs_set, constraint_dofs;
122 28 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> nonlocal_uncondensed_dofs;
123 :
124 : // Handle SCALAR dofs
125 1960 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
126 1470 : if (const auto & vg_description = _dof_map.variable_group(vg);
127 1470 : vg_description.type().family == SCALAR)
128 : {
129 8 : std::vector<dof_id_type> scalar_dof_indices;
130 140 : const processor_id_type last_pid = this->comm().size() - 1;
131 280 : for (const auto vg_vn : make_range(vg_description.n_variables()))
132 : {
133 8 : const auto vn = vg_description.number(vg_vn);
134 140 : _dof_map.SCALAR_dof_indices(scalar_dof_indices, vn);
135 140 : if (this->comm().rank() == last_pid)
136 20 : local_uncondensed_dofs_set.insert(scalar_dof_indices.begin(),
137 : scalar_dof_indices.end());
138 : else
139 116 : nonlocal_uncondensed_dofs[last_pid].insert(scalar_dof_indices.begin(),
140 : scalar_dof_indices.end());
141 : }
142 : }
143 :
144 : auto scalar_dofs_functor =
145 144 : [this,
146 : &uncondensed_global_to_local_map,
147 : &local_uncondensed_dofs_set,
148 : &nonlocal_uncondensed_dofs,
149 : &elem_uncondensed_dofs,
150 : &uncondensed_local_dof_number,
151 : &constraint_dofs](const Elem & /*elem*/,
152 : std::vector<dof_id_type> & dof_indices,
153 208 : const std::vector<dof_id_type> & scalar_dof_indices) {
154 192 : dof_indices.insert(dof_indices.end(), scalar_dof_indices.begin(), scalar_dof_indices.end());
155 352 : for (const auto global_dof : scalar_dof_indices)
156 176 : this->add_uncondensed_dof_plus_constraint_dofs(global_dof,
157 : false,
158 : *uncondensed_global_to_local_map,
159 : local_uncondensed_dofs_set,
160 : nonlocal_uncondensed_dofs,
161 : elem_uncondensed_dofs,
162 : uncondensed_local_dof_number,
163 : constraint_dofs);
164 652 : };
165 :
166 34074 : auto field_dofs_functor = [this,
167 : &condensed_local_dof_number,
168 : &condensed_global_to_local_map,
169 : &uncondensed_global_to_local_map,
170 : &local_uncondensed_dofs_set,
171 : &nonlocal_uncondensed_dofs,
172 : &elem_uncondensed_dofs,
173 : &uncondensed_local_dof_number,
174 : &constraint_dofs](const Elem & elem,
175 : const unsigned int node_num,
176 : const unsigned int var_num,
177 : std::vector<dof_id_type> & dof_indices,
178 50556 : const dof_id_type field_dof) {
179 41646 : dof_indices.push_back(field_dof);
180 :
181 3786 : bool uncondensed_dof = false;
182 7572 : if (_uncondensed_vars.count(var_num))
183 : {
184 192 : libmesh_assert_msg(
185 : node_num == invalid_uint,
186 : "Users should not be providing continuous FEM variables to the uncondensed vars API");
187 192 : uncondensed_dof = true;
188 : }
189 :
190 41646 : if (node_num != invalid_uint && !elem.is_internal(node_num))
191 2256 : uncondensed_dof = true;
192 :
193 19086 : if (uncondensed_dof)
194 26928 : this->add_uncondensed_dof_plus_constraint_dofs(field_dof,
195 : false,
196 : *uncondensed_global_to_local_map,
197 : local_uncondensed_dofs_set,
198 : nonlocal_uncondensed_dofs,
199 : elem_uncondensed_dofs,
200 : uncondensed_local_dof_number,
201 : constraint_dofs);
202 : else
203 14718 : (*condensed_global_to_local_map)[field_dof] = condensed_local_dof_number++;
204 42122 : };
205 :
206 8216 : for (auto elem : _mesh.active_local_element_ptr_range())
207 : {
208 2750 : auto & dof_data = _elem_to_dof_data[elem->id()];
209 2750 : condensed_local_dof_number = 0;
210 2750 : uncondensed_local_dof_number = 0;
211 2750 : condensed_global_to_local_map = &dof_data.condensed_global_to_local_map;
212 2750 : uncondensed_global_to_local_map = &dof_data.uncondensed_global_to_local_map;
213 :
214 2750 : const auto sub_id = elem->subdomain_id();
215 7788 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
216 : {
217 5038 : const auto & var_group = _dof_map.variable_group(vg);
218 12188 : for (const auto v : make_range(var_group.n_variables()))
219 : {
220 7150 : const auto var_num = var_group.number(v);
221 7150 : dof_data.reduced_space_indices.resize(var_num + 1);
222 6500 : if (!var_group.active_on_subdomain(sub_id))
223 0 : continue;
224 650 : elem_uncondensed_dofs.clear();
225 8450 : _dof_map.dof_indices(elem,
226 : elem_dofs,
227 : var_num,
228 : scalar_dofs_functor,
229 : field_dofs_functor,
230 1300 : elem->p_level());
231 7150 : if (!elem_uncondensed_dofs.empty())
232 : {
233 4334 : auto & var_reduced_space_indices = dof_data.reduced_space_indices[var_num];
234 3940 : var_reduced_space_indices.insert(var_reduced_space_indices.end(),
235 : elem_uncondensed_dofs.begin(),
236 1182 : elem_uncondensed_dofs.end());
237 3940 : full_vars_present_in_reduced_sys.insert(var_num);
238 : }
239 : }
240 : }
241 462 : }
242 :
243 : //
244 : // We've built our local uncondensed dofs container ... but only using local element dof_indices
245 : // calls. It can be the case that we own a degree of freedom that is not actually needed by our
246 : // local element assembly but is needed by other processes element assembly. One example we've run
247 : // into of this is a mid-edge coarse element node holding side hierarchic dofs which is also a
248 : // fine element's vertex node. This node may be owned by the process holding the fine element
249 : // which doesn't need those side hierarchic dofs for its assembly
250 : //
251 :
252 : // Build supported query type. Has to be map to contiguous data for calls to MPI
253 28 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> nonlocal_uncondensed_dofs_mapvec;
254 1239 : for (const auto & [pid, set] : nonlocal_uncondensed_dofs)
255 : {
256 10 : auto & vec = nonlocal_uncondensed_dofs_mapvec[pid];
257 739 : vec.assign(set.begin(), set.end());
258 : }
259 : // clear no longer needed memory
260 14 : nonlocal_uncondensed_dofs.clear();
261 :
262 : auto receive_needed_local_dofs =
263 729 : [&local_uncondensed_dofs_set](processor_id_type,
264 739 : const std::vector<dof_id_type> & local_dofs_to_insert) {
265 749 : local_uncondensed_dofs_set.insert(local_dofs_to_insert.begin(), local_dofs_to_insert.end());
266 1229 : };
267 :
268 490 : TIMPI::push_parallel_vector_data(
269 490 : _mesh.comm(), nonlocal_uncondensed_dofs_mapvec, receive_needed_local_dofs);
270 :
271 490 : _local_uncondensed_dofs.assign(local_uncondensed_dofs_set.begin(),
272 : local_uncondensed_dofs_set.end());
273 14 : local_uncondensed_dofs_set.clear();
274 :
275 : //
276 : // Build the reduced system data
277 : //
278 :
279 28 : const dof_id_type n_local = _local_uncondensed_dofs.size();
280 490 : dof_id_type n = n_local;
281 490 : this->comm().sum(n);
282 :
283 : // Get DOF counts on all processors
284 490 : this->compute_dof_info(n_local);
285 :
286 : // Build a map from the full size problem uncondensed dof indices to the reduced problem
287 : // (uncondensed) dof indices
288 28 : std::unordered_map<dof_id_type, dof_id_type> full_dof_to_reduced_dof;
289 518 : const auto local_start = _first_df[this->processor_id()];
290 15571 : for (const auto i : index_range(_local_uncondensed_dofs))
291 15081 : full_dof_to_reduced_dof[_local_uncondensed_dofs[i]] = i + local_start;
292 :
293 : // Build the condensed system sparsity pattern
294 490 : _reduced_sp = this->_dof_map.build_sparsity(
295 490 : this->_mesh, /*calculate_constrained=*/false, /*use_condensed_system=*/true);
296 14 : const auto & nnz = _reduced_sp->get_n_nz();
297 14 : const auto & noz = _reduced_sp->get_n_oz();
298 14 : libmesh_assert(nnz.size() == noz.size());
299 :
300 : // Optimization for PETSc. This is critical for problems in which there are SCALAR dofs that
301 : // introduce dense rows to avoid allocating a dense matrix
302 504 : _reduced_nnz.resize(_local_uncondensed_dofs.size());
303 504 : _reduced_noz.resize(_local_uncondensed_dofs.size());
304 15571 : for (const dof_id_type local_reduced_i : index_range(_local_uncondensed_dofs))
305 : {
306 15081 : const dof_id_type full_i = _local_uncondensed_dofs[local_reduced_i];
307 15081 : const dof_id_type local_full_i = full_i - _dof_map.first_dof();
308 1371 : libmesh_assert(local_full_i < nnz.size());
309 16452 : _reduced_nnz[local_reduced_i] = nnz[local_full_i];
310 17823 : _reduced_noz[local_reduced_i] = noz[local_full_i];
311 : }
312 :
313 : //
314 : // Now we need to pull our nonlocal data
315 : //
316 :
317 729 : auto gather_functor = [&full_dof_to_reduced_dof](processor_id_type,
318 : const std::vector<dof_id_type> & full_dof_ids,
319 3137 : std::vector<dof_id_type> & reduced_dof_ids) {
320 759 : reduced_dof_ids.resize(full_dof_ids.size());
321 3866 : for (const auto i : index_range(full_dof_ids))
322 3202 : reduced_dof_ids[i] = libmesh_map_find(full_dof_to_reduced_dof, full_dof_ids[i]);
323 1225 : };
324 :
325 : auto action_functor =
326 729 : [&full_dof_to_reduced_dof](processor_id_type,
327 : const std::vector<dof_id_type> & full_dof_ids,
328 3137 : const std::vector<dof_id_type> & reduced_dof_ids) {
329 3866 : for (const auto i : index_range(full_dof_ids))
330 : {
331 85 : libmesh_assert(!full_dof_to_reduced_dof.count(full_dof_ids[i]));
332 3202 : full_dof_to_reduced_dof[full_dof_ids[i]] = reduced_dof_ids[i];
333 : }
334 496 : };
335 :
336 490 : TIMPI::pull_parallel_vector_data(this->comm(),
337 : nonlocal_uncondensed_dofs_mapvec,
338 : gather_functor,
339 : action_functor,
340 : &DofObject::invalid_id);
341 14 : nonlocal_uncondensed_dofs_mapvec.clear();
342 :
343 : // Determine the variables with any degrees of freedom present in the reduced system
344 490 : _communicator.set_union(full_vars_present_in_reduced_sys);
345 504 : _reduced_vars.reserve(full_vars_present_in_reduced_sys.size());
346 14 : unsigned int first_local_number = 0;
347 1680 : for (const auto i : index_range(full_vars_present_in_reduced_sys))
348 : {
349 1190 : const auto full_var_num = *std::next(full_vars_present_in_reduced_sys.begin(), i);
350 1190 : const auto & full_var = _dof_map.variable(full_var_num);
351 3434 : _reduced_vars.push_back(Variable{nullptr,
352 34 : full_var.name(),
353 : cast_int<unsigned int>(i),
354 : first_local_number,
355 : full_var.type()});
356 1190 : first_local_number += _reduced_vars.back().n_components(_mesh);
357 : }
358 :
359 : // Now we can finally set our element reduced dof indices
360 28 : std::vector<dof_id_type> var_full_dof_indices;
361 3240 : for (auto & [elem, dof_data] : _elem_to_dof_data)
362 : {
363 250 : libmesh_ignore(elem);
364 2750 : auto & reduced_space_indices = dof_data.reduced_space_indices;
365 : // Keep around only those variables which are present in our reduced system
366 : {
367 2750 : std::size_t i = 0;
368 : reduced_space_indices.erase(
369 2250 : std::remove_if(
370 : reduced_space_indices.begin(),
371 : reduced_space_indices.end(),
372 6500 : [&full_vars_present_in_reduced_sys, &i](const std::vector<dof_id_type> &) {
373 6500 : return !full_vars_present_in_reduced_sys.count(i++);
374 500 : }),
375 500 : reduced_space_indices.end());
376 : }
377 250 : libmesh_assert(reduced_space_indices.size() == full_vars_present_in_reduced_sys.size());
378 :
379 7084 : for (auto & var_dof_indices : reduced_space_indices)
380 : {
381 4334 : var_full_dof_indices = var_dof_indices;
382 394 : var_dof_indices.clear();
383 32714 : for (const auto full_dof : var_full_dof_indices)
384 28380 : var_dof_indices.push_back(libmesh_map_find(full_dof_to_reduced_dof, full_dof));
385 : }
386 : }
387 :
388 : // Build our dof constraints map
389 2640 : for (const auto full_dof : constraint_dofs)
390 2150 : _full_to_reduced_constraint_dofs[full_dof] =
391 2150 : libmesh_map_find(full_dof_to_reduced_dof, full_dof);
392 14 : constraint_dofs.clear();
393 :
394 : // Prevent querying Nodes for dof indices
395 504 : std::vector<unsigned int> nvpg(_reduced_vars.size());
396 1680 : for (auto & elem : nvpg)
397 1190 : elem = 1;
398 :
399 : // add_system returns a system if it already exists instead of erroring so there's no harm if
400 : // we do this multiple times
401 504 : _reduced_system = &_system.get_equation_systems().add_system("Basic", "reduced");
402 490 : if (!_reduced_system->is_initialized())
403 : {
404 350 : _reduced_system->init();
405 20650 : for (auto * const nd : _mesh.active_node_ptr_range())
406 : {
407 10582 : nd->set_n_vars_per_group(_reduced_system->number(), nvpg);
408 36540 : for (const auto g : index_range(nvpg))
409 25958 : nd->set_n_comp_group(_reduced_system->number(), g, 0);
410 330 : }
411 : }
412 :
413 490 : _sc_is_initialized = true;
414 490 : }
415 :
416 0 : unsigned int StaticCondensationDofMap::n_variables() const { return _reduced_vars.size(); }
417 :
418 0 : const Variable & StaticCondensationDofMap::variable(const unsigned int c) const
419 : {
420 0 : return _reduced_vars[c];
421 : }
422 :
423 0 : void StaticCondensationDofMap::dof_indices(const Elem * const elem,
424 : std::vector<dof_id_type> & di,
425 : const unsigned int vn,
426 : int /*p_level*/) const
427 : {
428 0 : di.clear();
429 0 : di = libmesh_map_find(_elem_to_dof_data, elem->id()).reduced_space_indices[vn];
430 0 : }
431 :
432 0 : void StaticCondensationDofMap::dof_indices(const Node * const,
433 : std::vector<dof_id_type> &,
434 : const unsigned int) const
435 : {
436 0 : libmesh_error_msg("StaticCondensationDofMap dof indices are only meant to be queried with "
437 : "elements, not nodes");
438 : }
439 :
440 490 : void StaticCondensationDofMap::clear()
441 : {
442 490 : DofMapBase::clear();
443 14 : _elem_to_dof_data.clear();
444 14 : _uncondensed_vars.clear();
445 476 : _reduced_vars.clear();
446 14 : _reduced_sp = nullptr;
447 14 : _reduced_nnz.clear();
448 14 : _reduced_noz.clear();
449 490 : _sc_is_initialized = false;
450 490 : }
451 : }
|
