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 :
19 : // Local Includes
20 : #include "libmesh/default_coupling.h" // For downconversion
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/eigen_system.h"
23 : #include "libmesh/elem.h"
24 : #include "libmesh/explicit_system.h"
25 : #include "libmesh/fe_interface.h"
26 : #include "libmesh/frequency_system.h"
27 : #include "libmesh/int_range.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/linear_implicit_system.h"
30 : #include "libmesh/mesh_base.h"
31 : #include "libmesh/mesh_refinement.h"
32 : #include "libmesh/newmark_system.h"
33 : #include "libmesh/nonlinear_implicit_system.h"
34 : #include "libmesh/parallel.h"
35 : #include "libmesh/rb_construction.h"
36 : #include "libmesh/remote_elem.h"
37 : #include "libmesh/transient_rb_construction.h"
38 : #include "libmesh/transient_system.h"
39 :
40 : // System includes
41 : #include <functional> // std::plus
42 : #include <numeric> // std::iota
43 : #include <sstream>
44 :
45 : // Include the systems before this one to avoid
46 : // overlapping forward declarations.
47 : #include "libmesh/equation_systems.h"
48 :
49 : namespace libMesh
50 : {
51 :
52 238871 : EquationSystems::EquationSystems (MeshBase & m) :
53 : ParallelObject (m),
54 225279 : _mesh (m),
55 225279 : _refine_in_reinit(true),
56 245667 : _enable_default_ghosting(true)
57 : {
58 : // Set default parameters
59 238871 : this->parameters.set<Real> ("linear solver tolerance") = TOLERANCE * TOLERANCE;
60 238871 : this->parameters.set<unsigned int>("linear solver maximum iterations") = 5000;
61 238871 : }
62 :
63 :
64 :
65 411667 : EquationSystems::~EquationSystems () = default;
66 :
67 :
68 :
69 71 : void EquationSystems::clear ()
70 : {
71 : // Clear any additional parameters
72 2 : parameters.clear ();
73 :
74 : // Clear the systems.
75 2 : _systems.clear();
76 71 : }
77 :
78 :
79 :
80 238398 : void EquationSystems::init ()
81 : {
82 : #ifndef NDEBUG
83 13640 : for (auto i : make_range(this->n_systems()))
84 6856 : libmesh_assert(!this->get_system(i).is_initialized());
85 : #endif
86 :
87 238398 : this->reinit_mesh();
88 238304 : }
89 :
90 :
91 :
92 25688 : void EquationSystems::reinit ()
93 : {
94 25688 : const bool mesh_changed = this->reinit_solutions();
95 :
96 : // If the mesh has changed, systems will need to reinitialize their
97 : // own data on the new mesh.
98 25688 : if (mesh_changed)
99 25688 : this->reinit_systems();
100 25688 : }
101 :
102 238824 : void EquationSystems::reinit_mesh ()
103 : {
104 6796 : const unsigned int n_sys = this->n_systems();
105 :
106 6796 : libmesh_assert_not_equal_to (n_sys, 0);
107 :
108 : // Our DofMaps are going to expect a prepared mesh later, but our
109 : // older codes might have used some utilities that mark a mesh
110 : // unprepared without re-preparing it. At *this* point hopefully
111 : // everybody's done changing the mesh and we can make sure it's
112 : // prepared.
113 238824 : if (!_mesh.is_prepared())
114 0 : _mesh.complete_preparation();
115 :
116 : // Tell all the \p DofObject entities how many systems
117 : // there are.
118 27538916 : for (auto & node : _mesh.node_ptr_range())
119 14611100 : node->set_n_systems(n_sys);
120 :
121 : Threads::parallel_for
122 238824 : (_mesh.element_stored_range(),
123 840244 : [n_sys](const ElemRange & range)
124 : {
125 6635719 : for (Elem * elem : range)
126 6396797 : elem->set_n_systems(n_sys);
127 232094 : });
128 :
129 : //for (auto i : make_range(this->n_systems()))
130 : //this->get_system(i).init();
131 :
132 : #ifdef LIBMESH_ENABLE_AMR
133 252416 : MeshRefinement mesh_refine(_mesh);
134 238824 : mesh_refine.clean_refinement_flags();
135 : #endif
136 :
137 : // Now loop over all the systems belonging to this ES
138 : // and call reinit_mesh for each system
139 480170 : for (auto i : make_range(this->n_systems()))
140 241440 : this->get_system(i).reinit_mesh();
141 :
142 238816 : }
143 :
144 25688 : bool EquationSystems::reinit_solutions ()
145 : {
146 854 : parallel_object_only();
147 :
148 854 : const unsigned int n_sys = this->n_systems();
149 854 : libmesh_assert_not_equal_to (n_sys, 0);
150 :
151 : // And any new systems will need initialization
152 55623 : for (unsigned int i=0; i != n_sys; ++i)
153 29935 : if (!this->get_system(i).is_initialized())
154 143 : this->get_system(i).init();
155 :
156 : // We used to assert that all nodes and elements *already* had
157 : // n_systems() properly set; however this is false in the case where
158 : // user code has manually added nodes and/or elements to an
159 : // already-initialized system.
160 :
161 : // Make sure all the \p DofObject entities know how many systems
162 : // there are.
163 : {
164 : // All the nodes
165 24007926 : for (auto & node : _mesh.node_ptr_range())
166 12740412 : node->set_n_systems(n_sys);
167 :
168 : // All the elements
169 : Threads::parallel_for
170 25688 : (_mesh.element_stored_range(),
171 590632 : [n_sys](const ElemRange & range)
172 : {
173 10490319 : for (Elem * elem : range)
174 10464429 : elem->set_n_systems(n_sys);
175 24968 : });
176 : }
177 :
178 : // Localize each system's vectors
179 55623 : for (unsigned int i=0; i != n_sys; ++i)
180 29935 : this->get_system(i).re_update();
181 :
182 : #ifdef LIBMESH_ENABLE_AMR
183 :
184 854 : bool mesh_changed = false;
185 :
186 : // FIXME: For backwards compatibility, assume
187 : // refine_and_coarsen_elements or refine_uniformly have already
188 : // been called
189 : {
190 55623 : for (unsigned int i=0; i != n_sys; ++i)
191 : {
192 966 : System & sys = this->get_system(i);
193 :
194 : // Even if the system doesn't have any variables in it we want
195 : // consistent behavior; e.g. distribute_dofs should have the
196 : // opportunity to count up zero dofs on each processor.
197 : //
198 : // Who's been adding zero-var systems anyway, outside of my
199 : // unit tests? - RHS
200 : // if (!sys.n_vars())
201 : // continue;
202 :
203 29935 : sys.get_dof_map().distribute_dofs(_mesh);
204 :
205 : // Recreate any user or internal constraints
206 29935 : sys.reinit_constraints();
207 :
208 : // Even if there weren't any constraint changes,
209 : // reinit_constraints() did prepare_send_list() for us.
210 :
211 29935 : sys.prolong_vectors();
212 : }
213 854 : mesh_changed = true;
214 : }
215 :
216 25688 : if (this->_refine_in_reinit)
217 : {
218 : // Don't override any user refinement settings
219 27246 : MeshRefinement mesh_refine(_mesh);
220 25546 : mesh_refine.face_level_mismatch_limit() = 0; // unlimited
221 25546 : mesh_refine.overrefined_boundary_limit() = -1; // unlimited
222 25546 : mesh_refine.underrefined_boundary_limit() = -1; // unlimited
223 :
224 : // Try to coarsen the mesh, then restrict each system's vectors
225 : // if necessary
226 25546 : if (mesh_refine.coarsen_elements())
227 : {
228 0 : for (auto i : make_range(this->n_systems()))
229 : {
230 0 : System & sys = this->get_system(i);
231 0 : sys.get_dof_map().distribute_dofs(_mesh);
232 0 : sys.reinit_constraints();
233 :
234 : // Even if there weren't any constraint changes,
235 : // reinit_constraints() did prepare_send_list() for us.
236 :
237 0 : sys.restrict_vectors();
238 : }
239 0 : mesh_changed = true;
240 : }
241 :
242 : // Once vectors are all restricted, we can delete
243 : // children of coarsened elements
244 850 : if (mesh_changed)
245 25546 : this->get_mesh().contract();
246 :
247 : // Try to refine the mesh, then prolong each system's vectors
248 : // if necessary
249 25546 : if (mesh_refine.refine_elements())
250 : {
251 2414 : for (auto i : make_range(this->n_systems()))
252 : {
253 34 : System & sys = this->get_system(i);
254 1207 : sys.get_dof_map().distribute_dofs(_mesh);
255 1207 : sys.reinit_constraints();
256 :
257 : // Even if there weren't any constraint changes,
258 : // reinit_constraints() did prepare_send_list() for us.
259 :
260 1207 : sys.prolong_vectors();
261 : }
262 34 : mesh_changed = true;
263 : }
264 23846 : }
265 :
266 854 : return mesh_changed;
267 :
268 : #endif // #ifdef LIBMESH_ENABLE_AMR
269 :
270 : return false;
271 : }
272 :
273 :
274 :
275 25688 : void EquationSystems::reinit_systems()
276 : {
277 55623 : for (auto i : make_range(this->n_systems()))
278 29935 : this->get_system(i).reinit();
279 25688 : }
280 :
281 :
282 :
283 0 : void EquationSystems::allgather ()
284 : {
285 : // A serial mesh means nothing needs to be done
286 0 : if (_mesh.is_serial())
287 0 : return;
288 :
289 0 : const unsigned int n_sys = this->n_systems();
290 :
291 0 : libmesh_assert_not_equal_to (n_sys, 0);
292 :
293 : // Gather the mesh
294 0 : _mesh.allgather();
295 :
296 : // Tell all the \p DofObject entities how many systems
297 : // there are.
298 0 : for (auto & node : _mesh.node_ptr_range())
299 0 : node->set_n_systems(n_sys);
300 :
301 : Threads::parallel_for
302 0 : (_mesh.element_stored_range(),
303 0 : [n_sys](const ElemRange & range)
304 : {
305 0 : for (Elem * elem : range)
306 0 : elem->set_n_systems(n_sys);
307 0 : });
308 :
309 : // And distribute each system's dofs
310 0 : for (auto i : make_range(this->n_systems()))
311 : {
312 0 : System & sys = this->get_system(i);
313 0 : DofMap & dof_map = sys.get_dof_map();
314 0 : dof_map.distribute_dofs(_mesh);
315 :
316 : // The user probably won't need constraint equations or the
317 : // send_list after an allgather, but let's keep it in consistent
318 : // shape just in case.
319 0 : sys.reinit_constraints();
320 :
321 : // Even if there weren't any constraint changes,
322 : // reinit_constraints() did prepare_send_list() for us.
323 : }
324 : }
325 :
326 :
327 :
328 213 : void EquationSystems::enable_default_ghosting (bool enable)
329 : {
330 213 : _enable_default_ghosting = enable;
331 12 : MeshBase &mesh = this->get_mesh();
332 :
333 213 : if (enable)
334 71 : mesh.add_ghosting_functor(mesh.default_ghosting());
335 : else
336 142 : mesh.remove_ghosting_functor(mesh.default_ghosting());
337 :
338 426 : for (auto i : make_range(this->n_systems()))
339 : {
340 6 : DofMap & dof_map = this->get_system(i).get_dof_map();
341 213 : if (enable)
342 71 : dof_map.add_default_ghosting();
343 : else
344 142 : dof_map.remove_default_ghosting();
345 : }
346 213 : }
347 :
348 :
349 :
350 10928 : void EquationSystems::update ()
351 : {
352 628 : LOG_SCOPE("update()", "EquationSystems");
353 :
354 : // Localize each system's vectors
355 22557 : for (auto i : make_range(this->n_systems()))
356 11629 : this->get_system(i).update();
357 10928 : }
358 :
359 :
360 :
361 5049 : System & EquationSystems::add_system (std::string_view sys_type,
362 : std::string_view name)
363 : {
364 : // If the user already built a system with this name, we'll
365 : // trust them and we'll use it. That way they can pre-add
366 : // non-standard derived system classes, and if their restart file
367 : // has some non-standard sys_type we won't throw an error.
368 5049 : if (_systems.count(name))
369 : {
370 3645 : return this->get_system(name);
371 : }
372 : // Build a basic System
373 1380 : else if (sys_type == "Basic")
374 841 : this->add_system<System> (name);
375 :
376 : // Build a Newmark system
377 563 : else if (sys_type == "Newmark")
378 0 : this->add_system<NewmarkSystem> (name);
379 :
380 : // Build an Explicit system
381 561 : else if ((sys_type == "Explicit"))
382 71 : this->add_system<ExplicitSystem> (name);
383 :
384 : // Build an Implicit system
385 506 : else if ((sys_type == "Implicit") ||
386 506 : (sys_type == "Steady" ))
387 0 : this->add_system<ImplicitSystem> (name);
388 :
389 : // build a transient implicit linear system
390 1434 : else if ((sys_type == "Transient") ||
391 506 : (sys_type == "TransientImplicit") ||
392 506 : (sys_type == "TransientLinearImplicit"))
393 212 : this->add_system<TransientLinearImplicitSystem> (name);
394 :
395 : // build a transient implicit nonlinear system
396 280 : else if (sys_type == "TransientNonlinearImplicit")
397 0 : this->add_system<TransientNonlinearImplicitSystem> (name);
398 :
399 : // build a transient explicit system
400 280 : else if (sys_type == "TransientExplicit")
401 0 : this->add_system<TransientExplicitSystem> (name);
402 :
403 : // build a linear implicit system
404 280 : else if (sys_type == "LinearImplicit")
405 0 : this->add_system<LinearImplicitSystem> (name);
406 :
407 : // build a nonlinear implicit system
408 272 : else if (sys_type == "NonlinearImplicit")
409 280 : this->add_system<NonlinearImplicitSystem> (name);
410 :
411 : // build a Reduced Basis Construction system
412 0 : else if (sys_type == "RBConstruction")
413 0 : this->add_system<RBConstruction> (name);
414 :
415 : // build a transient Reduced Basis Construction system
416 0 : else if (sys_type == "TransientRBConstruction")
417 0 : this->add_system<TransientRBConstruction> (name);
418 :
419 : #ifdef LIBMESH_HAVE_SLEPC
420 : // build an eigen system
421 0 : else if (sys_type == "Eigen")
422 0 : this->add_system<EigenSystem> (name);
423 0 : else if (sys_type == "TransientEigenSystem")
424 0 : this->add_system<TransientEigenSystem> (name);
425 : #endif
426 :
427 : #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
428 : // build a frequency system
429 0 : else if (sys_type == "Frequency")
430 0 : this->add_system<FrequencySystem> (name);
431 : #endif
432 :
433 : else
434 0 : libmesh_error_msg("ERROR: Unknown system type: " << sys_type);
435 :
436 : // Return a reference to the new system
437 : //return (*this)(name);
438 1404 : return this->get_system(name);
439 : }
440 :
441 :
442 :
443 0 : void EquationSystems::solve ()
444 : {
445 0 : libmesh_assert (this->n_systems());
446 :
447 0 : for (auto i : make_range(this->n_systems()))
448 0 : this->get_system(i).solve();
449 0 : }
450 :
451 :
452 :
453 0 : void EquationSystems::sensitivity_solve (const ParameterVector & parameters_in)
454 : {
455 0 : libmesh_assert (this->n_systems());
456 :
457 0 : for (auto i : make_range(this->n_systems()))
458 0 : this->get_system(i).sensitivity_solve(parameters_in);
459 0 : }
460 :
461 :
462 :
463 0 : void EquationSystems::adjoint_solve (const QoISet & qoi_indices)
464 : {
465 0 : libmesh_assert (this->n_systems());
466 :
467 0 : for (unsigned int i=this->n_systems(); i != 0; --i)
468 0 : this->get_system(i-1).adjoint_solve(qoi_indices);
469 0 : }
470 :
471 :
472 :
473 86112 : void EquationSystems::build_variable_names (std::vector<std::string> & var_names,
474 : const FEType * type,
475 : const std::set<std::string> * system_names) const
476 : {
477 : // start indexing at end of possibly non-empty vector of variable names to avoid overwriting them
478 88812 : unsigned int var_num = var_names.size();
479 :
480 : // We'll want to double-check that we don't have any naming
481 : // conflicts; this API causes problems down the line if so.
482 5400 : std::unordered_multiset<std::string> seen_names;
483 :
484 : // Need to size var_names by scalar variables plus all the
485 : // vector components for all the vector variables
486 : //Could this be replaced by a/some convenience methods?[PB]
487 : {
488 2700 : unsigned int n_scalar_vars = 0;
489 2700 : unsigned int n_vector_vars = 0;
490 :
491 176010 : for (const auto & [sys_name, sys_ptr] : _systems)
492 : {
493 : // Check current system is listed in system_names, and skip pos if not
494 5656 : bool use_current_system = (system_names == nullptr);
495 89898 : if (!use_current_system)
496 846 : use_current_system = system_names->count(sys_name);
497 89898 : if (!use_current_system || sys_ptr->hide_output())
498 : {
499 4292 : for (auto vn : make_range(sys_ptr->n_vars()))
500 2246 : seen_names.insert(sys_ptr->variable_name(vn));
501 1994 : continue;
502 1942 : }
503 :
504 209280 : for (auto vn : make_range(sys_ptr->n_vars()))
505 : {
506 121428 : seen_names.insert(sys_ptr->variable_name(vn));
507 121428 : if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
508 6879 : n_vector_vars++;
509 : else
510 114549 : n_scalar_vars++;
511 : }
512 : }
513 :
514 : // Here, we're assuming the number of vector components is the same
515 : // as the mesh spatial dimension.
516 5400 : unsigned int dim = this->get_mesh().spatial_dimension();
517 86112 : unsigned int nv = n_scalar_vars + dim*n_vector_vars;
518 :
519 : // We'd better not have more than dim*this->n_vars() (all vector variables)
520 : // Treat the NodeElem-only mesh case as dim=1
521 2700 : libmesh_assert_less_equal ( nv, (dim > 0 ? dim : 1)*this->n_vars() );
522 :
523 : // 'nv' represents the max possible number of output variables, so allocate enough memory for
524 : // all variables in the system to be populated here. When this is called more than once on a
525 : // single 'var_names' vector, different filters should be used such that duplicates don't occur.
526 86112 : var_names.resize( nv );
527 : }
528 :
529 175939 : for (const auto & [sys_name, sys_ptr] : _systems)
530 : {
531 : // Check current system is listed in system_names, and skip pos if not
532 5656 : bool use_current_system = (system_names == nullptr);
533 89898 : if (!use_current_system)
534 846 : use_current_system = system_names->count(sys_name);
535 89898 : if (!use_current_system || sys_ptr->hide_output())
536 1994 : continue;
537 :
538 209209 : for (auto vn : make_range(sys_ptr->n_vars()))
539 : {
540 121428 : const std::string & var_name = sys_ptr->variable_name(vn);
541 121428 : const FEType & fe_type = sys_ptr->variable_type(vn);
542 :
543 121428 : unsigned int n_vec_dim = FEInterface::n_vec_dim( sys_ptr->get_mesh(), fe_type);
544 :
545 : // Filter on the type if requested
546 121428 : if (type == nullptr || (type && *type == fe_type))
547 : {
548 121215 : if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
549 : {
550 6879 : switch(n_vec_dim)
551 : {
552 0 : case 0:
553 : case 1:
554 0 : var_names[var_num++] = var_name;
555 0 : libmesh_error_msg_if(seen_names.count(var_name) > 1,
556 : "Duplicate variable name "+var_name);
557 0 : break;
558 5622 : case 2:
559 5622 : var_names[var_num++] = var_name+"_x";
560 5622 : var_names[var_num++] = var_name+"_y";
561 11307 : libmesh_error_msg_if(seen_names.count(var_name+"_x"),
562 : "Duplicate variable name "+var_name+"_x");
563 10960 : libmesh_error_msg_if(seen_names.count(var_name+"_y"),
564 : "Duplicate variable name "+var_name+"_y");
565 142 : break;
566 1257 : case 3:
567 1257 : var_names[var_num++] = var_name+"_x";
568 1257 : var_names[var_num++] = var_name+"_y";
569 1257 : var_names[var_num++] = var_name+"_z";
570 2486 : libmesh_error_msg_if(seen_names.count(var_name+"_x"),
571 : "Duplicate variable name "+var_name+"_x");
572 2486 : libmesh_error_msg_if(seen_names.count(var_name+"_y"),
573 : "Duplicate variable name "+var_name+"_y");
574 2553 : libmesh_error_msg_if(seen_names.count(var_name+"_z"),
575 : "Duplicate variable name "+var_name+"_z");
576 28 : break;
577 0 : default:
578 0 : libmesh_error_msg("Invalid dim in build_variable_names");
579 : }
580 : }
581 : else
582 114336 : var_names[var_num++] = var_name;
583 : }
584 : }
585 : }
586 : // Now resize again in case we filtered any names
587 86041 : var_names.resize(var_num);
588 86041 : }
589 :
590 11819 : bool EquationSystems::is_elemental_data_fe_type (const FEType & type)
591 : {
592 12151 : return type.order == CONSTANT &&
593 11749 : (type.family == MONOMIAL ||
594 658 : type.family == MONOMIAL_VEC ||
595 11837 : type.family == XYZ);
596 : }
597 :
598 :
599 :
600 70 : void EquationSystems::build_elemental_data_variable_names
601 : (std::vector<std::string> & var_names,
602 : const std::set<std::string> * system_names) const
603 : {
604 74 : const std::vector<std::string> name_filter = var_names;
605 2 : const bool is_names_empty = name_filter.empty();
606 2 : var_names.clear();
607 :
608 280 : const std::vector<std::string> component_suffix = {"_x", "_y", "_z"};
609 70 : const unsigned int dim = _mesh.spatial_dimension();
610 70 : libmesh_error_msg_if(dim > 3, "Invalid dim in build_elemental_data_variable_names");
611 :
612 140 : for (const auto & [sys_name, sys_ptr] : _systems)
613 : {
614 70 : const bool use_current_system = (system_names == nullptr) || system_names->count(sys_name);
615 70 : if (!use_current_system || sys_ptr->hide_output())
616 0 : continue;
617 :
618 350 : for (auto var : make_range(sys_ptr->n_vars()))
619 : {
620 280 : const FEType & var_type = sys_ptr->variable_type(var);
621 280 : if (!EquationSystems::is_elemental_data_fe_type(var_type))
622 68 : continue;
623 :
624 210 : if (FEInterface::field_type(var_type) == TYPE_VECTOR)
625 : {
626 280 : for (auto comp : make_range(dim))
627 : {
628 : const std::string name =
629 222 : sys_ptr->variable_name(var) + component_suffix[comp];
630 :
631 210 : if (is_names_empty ||
632 6 : std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end())
633 210 : var_names.push_back(name);
634 : }
635 : }
636 : else
637 : {
638 140 : const std::string & name = sys_ptr->variable_name(var);
639 :
640 140 : if (is_names_empty ||
641 4 : std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end())
642 140 : var_names.push_back(name);
643 : }
644 : }
645 : }
646 :
647 70 : std::sort(var_names.begin(), var_names.end());
648 202 : }
649 :
650 :
651 :
652 0 : void EquationSystems::build_solution_vector (std::vector<Number> &,
653 : std::string_view,
654 : std::string_view) const
655 : {
656 : // TODO:[BSK] re-implement this from the method below
657 0 : libmesh_not_implemented();
658 : }
659 :
660 :
661 :
662 :
663 : std::unique_ptr<NumericVector<Number>>
664 79685 : EquationSystems::build_parallel_solution_vector(const std::set<std::string> * system_names,
665 : bool add_sides) const
666 : {
667 4888 : LOG_SCOPE("build_parallel_solution_vector()", "EquationSystems");
668 :
669 : // This function must be run on all processors at once
670 2444 : parallel_object_only();
671 :
672 79685 : const unsigned int dim = _mesh.spatial_dimension();
673 79685 : const dof_id_type max_nn = _mesh.max_node_id();
674 :
675 : // allocate vector storage to hold
676 : // (max_node_id)*(number_of_variables) entries.
677 : //
678 : // If node renumbering is disabled and adaptive coarsening has
679 : // created gaps between node numbers, then this vector will be
680 : // sparse.
681 : //
682 : // We have to differentiate between between scalar and vector
683 : // variables. We intercept vector variables and treat each
684 : // component as a scalar variable (consistently with build_solution_names).
685 :
686 2444 : unsigned int nv = 0;
687 :
688 : //Could this be replaced by a/some convenience methods?[PB]
689 : {
690 2444 : unsigned int n_scalar_vars = 0;
691 2444 : unsigned int n_vector_vars = 0;
692 163014 : for (const auto & [sys_name, sys_ptr] : _systems)
693 : {
694 : // Check current system is listed in system_names, and skip pos if not
695 5136 : bool use_current_system = (system_names == nullptr);
696 83329 : if (!use_current_system)
697 352 : use_current_system = system_names->count(sys_name);
698 83329 : if (!use_current_system || sys_ptr->hide_output())
699 1582 : continue;
700 :
701 195360 : for (auto vn : make_range(sys_ptr->n_vars()))
702 : {
703 113653 : if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
704 6455 : n_vector_vars++;
705 : else
706 107198 : n_scalar_vars++;
707 : }
708 : }
709 : // Here, we're assuming the number of vector components is the same
710 : // as the mesh spatial dimension.
711 79685 : nv = n_scalar_vars + dim*n_vector_vars;
712 : }
713 :
714 : // Get the number of nodes to store locally.
715 : dof_id_type n_local_nodes = cast_int<dof_id_type>
716 156926 : (std::distance(_mesh.local_nodes_begin(),
717 159370 : _mesh.local_nodes_end()));
718 :
719 : // If node renumbering has been disabled, nodes may not be numbered
720 : // contiguously, and the number of nodes might not match the
721 : // max_node_id. In this case we just do our best.
722 79685 : dof_id_type n_total_nodes = n_local_nodes;
723 79685 : _mesh.comm().sum(n_total_nodes);
724 :
725 79685 : const processor_id_type n_proc = _mesh.comm().size();
726 4888 : const processor_id_type my_pid = _mesh.comm().rank();
727 79685 : const dof_id_type n_gaps = max_nn - n_total_nodes;
728 79685 : const dof_id_type gaps_per_processor = n_gaps / n_proc;
729 79685 : const dof_id_type remainder_gaps = n_gaps % n_proc;
730 :
731 84573 : n_local_nodes = n_local_nodes + // Actual nodes
732 79685 : gaps_per_processor + // Our even share of gaps
733 79685 : (my_pid < remainder_gaps); // Leftovers
734 :
735 : // If we've been asked to build added sides' data, we need space to
736 : // add it. Keep track of how much space.
737 79685 : dof_id_type local_added_side_nodes = 0,
738 2444 : added_side_nodes = 0;
739 :
740 : // others_added_side_nodes[p]: local_added_side_nodes on rank p
741 4888 : std::vector<dof_id_type> others_added_side_nodes;
742 :
743 : // A map of (element_id, side, side_node) pairs to the corresponding
744 : // added side node index.
745 : std::map<std::tuple<dof_id_type, unsigned short, unsigned short>,
746 4888 : dof_id_type> discontinuous_node_indices;
747 :
748 : // If we don't have any added side nodes, we'll have no offsets from
749 : // them, and we won't care about which offsets apply to which node
750 : // ids either.
751 :
752 : // Number of true nodes on processors [0,p]
753 4888 : std::vector<dof_id_type> true_node_offsets;
754 : // Number of added (fake) nodes on processors [0,p)
755 4888 : std::vector<dof_id_type> added_node_offsets;
756 :
757 : auto node_id_to_vec_id =
758 69683863 : [&true_node_offsets, &added_node_offsets]
759 84331379 : (const dof_id_type node_id)
760 : {
761 84321119 : if (true_node_offsets.empty())
762 76979919 : return node_id; // O(1) in the common !add_sides case
763 :
764 : // Find the processor id that has node_id in the parallel vec
765 20520 : const auto lb = std::upper_bound(true_node_offsets.begin(),
766 2052 : true_node_offsets.end(), node_id);
767 2052 : libmesh_assert(lb != true_node_offsets.end());
768 2052 : const processor_id_type p = lb - true_node_offsets.begin();
769 :
770 26676 : return node_id + added_node_offsets[p];
771 79685 : };
772 :
773 79685 : if (add_sides)
774 : {
775 1207 : true_node_offsets.resize(n_proc);
776 1207 : added_node_offsets.resize(n_proc);
777 :
778 : // One loop to count everyone's new side nodes
779 10940 : for (const auto & elem : _mesh.active_element_ptr_range())
780 : {
781 25064 : for (auto s : elem->side_index_range())
782 : {
783 20400 : if (redundant_added_side(*elem,s))
784 5460 : continue;
785 :
786 : const std::vector<unsigned int> side_nodes =
787 15548 : elem->nodes_on_side(s);
788 :
789 15548 : if (elem->processor_id() == this->processor_id())
790 3952 : local_added_side_nodes += side_nodes.size();
791 : }
792 1139 : }
793 :
794 1207 : others_added_side_nodes.resize(n_proc);
795 1207 : _mesh.comm().allgather(local_added_side_nodes,
796 : others_added_side_nodes);
797 :
798 1207 : added_side_nodes = std::accumulate(others_added_side_nodes.begin(),
799 : others_added_side_nodes.end(), 0,
800 : std::plus<>());
801 :
802 1207 : _mesh.comm().allgather(n_local_nodes, true_node_offsets);
803 11781 : for (auto p : make_range(n_proc-1))
804 10642 : true_node_offsets[p+1] += true_node_offsets[p];
805 34 : libmesh_assert_equal_to(true_node_offsets[n_proc-1], _mesh.max_node_id());
806 :
807 : // For nodes that exist in the mesh, we just need an offset to
808 : // tell where to put their solutions.
809 1207 : added_node_offsets[0] = 0;
810 11781 : for (auto p : make_range(n_proc-1))
811 10574 : added_node_offsets[p+1] =
812 10642 : added_node_offsets[p] + others_added_side_nodes[p];
813 :
814 : // For added side nodes, we need to fill a map. Start after all
815 : // the true node for our pid plus all the side nodes for
816 : // previous pids
817 1241 : dof_id_type node_counter = true_node_offsets[my_pid];
818 6494 : for (auto p : make_range(my_pid))
819 5304 : node_counter += others_added_side_nodes[p];
820 :
821 : // One loop to figure out whose added side nodes get which index
822 4492 : for (const auto & elem : _mesh.active_local_element_ptr_range())
823 : {
824 6240 : for (auto s : elem->side_index_range())
825 : {
826 4992 : if (redundant_added_side(*elem,s))
827 1344 : continue;
828 :
829 : const std::vector<unsigned int> side_nodes =
830 3952 : elem->nodes_on_side(s);
831 :
832 27264 : for (auto n : index_range(side_nodes))
833 : discontinuous_node_indices
834 25584 : [std::make_tuple(elem->id(),s,n)] = node_counter++;
835 : }
836 1139 : }
837 : }
838 :
839 : const dof_id_type
840 79685 : n_global_vals = (max_nn + added_side_nodes) * nv,
841 79685 : n_local_vals = (n_local_nodes + local_added_side_nodes) * nv;
842 :
843 : // Create a NumericVector to hold the parallel solution
844 79685 : std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
845 2444 : NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
846 79685 : parallel_soln.init(n_global_vals, n_local_vals, false, PARALLEL);
847 :
848 : // Create a NumericVector to hold the "repeat_count" for each node - this is essentially
849 : // the number of elements contributing to that node's value
850 82129 : std::unique_ptr<NumericVector<Number>> repeat_count_ptr = NumericVector<Number>::build(_communicator);
851 2444 : NumericVector<Number> & repeat_count = *repeat_count_ptr;
852 79685 : repeat_count.init(n_global_vals, n_local_vals, false, PARALLEL);
853 :
854 79685 : repeat_count.close();
855 :
856 2444 : unsigned int var_num=0;
857 :
858 : // For each system in this EquationSystems object,
859 : // update the global solution and if we are on processor 0,
860 : // loop over the elements and build the nodal solution
861 : // from the element solution. Then insert this nodal solution
862 : // into the vector passed to build_solution_vector.
863 163014 : for (const auto & [sys_name, sys_ptr] : _systems)
864 : {
865 : // Check current system is listed in system_names, and skip pos if not
866 5136 : bool use_current_system = (system_names == nullptr);
867 83329 : if (!use_current_system)
868 352 : use_current_system = system_names->count(sys_name);
869 83329 : if (!use_current_system || sys_ptr->hide_output())
870 1622 : continue;
871 :
872 2528 : const System & system = *sys_ptr;
873 81707 : const unsigned int nv_sys = system.n_vars();
874 5056 : const unsigned int sys_num = system.number();
875 :
876 : //Could this be replaced by a/some convenience methods?[PB]
877 2528 : unsigned int n_scalar_vars = 0;
878 2528 : unsigned int n_vector_vars = 0;
879 195360 : for (auto vn : make_range(sys_ptr->n_vars()))
880 : {
881 113653 : if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
882 6455 : n_vector_vars++;
883 : else
884 107198 : n_scalar_vars++;
885 : }
886 :
887 : // Here, we're assuming the number of vector components is the same
888 : // as the mesh spatial dimension.
889 81707 : unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
890 :
891 : // Update the current_local_solution
892 : {
893 2528 : System & non_const_sys = const_cast<System &>(system);
894 : // We used to simply call non_const_sys.solution->close()
895 : // here, but that is not allowed when the solution vector is
896 : // locked read-only, for example when printing the solution
897 : // during the middle of a solve... So try to be a bit
898 : // more careful about calling close() unnecessarily.
899 2528 : libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
900 81707 : if (!non_const_sys.solution->closed())
901 0 : non_const_sys.solution->close();
902 81707 : non_const_sys.update();
903 : }
904 :
905 2528 : NumericVector<Number> & sys_soln(*system.current_local_solution);
906 :
907 2528 : const DofMap & dof_map = system.get_dof_map();
908 :
909 5056 : std::vector<Number> elem_soln; // The finite element solution
910 5056 : std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
911 2528 : std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
912 :
913 2528 : unsigned var_inc = 0;
914 195360 : for (unsigned int var=0; var<nv_sys; var++)
915 : {
916 113653 : const FEType & fe_type = system.variable_type(var);
917 113653 : const Variable & var_description = system.variable(var);
918 113653 : unsigned int n_vec_dim = FEInterface::n_vec_dim( sys_ptr->get_mesh(), fe_type );
919 113653 : const bool add_p_level = fe_type.p_refinement;
920 :
921 28860254 : for (const auto & elem : _mesh.active_local_element_ptr_range())
922 : {
923 15738095 : if (var_description.active_on_subdomain(elem->subdomain_id()))
924 : {
925 15732563 : dof_map.dof_indices (elem, dof_indices, var);
926 15732563 : sys_soln.get(dof_indices, elem_soln);
927 :
928 15732563 : FEInterface::nodal_soln (elem->dim(),
929 : fe_type,
930 : elem,
931 : elem_soln,
932 : nodal_soln,
933 : add_p_level,
934 : n_vec_dim);
935 :
936 : // infinite elements should be skipped...
937 3404559 : if (!elem->infinite())
938 : {
939 1419394 : libmesh_assert_equal_to (nodal_soln.size(), n_vec_dim*elem->n_nodes());
940 :
941 101394896 : for (auto n : elem->node_index_range())
942 : {
943 84242939 : const Node & node = elem->node_ref(n);
944 :
945 : const dof_id_type node_idx =
946 84242939 : nv * node_id_to_vec_id(node.id());
947 :
948 183952975 : for (unsigned int d=0; d < n_vec_dim; d++)
949 : {
950 : // For vector-valued elements, all components are in nodal_soln. For each
951 : // node, the components are stored in order, i.e. node_0 -> s0_x, s0_y, s0_z
952 108010721 : parallel_soln.add(node_idx + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
953 :
954 : // Increment the repeat count for this position
955 99710036 : repeat_count.add(node_idx + (var_inc+d + var_num), 1);
956 : }
957 : }
958 :
959 15732563 : if (add_sides)
960 : {
961 9020 : for (auto s : elem->side_index_range())
962 : {
963 7200 : if (redundant_added_side(*elem,s))
964 1782 : continue;
965 :
966 : // Compute the FE solution at all the
967 : // side nodes
968 : FEInterface::side_nodal_soln
969 5256 : (fe_type, elem, s, elem_soln,
970 : nodal_soln, add_p_level, n_vec_dim);
971 :
972 : #ifdef DEBUG
973 : const std::vector<unsigned int> side_nodes =
974 876 : elem->nodes_on_side(s);
975 :
976 438 : libmesh_assert_equal_to
977 : (nodal_soln.size(),
978 : side_nodes.size());
979 : #endif
980 :
981 38736 : for (auto n : index_range(nodal_soln))
982 : {
983 : // Retrieve index into global solution vector.
984 : std::size_t node_index =
985 36270 : nv * libmesh_map_find(discontinuous_node_indices,
986 : std::make_tuple(elem->id(), s, n));
987 :
988 66960 : for (unsigned int d=0; d < n_vec_dim; d++)
989 : {
990 36270 : parallel_soln.add(node_index + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
991 33480 : repeat_count.add(node_index + (var_inc+d + var_num), 1);
992 : }
993 : }
994 : }
995 : }
996 : }
997 : }
998 : else // If this variable doesn't exist on this subdomain we have to still increment repeat_count so that we won't divide by 0 later:
999 100536 : for (auto n : elem->node_index_range())
1000 : {
1001 95004 : const Node & node = elem->node_ref(n);
1002 : // Only do this if this variable has NO DoFs at
1003 : // this node... it might have some from an
1004 : // adjoining element...
1005 95004 : if (!node.n_dofs(sys_num, var))
1006 : {
1007 : const dof_id_type node_idx =
1008 78180 : nv * node_id_to_vec_id(node.id());
1009 :
1010 156360 : for (unsigned int d=0; d < n_vec_dim; d++)
1011 78180 : repeat_count.add(node_idx + (var_inc+d + var_num), 1);
1012 : }
1013 : }
1014 :
1015 106589 : } // end loop over elements
1016 113653 : var_inc += n_vec_dim;
1017 : } // end loop on variables in this system
1018 :
1019 81707 : var_num += nv_sys_split;
1020 : } // end loop over systems
1021 :
1022 : // Sum the nodal solution values and repeat counts.
1023 79685 : parallel_soln.close();
1024 79685 : repeat_count.close();
1025 :
1026 : // If there were gaps in the node numbering, there will be
1027 : // corresponding zeros in the parallel_soln and repeat_count
1028 : // vectors. We need to set those repeat_count entries to 1
1029 : // in order to avoid dividing by zero.
1030 79685 : if (n_gaps)
1031 : {
1032 0 : for (numeric_index_type i=repeat_count.first_local_index();
1033 0 : i<repeat_count.last_local_index(); ++i)
1034 : {
1035 : // repeat_count entries are integral values but let's avoid a
1036 : // direct floating point comparison with 0 just in case some
1037 : // roundoff noise crept in during vector assembly?
1038 0 : if (std::abs(repeat_count(i)) < TOLERANCE)
1039 0 : repeat_count.set(i, 1.);
1040 : }
1041 :
1042 : // Make sure the repeat_count vector is up-to-date on all
1043 : // processors.
1044 0 : repeat_count.close();
1045 : }
1046 :
1047 : // Divide to get the average value at the nodes
1048 79685 : parallel_soln /= repeat_count;
1049 :
1050 82129 : return parallel_soln_ptr;
1051 74797 : }
1052 :
1053 :
1054 :
1055 79685 : void EquationSystems::build_solution_vector (std::vector<Number> & soln,
1056 : const std::set<std::string> * system_names,
1057 : bool add_sides) const
1058 : {
1059 4888 : LOG_SCOPE("build_solution_vector()", "EquationSystems");
1060 :
1061 : // Call the parallel implementation
1062 : std::unique_ptr<NumericVector<Number>> parallel_soln =
1063 82129 : this->build_parallel_solution_vector(system_names, add_sides);
1064 :
1065 : // Localize the NumericVector into the provided std::vector.
1066 79685 : parallel_soln->localize_to_one(soln);
1067 79685 : }
1068 :
1069 :
1070 :
1071 8017 : void EquationSystems::get_vars_active_subdomains(const std::vector<std::string> & names,
1072 : std::vector<std::set<subdomain_id_type>> & vars_active_subdomains) const
1073 : {
1074 8017 : vars_active_subdomains.clear();
1075 8243 : vars_active_subdomains.resize(names.size());
1076 :
1077 16034 : for (const auto & pr : _systems)
1078 : {
1079 226 : const auto & sys_ptr = pr.second;
1080 17092 : for (auto vn : make_range(sys_ptr->n_vars()))
1081 : {
1082 9075 : const std::string & var_name = sys_ptr->variable_name(vn);
1083 :
1084 9075 : auto names_it = std::find(names.begin(), names.end(), var_name);
1085 9331 : if(names_it != names.end())
1086 : {
1087 8375 : const Variable & variable = sys_ptr->variable(vn);
1088 236 : const std::set<subdomain_id_type> & active_subdomains = variable.active_subdomains();
1089 8611 : vars_active_subdomains[std::distance(names.begin(), names_it)] = active_subdomains;
1090 : }
1091 : }
1092 : }
1093 8017 : }
1094 :
1095 :
1096 :
1097 : void
1098 352 : EquationSystems::build_elemental_solution_vector (std::vector<Number> & soln,
1099 : std::vector<std::string> & names) const
1100 : {
1101 : // Call the parallel version of this function
1102 : std::unique_ptr<NumericVector<Number>> parallel_soln =
1103 362 : this->build_parallel_elemental_solution_vector(names);
1104 :
1105 : // Localize into 'soln', provided that parallel_soln is not empty.
1106 : // Note: parallel_soln will be empty in the event that none of the
1107 : // input names were elemental data variables, or there were simply none of these in
1108 : // the EquationSystems object.
1109 10 : soln.clear();
1110 352 : if (parallel_soln)
1111 352 : parallel_soln->localize_to_one(soln);
1112 352 : }
1113 :
1114 : std::vector<std::pair<unsigned int, unsigned int>>
1115 284 : EquationSystems::find_variable_numbers
1116 : (std::vector<std::string> & names, const FEType * type, const std::vector<FEType> * types) const
1117 : {
1118 : // Resolve class of type input and assert that at least one of them is null
1119 8 : libmesh_assert_msg(!type || !types,
1120 : "Input 'type', 'types', or neither in find_variable_numbers, but not both.");
1121 :
1122 8 : std::vector<FEType> type_filter;
1123 284 : if (type)
1124 0 : type_filter.push_back(*type);
1125 284 : else if (types)
1126 0 : type_filter = *types;
1127 :
1128 : return this->find_variable_numbers_by_predicate
1129 : (names,
1130 426 : [&type_filter](const FEType & var_type)
1131 : {
1132 426 : return type_filter.empty() ||
1133 12 : std::find(type_filter.begin(), type_filter.end(), var_type) != type_filter.end();
1134 568 : });
1135 : }
1136 :
1137 :
1138 :
1139 : std::vector<std::pair<unsigned int, unsigned int>>
1140 7947 : EquationSystems::find_elemental_data_variable_numbers (std::vector<std::string> & names) const
1141 : {
1142 : return this->find_variable_numbers_by_predicate
1143 15670 : (names, EquationSystems::is_elemental_data_fe_type);
1144 : }
1145 :
1146 :
1147 :
1148 : std::vector<std::pair<unsigned int, unsigned int>>
1149 8231 : EquationSystems::find_variable_numbers_by_predicate
1150 : (std::vector<std::string> & names,
1151 : const std::function<bool(const FEType &)> & type_filter) const
1152 : {
1153 : // This function must be run on all processors at once
1154 232 : parallel_object_only();
1155 :
1156 232 : libmesh_assert (this->n_systems());
1157 :
1158 : // Store a copy of the valid variable names, if any. The names vector will be repopulated with any
1159 : // valid names (or all if 'is_names_empty') in the system that passes through the type filter. If
1160 : // the variable is a vector, its name will be decomposed into its separate components in
1161 : // accordance with build_variable_names().
1162 8695 : std::vector<std::string> name_filter = names;
1163 232 : bool is_names_empty = name_filter.empty();
1164 232 : names.clear();
1165 :
1166 : // initialize convenience variables
1167 232 : FEType var_type;
1168 464 : std::string name;
1169 :
1170 33156 : const std::vector<std::string> component_suffix = {"_x", "_y", "_z"};
1171 8231 : unsigned int dim = _mesh.spatial_dimension();
1172 8231 : libmesh_error_msg_if(dim > 3, "Invalid dim in find_variable_numbers");
1173 :
1174 : // Now filter through the variables in each system and store the system index and their index
1175 : // within that system. This way, we know where to find their data even after we sort them.
1176 464 : std::vector<std::pair<unsigned int, unsigned int>> var_nums;
1177 :
1178 16462 : for (const auto & pr : _systems)
1179 : {
1180 232 : const System & system = *(pr.second);
1181 :
1182 17452 : for (auto var : make_range(system.n_vars()))
1183 : {
1184 : // apply the type filter
1185 9221 : var_type = system.variable_type(var);
1186 9221 : if (!type_filter(var_type))
1187 0 : continue;
1188 :
1189 : // apply the name filter (note that all variables pass if it is empty)
1190 9221 : if (FEInterface::field_type(var_type) == TYPE_VECTOR)
1191 : {
1192 44 : std::vector<std::string> component_names;
1193 2316 : for (unsigned int comp = 0; comp < dim; ++comp)
1194 : {
1195 1588 : name = system.variable_name(var) + component_suffix[comp];
1196 1568 : if (is_names_empty ||
1197 468 : (std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end()))
1198 1544 : component_names.push_back(name);
1199 : }
1200 :
1201 772 : if (! component_names.empty())
1202 772 : names.insert(names.end(), component_names.begin(), component_names.end());
1203 : else
1204 0 : continue;
1205 728 : }
1206 : else /*scalar-valued variable*/
1207 : {
1208 8449 : name = system.variable_name(var);
1209 8465 : if (is_names_empty ||
1210 522 : (std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end()))
1211 8449 : names.push_back(name);
1212 : else
1213 0 : continue;
1214 : }
1215 :
1216 : // if the variable made it through both filters get its system indices
1217 9221 : var_nums.emplace_back(system.number(), var);
1218 : }
1219 : }
1220 :
1221 : // Sort the var_nums vector pairs alphabetically based on the variable name
1222 8695 : std::vector<unsigned int> sort_index(var_nums.size());
1223 232 : std::iota(sort_index.begin(), sort_index.end(), 0);
1224 8231 : std::sort(sort_index.begin(), sort_index.end(),
1225 1980 : [&](const unsigned int & lhs, const unsigned int & rhs)
1226 2036 : {return this->get_system(var_nums[lhs].first).variable_name(var_nums[lhs].second) <
1227 2092 : this->get_system(var_nums[rhs].first).variable_name(var_nums[rhs].second);});
1228 :
1229 8463 : std::vector<std::pair<unsigned int, unsigned int>> var_nums_sorted(var_nums.size());
1230 17452 : for (auto i : index_range(var_nums_sorted))
1231 : {
1232 9741 : var_nums_sorted[i].first = var_nums[sort_index[i]].first;
1233 9481 : var_nums_sorted[i].second = var_nums[sort_index[i]].second;
1234 : }
1235 :
1236 : // Also sort the names vector
1237 8231 : std::sort(names.begin(), names.end());
1238 :
1239 : // Return the sorted vector pairs
1240 8463 : return var_nums_sorted;
1241 31068 : }
1242 :
1243 :
1244 : std::unique_ptr<NumericVector<Number>>
1245 352 : EquationSystems::build_parallel_elemental_solution_vector (std::vector<std::string> & names) const
1246 : {
1247 : // Filter any names that aren't elemental variables and get the system indices for those that are.
1248 : // Note that it's probably fine if the names vector is empty since we'll still filter out all
1249 : // non-elemental-data variables. If there are none, then nothing is output here.
1250 : std::vector<std::pair<unsigned int, unsigned int>> var_nums =
1251 362 : this->find_elemental_data_variable_numbers(names);
1252 :
1253 20 : const std::size_t nv = names.size(); /*total number of vars including vector components*/
1254 352 : const dof_id_type ne = _mesh.n_elem();
1255 10 : libmesh_assert_equal_to (ne, _mesh.max_elem_id());
1256 :
1257 : // If there are no variables to write out don't do anything...
1258 352 : if (!nv)
1259 0 : return std::unique_ptr<NumericVector<Number>>(nullptr);
1260 :
1261 : // We can handle the case where there are nullptrs in the Elem vector
1262 : // by just having extra zeros in the solution vector.
1263 352 : numeric_index_type parallel_soln_global_size = ne*nv;
1264 :
1265 352 : numeric_index_type div = parallel_soln_global_size / this->n_processors();
1266 352 : numeric_index_type mod = parallel_soln_global_size % this->n_processors();
1267 :
1268 : // Initialize all processors to the average size.
1269 10 : numeric_index_type parallel_soln_local_size = div;
1270 :
1271 : // The first "mod" processors get an extra entry.
1272 352 : if (this->processor_id() < mod)
1273 120 : parallel_soln_local_size = div+1;
1274 :
1275 : // Create a NumericVector to hold the parallel solution
1276 362 : std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
1277 10 : NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
1278 352 : parallel_soln.init(parallel_soln_global_size,
1279 : parallel_soln_local_size,
1280 : /*fast=*/false,
1281 20 : /*ParallelType=*/PARALLEL);
1282 :
1283 10 : unsigned int sys_ctr = 0;
1284 10 : unsigned int var_ctr = 0;
1285 1128 : for (auto i : index_range(var_nums))
1286 : {
1287 798 : std::pair<unsigned int, unsigned int> var_num = var_nums[i];
1288 22 : const System & system = this->get_system(var_num.first);
1289 :
1290 : // Update the current_local_solution if necessary
1291 776 : if (sys_ctr != var_num.first)
1292 : {
1293 0 : System & non_const_sys = const_cast<System &>(system);
1294 : // We used to simply call non_const_sys.solution->close()
1295 : // here, but that is not allowed when the solution vector is
1296 : // locked read-only, for example when printing the solution
1297 : // during during the middle of a solve... So try to be a bit
1298 : // more careful about calling close() unnecessarily.
1299 0 : libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
1300 0 : if (!non_const_sys.solution->closed())
1301 0 : non_const_sys.solution->close();
1302 0 : non_const_sys.update();
1303 0 : sys_ctr = var_num.first;
1304 : }
1305 :
1306 22 : NumericVector<Number> & sys_soln(*system.current_local_solution);
1307 :
1308 22 : const unsigned int var = var_num.second;
1309 :
1310 776 : const Variable & variable = system.variable(var);
1311 22 : const DofMap & dof_map = system.get_dof_map();
1312 :
1313 : // We need to check if the elemental data variable is a scalar or a vector and set the number of
1314 : // components for the latter as per es.find_variable_numbers().
1315 : // Even for the case where a variable is not active on any subdomain belonging to the
1316 : // processor, we still need to know this number to update 'var_ctr'.
1317 776 : const auto & var_type = system.variable_type(var);
1318 : const unsigned int n_comps =
1319 776 : (FEInterface::field_type(var_type) == TYPE_VECTOR) ?
1320 350 : FEInterface::n_vec_dim(_mesh, var_type) : 1;
1321 :
1322 : // Loop over all elements in the mesh and index all components of the variable if it's active
1323 : Threads::parallel_for
1324 776 : (_mesh.active_local_element_stored_range(),
1325 1464 : [&dof_map, &variable, ne, var, var_ctr, n_comps,
1326 29456 : ¶llel_soln, &sys_soln](const ConstElemRange & range)
1327 : {
1328 : // The DOF indices for the finite element
1329 44 : std::vector<dof_id_type> dof_indices;
1330 :
1331 4295 : for (const Elem * elem : range)
1332 : {
1333 3519 : if (variable.active_on_subdomain(elem->subdomain_id()))
1334 : {
1335 3519 : dof_map.dof_indices(elem, dof_indices, var);
1336 :
1337 : // The number of DOF components needs to be equal to the expected number so that we know
1338 : // where to store data to correctly correspond to variable names.
1339 315 : libmesh_assert_equal_to(dof_indices.size(), n_comps);
1340 :
1341 9909 : for (unsigned int comp = 0; comp < n_comps; comp++)
1342 6966 : parallel_soln.set(ne * (var_ctr + comp) + elem->id(), sys_soln(dof_indices[comp]));
1343 : }
1344 : }
1345 776 : });
1346 :
1347 776 : var_ctr += n_comps;
1348 : } // end loop over var_nums
1349 :
1350 : // NOTE: number of output names might not be equal to the number passed to this function. Any that
1351 : // aren't elemental data variables have been filtered out (see
1352 : // EquationSystems::find_variable_numbers).
1353 : //
1354 : // But, if everything is accounted for properly, then names.size() == var_ctr
1355 10 : libmesh_assert_equal_to(names.size(), var_ctr);
1356 :
1357 352 : parallel_soln.close();
1358 10 : return parallel_soln_ptr;
1359 332 : }
1360 :
1361 :
1362 :
1363 : void
1364 5507 : EquationSystems::build_discontinuous_solution_vector
1365 : (std::vector<Number> & soln,
1366 : const std::set<std::string> * system_names,
1367 : const std::vector<std::string> * var_names,
1368 : bool vertices_only,
1369 : bool add_sides) const
1370 : {
1371 460 : LOG_SCOPE("build_discontinuous_solution_vector()", "EquationSystems");
1372 :
1373 230 : libmesh_assert (this->n_systems());
1374 :
1375 22258 : const std::vector<std::string> component_suffix = {"_x", "_y", "_z"};
1376 : const auto requested_components =
1377 7337 : [this, var_names, &component_suffix](const System & system,
1378 1899 : const unsigned int var)
1379 : {
1380 405 : std::vector<unsigned int> components;
1381 :
1382 8147 : const std::string & var_name = system.variable_name(var);
1383 8147 : const FEType & fe_type = system.variable_type(var);
1384 8147 : const unsigned int n_vec_dim = FEInterface::n_vec_dim(_mesh, fe_type);
1385 :
1386 8147 : if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
1387 : {
1388 84 : libmesh_error_msg_if(n_vec_dim > component_suffix.size(),
1389 : "Invalid dim in build_discontinuous_solution_vector");
1390 :
1391 : const bool use_all_components =
1392 87 : (var_names == nullptr) ||
1393 81 : std::count(var_names->begin(), var_names->end(), var_name);
1394 :
1395 81 : if (n_vec_dim <= 1)
1396 : {
1397 0 : if (use_all_components)
1398 0 : components.push_back(0);
1399 : }
1400 : else
1401 324 : for (auto comp : make_range(n_vec_dim))
1402 : {
1403 261 : const std::string component_name = var_name + component_suffix[comp];
1404 261 : if (use_all_components ||
1405 243 : std::count(var_names->begin(), var_names->end(), component_name))
1406 243 : components.push_back(comp);
1407 : }
1408 : }
1409 8084 : else if (var_names == nullptr ||
1410 243 : std::count(var_names->begin(), var_names->end(), var_name))
1411 8066 : components.push_back(0);
1412 :
1413 8147 : return components;
1414 5507 : };
1415 :
1416 : // Get the number of variables (nv) by counting the number of variables
1417 : // in each system listed in system_names
1418 230 : unsigned int nv = 0;
1419 :
1420 11085 : for (const auto & [sys_name, sys_ptr] : _systems)
1421 : {
1422 : // Check current system is listed in system_names, and skip pos if not
1423 464 : bool use_current_system = (system_names == nullptr);
1424 5578 : if (!use_current_system)
1425 70 : use_current_system = system_names->count(sys_name);
1426 5578 : if (!use_current_system || sys_ptr->hide_output())
1427 0 : continue;
1428 :
1429 : // Loop over all variables in this System and check whether we
1430 : // are supposed to use each one.
1431 12502 : for (auto var_id : make_range(sys_ptr->n_vars()))
1432 7194 : nv += cast_int<unsigned int>(requested_components(*sys_ptr, var_id).size());
1433 : }
1434 :
1435 : // get the total "weight" - the number of nodal values to write for
1436 : // each variable.
1437 230 : unsigned int tw=0;
1438 408318 : for (const auto & elem : _mesh.active_element_ptr_range())
1439 : {
1440 222787 : tw += vertices_only ? elem->n_vertices() : elem->n_nodes();
1441 :
1442 222787 : if (add_sides)
1443 : {
1444 25064 : for (auto s : elem->side_index_range())
1445 : {
1446 20400 : if (redundant_added_side(*elem,s))
1447 5460 : continue;
1448 :
1449 : const std::vector<unsigned int> side_nodes =
1450 15548 : elem->nodes_on_side(s);
1451 :
1452 14940 : if (!vertices_only)
1453 15548 : tw += side_nodes.size();
1454 : else
1455 0 : for (auto n : index_range(side_nodes))
1456 0 : if (elem->is_vertex(side_nodes[n]))
1457 0 : ++tw;
1458 : }
1459 : }
1460 5047 : }
1461 :
1462 : // Only if we are on processor zero, allocate the storage
1463 : // to hold (number_of_nodes)*(number_of_variables) entries.
1464 5737 : if (_mesh.processor_id() == 0)
1465 986 : soln.resize(tw*nv);
1466 :
1467 460 : std::vector<Number> sys_soln;
1468 :
1469 : // Keep track of the variable "offset". This is used for indexing
1470 : // into the global solution vector.
1471 230 : unsigned int var_offset = 0;
1472 :
1473 : // For each system in this EquationSystems object,
1474 : // update the global solution and if we are on processor 0,
1475 : // loop over the elements and build the nodal solution
1476 : // from the element solution. Then insert this nodal solution
1477 : // into the vector passed to build_solution_vector.
1478 11085 : for (const auto & [sys_name, system] : _systems)
1479 : {
1480 : // Check current system is listed in system_names, and skip pos if not
1481 464 : bool use_current_system = (system_names == nullptr);
1482 5578 : if (!use_current_system)
1483 70 : use_current_system = system_names->count(sys_name);
1484 5578 : if (!use_current_system || system->hide_output())
1485 0 : continue;
1486 :
1487 5578 : const unsigned int nv_sys = system->n_vars();
1488 232 : const auto & dof_map = system->get_dof_map();
1489 :
1490 5578 : system->update_global_solution (sys_soln, 0);
1491 :
1492 : // Keep track of the number of vars actually written.
1493 232 : unsigned int n_vars_written_current_system = 0;
1494 :
1495 5810 : if (_mesh.processor_id() == 0)
1496 : {
1497 232 : std::vector<Number> soln_coeffs; // The finite element solution coeffs
1498 232 : std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
1499 232 : std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
1500 :
1501 : // For each variable, determine if we are supposed to
1502 : // write it, then loop over the active elements, compute
1503 : // the nodal_soln and store it to the "soln" vector. We
1504 : // store zeros for subdomain-restricted variables on
1505 : // elements where they are not active.
1506 2221 : for (auto var : make_range(nv_sys))
1507 : {
1508 : const std::vector<unsigned int> components_to_write =
1509 1223 : requested_components(*system, var);
1510 :
1511 : // If we aren't supposed to write this var, go to the
1512 : // next loop iteration.
1513 1223 : if (components_to_write.empty())
1514 0 : continue;
1515 :
1516 1223 : const FEType & fe_type = system->variable_type(var);
1517 1223 : const Variable & var_description = system->variable(var);
1518 1223 : const bool add_p_level = fe_type.p_refinement;
1519 1223 : const unsigned int n_vec_dim = FEInterface::n_vec_dim(_mesh, fe_type);
1520 :
1521 135 : unsigned int nn=0;
1522 :
1523 232321 : for (auto & elem : _mesh.active_element_ptr_range())
1524 : {
1525 136299 : if (var_description.active_on_subdomain(elem->subdomain_id()))
1526 : {
1527 136299 : dof_map.dof_indices (elem, dof_indices, var);
1528 :
1529 157593 : soln_coeffs.resize(dof_indices.size());
1530 :
1531 977966 : for (auto i : index_range(dof_indices))
1532 1082375 : soln_coeffs[i] = sys_soln[dof_indices[i]];
1533 :
1534 : // Compute the FE solution at all the nodes, but
1535 : // only use the first n_vertices() entries if
1536 : // vertices_only == true.
1537 136299 : FEInterface::nodal_soln (elem->dim(),
1538 : fe_type,
1539 : elem,
1540 : soln_coeffs,
1541 : nodal_soln,
1542 : add_p_level,
1543 : n_vec_dim);
1544 :
1545 : // infinite elements should be skipped...
1546 61566 : if (!elem->infinite())
1547 : {
1548 21294 : libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes()*n_vec_dim);
1549 :
1550 : const unsigned int n_vals =
1551 136299 : vertices_only ? elem->n_vertices() : elem->n_nodes();
1552 :
1553 1065615 : for (auto n : make_range(n_vals))
1554 : {
1555 : // Compute index into global solution vector.
1556 : std::size_t index =
1557 929316 : nv * (nn++) + (n_vars_written_current_system + var_offset);
1558 :
1559 1860040 : for (auto component_index : index_range(components_to_write))
1560 1070656 : soln[index + component_index] +=
1561 1070656 : nodal_soln[n_vec_dim*n + components_to_write[component_index]];
1562 : }
1563 : }
1564 : }
1565 : else
1566 0 : nn += vertices_only ? elem->n_vertices() : elem->n_nodes();
1567 953 : } // end loop over active elements writing interiors
1568 :
1569 : // Loop writing "fake" sides, if requested
1570 1223 : if (add_sides)
1571 : {
1572 : // We don't build discontinuous solution vectors in
1573 : // parallel yet, but we'll do ordering of fake side
1574 : // values as if we did, for consistency with the
1575 : // parallel continuous ordering and for future
1576 : // compatibility.
1577 : std::vector<std::vector<const Elem *>>
1578 375 : elems_by_pid(_mesh.n_processors());
1579 :
1580 3655 : for (const auto & elem : _mesh.active_element_ptr_range())
1581 2070 : elems_by_pid[elem->processor_id()].push_back(elem);
1582 :
1583 2075 : for (auto p : index_range(elems_by_pid))
1584 3455 : for (const Elem * elem : elems_by_pid[p])
1585 : {
1586 1680 : if (var_description.active_on_subdomain(elem->subdomain_id()))
1587 : {
1588 1680 : dof_map.dof_indices (elem, dof_indices, var);
1589 :
1590 1820 : soln_coeffs.resize(dof_indices.size());
1591 :
1592 32064 : for (auto i : index_range(dof_indices))
1593 35448 : soln_coeffs[i] = sys_soln[dof_indices[i]];
1594 :
1595 8880 : for (auto s : elem->side_index_range())
1596 : {
1597 7200 : if (redundant_added_side(*elem,s))
1598 1944 : continue;
1599 :
1600 : const std::vector<unsigned int> side_nodes =
1601 5694 : elem->nodes_on_side(s);
1602 :
1603 : // Compute the FE solution at all the
1604 : // side nodes, but only use those for
1605 : // which is_vertex() == true if
1606 : // vertices_only == true.
1607 : FEInterface::side_nodal_soln
1608 5256 : (fe_type, elem, s, soln_coeffs,
1609 : nodal_soln, add_p_level,
1610 : n_vec_dim);
1611 :
1612 438 : libmesh_assert_equal_to
1613 : (nodal_soln.size(),
1614 : side_nodes.size()*n_vec_dim);
1615 :
1616 : // If we don't have a continuous FE
1617 : // then we want to average between
1618 : // sides, at least in the equal-level
1619 : // case where it's easy. This is
1620 : // analogous to our repeat_count
1621 : // behavior elsewhere.
1622 : const FEContinuity cont =
1623 5256 : FEInterface::get_continuity(fe_type);
1624 876 : const Elem * const neigh = elem->neighbor_ptr(s);
1625 :
1626 5256 : if ((cont == DISCONTINUOUS || cont == H_CURL || cont == H_DIV) &&
1627 0 : neigh &&
1628 5694 : neigh->level() == elem->level() &&
1629 0 : var_description.active_on_subdomain(neigh->subdomain_id()))
1630 : {
1631 0 : std::vector<dof_id_type> neigh_indices;
1632 0 : dof_map.dof_indices (neigh, neigh_indices, var);
1633 0 : std::vector<Number> neigh_coeffs(neigh_indices.size());
1634 :
1635 0 : for (auto i : index_range(neigh_indices))
1636 0 : neigh_coeffs[i] = sys_soln[neigh_indices[i]];
1637 :
1638 : const unsigned int s_neigh =
1639 0 : neigh->which_neighbor_am_i(elem);
1640 0 : std::vector<Number> neigh_soln;
1641 : FEInterface::side_nodal_soln
1642 0 : (fe_type, neigh, s_neigh,
1643 : neigh_coeffs, neigh_soln, add_p_level,
1644 : n_vec_dim);
1645 :
1646 : const std::vector<unsigned int> neigh_nodes =
1647 0 : neigh->nodes_on_side(s_neigh);
1648 0 : for (auto n : index_range(side_nodes))
1649 0 : for (auto neigh_n : index_range(neigh_nodes))
1650 0 : if (neigh->node_ptr(neigh_nodes[neigh_n])
1651 0 : == elem->node_ptr(side_nodes[n]))
1652 0 : for (auto comp : make_range(n_vec_dim))
1653 : {
1654 0 : const auto nodal_index = n_vec_dim*n + comp;
1655 0 : nodal_soln[nodal_index] +=
1656 0 : neigh_soln[n_vec_dim*neigh_n + comp];
1657 0 : nodal_soln[nodal_index] /= 2;
1658 : }
1659 : }
1660 :
1661 38736 : for (auto n : index_range(side_nodes))
1662 : {
1663 33480 : if (vertices_only &&
1664 0 : !elem->is_vertex(n))
1665 0 : continue;
1666 :
1667 : // Compute index into global solution vector.
1668 : std::size_t index =
1669 33480 : nv * (nn++) + (n_vars_written_current_system + var_offset);
1670 :
1671 66960 : for (auto component_index : index_range(components_to_write))
1672 36270 : soln[index + component_index] +=
1673 36270 : nodal_soln[n_vec_dim*n + components_to_write[component_index]];
1674 : }
1675 : }
1676 : }
1677 : else
1678 : {
1679 0 : nn += vertices_only ? elem->n_vertices() : elem->n_nodes();
1680 :
1681 0 : for (auto s : elem->side_index_range())
1682 : {
1683 0 : if (redundant_added_side(*elem,s))
1684 0 : continue;
1685 :
1686 : const std::vector<unsigned int> side_nodes =
1687 0 : elem->nodes_on_side(s);
1688 :
1689 0 : for (auto n : index_range(side_nodes))
1690 : {
1691 0 : if (vertices_only &&
1692 0 : !elem->is_vertex(n))
1693 0 : continue;
1694 0 : nn++;
1695 : }
1696 : }
1697 : }
1698 : } // end loop over active elements, writing "fake" sides
1699 250 : }
1700 : // If we made it here, we actually wrote a variable, so increment
1701 : // the number of variables actually written for the current system.
1702 1358 : n_vars_written_current_system += cast_int<unsigned int>(components_to_write.size());
1703 :
1704 : } // end loop over vars
1705 : } // end if proc 0
1706 :
1707 : // Update offset for next loop iteration.
1708 5578 : var_offset += n_vars_written_current_system;
1709 : } // end loop over systems
1710 15601 : }
1711 :
1712 :
1713 :
1714 188928 : bool EquationSystems::redundant_added_side(const Elem & elem, unsigned int side)
1715 : {
1716 10536 : libmesh_assert(elem.active());
1717 :
1718 21072 : const Elem * neigh = elem.neighbor_ptr(side);
1719 :
1720 : // Write boundary sides.
1721 188928 : if (!neigh)
1722 4860 : return false;
1723 :
1724 : // Write ghost sides in Nemesis
1725 102432 : if (neigh == remote_elem)
1726 0 : return false;
1727 :
1728 : // Don't write a coarser side if a finer side exists
1729 5676 : if (!neigh->active())
1730 0 : return true;
1731 :
1732 : // Don't write a side redundantly from both of the
1733 : // elements sharing it. We'll disambiguate with id().
1734 101376 : return (neigh->id() < elem.id());
1735 : }
1736 :
1737 :
1738 :
1739 0 : bool EquationSystems::compare (const EquationSystems & other_es,
1740 : const Real threshold,
1741 : const bool verbose) const
1742 : {
1743 : // safety check, whether we handle at least the same number
1744 : // of systems
1745 0 : std::vector<bool> os_result;
1746 :
1747 0 : if (this->n_systems() != other_es.n_systems())
1748 : {
1749 0 : if (verbose)
1750 : {
1751 0 : libMesh::out << " Fatal difference. This system handles "
1752 0 : << this->n_systems() << " systems," << std::endl
1753 0 : << " while the other system handles "
1754 0 : << other_es.n_systems()
1755 0 : << " systems." << std::endl
1756 0 : << " Aborting comparison." << std::endl;
1757 : }
1758 0 : return false;
1759 : }
1760 : else
1761 : {
1762 : // start comparing each system
1763 0 : for (const auto & [sys_name, sys_ptr] : _systems)
1764 : {
1765 : // get the other system
1766 0 : const System & other_system = other_es.get_system (sys_name);
1767 :
1768 0 : os_result.push_back (sys_ptr->compare (other_system, threshold, verbose));
1769 :
1770 : }
1771 :
1772 : }
1773 :
1774 :
1775 : // sum up the results
1776 0 : if (os_result.size()==0)
1777 0 : return true;
1778 : else
1779 : {
1780 : bool os_identical;
1781 0 : unsigned int n = 0;
1782 0 : do
1783 : {
1784 0 : os_identical = os_result[n];
1785 0 : n++;
1786 : }
1787 0 : while (os_identical && n<os_result.size());
1788 0 : return os_identical;
1789 : }
1790 : }
1791 :
1792 :
1793 :
1794 16683 : std::string EquationSystems::get_info () const
1795 : {
1796 17607 : std::ostringstream oss;
1797 :
1798 462 : unsigned int n_hidden_sys = 0;
1799 35946 : for (auto & pr : _systems)
1800 19263 : n_hidden_sys += pr.second->hide_output();
1801 :
1802 : oss << " EquationSystems\n"
1803 16221 : << " n_systems()=" << this->n_systems()
1804 36090 : << (n_hidden_sys ? " (hidden: " + std::to_string(n_hidden_sys) + ")" : "")
1805 16683 : << "\n";
1806 :
1807 : // Print the info for the individual systems
1808 35946 : for (const auto & pr : _systems)
1809 19263 : if (!pr.second->hide_output())
1810 35278 : oss << pr.second->get_info();
1811 :
1812 :
1813 : // // Possibly print the parameters
1814 : // if (!this->parameters.empty())
1815 : // {
1816 : // oss << " n_parameters()=" << this->n_parameters() << '\n';
1817 : // oss << " Parameters:\n";
1818 :
1819 : // for (const auto & [key, val] : _parameters)
1820 : // oss << " "
1821 : // << "\""
1822 : // << key
1823 : // << "\""
1824 : // << "="
1825 : // << val
1826 : // << '\n';
1827 : // }
1828 :
1829 17145 : return oss.str();
1830 15759 : }
1831 :
1832 :
1833 :
1834 16683 : void EquationSystems::print_info (std::ostream & os) const
1835 : {
1836 17145 : os << this->get_info()
1837 462 : << std::endl;
1838 16683 : }
1839 :
1840 :
1841 :
1842 0 : std::ostream & operator << (std::ostream & os,
1843 : const EquationSystems & es)
1844 : {
1845 0 : es.print_info(os);
1846 0 : return os;
1847 : }
1848 :
1849 :
1850 :
1851 2700 : unsigned int EquationSystems::n_vars () const
1852 : {
1853 2700 : unsigned int tot=0;
1854 :
1855 5528 : for (const auto & pr : _systems)
1856 2828 : tot += pr.second->n_vars();
1857 :
1858 2700 : return tot;
1859 : }
1860 :
1861 :
1862 :
1863 0 : std::size_t EquationSystems::n_dofs () const
1864 : {
1865 0 : std::size_t tot=0;
1866 :
1867 0 : for (const auto & pr : _systems)
1868 0 : tot += pr.second->n_dofs();
1869 :
1870 0 : return tot;
1871 : }
1872 :
1873 :
1874 :
1875 :
1876 22948 : std::size_t EquationSystems::n_active_dofs () const
1877 : {
1878 730 : std::size_t tot=0;
1879 :
1880 45896 : for (const auto & pr : _systems)
1881 22948 : tot += pr.second->n_active_dofs();
1882 :
1883 22948 : return tot;
1884 : }
1885 :
1886 :
1887 242602 : void EquationSystems::_add_system_to_nodes_and_elems()
1888 : {
1889 : // All the nodes
1890 31115768 : for (auto & node : _mesh.node_ptr_range())
1891 16488544 : node->add_system();
1892 :
1893 : // All the elements
1894 : Threads::parallel_for
1895 242602 : (_mesh.element_stored_range(),
1896 235786 : [](const ElemRange & range)
1897 : {
1898 7317468 : for (Elem * elem : range)
1899 7074824 : elem->add_system();
1900 235786 : });
1901 242602 : }
1902 :
1903 71 : void EquationSystems::_remove_default_ghosting(unsigned int sys_num)
1904 : {
1905 71 : this->get_system(sys_num).get_dof_map().remove_default_ghosting();
1906 71 : }
1907 :
1908 : } // namespace libMesh
|
