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 490 : _local_uncondensed_dofs.assign(local_uncondensed_dofs_set.begin(),
244 : local_uncondensed_dofs_set.end());
245 14 : local_uncondensed_dofs_set.clear();
246 :
247 : //
248 : // Build the reduced system data
249 : //
250 :
251 28 : const dof_id_type n_local = _local_uncondensed_dofs.size();
252 490 : dof_id_type n = n_local;
253 490 : this->comm().sum(n);
254 :
255 : // Get DOF counts on all processors
256 490 : this->compute_dof_info(n_local);
257 :
258 : // Build a map from the full size problem uncondensed dof indices to the reduced problem
259 : // (uncondensed) dof indices
260 28 : std::unordered_map<dof_id_type, dof_id_type> full_dof_to_reduced_dof;
261 518 : const auto local_start = _first_df[this->processor_id()];
262 15571 : for (const auto i : index_range(_local_uncondensed_dofs))
263 15081 : full_dof_to_reduced_dof[_local_uncondensed_dofs[i]] = i + local_start;
264 :
265 : // Build the condensed system sparsity pattern
266 490 : _reduced_sp = this->_dof_map.build_sparsity(
267 490 : this->_mesh, /*calculate_constrained=*/false, /*use_condensed_system=*/true);
268 14 : const auto & nnz = _reduced_sp->get_n_nz();
269 14 : const auto & noz = _reduced_sp->get_n_oz();
270 14 : libmesh_assert(nnz.size() == noz.size());
271 :
272 : // Optimization for PETSc. This is critical for problems in which there are SCALAR dofs that
273 : // introduce dense rows to avoid allocating a dense matrix
274 504 : _reduced_nnz.resize(_local_uncondensed_dofs.size());
275 504 : _reduced_noz.resize(_local_uncondensed_dofs.size());
276 15571 : for (const dof_id_type local_reduced_i : index_range(_local_uncondensed_dofs))
277 : {
278 15081 : const dof_id_type full_i = _local_uncondensed_dofs[local_reduced_i];
279 15081 : const dof_id_type local_full_i = full_i - _dof_map.first_dof();
280 1371 : libmesh_assert(local_full_i < nnz.size());
281 16452 : _reduced_nnz[local_reduced_i] = nnz[local_full_i];
282 17823 : _reduced_noz[local_reduced_i] = noz[local_full_i];
283 : }
284 :
285 : //
286 : // Now we need to pull our nonlocal data
287 : //
288 :
289 : // build our queries
290 28 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> nonlocal_uncondensed_dofs_mapvec;
291 1239 : for (const auto & [pid, set] : nonlocal_uncondensed_dofs)
292 : {
293 10 : auto & vec = nonlocal_uncondensed_dofs_mapvec[pid];
294 739 : vec.assign(set.begin(), set.end());
295 : }
296 : // clear no longer needed memory
297 14 : nonlocal_uncondensed_dofs.clear();
298 :
299 729 : auto gather_functor = [&full_dof_to_reduced_dof](processor_id_type,
300 : const std::vector<dof_id_type> & full_dof_ids,
301 3137 : std::vector<dof_id_type> & reduced_dof_ids) {
302 759 : reduced_dof_ids.resize(full_dof_ids.size());
303 3866 : for (const auto i : index_range(full_dof_ids))
304 3202 : reduced_dof_ids[i] = libmesh_map_find(full_dof_to_reduced_dof, full_dof_ids[i]);
305 1225 : };
306 :
307 : auto action_functor =
308 729 : [&full_dof_to_reduced_dof](processor_id_type,
309 : const std::vector<dof_id_type> & full_dof_ids,
310 3137 : const std::vector<dof_id_type> & reduced_dof_ids) {
311 3866 : for (const auto i : index_range(full_dof_ids))
312 : {
313 85 : libmesh_assert(!full_dof_to_reduced_dof.count(full_dof_ids[i]));
314 3202 : full_dof_to_reduced_dof[full_dof_ids[i]] = reduced_dof_ids[i];
315 : }
316 496 : };
317 :
318 490 : TIMPI::pull_parallel_vector_data(this->comm(),
319 : nonlocal_uncondensed_dofs_mapvec,
320 : gather_functor,
321 : action_functor,
322 : &DofObject::invalid_id);
323 14 : nonlocal_uncondensed_dofs_mapvec.clear();
324 :
325 : // Determine the variables with any degrees of freedom present in the reduced system
326 490 : _communicator.set_union(full_vars_present_in_reduced_sys);
327 504 : _reduced_vars.reserve(full_vars_present_in_reduced_sys.size());
328 14 : unsigned int first_local_number = 0;
329 1680 : for (const auto i : index_range(full_vars_present_in_reduced_sys))
330 : {
331 1190 : const auto full_var_num = *std::next(full_vars_present_in_reduced_sys.begin(), i);
332 1190 : const auto & full_var = _dof_map.variable(full_var_num);
333 3434 : _reduced_vars.push_back(Variable{nullptr,
334 34 : full_var.name(),
335 : cast_int<unsigned int>(i),
336 : first_local_number,
337 : full_var.type()});
338 1190 : first_local_number += _reduced_vars.back().n_components(_mesh);
339 : }
340 :
341 : // Now we can finally set our element reduced dof indices
342 28 : std::vector<dof_id_type> var_full_dof_indices;
343 3240 : for (auto & [elem, dof_data] : _elem_to_dof_data)
344 : {
345 250 : libmesh_ignore(elem);
346 2750 : auto & reduced_space_indices = dof_data.reduced_space_indices;
347 : // Start erasing from the back to reduce the number of copy assignment operations
348 3000 : if (reduced_space_indices.size())
349 7150 : for (typename std::vector<dof_id_type>::difference_type i =
350 2750 : reduced_space_indices.size() - 1;
351 9900 : i >= 0;
352 : --i)
353 7150 : if (!full_vars_present_in_reduced_sys.count(i))
354 256 : reduced_space_indices.erase(reduced_space_indices.begin() + i);
355 : // It is theoretically possible that we have an element that doesn't have dofs for one of the
356 : // variables present in our reduced system, which is why the assertion below is not an
357 : // equality assertion
358 250 : libmesh_assert(reduced_space_indices.size() <= full_vars_present_in_reduced_sys.size());
359 :
360 7084 : for (auto & var_dof_indices : reduced_space_indices)
361 : {
362 4334 : var_full_dof_indices = var_dof_indices;
363 394 : var_dof_indices.clear();
364 32714 : for (const auto full_dof : var_full_dof_indices)
365 28380 : var_dof_indices.push_back(libmesh_map_find(full_dof_to_reduced_dof, full_dof));
366 : }
367 : }
368 :
369 : // Build our dof constraints map
370 2640 : for (const auto full_dof : constraint_dofs)
371 2150 : _full_to_reduced_constraint_dofs[full_dof] =
372 2150 : libmesh_map_find(full_dof_to_reduced_dof, full_dof);
373 14 : constraint_dofs.clear();
374 :
375 : // Prevent querying Nodes for dof indices
376 504 : std::vector<unsigned int> nvpg(_reduced_vars.size());
377 1680 : for (auto & elem : nvpg)
378 1190 : elem = 1;
379 :
380 : // add_system returns a system if it already exists instead of erroring so there's no harm if
381 : // we do this multiple times
382 504 : _reduced_system = &_system.get_equation_systems().add_system("Basic", "reduced");
383 490 : if (!_reduced_system->is_initialized())
384 : {
385 350 : _reduced_system->init();
386 20650 : for (auto * const nd : _mesh.active_node_ptr_range())
387 : {
388 10582 : nd->set_n_vars_per_group(_reduced_system->number(), nvpg);
389 36540 : for (const auto g : index_range(nvpg))
390 25958 : nd->set_n_comp_group(_reduced_system->number(), g, 0);
391 330 : }
392 : }
393 :
394 490 : _sc_is_initialized = true;
395 490 : }
396 :
397 0 : unsigned int StaticCondensationDofMap::n_variables() const { return _reduced_vars.size(); }
398 :
399 0 : const Variable & StaticCondensationDofMap::variable(const unsigned int c) const
400 : {
401 0 : return _reduced_vars[c];
402 : }
403 :
404 0 : void StaticCondensationDofMap::dof_indices(const Elem * const elem,
405 : std::vector<dof_id_type> & di,
406 : const unsigned int vn,
407 : int /*p_level*/) const
408 : {
409 0 : di.clear();
410 0 : di = libmesh_map_find(_elem_to_dof_data, elem->id()).reduced_space_indices[vn];
411 0 : }
412 :
413 0 : void StaticCondensationDofMap::dof_indices(const Node * const,
414 : std::vector<dof_id_type> &,
415 : const unsigned int) const
416 : {
417 0 : libmesh_error_msg("StaticCondensationDofMap dof indices are only meant to be queried with "
418 : "elements, not nodes");
419 : }
420 :
421 490 : void StaticCondensationDofMap::clear()
422 : {
423 490 : DofMapBase::clear();
424 14 : _elem_to_dof_data.clear();
425 14 : _uncondensed_vars.clear();
426 476 : _reduced_vars.clear();
427 14 : _reduced_sp = nullptr;
428 14 : _reduced_nnz.clear();
429 14 : _reduced_noz.clear();
430 490 : _sc_is_initialized = false;
431 490 : }
432 : }
|
