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