equation_systems.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // System includes
20 #include <sstream>
21 
22 // Local Includes
23 #include "libmesh/default_coupling.h" // For downconversion
25 #include "libmesh/fe_interface.h"
29 #include "libmesh/newmark_system.h"
33 #include "libmesh/eigen_system.h"
34 #include "libmesh/parallel.h"
36 #include "libmesh/dof_map.h"
37 #include "libmesh/mesh_base.h"
38 #include "libmesh/elem.h"
40 
41 // Include the systems before this one to avoid
42 // overlapping forward declarations.
44 
45 namespace libMesh
46 {
47 
48 // Forward Declarations
49 
50 
51 
52 
53 // ------------------------------------------------------------
54 // EquationSystems class implementation
56  ParallelObject (m),
57  _mesh (m),
58  _refine_in_reinit(true),
59  _enable_default_ghosting(true)
60 {
61  // Set default parameters
62  this->parameters.set<Real> ("linear solver tolerance") = TOLERANCE * TOLERANCE;
63  this->parameters.set<unsigned int>("linear solver maximum iterations") = 5000;
64 }
65 
66 
67 
69 {
70  this->clear ();
71 }
72 
73 
74 
76 {
77  // Clear any additional parameters
78  parameters.clear ();
79 
80  // clear the systems. We must delete them
81  // since we newed them!
82  while (!_systems.empty())
83  {
84  system_iterator pos = _systems.begin();
85 
86  System * sys = pos->second;
87  delete sys;
88  sys = nullptr;
89 
90  _systems.erase (pos);
91  }
92 }
93 
94 
95 
97 {
98  const unsigned int n_sys = this->n_systems();
99 
100  libmesh_assert_not_equal_to (n_sys, 0);
101 
102  // Tell all the \p DofObject entities how many systems
103  // there are.
104  for (auto & node : _mesh.node_ptr_range())
105  node->set_n_systems(n_sys);
106 
107  for (auto & elem : _mesh.element_ptr_range())
108  elem->set_n_systems(n_sys);
109 
110  for (unsigned int i=0; i != this->n_systems(); ++i)
111  this->get_system(i).init();
112 
113 #ifdef LIBMESH_ENABLE_AMR
114  MeshRefinement mesh_refine(_mesh);
115  mesh_refine.clean_refinement_flags();
116 #endif
117 }
118 
119 
120 
122 {
123  const bool mesh_changed = this->reinit_solutions();
124 
125  // If the mesh has changed, systems will need to reinitialize their
126  // own data on the new mesh.
127  if (mesh_changed)
128  this->reinit_systems();
129 }
130 
131 
132 
134 {
135  parallel_object_only();
136 
137  const unsigned int n_sys = this->n_systems();
138  libmesh_assert_not_equal_to (n_sys, 0);
139 
140  // We may have added new systems since our last
141  // EquationSystems::(re)init call
142  bool _added_new_systems = false;
143  for (unsigned int i=0; i != n_sys; ++i)
144  if (!this->get_system(i).is_initialized())
145  _added_new_systems = true;
146 
147  if (_added_new_systems)
148  {
149  // Our DofObjects will need space for the additional systems
150  for (auto & node : _mesh.node_ptr_range())
151  node->set_n_systems(n_sys);
152 
153  for (auto & elem : _mesh.element_ptr_range())
154  elem->set_n_systems(n_sys);
155 
156  // And any new systems will need initialization
157  for (unsigned int i=0; i != n_sys; ++i)
158  if (!this->get_system(i).is_initialized())
159  this->get_system(i).init();
160  }
161 
162 
163  // We used to assert that all nodes and elements *already* had
164  // n_systems() properly set; however this is false in the case where
165  // user code has manually added nodes and/or elements to an
166  // already-initialized system.
167 
168  // Make sure all the \p DofObject entities know how many systems
169  // there are.
170  {
171  // All the nodes
172  for (auto & node : _mesh.node_ptr_range())
173  node->set_n_systems(this->n_systems());
174 
175  // All the elements
176  for (auto & elem : _mesh.element_ptr_range())
177  elem->set_n_systems(this->n_systems());
178  }
179 
180  // Localize each system's vectors
181  for (unsigned int i=0; i != this->n_systems(); ++i)
182  this->get_system(i).re_update();
183 
184 #ifdef LIBMESH_ENABLE_AMR
185 
186  bool dof_constraints_created = false;
187  bool mesh_changed = false;
188 
189  // FIXME: For backwards compatibility, assume
190  // refine_and_coarsen_elements or refine_uniformly have already
191  // been called
192  {
193  for (unsigned int i=0; i != this->n_systems(); ++i)
194  {
195  System & sys = this->get_system(i);
196 
197  // Even if the system doesn't have any variables in it we want
198  // consistent behavior; e.g. distribute_dofs should have the
199  // opportunity to count up zero dofs on each processor.
200  //
201  // Who's been adding zero-var systems anyway, outside of my
202  // unit tests? - RHS
203  // if (!sys.n_vars())
204  // continue;
205 
207 
208  // Recreate any user or internal constraints
209  sys.reinit_constraints();
210 
211  sys.prolong_vectors();
212  }
213  mesh_changed = true;
214  dof_constraints_created = true;
215  }
216 
217  if (this->_refine_in_reinit)
218  {
219  // Don't override any user refinement settings
220  MeshRefinement mesh_refine(_mesh);
221  mesh_refine.face_level_mismatch_limit() = 0; // unlimited
222  mesh_refine.overrefined_boundary_limit() = -1; // unlimited
223  mesh_refine.underrefined_boundary_limit() = -1; // unlimited
224 
225  // Try to coarsen the mesh, then restrict each system's vectors
226  // if necessary
227  if (mesh_refine.coarsen_elements())
228  {
229  for (unsigned int i=0; i != this->n_systems(); ++i)
230  {
231  System & sys = this->get_system(i);
232  if (!dof_constraints_created)
233  {
235  sys.reinit_constraints();
236  }
237  sys.restrict_vectors();
238  }
239  mesh_changed = true;
240  dof_constraints_created = true;
241  }
242 
243  // Once vectors are all restricted, we can delete
244  // children of coarsened elements
245  if (mesh_changed)
246  this->get_mesh().contract();
247 
248  // Try to refine the mesh, then prolong each system's vectors
249  // if necessary
250  if (mesh_refine.refine_elements())
251  {
252  for (unsigned int i=0; i != this->n_systems(); ++i)
253  {
254  System & sys = this->get_system(i);
255  if (!dof_constraints_created)
256  {
258  sys.reinit_constraints();
259  }
260  sys.prolong_vectors();
261  }
262  mesh_changed = true;
263  // dof_constraints_created = true;
264  }
265  }
266 
267  return mesh_changed;
268 
269 #endif // #ifdef LIBMESH_ENABLE_AMR
270 
271  return false;
272 }
273 
274 
275 
277 {
278  for (unsigned int i=0; i != this->n_systems(); ++i)
279  this->get_system(i).reinit();
280 }
281 
282 
283 
285 {
286  // A serial mesh means nothing needs to be done
287  if (_mesh.is_serial())
288  return;
289 
290  const unsigned int n_sys = this->n_systems();
291 
292  libmesh_assert_not_equal_to (n_sys, 0);
293 
294  // Gather the mesh
295  _mesh.allgather();
296 
297  // Tell all the \p DofObject entities how many systems
298  // there are.
299  for (auto & node : _mesh.node_ptr_range())
300  node->set_n_systems(n_sys);
301 
302  for (auto & elem : _mesh.element_ptr_range())
303  elem->set_n_systems(n_sys);
304 
305  // And distribute each system's dofs
306  for (unsigned int i=0; i != this->n_systems(); ++i)
307  {
308  System & sys = this->get_system(i);
309  DofMap & dof_map = sys.get_dof_map();
310  dof_map.distribute_dofs(_mesh);
311 
312  // The user probably won't need constraint equations or the
313  // send_list after an allgather, but let's keep it in consistent
314  // shape just in case.
315  sys.reinit_constraints();
316  dof_map.prepare_send_list();
317  }
318 }
319 
320 
321 
323 {
324  _enable_default_ghosting = enable;
325  MeshBase &mesh = this->get_mesh();
326 
327  if (enable)
328  mesh.add_ghosting_functor(mesh.default_ghosting());
329  else
330  mesh.remove_ghosting_functor(mesh.default_ghosting());
331 
332  for (unsigned int i=0; i != this->n_systems(); ++i)
333  {
334  DofMap & dof_map = this->get_system(i).get_dof_map();
335  if (enable)
336  dof_map.add_default_ghosting();
337  else
338  dof_map.remove_default_ghosting();
339  }
340 }
341 
342 
343 
345 {
346  LOG_SCOPE("update()", "EquationSystems");
347 
348  // Localize each system's vectors
349  for (unsigned int i=0; i != this->n_systems(); ++i)
350  this->get_system(i).update();
351 }
352 
353 
354 
355 System & EquationSystems::add_system (const std::string & sys_type,
356  const std::string & name)
357 {
358  // If the user already built a system with this name, we'll
359  // trust them and we'll use it. That way they can pre-add
360  // non-standard derived system classes, and if their restart file
361  // has some non-standard sys_type we won't throw an error.
362  if (_systems.count(name))
363  {
364  return this->get_system(name);
365  }
366  // Build a basic System
367  else if (sys_type == "Basic")
368  this->add_system<System> (name);
369 
370  // Build a Newmark system
371  else if (sys_type == "Newmark")
372  this->add_system<NewmarkSystem> (name);
373 
374  // Build an Explicit system
375  else if ((sys_type == "Explicit"))
376  this->add_system<ExplicitSystem> (name);
377 
378  // Build an Implicit system
379  else if ((sys_type == "Implicit") ||
380  (sys_type == "Steady" ))
381  this->add_system<ImplicitSystem> (name);
382 
383  // build a transient implicit linear system
384  else if ((sys_type == "Transient") ||
385  (sys_type == "TransientImplicit") ||
386  (sys_type == "TransientLinearImplicit"))
387  this->add_system<TransientLinearImplicitSystem> (name);
388 
389  // build a transient implicit nonlinear system
390  else if (sys_type == "TransientNonlinearImplicit")
391  this->add_system<TransientNonlinearImplicitSystem> (name);
392 
393  // build a transient explicit system
394  else if (sys_type == "TransientExplicit")
395  this->add_system<TransientExplicitSystem> (name);
396 
397  // build a linear implicit system
398  else if (sys_type == "LinearImplicit")
399  this->add_system<LinearImplicitSystem> (name);
400 
401  // build a nonlinear implicit system
402  else if (sys_type == "NonlinearImplicit")
403  this->add_system<NonlinearImplicitSystem> (name);
404 
405  // build a Reduced Basis Construction system
406  else if (sys_type == "RBConstruction")
407  this->add_system<RBConstruction> (name);
408 
409  // build a transient Reduced Basis Construction system
410  else if (sys_type == "TransientRBConstruction")
411  this->add_system<TransientRBConstruction> (name);
412 
413 #ifdef LIBMESH_HAVE_SLEPC
414  // build an eigen system
415  else if (sys_type == "Eigen")
416  this->add_system<EigenSystem> (name);
417  else if (sys_type == "TransientEigenSystem")
418  this->add_system<TransientEigenSystem> (name);
419 #endif
420 
421 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
422  // build a frequency system
423  else if (sys_type == "Frequency")
424  this->add_system<FrequencySystem> (name);
425 #endif
426 
427  else
428  libmesh_error_msg("ERROR: Unknown system type: " << sys_type);
429 
430  // Return a reference to the new system
431  //return (*this)(name);
432  return this->get_system(name);
433 }
434 
435 
436 
437 
438 
439 
440 #ifdef LIBMESH_ENABLE_DEPRECATED
441 void EquationSystems::delete_system (const std::string & name)
442 {
443  libmesh_deprecated();
444 
445  if (!_systems.count(name))
446  libmesh_error_msg("ERROR: no system named " << name);
447 
448  delete _systems[name];
449 
450  _systems.erase (name);
451 }
452 #endif
453 
454 
455 
457 {
458  libmesh_assert (this->n_systems());
459 
460  for (unsigned int i=0; i != this->n_systems(); ++i)
461  this->get_system(i).solve();
462 }
463 
464 
465 
467 {
468  libmesh_assert (this->n_systems());
469 
470  for (unsigned int i=0; i != this->n_systems(); ++i)
471  this->get_system(i).sensitivity_solve(parameters_in);
472 }
473 
474 
475 
476 void EquationSystems::adjoint_solve (const QoISet & qoi_indices)
477 {
478  libmesh_assert (this->n_systems());
479 
480  for (unsigned int i=this->n_systems(); i != 0; --i)
481  this->get_system(i-1).adjoint_solve(qoi_indices);
482 }
483 
484 
485 
486 void EquationSystems::build_variable_names (std::vector<std::string> & var_names,
487  const FEType * type,
488  const std::set<std::string> * system_names) const
489 {
490  unsigned int var_num=0;
491 
492  const_system_iterator pos = _systems.begin();
493  const const_system_iterator end = _systems.end();
494 
495  // Need to size var_names by scalar variables plus all the
496  // vector components for all the vector variables
497  //Could this be replaced by a/some convenience methods?[PB]
498  {
499  unsigned int n_scalar_vars = 0;
500  unsigned int n_vector_vars = 0;
501 
502  for (; pos != end; ++pos)
503  {
504  // Check current system is listed in system_names, and skip pos if not
505  bool use_current_system = (system_names == nullptr);
506  if (!use_current_system)
507  use_current_system = system_names->count(pos->first);
508  if (!use_current_system || pos->second->hide_output())
509  continue;
510 
511  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
512  {
513  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
514  n_vector_vars++;
515  else
516  n_scalar_vars++;
517  }
518  }
519 
520  // Here, we're assuming the number of vector components is the same
521  // as the mesh dimension. Will break for mixed dimension meshes.
522  unsigned int dim = this->get_mesh().mesh_dimension();
523  unsigned int nv = n_scalar_vars + dim*n_vector_vars;
524 
525  // We'd better not have more than dim*his->n_vars() (all vector variables)
526  libmesh_assert_less_equal ( nv, dim*this->n_vars() );
527 
528  // Here, we're assuming the number of vector components is the same
529  // as the mesh dimension. Will break for mixed dimension meshes.
530 
531  var_names.resize( nv );
532  }
533 
534  // reset
535  pos = _systems.begin();
536 
537  for (; pos != end; ++pos)
538  {
539  // Check current system is listed in system_names, and skip pos if not
540  bool use_current_system = (system_names == nullptr);
541  if (!use_current_system)
542  use_current_system = system_names->count(pos->first);
543  if (!use_current_system || pos->second->hide_output())
544  continue;
545 
546  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
547  {
548  std::string var_name = pos->second->variable_name(vn);
549  FEType fe_type = pos->second->variable_type(vn);
550 
551  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type);
552 
553  // Filter on the type if requested
554  if (type == nullptr || (type && *type == fe_type))
555  {
556  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
557  {
558  switch(n_vec_dim)
559  {
560  case 0:
561  case 1:
562  var_names[var_num++] = var_name;
563  break;
564  case 2:
565  var_names[var_num++] = var_name+"_x";
566  var_names[var_num++] = var_name+"_y";
567  break;
568  case 3:
569  var_names[var_num++] = var_name+"_x";
570  var_names[var_num++] = var_name+"_y";
571  var_names[var_num++] = var_name+"_z";
572  break;
573  default:
574  libmesh_error_msg("Invalid dim in build_variable_names");
575  }
576  }
577  else
578  var_names[var_num++] = var_name;
579  }
580  }
581  }
582  // Now resize again in case we filtered any names
583  var_names.resize(var_num);
584 }
585 
586 
587 
588 void EquationSystems::build_solution_vector (std::vector<Number> &,
589  const std::string &,
590  const std::string &) const
591 {
592  // TODO:[BSK] re-implement this from the method below
593  libmesh_not_implemented();
594 }
595 
596 
597 
598 
599 std::unique_ptr<NumericVector<Number>>
600 EquationSystems::build_parallel_solution_vector(const std::set<std::string> * system_names) const
601 {
602  LOG_SCOPE("build_parallel_solution_vector()", "EquationSystems");
603 
604  // This function must be run on all processors at once
605  parallel_object_only();
606 
607  const unsigned int dim = _mesh.mesh_dimension();
608  const dof_id_type max_nn = _mesh.max_node_id();
609 
610  // allocate vector storage to hold
611  // (max_node_id)*(number_of_variables) entries.
612  //
613  // If node renumbering is disabled and adaptive coarsening has
614  // created gaps between node numbers, then this vector will be
615  // sparse.
616  //
617  // We have to differentiate between between scalar and vector
618  // variables. We intercept vector variables and treat each
619  // component as a scalar variable (consistently with build_solution_names).
620 
621  unsigned int nv = 0;
622 
623  //Could this be replaced by a/some convenience methods?[PB]
624  {
625  unsigned int n_scalar_vars = 0;
626  unsigned int n_vector_vars = 0;
627  const_system_iterator pos = _systems.begin();
628  const const_system_iterator end = _systems.end();
629 
630  for (; pos != end; ++pos)
631  {
632  // Check current system is listed in system_names, and skip pos if not
633  bool use_current_system = (system_names == nullptr);
634  if (!use_current_system)
635  use_current_system = system_names->count(pos->first);
636  if (!use_current_system || pos->second->hide_output())
637  continue;
638 
639  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
640  {
641  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
642  n_vector_vars++;
643  else
644  n_scalar_vars++;
645  }
646  }
647  // Here, we're assuming the number of vector components is the same
648  // as the mesh dimension. Will break for mixed dimension meshes.
649  nv = n_scalar_vars + dim*n_vector_vars;
650  }
651 
652  // Get the number of nodes to store locally.
653  dof_id_type n_local_nodes = cast_int<dof_id_type>
654  (std::distance(_mesh.local_nodes_begin(),
656 
657  // If node renumbering has been disabled, nodes may not be numbered
658  // contiguously, and the number of nodes might not match the
659  // max_node_id. In this case we just do our best.
660  dof_id_type n_total_nodes = n_local_nodes;
661  _mesh.comm().sum(n_total_nodes);
662 
663  const dof_id_type n_gaps = max_nn - n_total_nodes;
664  const dof_id_type gaps_per_processor = n_gaps / _mesh.comm().size();
665  const dof_id_type remainder_gaps = n_gaps % _mesh.comm().size();
666 
667  n_local_nodes = n_local_nodes + // Actual nodes
668  gaps_per_processor + // Our even share of gaps
669  (_mesh.comm().rank() < remainder_gaps); // Leftovers
670 
671  // Create a NumericVector to hold the parallel solution
672  std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
673  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
674  parallel_soln.init(max_nn*nv, n_local_nodes*nv, false, PARALLEL);
675 
676  // Create a NumericVector to hold the "repeat_count" for each node - this is essentially
677  // the number of elements contributing to that node's value
678  std::unique_ptr<NumericVector<Number>> repeat_count_ptr = NumericVector<Number>::build(_communicator);
679  NumericVector<Number> & repeat_count = *repeat_count_ptr;
680  repeat_count.init(max_nn*nv, n_local_nodes*nv, false, PARALLEL);
681 
682  repeat_count.close();
683 
684  unsigned int var_num=0;
685 
686  // For each system in this EquationSystems object,
687  // update the global solution and if we are on processor 0,
688  // loop over the elements and build the nodal solution
689  // from the element solution. Then insert this nodal solution
690  // into the vector passed to build_solution_vector.
691  const_system_iterator pos = _systems.begin();
692  const const_system_iterator end = _systems.end();
693 
694  for (; pos != end; ++pos)
695  {
696  // Check current system is listed in system_names, and skip pos if not
697  bool use_current_system = (system_names == nullptr);
698  if (!use_current_system)
699  use_current_system = system_names->count(pos->first);
700  if (!use_current_system || pos->second->hide_output())
701  continue;
702 
703  const System & system = *(pos->second);
704  const unsigned int nv_sys = system.n_vars();
705  const unsigned int sys_num = system.number();
706 
707  //Could this be replaced by a/some convenience methods?[PB]
708  unsigned int n_scalar_vars = 0;
709  unsigned int n_vector_vars = 0;
710  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
711  {
712  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
713  n_vector_vars++;
714  else
715  n_scalar_vars++;
716  }
717 
718  // Here, we're assuming the number of vector components is the same
719  // as the mesh dimension. Will break for mixed dimension meshes.
720  unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
721 
722  // Update the current_local_solution
723  {
724  System & non_const_sys = const_cast<System &>(system);
725  // We used to simply call non_const_sys.solution->close()
726  // here, but that is not allowed when the solution vector is
727  // locked read-only, for example when printing the solution
728  // during the middle of a solve... So try to be a bit
729  // more careful about calling close() unnecessarily.
730  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
731  if (!non_const_sys.solution->closed())
732  non_const_sys.solution->close();
733  non_const_sys.update();
734  }
735 
736  NumericVector<Number> & sys_soln(*system.current_local_solution);
737 
738  std::vector<Number> elem_soln; // The finite element solution
739  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
740  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
741 
742  unsigned var_inc = 0;
743  for (unsigned int var=0; var<nv_sys; var++)
744  {
745  const FEType & fe_type = system.variable_type(var);
746  const Variable & var_description = system.variable(var);
747  const DofMap & dof_map = system.get_dof_map();
748 
749  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type );
750 
751  for (const auto & elem : _mesh.active_local_element_ptr_range())
752  {
753  if (var_description.active_on_subdomain(elem->subdomain_id()))
754  {
755  dof_map.dof_indices (elem, dof_indices, var);
756 
757  elem_soln.resize(dof_indices.size());
758 
759  for (std::size_t i=0; i<dof_indices.size(); i++)
760  elem_soln[i] = sys_soln(dof_indices[i]);
761 
763  fe_type,
764  elem,
765  elem_soln,
766  nodal_soln);
767 
768 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
769  // infinite elements should be skipped...
770  if (!elem->infinite())
771 #endif
772  {
773  libmesh_assert_equal_to (nodal_soln.size(), n_vec_dim*elem->n_nodes());
774 
775  for (unsigned int n=0; n<elem->n_nodes(); n++)
776  {
777  for (unsigned int d=0; d < n_vec_dim; d++)
778  {
779  // For vector-valued elements, all components are in nodal_soln. For each
780  // node, the components are stored in order, i.e. node_0 -> s0_x, s0_y, s0_z
781  parallel_soln.add(nv*(elem->node_id(n)) + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
782 
783  // Increment the repeat count for this position
784  repeat_count.add(nv*(elem->node_id(n)) + (var_inc+d + var_num), 1);
785  }
786  }
787  }
788  }
789  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:
790  for (unsigned int n=0; n<elem->n_nodes(); n++)
791  // Only do this if this variable has NO DoFs at this node... it might have some from an adjoining element...
792  if (!elem->node_ptr(n)->n_dofs(sys_num, var))
793  for (unsigned int d=0; d < n_vec_dim; d++)
794  repeat_count.add(nv*(elem->node_id(n)) + (var+d + var_num), 1);
795 
796  } // end loop over elements
797  var_inc += n_vec_dim;
798  } // end loop on variables in this system
799 
800  var_num += nv_sys_split;
801  } // end loop over systems
802 
803  // Sum the nodal solution values and repeat counts.
804  parallel_soln.close();
805  repeat_count.close();
806 
807  // If there were gaps in the node numbering, there will be
808  // corresponding zeros in the parallel_soln and repeat_count
809  // vectors. We need to set those repeat_count entries to 1
810  // in order to avoid dividing by zero.
811  if (n_gaps)
812  {
813  for (numeric_index_type i=repeat_count.first_local_index();
814  i<repeat_count.last_local_index(); ++i)
815  {
816  // repeat_count entries are integral values but let's avoid a
817  // direct floating point comparison with 0 just in case some
818  // roundoff noise crept in during vector assembly?
819  if (std::abs(repeat_count(i)) < TOLERANCE)
820  repeat_count.set(i, 1.);
821  }
822 
823  // Make sure the repeat_count vector is up-to-date on all
824  // processors.
825  repeat_count.close();
826  }
827 
828  // Divide to get the average value at the nodes
829  parallel_soln /= repeat_count;
830 
831  return std::unique_ptr<NumericVector<Number>>(parallel_soln_ptr.release());
832 }
833 
834 
835 
836 void EquationSystems::build_solution_vector (std::vector<Number> & soln,
837  const std::set<std::string> * system_names) const
838 {
839  LOG_SCOPE("build_solution_vector()", "EquationSystems");
840 
841  // Call the parallel implementation
842  std::unique_ptr<NumericVector<Number>> parallel_soln =
843  this->build_parallel_solution_vector(system_names);
844 
845  // Localize the NumericVector into the provided std::vector.
846  parallel_soln->localize_to_one(soln);
847 }
848 
849 
850 
851 void EquationSystems::get_vars_active_subdomains(const std::vector<std::string> & names,
852  std::vector<std::set<subdomain_id_type>> & vars_active_subdomains) const
853 {
854  unsigned int var_num=0;
855 
856  vars_active_subdomains.clear();
857  vars_active_subdomains.resize(names.size());
858 
859  const_system_iterator pos = _systems.begin();
860  const const_system_iterator end = _systems.end();
861 
862  for (; pos != end; ++pos)
863  {
864  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
865  {
866  std::string var_name = pos->second->variable_name(vn);
867 
868  auto names_it = std::find(names.begin(), names.end(), var_name);
869  if(names_it != names.end())
870  {
871  const Variable & variable = pos->second->variable(vn);
872  const std::set<subdomain_id_type> & active_subdomains = variable.active_subdomains();
873  vars_active_subdomains[var_num++] = active_subdomains;
874  }
875  }
876  }
877 
878  libmesh_assert_equal_to(var_num, names.size());
879 }
880 
881 
882 
883 void EquationSystems::get_solution (std::vector<Number> & soln,
884  std::vector<std::string> & names) const
885 {
886  libmesh_deprecated();
887  this->build_elemental_solution_vector(soln, names);
888 }
889 
890 
891 
892 void
894  std::vector<std::string> & names) const
895 {
896  // Call the parallel version of this function
897  std::unique_ptr<NumericVector<Number>> parallel_soln =
899 
900  // Localize into 'soln', provided that parallel_soln is not empty.
901  // Note: parallel_soln will be empty in the event that none of the
902  // input names were CONSTANT, MONOMIAL variables or there were
903  // simply no CONSTANT, MONOMIAL variables in the EquationSystems
904  // object.
905  soln.clear();
906  if (parallel_soln)
907  parallel_soln->localize_to_one(soln);
908 }
909 
910 
911 
912 std::unique_ptr<NumericVector<Number>>
913 EquationSystems::build_parallel_elemental_solution_vector (std::vector<std::string> & names) const
914 {
915  // This function must be run on all processors at once
916  parallel_object_only();
917 
918  libmesh_assert (this->n_systems());
919 
920  const dof_id_type ne = _mesh.n_elem();
921 
922  libmesh_assert_equal_to (ne, _mesh.max_elem_id());
923 
924  // If the names vector has entries, we will only populate the soln vector
925  // with names included in that list. Note: The names vector may be
926  // reordered upon exiting this function
927  std::vector<std::string> filter_names = names;
928  bool is_filter_names = !filter_names.empty();
929 
930  names.clear();
931 
932  const FEType type(CONSTANT, MONOMIAL);
933 
934  dof_id_type nv = 0;
935 
936  // Find the total number of variables to output
937  std::vector<std::vector<unsigned>> do_output(_systems.size());
938  {
939  const_system_iterator pos = _systems.begin();
940  const const_system_iterator end = _systems.end();
941  unsigned sys_ctr = 0;
942 
943  for (; pos != end; ++pos, ++sys_ctr)
944  {
945  const System & system = *(pos->second);
946  const unsigned int nv_sys = system.n_vars();
947 
948  do_output[sys_ctr].resize(nv_sys);
949 
950  for (unsigned int var=0; var < nv_sys; ++var)
951  {
952  if (system.variable_type(var) != type ||
953  (is_filter_names && std::find(filter_names.begin(), filter_names.end(), system.variable_name(var)) == filter_names.end()))
954  continue;
955 
956  // Otherwise, this variable should be output
957  nv++;
958  do_output[sys_ctr][var] = 1;
959  }
960  }
961  }
962 
963  // If there are no variables to write out don't do anything...
964  if (!nv)
965  return std::unique_ptr<NumericVector<Number>>(nullptr);
966 
967  // We can handle the case where there are nullptrs in the Elem vector
968  // by just having extra zeros in the solution vector.
969  numeric_index_type parallel_soln_global_size = ne*nv;
970 
971  numeric_index_type div = parallel_soln_global_size / this->n_processors();
972  numeric_index_type mod = parallel_soln_global_size % this->n_processors();
973 
974  // Initialize all processors to the average size.
975  numeric_index_type parallel_soln_local_size = div;
976 
977  // The first "mod" processors get an extra entry.
978  if (this->processor_id() < mod)
979  parallel_soln_local_size = div+1;
980 
981  // Create a NumericVector to hold the parallel solution
982  std::unique_ptr<NumericVector<Number>> parallel_soln_ptr = NumericVector<Number>::build(_communicator);
983  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
984  parallel_soln.init(parallel_soln_global_size,
985  parallel_soln_local_size,
986  /*fast=*/false,
987  /*ParallelType=*/PARALLEL);
988 
989  dof_id_type var_num = 0;
990 
991  // For each system in this EquationSystems object,
992  // update the global solution and collect the
993  // CONSTANT MONOMIALs. The entries are in variable-major
994  // format.
995  const_system_iterator pos = _systems.begin();
996  const const_system_iterator end = _systems.end();
997  unsigned sys_ctr = 0;
998 
999  for (; pos != end; ++pos, ++sys_ctr)
1000  {
1001  const System & system = *(pos->second);
1002  const unsigned int nv_sys = system.n_vars();
1003 
1004  // Update the current_local_solution
1005  {
1006  System & non_const_sys = const_cast<System &>(system);
1007  // We used to simply call non_const_sys.solution->close()
1008  // here, but that is not allowed when the solution vector is
1009  // locked read-only, for example when printing the solution
1010  // during during the middle of a solve... So try to be a bit
1011  // more careful about calling close() unnecessarily.
1012  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
1013  if (!non_const_sys.solution->closed())
1014  non_const_sys.solution->close();
1015  non_const_sys.update();
1016  }
1017 
1018  NumericVector<Number> & sys_soln(*system.current_local_solution);
1019 
1020  // The DOF indices for the finite element
1021  std::vector<dof_id_type> dof_indices;
1022 
1023  // Loop over the variable names and load them in order
1024  for (unsigned int var=0; var < nv_sys; ++var)
1025  {
1026  // Skip this variable if we are not outputting it.
1027  if (!do_output[sys_ctr][var])
1028  continue;
1029 
1030  names.push_back(system.variable_name(var));
1031 
1032  const Variable & variable = system.variable(var);
1033  const DofMap & dof_map = system.get_dof_map();
1034 
1035  for (const auto & elem : _mesh.active_local_element_ptr_range())
1036  if (variable.active_on_subdomain(elem->subdomain_id()))
1037  {
1038  dof_map.dof_indices (elem, dof_indices, var);
1039 
1040  libmesh_assert_equal_to (1, dof_indices.size());
1041 
1042  parallel_soln.set((ne*var_num)+elem->id(), sys_soln(dof_indices[0]));
1043  }
1044 
1045  var_num++;
1046  } // end loop on variables in this system
1047  } // end loop over systems
1048 
1049  parallel_soln.close();
1050  return std::unique_ptr<NumericVector<Number>>(parallel_soln_ptr.release());
1051 }
1052 
1053 
1054 
1056  const std::set<std::string> * system_names) const
1057 {
1058  LOG_SCOPE("build_discontinuous_solution_vector()", "EquationSystems");
1059 
1060  libmesh_assert (this->n_systems());
1061 
1062  const unsigned int dim = _mesh.mesh_dimension();
1063 
1064  // Get the number of variables (nv) by counting the number of variables
1065  // in each system listed in system_names
1066  unsigned int nv = 0;
1067 
1068  {
1069  const_system_iterator pos = _systems.begin();
1070  const const_system_iterator end = _systems.end();
1071 
1072  for (; pos != end; ++pos)
1073  {
1074  // Check current system is listed in system_names, and skip pos if not
1075  bool use_current_system = (system_names == nullptr);
1076  if (!use_current_system)
1077  use_current_system = system_names->count(pos->first);
1078  if (!use_current_system || pos->second->hide_output())
1079  continue;
1080 
1081  const System & system = *(pos->second);
1082  nv += system.n_vars();
1083  }
1084  }
1085 
1086  unsigned int tw=0;
1087 
1088  // get the total weight
1089  for (const auto & elem : _mesh.active_element_ptr_range())
1090  tw += elem->n_nodes();
1091 
1092  // Only if we are on processor zero, allocate the storage
1093  // to hold (number_of_nodes)*(number_of_variables) entries.
1094  if (_mesh.processor_id() == 0)
1095  soln.resize(tw*nv);
1096 
1097  std::vector<Number> sys_soln;
1098 
1099 
1100  unsigned int var_num=0;
1101 
1102  // For each system in this EquationSystems object,
1103  // update the global solution and if we are on processor 0,
1104  // loop over the elements and build the nodal solution
1105  // from the element solution. Then insert this nodal solution
1106  // into the vector passed to build_solution_vector.
1107  {
1108  const_system_iterator pos = _systems.begin();
1109  const const_system_iterator end = _systems.end();
1110 
1111  for (; pos != end; ++pos)
1112  {
1113  // Check current system is listed in system_names, and skip pos if not
1114  bool use_current_system = (system_names == nullptr);
1115  if (!use_current_system)
1116  use_current_system = system_names->count(pos->first);
1117  if (!use_current_system || pos->second->hide_output())
1118  continue;
1119 
1120  const System & system = *(pos->second);
1121  const unsigned int nv_sys = system.n_vars();
1122 
1123  system.update_global_solution (sys_soln, 0);
1124 
1125  if (_mesh.processor_id() == 0)
1126  {
1127  std::vector<Number> elem_soln; // The finite element solution
1128  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
1129  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
1130 
1131  for (unsigned int var=0; var<nv_sys; var++)
1132  {
1133  const FEType & fe_type = system.variable_type(var);
1134  const Variable & var_description = system.variable(var);
1135 
1136  unsigned int nn=0;
1137 
1138  for (auto & elem : _mesh.active_element_ptr_range())
1139  {
1140  if (var_description.active_on_subdomain(elem->subdomain_id()))
1141  {
1142  system.get_dof_map().dof_indices (elem, dof_indices, var);
1143 
1144  elem_soln.resize(dof_indices.size());
1145 
1146  for (std::size_t i=0; i<dof_indices.size(); i++)
1147  elem_soln[i] = sys_soln[dof_indices[i]];
1148 
1150  fe_type,
1151  elem,
1152  elem_soln,
1153  nodal_soln);
1154 
1155 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1156  // infinite elements should be skipped...
1157  if (!elem->infinite())
1158 #endif
1159  {
1160  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
1161 
1162  for (unsigned int n=0; n<elem->n_nodes(); n++)
1163  {
1164  soln[nv*(nn++) + (var + var_num)] +=
1165  nodal_soln[n];
1166  }
1167  }
1168  }
1169  else
1170  nn += elem->n_nodes();
1171  }
1172  }
1173  }
1174 
1175  var_num += nv_sys;
1176  }
1177  }
1178 }
1179 
1180 
1181 
1183  const Real threshold,
1184  const bool verbose) const
1185 {
1186  // safety check, whether we handle at least the same number
1187  // of systems
1188  std::vector<bool> os_result;
1189 
1190  if (this->n_systems() != other_es.n_systems())
1191  {
1192  if (verbose)
1193  {
1194  libMesh::out << " Fatal difference. This system handles "
1195  << this->n_systems() << " systems," << std::endl
1196  << " while the other system handles "
1197  << other_es.n_systems()
1198  << " systems." << std::endl
1199  << " Aborting comparison." << std::endl;
1200  }
1201  return false;
1202  }
1203  else
1204  {
1205  // start comparing each system
1206  const_system_iterator pos = _systems.begin();
1207  const const_system_iterator end = _systems.end();
1208 
1209  for (; pos != end; ++pos)
1210  {
1211  const std::string & sys_name = pos->first;
1212  const System & system = *(pos->second);
1213 
1214  // get the other system
1215  const System & other_system = other_es.get_system (sys_name);
1216 
1217  os_result.push_back (system.compare (other_system, threshold, verbose));
1218 
1219  }
1220 
1221  }
1222 
1223 
1224  // sum up the results
1225  if (os_result.size()==0)
1226  return true;
1227  else
1228  {
1229  bool os_identical;
1230  unsigned int n = 0;
1231  do
1232  {
1233  os_identical = os_result[n];
1234  n++;
1235  }
1236  while (os_identical && n<os_result.size());
1237  return os_identical;
1238  }
1239 }
1240 
1241 
1242 
1243 std::string EquationSystems::get_info () const
1244 {
1245  std::ostringstream oss;
1246 
1247  oss << " EquationSystems\n"
1248  << " n_systems()=" << this->n_systems() << '\n';
1249 
1250  // Print the info for the individual systems
1251  const_system_iterator pos = _systems.begin();
1252  const const_system_iterator end = _systems.end();
1253 
1254  for (; pos != end; ++pos)
1255  oss << pos->second->get_info();
1256 
1257 
1258  // // Possibly print the parameters
1259  // if (!this->parameters.empty())
1260  // {
1261  // oss << " n_parameters()=" << this->n_parameters() << '\n';
1262  // oss << " Parameters:\n";
1263 
1264  // for (const auto & pr : _parameters)
1265  // oss << " "
1266  // << "\""
1267  // << pr.first
1268  // << "\""
1269  // << "="
1270  // << pr.second
1271  // << '\n';
1272  // }
1273 
1274  return oss.str();
1275 }
1276 
1277 
1278 
1279 void EquationSystems::print_info (std::ostream & os) const
1280 {
1281  os << this->get_info()
1282  << std::endl;
1283 }
1284 
1285 
1286 
1287 std::ostream & operator << (std::ostream & os,
1288  const EquationSystems & es)
1289 {
1290  es.print_info(os);
1291  return os;
1292 }
1293 
1294 
1295 
1296 unsigned int EquationSystems::n_vars () const
1297 {
1298  unsigned int tot=0;
1299 
1300  const_system_iterator pos = _systems.begin();
1301  const const_system_iterator end = _systems.end();
1302 
1303  for (; pos != end; ++pos)
1304  tot += pos->second->n_vars();
1305 
1306  return tot;
1307 }
1308 
1309 
1310 
1311 std::size_t EquationSystems::n_dofs () const
1312 {
1313  std::size_t tot=0;
1314 
1315  const_system_iterator pos = _systems.begin();
1316  const const_system_iterator end = _systems.end();
1317 
1318  for (; pos != end; ++pos)
1319  tot += pos->second->n_dofs();
1320 
1321  return tot;
1322 }
1323 
1324 
1325 
1326 
1328 {
1329  std::size_t tot=0;
1330 
1331  const_system_iterator pos = _systems.begin();
1332  const const_system_iterator end = _systems.end();
1333 
1334  for (; pos != end; ++pos)
1335  tot += pos->second->n_active_dofs();
1336 
1337  return tot;
1338 }
1339 
1340 
1342 {
1343  // All the nodes
1344  for (auto & node : _mesh.node_ptr_range())
1345  node->add_system();
1346 
1347  // All the elements
1348  for (auto & elem : _mesh.element_ptr_range())
1349  elem->add_system();
1350 }
1351 
1353 {
1354  this->get_system(sys_num).get_dof_map().remove_default_ghosting();
1355 }
1356 
1357 } // namespace libMesh
EquationSystems(MeshBase &mesh)
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
double abs(double a)
Manages multiples systems of equations.
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
void update_global_solution(std::vector< Number > &global_soln) const
Definition: system.C:642
unsigned int n_systems() const
const Variable & variable(unsigned int var) const
Definition: system.h:2133
Specifies parameters for parameter sensitivity calculations.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1930
std::unique_ptr< NumericVector< Number > > build_parallel_elemental_solution_vector(std::vector< std::string > &names) const
std::size_t n_dofs() const
Used to specify quantities of interest in a simulation.
Definition: qoi_set.h:45
processor_id_type size() const
Definition: communicator.h:175
void add_default_ghosting()
Definition: dof_map.C:1800
const T_sys & get_system(const std::string &name) const
virtual void allgather()
Definition: mesh_base.h:183
void print_info(std::ostream &os=libMesh::out) const
MeshBase & mesh
static FEFieldType field_type(const FEType &fe_type)
std::unique_ptr< NumericVector< Number > > build_parallel_solution_vector(const std::set< std::string > *system_names=nullptr) const
const Parallel::Communicator & comm() const
virtual void enable_default_ghosting(bool enable)
static const Real TOLERANCE
IterBase * end
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
const Parallel::Communicator & _communicator
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Base class for Mesh.
Definition: mesh_base.h:77
virtual std::string get_info() const
virtual void prolong_vectors()
Definition: system.C:380
signed char & underrefined_boundary_limit()
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
void get_vars_active_subdomains(const std::vector< std::string > &names, std::vector< std::set< subdomain_id_type >> &vars_active_subdomains) const
virtual void sensitivity_solve(const ParameterVector &parameters)
std::map< std::string, System * >::const_iterator const_system_iterator
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:154
unsigned int number() const
Definition: system.h:2025
A variable which is solved for in a System of equations.
Definition: variable.h:49
const std::set< subdomain_id_type > & active_subdomains() const
Definition: variable.h:150
Responsible for mesh refinement algorithms and data.
virtual SimpleRange< element_iterator > element_ptr_range()=0
dof_id_type numeric_index_type
Definition: id_types.h:92
processor_id_type rank() const
Definition: communicator.h:173
virtual node_iterator local_nodes_end()=0
std::map< std::string, System * >::iterator system_iterator
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
void delete_system(const std::string &name)
virtual SimpleRange< node_iterator > node_ptr_range()=0
virtual dof_id_type max_elem_id() const =0
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
void distribute_dofs(MeshBase &)
Definition: dof_map.C:902
virtual System & add_system(const std::string &system_type, const std::string &name)
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2153
virtual void restrict_vectors()
Definition: system.C:324
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
unsigned char & face_level_mismatch_limit()
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()
Definition: system.C:397
virtual void close()=0
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
virtual numeric_index_type first_local_index() const =0
virtual node_iterator local_nodes_begin()=0
void get_solution(std::vector< Number > &soln, std::vector< std::string > &names) const
virtual void update()
Definition: system.C:408
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2183
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool contract()=0
T & set(const std::string &)
Definition: parameters.h:464
void build_solution_vector(std::vector< Number > &soln, const std::string &system_name, const std::string &variable_name="all_vars") const
virtual void clear()
Definition: parameters.h:317
signed char & overrefined_boundary_limit()
void remove_default_ghosting()
Definition: dof_map.C:1792
void _remove_default_ghosting(unsigned int sys_num)
const MeshBase & get_mesh() const
std::size_t n_active_dofs() const
unsigned int mesh_dimension() const
Definition: mesh_base.C:126
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1535
virtual bool compare(const System &other_system, const Real threshold, const bool verbose) const
Definition: system.C:514
virtual void set(const numeric_index_type i, const T value)=0
void prepare_send_list()
Definition: dof_map.C:1639
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)
Definition: fe_interface.C:550
virtual bool compare(const EquationSystems &other_es, const Real threshold, const bool verbose) const
unsigned int n_vars() const
Definition: system.h:2105
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
processor_id_type processor_id() const
OStreamProxy out(std::cout)
void build_discontinuous_solution_vector(std::vector< Number > &soln, const std::set< std::string > *system_names=nullptr) const
const DofMap & get_dof_map() const
Definition: system.h:2049
std::ostream & operator<<(std::ostream &os, const FEAbstract &fe)
Definition: fe_abstract.C:809
virtual numeric_index_type last_local_index() const =0
void build_elemental_solution_vector(std::vector< Number > &soln, std::vector< std::string > &names) const
uint8_t dof_id_type
Definition: id_types.h:64
std::map< std::string, System * > _systems