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