libMesh
equation_systems.C
Go to the documentation of this file.
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 
53  ParallelObject (m),
54  _mesh (m),
55  _refine_in_reinit(true),
56  _enable_default_ghosting(true)
57 {
58  // Set default parameters
59  this->parameters.set<Real> ("linear solver tolerance") = TOLERANCE * TOLERANCE;
60  this->parameters.set<unsigned int>("linear solver maximum iterations") = 5000;
61 }
62 
63 
64 
66 
67 
68 
70 {
71  // Clear any additional parameters
72  parameters.clear ();
73 
74  // Clear the systems.
75  _systems.clear();
76 }
77 
78 
79 
81 {
82 #ifndef NDEBUG
83  for (auto i : make_range(this->n_systems()))
84  libmesh_assert(!this->get_system(i).is_initialized());
85 #endif
86 
87  this->reinit_mesh();
88 }
89 
90 
91 
93 {
94  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  if (mesh_changed)
99  this->reinit_systems();
100 }
101 
103 {
104  const unsigned int n_sys = this->n_systems();
105 
106  libmesh_assert_not_equal_to (n_sys, 0);
107 
108  // Tell all the \p DofObject entities how many systems
109  // there are.
110  for (auto & node : _mesh.node_ptr_range())
111  node->set_n_systems(n_sys);
112 
113  for (auto & elem : _mesh.element_ptr_range())
114  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  MeshRefinement mesh_refine(_mesh);
121  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  for (auto i : make_range(this->n_systems()))
127  this->get_system(i).reinit_mesh();
128 
129 }
130 
132 {
133  parallel_object_only();
134 
135  const unsigned int n_sys = this->n_systems();
136  libmesh_assert_not_equal_to (n_sys, 0);
137 
138  // And any new systems will need initialization
139  for (unsigned int i=0; i != n_sys; ++i)
140  if (!this->get_system(i).is_initialized())
141  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  for (auto & node : _mesh.node_ptr_range())
153  node->set_n_systems(n_sys);
154 
155  // All the elements
156  for (auto & elem : _mesh.element_ptr_range())
157  elem->set_n_systems(n_sys);
158  }
159 
160  // Localize each system's vectors
161  for (unsigned int i=0; i != n_sys; ++i)
162  this->get_system(i).re_update();
163 
164 #ifdef LIBMESH_ENABLE_AMR
165 
166  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  for (unsigned int i=0; i != n_sys; ++i)
173  {
174  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 
186 
187  // Recreate any user or internal constraints
188  sys.reinit_constraints();
189 
190  // Even if there weren't any constraint changes,
191  // reinit_constraints() did prepare_send_list() for us.
192 
193  sys.prolong_vectors();
194  }
195  mesh_changed = true;
196  }
197 
198  if (this->_refine_in_reinit)
199  {
200  // Don't override any user refinement settings
201  MeshRefinement mesh_refine(_mesh);
202  mesh_refine.face_level_mismatch_limit() = 0; // unlimited
203  mesh_refine.overrefined_boundary_limit() = -1; // unlimited
204  mesh_refine.underrefined_boundary_limit() = -1; // unlimited
205 
206  // Try to coarsen the mesh, then restrict each system's vectors
207  // if necessary
208  if (mesh_refine.coarsen_elements())
209  {
210  for (auto i : make_range(this->n_systems()))
211  {
212  System & sys = this->get_system(i);
214  sys.reinit_constraints();
215 
216  // Even if there weren't any constraint changes,
217  // reinit_constraints() did prepare_send_list() for us.
218 
219  sys.restrict_vectors();
220  }
221  mesh_changed = true;
222  }
223 
224  // Once vectors are all restricted, we can delete
225  // children of coarsened elements
226  if (mesh_changed)
227  this->get_mesh().contract();
228 
229  // Try to refine the mesh, then prolong each system's vectors
230  // if necessary
231  if (mesh_refine.refine_elements())
232  {
233  for (auto i : make_range(this->n_systems()))
234  {
235  System & sys = this->get_system(i);
237  sys.reinit_constraints();
238 
239  // Even if there weren't any constraint changes,
240  // reinit_constraints() did prepare_send_list() for us.
241 
242  sys.prolong_vectors();
243  }
244  mesh_changed = true;
245  }
246  }
247 
248  return mesh_changed;
249 
250 #endif // #ifdef LIBMESH_ENABLE_AMR
251 
252  return false;
253 }
254 
255 
256 
258 {
259  for (auto i : make_range(this->n_systems()))
260  this->get_system(i).reinit();
261 }
262 
263 
264 
266 {
267  // A serial mesh means nothing needs to be done
268  if (_mesh.is_serial())
269  return;
270 
271  const unsigned int n_sys = this->n_systems();
272 
273  libmesh_assert_not_equal_to (n_sys, 0);
274 
275  // Gather the mesh
276  _mesh.allgather();
277 
278  // Tell all the \p DofObject entities how many systems
279  // there are.
280  for (auto & node : _mesh.node_ptr_range())
281  node->set_n_systems(n_sys);
282 
283  for (auto & elem : _mesh.element_ptr_range())
284  elem->set_n_systems(n_sys);
285 
286  // And distribute each system's dofs
287  for (auto i : make_range(this->n_systems()))
288  {
289  System & sys = this->get_system(i);
290  DofMap & dof_map = sys.get_dof_map();
291  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  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 
306 {
307  _enable_default_ghosting = enable;
308  MeshBase &mesh = this->get_mesh();
309 
310  if (enable)
311  mesh.add_ghosting_functor(mesh.default_ghosting());
312  else
313  mesh.remove_ghosting_functor(mesh.default_ghosting());
314 
315  for (auto i : make_range(this->n_systems()))
316  {
317  DofMap & dof_map = this->get_system(i).get_dof_map();
318  if (enable)
319  dof_map.add_default_ghosting();
320  else
321  dof_map.remove_default_ghosting();
322  }
323 }
324 
325 
326 
328 {
329  LOG_SCOPE("update()", "EquationSystems");
330 
331  // Localize each system's vectors
332  for (auto i : make_range(this->n_systems()))
333  this->get_system(i).update();
334 }
335 
336 
337 
338 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  if (_systems.count(name))
346  {
347  return this->get_system(name);
348  }
349  // Build a basic System
350  else if (sys_type == "Basic")
351  this->add_system<System> (name);
352 
353  // Build a Newmark system
354  else if (sys_type == "Newmark")
355  this->add_system<NewmarkSystem> (name);
356 
357  // Build an Explicit system
358  else if ((sys_type == "Explicit"))
359  this->add_system<ExplicitSystem> (name);
360 
361  // Build an Implicit system
362  else if ((sys_type == "Implicit") ||
363  (sys_type == "Steady" ))
364  this->add_system<ImplicitSystem> (name);
365 
366  // build a transient implicit linear system
367  else if ((sys_type == "Transient") ||
368  (sys_type == "TransientImplicit") ||
369  (sys_type == "TransientLinearImplicit"))
370  this->add_system<TransientLinearImplicitSystem> (name);
371 
372  // build a transient implicit nonlinear system
373  else if (sys_type == "TransientNonlinearImplicit")
374  this->add_system<TransientNonlinearImplicitSystem> (name);
375 
376  // build a transient explicit system
377  else if (sys_type == "TransientExplicit")
378  this->add_system<TransientExplicitSystem> (name);
379 
380  // build a linear implicit system
381  else if (sys_type == "LinearImplicit")
382  this->add_system<LinearImplicitSystem> (name);
383 
384  // build a nonlinear implicit system
385  else if (sys_type == "NonlinearImplicit")
386  this->add_system<NonlinearImplicitSystem> (name);
387 
388  // build a Reduced Basis Construction system
389  else if (sys_type == "RBConstruction")
390  this->add_system<RBConstruction> (name);
391 
392  // build a transient Reduced Basis Construction system
393  else if (sys_type == "TransientRBConstruction")
394  this->add_system<TransientRBConstruction> (name);
395 
396 #ifdef LIBMESH_HAVE_SLEPC
397  // build an eigen system
398  else if (sys_type == "Eigen")
399  this->add_system<EigenSystem> (name);
400  else if (sys_type == "TransientEigenSystem")
401  this->add_system<TransientEigenSystem> (name);
402 #endif
403 
404 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
405  // build a frequency system
406  else if (sys_type == "Frequency")
407  this->add_system<FrequencySystem> (name);
408 #endif
409 
410  else
411  libmesh_error_msg("ERROR: Unknown system type: " << sys_type);
412 
413  // Return a reference to the new system
414  //return (*this)(name);
415  return this->get_system(name);
416 }
417 
418 
419 
421 {
422  libmesh_assert (this->n_systems());
423 
424  for (auto i : make_range(this->n_systems()))
425  this->get_system(i).solve();
426 }
427 
428 
429 
431 {
432  libmesh_assert (this->n_systems());
433 
434  for (auto i : make_range(this->n_systems()))
435  this->get_system(i).sensitivity_solve(parameters_in);
436 }
437 
438 
439 
440 void EquationSystems::adjoint_solve (const QoISet & qoi_indices)
441 {
442  libmesh_assert (this->n_systems());
443 
444  for (unsigned int i=this->n_systems(); i != 0; --i)
445  this->get_system(i-1).adjoint_solve(qoi_indices);
446 }
447 
448 
449 
450 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  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  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  unsigned int n_scalar_vars = 0;
466  unsigned int n_vector_vars = 0;
467 
468  for (const auto & [sys_name, sys_ptr] : _systems)
469  {
470  // Check current system is listed in system_names, and skip pos if not
471  bool use_current_system = (system_names == nullptr);
472  if (!use_current_system)
473  use_current_system = system_names->count(sys_name);
474  if (!use_current_system || sys_ptr->hide_output())
475  {
476  for (auto vn : make_range(sys_ptr->n_vars()))
477  seen_names.insert(sys_ptr->variable_name(vn));
478  continue;
479  }
480 
481  for (auto vn : make_range(sys_ptr->n_vars()))
482  {
483  seen_names.insert(sys_ptr->variable_name(vn));
484  if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
485  n_vector_vars++;
486  else
487  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  unsigned int dim = this->get_mesh().spatial_dimension();
494  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  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  var_names.resize( nv );
504  }
505 
506  for (const auto & [sys_name, sys_ptr] : _systems)
507  {
508  // Check current system is listed in system_names, and skip pos if not
509  bool use_current_system = (system_names == nullptr);
510  if (!use_current_system)
511  use_current_system = system_names->count(sys_name);
512  if (!use_current_system || sys_ptr->hide_output())
513  continue;
514 
515  for (auto vn : make_range(sys_ptr->n_vars()))
516  {
517  const std::string & var_name = sys_ptr->variable_name(vn);
518  const FEType & fe_type = sys_ptr->variable_type(vn);
519 
520  unsigned int n_vec_dim = FEInterface::n_vec_dim( sys_ptr->get_mesh(), fe_type);
521 
522  // Filter on the type if requested
523  if (type == nullptr || (type && *type == fe_type))
524  {
525  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
526  {
527  switch(n_vec_dim)
528  {
529  case 0:
530  case 1:
531  var_names[var_num++] = var_name;
532  libmesh_error_msg_if(seen_names.count(var_name) > 1,
533  "Duplicate variable name "+var_name);
534  break;
535  case 2:
536  var_names[var_num++] = var_name+"_x";
537  var_names[var_num++] = var_name+"_y";
538  libmesh_error_msg_if(seen_names.count(var_name+"_x"),
539  "Duplicate variable name "+var_name+"_x");
540  libmesh_error_msg_if(seen_names.count(var_name+"_y"),
541  "Duplicate variable name "+var_name+"_y");
542  break;
543  case 3:
544  var_names[var_num++] = var_name+"_x";
545  var_names[var_num++] = var_name+"_y";
546  var_names[var_num++] = var_name+"_z";
547  libmesh_error_msg_if(seen_names.count(var_name+"_x"),
548  "Duplicate variable name "+var_name+"_x");
549  libmesh_error_msg_if(seen_names.count(var_name+"_y"),
550  "Duplicate variable name "+var_name+"_y");
551  libmesh_error_msg_if(seen_names.count(var_name+"_z"),
552  "Duplicate variable name "+var_name+"_z");
553  break;
554  default:
555  libmesh_error_msg("Invalid dim in build_variable_names");
556  }
557  }
558  else
559  var_names[var_num++] = var_name;
560  }
561  }
562  }
563  // Now resize again in case we filtered any names
564  var_names.resize(var_num);
565 }
566 
567 
568 
569 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  libmesh_not_implemented();
575 }
576 
577 
578 
579 
580 std::unique_ptr<NumericVector<Number>>
581 EquationSystems::build_parallel_solution_vector(const std::set<std::string> * system_names,
582  bool add_sides) const
583 {
584  LOG_SCOPE("build_parallel_solution_vector()", "EquationSystems");
585 
586  // This function must be run on all processors at once
587  parallel_object_only();
588 
589  const unsigned int dim = _mesh.spatial_dimension();
590  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  unsigned int nv = 0;
604 
605  //Could this be replaced by a/some convenience methods?[PB]
606  {
607  unsigned int n_scalar_vars = 0;
608  unsigned int n_vector_vars = 0;
609  for (const auto & [sys_name, sys_ptr] : _systems)
610  {
611  // Check current system is listed in system_names, and skip pos if not
612  bool use_current_system = (system_names == nullptr);
613  if (!use_current_system)
614  use_current_system = system_names->count(sys_name);
615  if (!use_current_system || sys_ptr->hide_output())
616  continue;
617 
618  for (auto vn : make_range(sys_ptr->n_vars()))
619  {
620  if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
621  n_vector_vars++;
622  else
623  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  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  (std::distance(_mesh.local_nodes_begin(),
634  _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  dof_id_type n_total_nodes = n_local_nodes;
640  _mesh.comm().sum(n_total_nodes);
641 
642  const processor_id_type n_proc = _mesh.comm().size();
643  const processor_id_type my_pid = _mesh.comm().rank();
644  const dof_id_type n_gaps = max_nn - n_total_nodes;
645  const dof_id_type gaps_per_processor = n_gaps / n_proc;
646  const dof_id_type remainder_gaps = n_gaps % n_proc;
647 
648  n_local_nodes = n_local_nodes + // Actual nodes
649  gaps_per_processor + // Our even share of gaps
650  (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  dof_id_type local_added_side_nodes = 0,
655  added_side_nodes = 0;
656 
657  // others_added_side_nodes[p]: local_added_side_nodes on rank p
658  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  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  std::vector<dof_id_type> true_node_offsets;
671  // Number of added (fake) nodes on processors [0,p)
672  std::vector<dof_id_type> added_node_offsets;
673 
674  auto node_id_to_vec_id =
675  [&true_node_offsets, &added_node_offsets]
676  (const dof_id_type node_id)
677  {
678  if (true_node_offsets.empty())
679  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  const auto lb = std::upper_bound(true_node_offsets.begin(),
683  true_node_offsets.end(), node_id);
684  libmesh_assert(lb != true_node_offsets.end());
685  const processor_id_type p = lb - true_node_offsets.begin();
686 
687  return node_id + added_node_offsets[p];
688  };
689 
690  if (add_sides)
691  {
692  true_node_offsets.resize(n_proc);
693  added_node_offsets.resize(n_proc);
694 
695  // One loop to count everyone's new side nodes
696  for (const auto & elem : _mesh.active_element_ptr_range())
697  {
698  for (auto s : elem->side_index_range())
699  {
700  if (redundant_added_side(*elem,s))
701  continue;
702 
703  const std::vector<unsigned int> side_nodes =
704  elem->nodes_on_side(s);
705 
706  if (elem->processor_id() == this->processor_id())
707  local_added_side_nodes += side_nodes.size();
708  }
709  }
710 
711  others_added_side_nodes.resize(n_proc);
712  _mesh.comm().allgather(local_added_side_nodes,
713  others_added_side_nodes);
714 
715  added_side_nodes = std::accumulate(others_added_side_nodes.begin(),
716  others_added_side_nodes.end(), 0,
717  std::plus<>());
718 
719  _mesh.comm().allgather(n_local_nodes, true_node_offsets);
720  for (auto p : make_range(n_proc-1))
721  true_node_offsets[p+1] += true_node_offsets[p];
722  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  added_node_offsets[0] = 0;
727  for (auto p : make_range(n_proc-1))
728  added_node_offsets[p+1] =
729  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  dof_id_type node_counter = true_node_offsets[my_pid];
735  for (auto p : make_range(my_pid))
736  node_counter += others_added_side_nodes[p];
737 
738  // One loop to figure out whose added side nodes get which index
739  for (const auto & elem : _mesh.active_local_element_ptr_range())
740  {
741  for (auto s : elem->side_index_range())
742  {
743  if (redundant_added_side(*elem,s))
744  continue;
745 
746  const std::vector<unsigned int> side_nodes =
747  elem->nodes_on_side(s);
748 
749  for (auto n : index_range(side_nodes))
750  discontinuous_node_indices
751  [std::make_tuple(elem->id(),s,n)] = node_counter++;
752  }
753  }
754  }
755 
756  const dof_id_type
757  n_global_vals = (max_nn + added_side_nodes) * nv,
758  n_local_vals = (n_local_nodes + local_added_side_nodes) * nv;
759 
760  // Create a NumericVector to hold the parallel solution
761  std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
762  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
763  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  std::unique_ptr<NumericVector<Number>> repeat_count_ptr = NumericVector<Number>::build(_communicator);
768  NumericVector<Number> & repeat_count = *repeat_count_ptr;
769  repeat_count.init(n_global_vals, n_local_vals, false, PARALLEL);
770 
771  repeat_count.close();
772 
773  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  for (const auto & [sys_name, sys_ptr] : _systems)
781  {
782  // Check current system is listed in system_names, and skip pos if not
783  bool use_current_system = (system_names == nullptr);
784  if (!use_current_system)
785  use_current_system = system_names->count(sys_name);
786  if (!use_current_system || sys_ptr->hide_output())
787  continue;
788 
789  const System & system = *sys_ptr;
790  const unsigned int nv_sys = system.n_vars();
791  const unsigned int sys_num = system.number();
792 
793  //Could this be replaced by a/some convenience methods?[PB]
794  unsigned int n_scalar_vars = 0;
795  unsigned int n_vector_vars = 0;
796  for (auto vn : make_range(sys_ptr->n_vars()))
797  {
798  if (FEInterface::field_type(sys_ptr->variable_type(vn)) == TYPE_VECTOR)
799  n_vector_vars++;
800  else
801  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  unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
807 
808  // Update the current_local_solution
809  {
810  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  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
817  if (!non_const_sys.solution->closed())
818  non_const_sys.solution->close();
819  non_const_sys.update();
820  }
821 
822  NumericVector<Number> & sys_soln(*system.current_local_solution);
823 
824  const DofMap & dof_map = system.get_dof_map();
825 
826  std::vector<Number> elem_soln; // The finite element solution
827  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
828  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
829 
830  unsigned var_inc = 0;
831  for (unsigned int var=0; var<nv_sys; var++)
832  {
833  const FEType & fe_type = system.variable_type(var);
834  const Variable & var_description = system.variable(var);
835  unsigned int n_vec_dim = FEInterface::n_vec_dim( sys_ptr->get_mesh(), fe_type );
836  const auto vg = dof_map.var_group_from_var_number(var);
837  const bool add_p_level = dof_map.should_p_refine(vg);
838 
839  for (const auto & elem : _mesh.active_local_element_ptr_range())
840  {
841  if (var_description.active_on_subdomain(elem->subdomain_id()))
842  {
843  dof_map.dof_indices (elem, dof_indices, var);
844  sys_soln.get(dof_indices, elem_soln);
845 
846  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  if (!elem->infinite())
856  {
857  libmesh_assert_equal_to (nodal_soln.size(), n_vec_dim*elem->n_nodes());
858 
859  for (auto n : elem->node_index_range())
860  {
861  const Node & node = elem->node_ref(n);
862 
863  const dof_id_type node_idx =
864  nv * node_id_to_vec_id(node.id());
865 
866  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  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  repeat_count.add(node_idx + (var_inc+d + var_num), 1);
874  }
875  }
876 
877  if (add_sides)
878  {
879  for (auto s : elem->side_index_range())
880  {
881  if (redundant_added_side(*elem,s))
882  continue;
883 
884  // Compute the FE solution at all the
885  // side nodes
887  (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  elem->nodes_on_side(s);
893 
894  libmesh_assert_equal_to
895  (nodal_soln.size(),
896  side_nodes.size());
897 #endif
898 
899  for (auto n : index_range(nodal_soln))
900  {
901  // Retrieve index into global solution vector.
902  std::size_t node_index =
903  nv * libmesh_map_find(discontinuous_node_indices,
904  std::make_tuple(elem->id(), s, n));
905 
906  for (unsigned int d=0; d < n_vec_dim; d++)
907  {
908  parallel_soln.add(node_index + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
909  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  for (auto n : elem->node_index_range())
918  {
919  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  if (!node.n_dofs(sys_num, var))
924  {
925  const dof_id_type node_idx =
926  nv * node_id_to_vec_id(node.id());
927 
928  for (unsigned int d=0; d < n_vec_dim; d++)
929  repeat_count.add(node_idx + (var_inc+d + var_num), 1);
930  }
931  }
932 
933  } // end loop over elements
934  var_inc += n_vec_dim;
935  } // end loop on variables in this system
936 
937  var_num += nv_sys_split;
938  } // end loop over systems
939 
940  // Sum the nodal solution values and repeat counts.
941  parallel_soln.close();
942  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  if (n_gaps)
949  {
950  for (numeric_index_type i=repeat_count.first_local_index();
951  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  if (std::abs(repeat_count(i)) < TOLERANCE)
957  repeat_count.set(i, 1.);
958  }
959 
960  // Make sure the repeat_count vector is up-to-date on all
961  // processors.
962  repeat_count.close();
963  }
964 
965  // Divide to get the average value at the nodes
966  parallel_soln /= repeat_count;
967 
968  return parallel_soln_ptr;
969 }
970 
971 
972 
973 void EquationSystems::build_solution_vector (std::vector<Number> & soln,
974  const std::set<std::string> * system_names,
975  bool add_sides) const
976 {
977  LOG_SCOPE("build_solution_vector()", "EquationSystems");
978 
979  // Call the parallel implementation
980  std::unique_ptr<NumericVector<Number>> parallel_soln =
981  this->build_parallel_solution_vector(system_names, add_sides);
982 
983  // Localize the NumericVector into the provided std::vector.
984  parallel_soln->localize_to_one(soln);
985 }
986 
987 
988 
989 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  vars_active_subdomains.clear();
993  vars_active_subdomains.resize(names.size());
994 
995  for (const auto & pr : _systems)
996  {
997  const auto & sys_ptr = pr.second;
998  for (auto vn : make_range(sys_ptr->n_vars()))
999  {
1000  const std::string & var_name = sys_ptr->variable_name(vn);
1001 
1002  auto names_it = std::find(names.begin(), names.end(), var_name);
1003  if(names_it != names.end())
1004  {
1005  const Variable & variable = sys_ptr->variable(vn);
1006  const std::set<subdomain_id_type> & active_subdomains = variable.active_subdomains();
1007  vars_active_subdomains[std::distance(names.begin(), names_it)] = active_subdomains;
1008  }
1009  }
1010  }
1011 }
1012 
1013 
1014 
1015 #ifdef LIBMESH_ENABLE_DEPRECATED
1016 void EquationSystems::get_solution (std::vector<Number> & soln,
1017  std::vector<std::string> & names) const
1018 {
1019  libmesh_deprecated();
1020  this->build_elemental_solution_vector(soln, names);
1021 }
1022 #endif // LIBMESH_ENABLE_DEPRECATED
1023 
1024 
1025 
1026 void
1028  std::vector<std::string> & names) const
1029 {
1030  // Call the parallel version of this function
1031  std::unique_ptr<NumericVector<Number>> parallel_soln =
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  soln.clear();
1040  if (parallel_soln)
1041  parallel_soln->localize_to_one(soln);
1042 }
1043 
1044 std::vector<std::pair<unsigned int, unsigned int>>
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  parallel_object_only();
1050 
1051  libmesh_assert (this->n_systems());
1052 
1053  // Resolve class of type input and assert that at least one of them is null
1054  libmesh_assert_msg(!type || !types,
1055  "Input 'type', 'types', or neither in find_variable_numbers, but not both.");
1056 
1057  std::vector<FEType> type_filter;
1058  if (type)
1059  type_filter.push_back(*type);
1060  else if (types)
1061  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  std::vector<std::string> name_filter = names;
1068  bool is_names_empty = name_filter.empty();
1069  names.clear();
1070 
1071  // initialize convenience variables
1072  FEType var_type;
1073  std::string name;
1074 
1075  const std::vector<std::string> component_suffix = {"_x", "_y", "_z"};
1076  unsigned int dim = _mesh.spatial_dimension();
1077  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  std::vector<std::pair<unsigned int, unsigned int>> var_nums;
1082 
1083  for (const auto & pr : _systems)
1084  {
1085  const System & system = *(pr.second);
1086 
1087  for (auto var : make_range(system.n_vars()))
1088  {
1089  // apply the type filter
1090  var_type = system.variable_type(var);
1091  if (type_filter.size() &&
1092  std::find(type_filter.begin(), type_filter.end(), var_type) == type_filter.end())
1093  continue;
1094 
1095  // apply the name filter (note that all variables pass if it is empty)
1096  if (FEInterface::field_type(var_type) == TYPE_VECTOR)
1097  {
1098  std::vector<std::string> component_names;
1099  for (unsigned int comp = 0; comp < dim; ++comp)
1100  {
1101  name = system.variable_name(var) + component_suffix[comp];
1102  if (is_names_empty ||
1103  (std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end()))
1104  component_names.push_back(name);
1105  }
1106 
1107  if (! component_names.empty())
1108  names.insert(names.end(), component_names.begin(), component_names.end());
1109  else
1110  continue;
1111  }
1112  else /*scalar-valued variable*/
1113  {
1114  name = system.variable_name(var);
1115  if (is_names_empty ||
1116  (std::find(name_filter.begin(), name_filter.end(), name) != name_filter.end()))
1117  names.push_back(name);
1118  else
1119  continue;
1120  }
1121 
1122  // if the variable made it through both filters get its system indices
1123  var_nums.emplace_back(system.number(), var);
1124  }
1125  }
1126 
1127  // Sort the var_nums vector pairs alphabetically based on the variable name
1128  std::vector<unsigned int> sort_index(var_nums.size());
1129  std::iota(sort_index.begin(), sort_index.end(), 0);
1130  std::sort(sort_index.begin(), sort_index.end(),
1131  [&](const unsigned int & lhs, const unsigned int & rhs)
1132  {return this->get_system(var_nums[lhs].first).variable_name(var_nums[lhs].second) <
1133  this->get_system(var_nums[rhs].first).variable_name(var_nums[rhs].second);});
1134 
1135  std::vector<std::pair<unsigned int, unsigned int>> var_nums_sorted(var_nums.size());
1136  for (auto i : index_range(var_nums_sorted))
1137  {
1138  var_nums_sorted[i].first = var_nums[sort_index[i]].first;
1139  var_nums_sorted[i].second = var_nums[sort_index[i]].second;
1140  }
1141 
1142  // Also sort the names vector
1143  std::sort(names.begin(), names.end());
1144 
1145  // Return the sorted vector pairs
1146  return var_nums_sorted;
1147 }
1148 
1149 
1150 std::unique_ptr<NumericVector<Number>>
1151 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  std::vector<FEType> type = {FEType(CONSTANT, MONOMIAL), FEType(CONSTANT, MONOMIAL_VEC)};
1157  std::vector<std::pair<unsigned int, unsigned int>> var_nums =
1158  this->find_variable_numbers(names, /*type=*/nullptr, &type);
1159 
1160  const std::size_t nv = names.size(); /*total number of vars including vector components*/
1161  const dof_id_type ne = _mesh.n_elem();
1162  libmesh_assert_equal_to (ne, _mesh.max_elem_id());
1163 
1164  // If there are no variables to write out don't do anything...
1165  if (!nv)
1166  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  numeric_index_type parallel_soln_global_size = ne*nv;
1171 
1172  numeric_index_type div = parallel_soln_global_size / this->n_processors();
1173  numeric_index_type mod = parallel_soln_global_size % this->n_processors();
1174 
1175  // Initialize all processors to the average size.
1176  numeric_index_type parallel_soln_local_size = div;
1177 
1178  // The first "mod" processors get an extra entry.
1179  if (this->processor_id() < mod)
1180  parallel_soln_local_size = div+1;
1181 
1182  // Create a NumericVector to hold the parallel solution
1183  std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
1184  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
1185  parallel_soln.init(parallel_soln_global_size,
1186  parallel_soln_local_size,
1187  /*fast=*/false,
1188  /*ParallelType=*/PARALLEL);
1189 
1190  unsigned int sys_ctr = 0;
1191  unsigned int var_ctr = 0;
1192  for (auto i : index_range(var_nums))
1193  {
1194  std::pair<unsigned int, unsigned int> var_num = var_nums[i];
1195  const System & system = this->get_system(var_num.first);
1196 
1197  // Update the current_local_solution if necessary
1198  if (sys_ctr != var_num.first)
1199  {
1200  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  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
1207  if (!non_const_sys.solution->closed())
1208  non_const_sys.solution->close();
1209  non_const_sys.update();
1210  sys_ctr = var_num.first;
1211  }
1212 
1213  NumericVector<Number> & sys_soln(*system.current_local_solution);
1214 
1215  // The DOF indices for the finite element
1216  std::vector<dof_id_type> dof_indices;
1217 
1218  const unsigned int var = var_num.second;
1219 
1220  const Variable & variable = system.variable(var);
1221  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  (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  for (const auto & elem : _mesh.active_local_element_ptr_range())
1232  if (variable.active_on_subdomain(elem->subdomain_id()))
1233  {
1234  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  libmesh_assert_equal_to(dof_indices.size(), n_comps);
1239 
1240  for (unsigned int comp = 0; comp < n_comps; comp++)
1241  parallel_soln.set(ne * (var_ctr + comp) + elem->id(), sys_soln(dof_indices[comp]));
1242  }
1243 
1244  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  libmesh_assert_equal_to(names.size(), var_ctr);
1253 
1254  parallel_soln.close();
1255  return parallel_soln_ptr;
1256 }
1257 
1258 
1259 
1260 void
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  LOG_SCOPE("build_discontinuous_solution_vector()", "EquationSystems");
1269 
1270  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  unsigned int nv = 0;
1275 
1276  for (const auto & [sys_name, sys_ptr] : _systems)
1277  {
1278  // Check current system is listed in system_names, and skip pos if not
1279  bool use_current_system = (system_names == nullptr);
1280  if (!use_current_system)
1281  use_current_system = system_names->count(sys_name);
1282  if (!use_current_system || sys_ptr->hide_output())
1283  continue;
1284 
1285  // Loop over all variables in this System and check whether we
1286  // are supposed to use each one.
1287  for (auto var_id : make_range(sys_ptr->n_vars()))
1288  {
1289  bool use_current_var = (var_names == nullptr);
1290  if (!use_current_var)
1291  use_current_var = std::count(var_names->begin(),
1292  var_names->end(),
1293  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  if (use_current_var)
1298  nv++;
1299  }
1300  }
1301 
1302  // get the total "weight" - the number of nodal values to write for
1303  // each variable.
1304  unsigned int tw=0;
1305  for (const auto & elem : _mesh.active_element_ptr_range())
1306  {
1307  tw += vertices_only ? elem->n_vertices() : elem->n_nodes();
1308 
1309  if (add_sides)
1310  {
1311  for (auto s : elem->side_index_range())
1312  {
1313  if (redundant_added_side(*elem,s))
1314  continue;
1315 
1316  const std::vector<unsigned int> side_nodes =
1317  elem->nodes_on_side(s);
1318 
1319  if (!vertices_only)
1320  tw += side_nodes.size();
1321  else
1322  for (auto n : index_range(side_nodes))
1323  if (elem->is_vertex(side_nodes[n]))
1324  ++tw;
1325  }
1326  }
1327  }
1328 
1329  // Only if we are on processor zero, allocate the storage
1330  // to hold (number_of_nodes)*(number_of_variables) entries.
1331  if (_mesh.processor_id() == 0)
1332  soln.resize(tw*nv);
1333 
1334  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  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  for (const auto & [sys_name, system] : _systems)
1346  {
1347  // Check current system is listed in system_names, and skip pos if not
1348  bool use_current_system = (system_names == nullptr);
1349  if (!use_current_system)
1350  use_current_system = system_names->count(sys_name);
1351  if (!use_current_system || system->hide_output())
1352  continue;
1353 
1354  const unsigned int nv_sys = system->n_vars();
1355  const auto & dof_map = system->get_dof_map();
1356 
1357  system->update_global_solution (sys_soln, 0);
1358 
1359  // Keep track of the number of vars actually written.
1360  unsigned int n_vars_written_current_system = 0;
1361 
1362  if (_mesh.processor_id() == 0)
1363  {
1364  std::vector<Number> soln_coeffs; // The finite element solution coeffs
1365  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
1366  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  for (unsigned int var=0; var<nv_sys; var++)
1374  {
1375  bool use_current_var = (var_names == nullptr);
1376  if (!use_current_var)
1377  use_current_var = std::count(var_names->begin(),
1378  var_names->end(),
1379  system->variable_name(var));
1380 
1381  // If we aren't supposed to write this var, go to the
1382  // next loop iteration.
1383  if (!use_current_var)
1384  continue;
1385 
1386  const FEType & fe_type = system->variable_type(var);
1387  const Variable & var_description = system->variable(var);
1388  const auto vg = dof_map.var_group_from_var_number(var);
1389  const bool add_p_level = dof_map.should_p_refine(vg);
1390 
1391  unsigned int nn=0;
1392 
1393  for (auto & elem : _mesh.active_element_ptr_range())
1394  {
1395  if (var_description.active_on_subdomain(elem->subdomain_id()))
1396  {
1397  system->get_dof_map().dof_indices (elem, dof_indices, var);
1398 
1399  soln_coeffs.resize(dof_indices.size());
1400 
1401  for (auto i : index_range(dof_indices))
1402  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  FEInterface::nodal_soln (elem->dim(),
1408  fe_type,
1409  elem,
1410  soln_coeffs,
1411  nodal_soln,
1412  add_p_level,
1413  FEInterface::n_vec_dim(_mesh, fe_type));
1414 
1415  // infinite elements should be skipped...
1416  if (!elem->infinite())
1417  {
1418  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
1419 
1420  const unsigned int n_vals =
1421  vertices_only ? elem->n_vertices() : elem->n_nodes();
1422 
1423  for (unsigned int n=0; n<n_vals; n++)
1424  {
1425  // Compute index into global solution vector.
1426  std::size_t index =
1427  nv * (nn++) + (n_vars_written_current_system + var_offset);
1428 
1429  soln[index] += nodal_soln[n];
1430  }
1431  }
1432  }
1433  else
1434  nn += vertices_only ? elem->n_vertices() : elem->n_nodes();
1435  } // end loop over active elements writing interiors
1436 
1437  // Loop writing "fake" sides, if requested
1438  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  elems_by_pid(_mesh.n_processors());
1447 
1448  for (const auto & elem : _mesh.active_element_ptr_range())
1449  elems_by_pid[elem->processor_id()].push_back(elem);
1450 
1451  for (auto p : index_range(elems_by_pid))
1452  for (const Elem * elem : elems_by_pid[p])
1453  {
1454  if (var_description.active_on_subdomain(elem->subdomain_id()))
1455  {
1456  system->get_dof_map().dof_indices (elem, dof_indices, var);
1457 
1458  soln_coeffs.resize(dof_indices.size());
1459 
1460  for (auto i : index_range(dof_indices))
1461  soln_coeffs[i] = sys_soln[dof_indices[i]];
1462 
1463  for (auto s : elem->side_index_range())
1464  {
1465  if (redundant_added_side(*elem,s))
1466  continue;
1467 
1468  const std::vector<unsigned int> side_nodes =
1469  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.
1476  (fe_type, elem, s, soln_coeffs,
1477  nodal_soln, add_p_level,
1478  FEInterface::n_vec_dim(_mesh, fe_type));
1479 
1480  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  FEInterface::get_continuity(fe_type);
1492  const Elem * const neigh = elem->neighbor_ptr(s);
1493 
1494  if ((cont == DISCONTINUOUS || cont == H_CURL || cont == H_DIV) &&
1495  neigh &&
1496  neigh->level() == elem->level() &&
1497  var_description.active_on_subdomain(neigh->subdomain_id()))
1498  {
1499  std::vector<dof_id_type> neigh_indices;
1500  system->get_dof_map().dof_indices (neigh, neigh_indices, var);
1501  std::vector<Number> neigh_coeffs(neigh_indices.size());
1502 
1503  for (auto i : index_range(neigh_indices))
1504  neigh_coeffs[i] = sys_soln[neigh_indices[i]];
1505 
1506  const unsigned int s_neigh =
1507  neigh->which_neighbor_am_i(elem);
1508  std::vector<Number> neigh_soln;
1510  (fe_type, neigh, s_neigh,
1511  neigh_coeffs, neigh_soln, add_p_level,
1512  FEInterface::n_vec_dim(_mesh, fe_type));
1513 
1514  const std::vector<unsigned int> neigh_nodes =
1515  neigh->nodes_on_side(s_neigh);
1516  for (auto n : index_range(side_nodes))
1517  for (auto neigh_n : index_range(neigh_nodes))
1518  if (neigh->node_ptr(neigh_nodes[neigh_n])
1519  == elem->node_ptr(side_nodes[n]))
1520  {
1521  nodal_soln[n] += neigh_soln[neigh_n];
1522  nodal_soln[n] /= 2;
1523  }
1524  }
1525 
1526  for (auto n : index_range(side_nodes))
1527  {
1528  if (vertices_only &&
1529  !elem->is_vertex(n))
1530  continue;
1531 
1532  // Compute index into global solution vector.
1533  std::size_t index =
1534  nv * (nn++) + (n_vars_written_current_system + var_offset);
1535 
1536  soln[index] += nodal_soln[n];
1537  }
1538  }
1539  }
1540  else
1541  {
1542  nn += vertices_only ? elem->n_vertices() : elem->n_nodes();
1543 
1544  for (auto s : elem->side_index_range())
1545  {
1546  if (redundant_added_side(*elem,s))
1547  continue;
1548 
1549  const std::vector<unsigned int> side_nodes =
1550  elem->nodes_on_side(s);
1551 
1552  for (auto n : index_range(side_nodes))
1553  {
1554  if (vertices_only &&
1555  !elem->is_vertex(n))
1556  continue;
1557  nn++;
1558  }
1559  }
1560  }
1561  } // end loop over active elements, writing "fake" sides
1562  }
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  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  var_offset += n_vars_written_current_system;
1572  } // end loop over systems
1573 }
1574 
1575 
1576 
1577 bool EquationSystems::redundant_added_side(const Elem & elem, unsigned int side)
1578 {
1579  libmesh_assert(elem.active());
1580 
1581  const Elem * neigh = elem.neighbor_ptr(side);
1582 
1583  // Write boundary sides.
1584  if (!neigh)
1585  return false;
1586 
1587  // Write ghost sides in Nemesis
1588  if (neigh == remote_elem)
1589  return false;
1590 
1591  // Don't write a coarser side if a finer side exists
1592  if (!neigh->active())
1593  return true;
1594 
1595  // Don't write a side redundantly from both of the
1596  // elements sharing it. We'll disambiguate with id().
1597  return (neigh->id() < elem.id());
1598 }
1599 
1600 
1601 
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  std::vector<bool> os_result;
1609 
1610  if (this->n_systems() != other_es.n_systems())
1611  {
1612  if (verbose)
1613  {
1614  libMesh::out << " Fatal difference. This system handles "
1615  << this->n_systems() << " systems," << std::endl
1616  << " while the other system handles "
1617  << other_es.n_systems()
1618  << " systems." << std::endl
1619  << " Aborting comparison." << std::endl;
1620  }
1621  return false;
1622  }
1623  else
1624  {
1625  // start comparing each system
1626  for (const auto & [sys_name, sys_ptr] : _systems)
1627  {
1628  // get the other system
1629  const System & other_system = other_es.get_system (sys_name);
1630 
1631  os_result.push_back (sys_ptr->compare (other_system, threshold, verbose));
1632 
1633  }
1634 
1635  }
1636 
1637 
1638  // sum up the results
1639  if (os_result.size()==0)
1640  return true;
1641  else
1642  {
1643  bool os_identical;
1644  unsigned int n = 0;
1645  do
1646  {
1647  os_identical = os_result[n];
1648  n++;
1649  }
1650  while (os_identical && n<os_result.size());
1651  return os_identical;
1652  }
1653 }
1654 
1655 
1656 
1657 std::string EquationSystems::get_info () const
1658 {
1659  std::ostringstream oss;
1660 
1661  oss << " EquationSystems\n"
1662  << " n_systems()=" << this->n_systems() << '\n';
1663 
1664  // Print the info for the individual systems
1665  for (const auto & pr : _systems)
1666  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  return oss.str();
1686 }
1687 
1688 
1689 
1690 void EquationSystems::print_info (std::ostream & os) const
1691 {
1692  os << this->get_info()
1693  << std::endl;
1694 }
1695 
1696 
1697 
1698 std::ostream & operator << (std::ostream & os,
1699  const EquationSystems & es)
1700 {
1701  es.print_info(os);
1702  return os;
1703 }
1704 
1705 
1706 
1707 unsigned int EquationSystems::n_vars () const
1708 {
1709  unsigned int tot=0;
1710 
1711  for (const auto & pr : _systems)
1712  tot += pr.second->n_vars();
1713 
1714  return tot;
1715 }
1716 
1717 
1718 
1719 std::size_t EquationSystems::n_dofs () const
1720 {
1721  std::size_t tot=0;
1722 
1723  for (const auto & pr : _systems)
1724  tot += pr.second->n_dofs();
1725 
1726  return tot;
1727 }
1728 
1729 
1730 
1731 
1733 {
1734  std::size_t tot=0;
1735 
1736  for (const auto & pr : _systems)
1737  tot += pr.second->n_active_dofs();
1738 
1739  return tot;
1740 }
1741 
1742 
1744 {
1745  // All the nodes
1746  for (auto & node : _mesh.node_ptr_range())
1747  node->add_system();
1748 
1749  // All the elements
1750  for (auto & elem : _mesh.element_ptr_range())
1751  elem->add_system();
1752 }
1753 
1755 {
1756  this->get_system(sys_num).get_dof_map().remove_default_ghosting();
1757 }
1758 
1759 } // namespace libMesh
EquationSystems(MeshBase &mesh)
Constructor.
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
This is the EquationSystems class.
unsigned int n_vars() const
void build_variable_names(std::vector< std::string > &var_names, const FEType *type=nullptr, const std::set< std::string > *system_names=nullptr) const
Fill the input vector var_names with the names of the variables for each system.
A Node is like a Point, but with more information.
Definition: node.h:52
unsigned int n_systems() const
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2458
bool _enable_default_ghosting
Flag for whether to enable default ghosting on newly added Systems.
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
Data structure for specifying which Parameters should be independent variables in a parameter sensiti...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
std::unique_ptr< NumericVector< Number > > build_parallel_elemental_solution_vector(std::vector< std::string > &names) const
Builds a parallel vector of CONSTANT MONOMIAL and/or components of CONSTANT MONOMIAL_VEC solution val...
std::size_t distribute_dofs(MeshBase &)
Distribute dofs on the current mesh.
Definition: dof_map.C:978
static constexpr Real TOLERANCE
std::size_t n_dofs() const
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
unsigned int dim
void build_discontinuous_solution_vector(std::vector< Number > &soln, const std::set< std::string > *system_names=nullptr, const std::vector< std::string > *var_names=nullptr, bool vertices_only=false, bool add_sides=false) const
Fill the input vector soln with solution values.
void add_default_ghosting()
Add the default functor(s) for coupling and algebraic ghosting.
Definition: dof_map.C:2024
bool refine_elements()
Only refines the user-requested elements.
bool reinit_solutions()
Handle any mesh changes and project any solutions onto the updated mesh.
virtual void clear()
Restores the data structure to a pristine state.
virtual void allgather()
Gathers all elements and nodes of the mesh onto every processor.
Definition: mesh_base.h:240
virtual ~EquationSystems()
Destructor.
void sum(T &r) const
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
MeshBase & mesh
static FEFieldType field_type(const FEType &fe_type)
std::ostream & operator<<(std::ostream &os, const OrderWrapper &order)
Overload stream operators.
Definition: fe_type.h:182
processor_id_type rank() const
const Parallel::Communicator & comm() const
virtual void enable_default_ghosting(bool enable)
Enable or disable default ghosting functors on the Mesh and on all Systems.
void clean_refinement_flags()
Sets the refinement flag to Elem::DO_NOTHING for each element in the mesh.
void should_p_refine(unsigned int g, bool p_refine)
Describe whether the given variable group should be p-refined.
Definition: dof_map.h:2402
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
Call adjoint_solve on all the individual equation systems.
const Parallel::Communicator & _communicator
The libMesh namespace provides an interface to certain functionality in the library.
const T_sys & get_system(std::string_view name) const
Real distance(const Point &p)
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Change the dimension of the vector to n.
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
static void side_nodal_soln(const FEType &fe_t, const Elem *elem, const unsigned int side, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, const bool add_p_level=true, const unsigned int vdim=1)
Build the nodal soln on one side from the (full) element soln.
Definition: fe_interface.C:650
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
virtual std::string get_info() const
virtual void prolong_vectors()
Prolong vectors after the mesh has refined.
Definition: system.C:444
MeshBase & _mesh
The mesh data structure.
bool coarsen_elements()
Only coarsens the user-requested elements.
unsigned int var_group_from_var_number(unsigned int var_num) const
Definition: dof_map.h:2430
signed char & underrefined_boundary_limit()
If underrefined_boundary_limit is set to a nonnegative value, then refinement and coarsening will pro...
unsigned int n_dofs(const unsigned int s, const unsigned int var=libMesh::invalid_uint) const
Definition: dof_object.h:806
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
void allgather()
Serializes a distributed mesh and its associated degree of freedom numbering for all systems...
void get_vars_active_subdomains(const std::vector< std::string > &names, std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) const
Retrieve vars_active_subdomains, which indicates the active subdomains for each variable in names...
processor_id_type size() const
virtual void sensitivity_solve(const ParameterVector &parameters)
Call sensitivity_solve on all the individual equation systems.
virtual void reinit_systems()
Reinitialize all systems on the current mesh.
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
unsigned int number() const
Definition: system.h:2350
This class defines the notion of a variable in the system.
Definition: variable.h:50
const std::set< subdomain_id_type > & active_subdomains() const
Definition: variable.h:181
Implements (adaptive) mesh refinement algorithms for a MeshBase.
dof_id_type id() const
Definition: dof_object.h:828
dof_id_type numeric_index_type
Definition: id_types.h:99
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2919
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
virtual dof_id_type max_elem_id() const =0
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
libmesh_assert(ctx)
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2478
void build_solution_vector(std::vector< Number > &soln, std::string_view system_name, std::string_view variable_name="all_vars") const
Fill the input vector soln with the solution values for the system named name.
virtual void restrict_vectors()
Restrict vectors after the mesh has coarsened.
Definition: system.C:386
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
bool _refine_in_reinit
Flag for whether to call coarsen/refine in reinit().
virtual void reinit_mesh()
Handle the association of a completely new mesh with the EquationSystem and all the Systems assigned ...
unsigned char & face_level_mismatch_limit()
If face_level_mismatch_limit is set to a nonzero value, then refinement and coarsening will produce m...
An object whose state is distributed along a set of processors.
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
virtual void reinit_constraints()
Reinitializes the constraints for this system.
Definition: system.C:495
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType&#39;s FEContinuity based on the underlying FEFamily and potentially the Order...
virtual void solve()
Call solve on all the individual equation systems.
virtual numeric_index_type first_local_index() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
void get_solution(std::vector< Number > &soln, std::vector< std::string > &names) const
Retrieve the solution data for CONSTANT MONOMIALs and/or components of CONSTANT MONOMIAL_VECs.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
unsigned int level() const
Definition: elem.h:3074
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2508
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool contract()=0
Delete subactive (i.e.
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
unsigned int spatial_dimension() const
Definition: mesh_base.C:543
T & set(const std::string &)
Definition: parameters.h:469
std::unique_ptr< NumericVector< Number > > build_parallel_solution_vector(const std::set< std::string > *system_names=nullptr, bool add_sides=false) const
A version of build_solution_vector which is appropriate for "parallel" output formats like Nemesis...
OStreamProxy out
virtual void clear()
Clears internal data structures & frees any allocated memory.
Definition: parameters.h:324
signed char & overrefined_boundary_limit()
If overrefined_boundary_limit is set to a nonnegative value, then refinement and coarsening will prod...
void remove_default_ghosting()
Remove any default ghosting functor(s).
Definition: dof_map.C:2016
void _remove_default_ghosting(unsigned int sys_num)
This just calls DofMap::remove_default_ghosting() but using a shim lets us forward-declare DofMap...
const MeshBase & get_mesh() const
std::size_t n_active_dofs() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
std::vector< std::pair< unsigned int, unsigned int > > find_variable_numbers(std::vector< std::string > &names, const FEType *type=nullptr, const std::vector< FEType > *types=nullptr) const
Finds system and variable numbers for any variables of &#39;type&#39; or of &#39;types&#39; corresponding to the entr...
Parameters parameters
Data structure holding arbitrary parameters.
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
virtual void init()
Initialize all the systems.
virtual bool compare(const EquationSystems &other_es, const Real threshold, const bool verbose) const
unsigned int n_vars() const
Definition: system.h:2430
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
virtual void add(const numeric_index_type i, const T value)=0
Adds value to the vector entry specified by i.
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
processor_id_type processor_id() const
static void nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, const bool add_p_level=true, const unsigned int vdim=1)
Build the nodal soln from the element soln.
Definition: fe_interface.C:625
bool active() const
Definition: elem.h:2941
void _add_system_to_nodes_and_elems()
This function is used in the implementation of add_system, it loops over the nodes and elements of th...
const DofMap & get_dof_map() const
Definition: system.h:2374
virtual numeric_index_type last_local_index() const =0
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
void update()
Updates local values for all the systems.
virtual std::vector< unsigned int > nodes_on_side(const unsigned int) const =0
void build_elemental_solution_vector(std::vector< Number > &soln, std::vector< std::string > &names) const
Retrieve the solution data for CONSTANT MONOMIALs and/or components of CONSTANT MONOMIAL_VECs.
uint8_t dof_id_type
Definition: id_types.h:67
std::map< std::string, std::unique_ptr< System >, std::less<> > _systems
Data structure holding the systems.
static bool redundant_added_side(const Elem &elem, unsigned int side)
const RemoteElem * remote_elem
Definition: remote_elem.C:57