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