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/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 560 : StaticCondensationDofMap::StaticCondensationDofMap(MeshBase & mesh,
32 : System & system,
33 560 : const DofMap & dof_map)
34 : : DofMapBase(dof_map.comm()),
35 528 : _mesh(mesh),
36 528 : _system(system),
37 528 : _dof_map(dof_map),
38 560 : _sc_is_initialized(false)
39 : {
40 : libmesh_experimental();
41 560 : }
42 :
43 2112 : StaticCondensationDofMap::~StaticCondensationDofMap() = default;
44 :
45 2607326 : 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 443714 : if (uncondensed_global_to_local_map.count(full_dof_number))
57 : // We've already seen this dof on this element
58 267246 : return;
59 :
60 2315086 : if (_dof_map.local_index(full_dof_number))
61 191825 : local_uncondensed_dofs_set.insert(full_dof_number);
62 : else
63 337506 : nonlocal_uncondensed_dofs[_dof_map.dof_owner(full_dof_number)].insert(full_dof_number);
64 :
65 2315086 : elem_uncondensed_dofs.push_back(full_dof_number);
66 2315086 : (uncondensed_global_to_local_map)[full_dof_number] = uncondensed_local_dof_number++;
67 2315086 : if (involved_in_constraints)
68 57042 : constraint_dofs.insert(full_dof_number);
69 : }
70 :
71 2607326 : 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 2607326 : const auto & full_dof_constraints = _dof_map.get_dof_constraints();
83 221857 : auto it = full_dof_constraints.find(full_dof_number);
84 221857 : const bool is_constrained = it != full_dof_constraints.end();
85 2607326 : involved_in_constraints = involved_in_constraints || is_constrained;
86 :
87 2607326 : 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 2607326 : if (is_constrained)
96 920266 : for (const auto & [full_constraining_dof, weight] : it->second)
97 : {
98 52887 : libmesh_ignore(weight);
99 : // Our constraining dofs may themselves be constrained
100 620595 : 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 2607326 : }
110 :
111 4130 : void StaticCondensationDofMap::reinit()
112 : {
113 4130 : if (this->initialized())
114 3640 : this->clear();
115 :
116 236 : std::vector<dof_id_type> elem_dofs, elem_uncondensed_dofs; // only used to satisfy API
117 4130 : dof_id_type condensed_local_dof_number = 0, uncondensed_local_dof_number = 0;
118 4130 : std::unordered_map<dof_id_type, dof_id_type> *condensed_global_to_local_map = nullptr,
119 4130 : *uncondensed_global_to_local_map = nullptr;
120 236 : std::set<unsigned int> full_vars_present_in_reduced_sys;
121 236 : std::unordered_set<dof_id_type> local_uncondensed_dofs_set, constraint_dofs;
122 236 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> nonlocal_uncondensed_dofs;
123 :
124 : // Handle SCALAR dofs
125 9240 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
126 5110 : if (const auto & vg_description = _dof_map.variable_group(vg);
127 5110 : 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 4188 : };
165 :
166 2069760 : 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 2965076 : const dof_id_type field_dof) {
179 2495468 : dof_indices.push_back(field_dof);
180 :
181 212854 : bool uncondensed_dof = false;
182 425708 : 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 2495468 : if (node_num != invalid_uint && !elem.is_internal(node_num))
191 168762 : uncondensed_dof = true;
192 :
193 679787 : if (uncondensed_dof)
194 1986555 : 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 508913 : (*condensed_global_to_local_map)[field_dof] = condensed_local_dof_number++;
204 2499480 : };
205 :
206 502957 : for (auto elem : _mesh.active_local_element_ptr_range())
207 : {
208 185929 : auto & dof_data = _elem_to_dof_data[elem->id()];
209 185929 : condensed_local_dof_number = 0;
210 185929 : uncondensed_local_dof_number = 0;
211 185929 : condensed_global_to_local_map = &dof_data.condensed_global_to_local_map;
212 185929 : uncondensed_global_to_local_map = &dof_data.uncondensed_global_to_local_map;
213 :
214 185929 : const auto sub_id = elem->subdomain_id();
215 374146 : for (const auto vg : make_range(_dof_map.n_variable_groups()))
216 : {
217 188217 : const auto & var_group = _dof_map.variable_group(vg);
218 378546 : for (const auto v : make_range(var_group.n_variables()))
219 : {
220 190329 : const auto var_num = var_group.number(v);
221 190329 : dof_data.reduced_space_indices.resize(var_num + 1);
222 174186 : if (!var_group.active_on_subdomain(sub_id))
223 0 : continue;
224 16143 : elem_uncondensed_dofs.clear();
225 222615 : _dof_map.dof_indices(elem,
226 : elem_dofs,
227 : var_num,
228 : scalar_dofs_functor,
229 : field_dofs_functor,
230 32286 : elem->p_level());
231 190329 : if (!elem_uncondensed_dofs.empty())
232 : {
233 187513 : auto & var_reduced_space_indices = dof_data.reduced_space_indices[var_num];
234 171626 : var_reduced_space_indices.insert(var_reduced_space_indices.end(),
235 : elem_uncondensed_dofs.begin(),
236 47661 : elem_uncondensed_dofs.end());
237 171626 : full_vars_present_in_reduced_sys.insert(var_num);
238 : }
239 : }
240 : }
241 3894 : }
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 236 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> nonlocal_uncondensed_dofs_mapvec;
254 12887 : for (const auto & [pid, set] : nonlocal_uncondensed_dofs)
255 : {
256 87 : auto & vec = nonlocal_uncondensed_dofs_mapvec[pid];
257 8670 : vec.assign(set.begin(), set.end());
258 : }
259 : // clear no longer needed memory
260 118 : nonlocal_uncondensed_dofs.clear();
261 :
262 : auto receive_needed_local_dofs =
263 8583 : [&local_uncondensed_dofs_set](processor_id_type,
264 8670 : const std::vector<dof_id_type> & local_dofs_to_insert) {
265 8757 : local_uncondensed_dofs_set.insert(local_dofs_to_insert.begin(), local_dofs_to_insert.end());
266 12800 : };
267 :
268 4130 : TIMPI::push_parallel_vector_data(
269 4130 : _mesh.comm(), nonlocal_uncondensed_dofs_mapvec, receive_needed_local_dofs);
270 :
271 4130 : _local_uncondensed_dofs.assign(local_uncondensed_dofs_set.begin(),
272 : local_uncondensed_dofs_set.end());
273 118 : local_uncondensed_dofs_set.clear();
274 :
275 : //
276 : // Build the reduced system data
277 : //
278 :
279 236 : const dof_id_type n_local = _local_uncondensed_dofs.size();
280 4130 : dof_id_type n = n_local;
281 4130 : this->comm().sum(n);
282 :
283 : // Get DOF counts on all processors
284 4130 : 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 236 : std::unordered_map<dof_id_type, dof_id_type> full_dof_to_reduced_dof;
289 4366 : const auto local_start = _first_df[this->processor_id()];
290 1003940 : for (const auto i : index_range(_local_uncondensed_dofs))
291 999810 : full_dof_to_reduced_dof[_local_uncondensed_dofs[i]] = i + local_start;
292 :
293 : // Build the condensed system sparsity pattern
294 4130 : _reduced_sp = this->_dof_map.build_sparsity(
295 4130 : this->_mesh, /*calculate_constrained=*/false, /*use_condensed_system=*/true);
296 118 : const auto & nnz = _reduced_sp->get_n_nz();
297 118 : const auto & noz = _reduced_sp->get_n_oz();
298 118 : 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 4248 : _reduced_nnz.resize(_local_uncondensed_dofs.size());
303 4248 : _reduced_noz.resize(_local_uncondensed_dofs.size());
304 1003940 : for (const dof_id_type local_reduced_i : index_range(_local_uncondensed_dofs))
305 : {
306 999810 : const dof_id_type full_i = _local_uncondensed_dofs[local_reduced_i];
307 999810 : const dof_id_type local_full_i = full_i - _dof_map.first_dof();
308 85309 : libmesh_assert(local_full_i < nnz.size());
309 1085119 : _reduced_nnz[local_reduced_i] = nnz[local_full_i];
310 1170428 : _reduced_noz[local_reduced_i] = noz[local_full_i];
311 : }
312 :
313 : //
314 : // Now we need to pull our nonlocal data
315 : //
316 :
317 8583 : auto gather_functor = [&full_dof_to_reduced_dof](processor_id_type,
318 : const std::vector<dof_id_type> & full_dof_ids,
319 112856 : std::vector<dof_id_type> & reduced_dof_ids) {
320 8844 : reduced_dof_ids.resize(full_dof_ids.size());
321 121439 : for (const auto i : index_range(full_dof_ids))
322 115871 : reduced_dof_ids[i] = libmesh_map_find(full_dof_to_reduced_dof, full_dof_ids[i]);
323 12769 : };
324 :
325 : auto action_functor =
326 8583 : [&full_dof_to_reduced_dof](processor_id_type,
327 : const std::vector<dof_id_type> & full_dof_ids,
328 112856 : const std::vector<dof_id_type> & reduced_dof_ids) {
329 121439 : for (const auto i : index_range(full_dof_ids))
330 : {
331 3189 : libmesh_assert(!full_dof_to_reduced_dof.count(full_dof_ids[i]));
332 115871 : full_dof_to_reduced_dof[full_dof_ids[i]] = reduced_dof_ids[i];
333 : }
334 4186 : };
335 :
336 4130 : TIMPI::pull_parallel_vector_data(this->comm(),
337 : nonlocal_uncondensed_dofs_mapvec,
338 : gather_functor,
339 : action_functor,
340 : &DofObject::invalid_id);
341 118 : nonlocal_uncondensed_dofs_mapvec.clear();
342 :
343 : // Determine the variables with any degrees of freedom present in the reduced system
344 4130 : _communicator.set_union(full_vars_present_in_reduced_sys);
345 4248 : _reduced_vars.reserve(full_vars_present_in_reduced_sys.size());
346 118 : unsigned int first_local_number = 0;
347 8960 : for (const auto i : index_range(full_vars_present_in_reduced_sys))
348 : {
349 4830 : const auto full_var_num = *std::next(full_vars_present_in_reduced_sys.begin(), i);
350 4830 : const auto & full_var = _dof_map.variable(full_var_num);
351 13938 : _reduced_vars.push_back(Variable{nullptr,
352 138 : full_var.name(),
353 : cast_int<unsigned int>(i),
354 : first_local_number,
355 : full_var.type()});
356 4830 : first_local_number += _reduced_vars.back().n_components(_mesh);
357 : }
358 :
359 : // Now we can finally set our element reduced dof indices
360 236 : std::vector<dof_id_type> var_full_dof_indices;
361 190059 : for (auto & [elem, dof_data] : _elem_to_dof_data)
362 : {
363 15743 : libmesh_ignore(elem);
364 185929 : auto & reduced_space_indices = dof_data.reduced_space_indices;
365 : // Keep around only those variables which are present in our reduced system
366 : {
367 185929 : std::size_t i = 0;
368 : reduced_space_indices.erase(
369 154443 : std::remove_if(
370 : reduced_space_indices.begin(),
371 : reduced_space_indices.end(),
372 174186 : [&full_vars_present_in_reduced_sys, &i](const std::vector<dof_id_type> &) {
373 174186 : return !full_vars_present_in_reduced_sys.count(i++);
374 31486 : }),
375 31486 : reduced_space_indices.end());
376 : }
377 15743 : libmesh_assert(reduced_space_indices.size() == full_vars_present_in_reduced_sys.size());
378 :
379 373442 : for (auto & var_dof_indices : reduced_space_indices)
380 : {
381 187513 : var_full_dof_indices = var_dof_indices;
382 15887 : var_dof_indices.clear();
383 2502599 : for (const auto full_dof : var_full_dof_indices)
384 2315086 : 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 398068 : for (const auto full_dof : constraint_dofs)
390 393938 : _full_to_reduced_constraint_dofs[full_dof] =
391 393938 : libmesh_map_find(full_dof_to_reduced_dof, full_dof);
392 118 : constraint_dofs.clear();
393 :
394 : // Prevent querying Nodes for dof indices
395 4248 : std::vector<unsigned int> nvpg(_reduced_vars.size());
396 8960 : for (auto & elem : nvpg)
397 4830 : 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 4248 : _reduced_system = &_system.get_equation_systems().add_system("Basic", "reduced");
402 4130 : if (!_reduced_system->is_initialized())
403 : {
404 490 : _reduced_system->init();
405 464078 : for (auto * const nd : _mesh.active_node_ptr_range())
406 : {
407 238682 : nd->set_n_vars_per_group(_reduced_system->number(), nvpg);
408 492740 : for (const auto g : index_range(nvpg))
409 254058 : nd->set_n_comp_group(_reduced_system->number(), g, 0);
410 462 : }
411 : }
412 :
413 : // We don't want to write the reduced system
414 4130 : _reduced_system->hide_output() = true;
415 :
416 4130 : _sc_is_initialized = true;
417 4130 : }
418 :
419 0 : unsigned int StaticCondensationDofMap::n_variables() const { return _reduced_vars.size(); }
420 :
421 0 : const Variable & StaticCondensationDofMap::variable(const unsigned int c) const
422 : {
423 0 : return _reduced_vars[c];
424 : }
425 :
426 0 : void StaticCondensationDofMap::dof_indices(const Elem * const elem,
427 : std::vector<dof_id_type> & di,
428 : const unsigned int vn,
429 : int /*p_level*/) const
430 : {
431 0 : di.clear();
432 0 : di = libmesh_map_find(_elem_to_dof_data, elem->id()).reduced_space_indices[vn];
433 0 : }
434 :
435 0 : void StaticCondensationDofMap::dof_indices(const Node * const,
436 : std::vector<dof_id_type> &,
437 : const unsigned int) const
438 : {
439 0 : libmesh_error_msg("StaticCondensationDofMap dof indices are only meant to be queried with "
440 : "elements, not nodes");
441 : }
442 :
443 4200 : void StaticCondensationDofMap::clear()
444 : {
445 4200 : DofMapBase::clear();
446 120 : _elem_to_dof_data.clear();
447 120 : _uncondensed_vars.clear();
448 4080 : _reduced_vars.clear();
449 120 : _reduced_sp = nullptr;
450 120 : _reduced_nnz.clear();
451 120 : _reduced_noz.clear();
452 4200 : _sc_is_initialized = false;
453 4200 : }
454 : }
|
