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