libMesh
equation_systems.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/explicit_system.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/frequency_system.h"
26 #include "libmesh/linear_implicit_system.h"
27 #include "libmesh/mesh_refinement.h"
28 #include "libmesh/newmark_system.h"
29 #include "libmesh/nonlinear_implicit_system.h"
30 #include "libmesh/rb_construction.h"
31 #include "libmesh/transient_rb_construction.h"
32 #include "libmesh/eigen_system.h"
33 #include "libmesh/parallel.h"
34 #include "libmesh/transient_system.h"
35 #include "libmesh/dof_map.h"
36 #include "libmesh/mesh_base.h"
37 #include "libmesh/elem.h"
38 #include "libmesh/libmesh_logging.h"
39 
40 // Include the systems before this one to avoid
41 // overlapping forward declarations.
42 #include "libmesh/equation_systems.h"
43 
44 namespace libMesh
45 {
46 
47 // Forward Declarations
48 
49 
50 
51 
52 // ------------------------------------------------------------
53 // EquationSystems class implementation
55  ParallelObject (m),
56  _mesh (m)
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  this->_refine_in_reinit = true; // default value
62 }
63 
64 
65 
67 {
68  this->clear ();
69 }
70 
71 
72 
74 {
75  // Clear any additional parameters
76  parameters.clear ();
77 
78  // clear the systems. We must delete them
79  // since we newed them!
80  while (!_systems.empty())
81  {
82  system_iterator pos = _systems.begin();
83 
84  System * sys = pos->second;
85  delete sys;
86  sys = libmesh_nullptr;
87 
88  _systems.erase (pos);
89  }
90 }
91 
92 
93 
95 {
96  const unsigned int n_sys = this->n_systems();
97 
98  libmesh_assert_not_equal_to (n_sys, 0);
99 
100  // Tell all the \p DofObject entities how many systems
101  // there are.
102  for (auto & node : _mesh.node_ptr_range())
103  node->set_n_systems(n_sys);
104 
105  for (auto & elem : _mesh.element_ptr_range())
106  elem->set_n_systems(n_sys);
107 
108  for (unsigned int i=0; i != this->n_systems(); ++i)
109  this->get_system(i).init();
110 
111 #ifdef LIBMESH_ENABLE_AMR
112  MeshRefinement mesh_refine(_mesh);
113  mesh_refine.clean_refinement_flags();
114 #endif
115 }
116 
117 
118 
120 {
121  parallel_object_only();
122 
123  const unsigned int n_sys = this->n_systems();
124  libmesh_assert_not_equal_to (n_sys, 0);
125 
126  // We may have added new systems since our last
127  // EquationSystems::(re)init call
128  bool _added_new_systems = false;
129  for (unsigned int i=0; i != n_sys; ++i)
130  if (!this->get_system(i).is_initialized())
131  _added_new_systems = true;
132 
133  if (_added_new_systems)
134  {
135  // Our DofObjects will need space for the additional systems
136  for (auto & node : _mesh.node_ptr_range())
137  node->set_n_systems(n_sys);
138 
139  for (auto & elem : _mesh.element_ptr_range())
140  elem->set_n_systems(n_sys);
141 
142  // And any new systems will need initialization
143  for (unsigned int i=0; i != n_sys; ++i)
144  if (!this->get_system(i).is_initialized())
145  this->get_system(i).init();
146  }
147 
148 
149  // We used to assert that all nodes and elements *already* had
150  // n_systems() properly set; however this is false in the case where
151  // user code has manually added nodes and/or elements to an
152  // already-initialized system.
153 
154  // Make sure all the \p DofObject entities know how many systems
155  // there are.
156  {
157  // All the nodes
158  for (auto & node : _mesh.node_ptr_range())
159  node->set_n_systems(this->n_systems());
160 
161  // All the elements
162  for (auto & elem : _mesh.element_ptr_range())
163  elem->set_n_systems(this->n_systems());
164  }
165 
166  // Localize each system's vectors
167  for (unsigned int i=0; i != this->n_systems(); ++i)
168  this->get_system(i).re_update();
169 
170 #ifdef LIBMESH_ENABLE_AMR
171 
172  bool dof_constraints_created = false;
173  bool mesh_changed = false;
174 
175  // FIXME: For backwards compatibility, assume
176  // refine_and_coarsen_elements or refine_uniformly have already
177  // been called
178  {
179  for (unsigned int i=0; i != this->n_systems(); ++i)
180  {
181  System & sys = this->get_system(i);
182 
183  // Even if the system doesn't have any variables in it we want
184  // consistent behavior; e.g. distribute_dofs should have the
185  // opportunity to count up zero dofs on each processor.
186  //
187  // Who's been adding zero-var systems anyway, outside of my
188  // unit tests? - RHS
189  // if (!sys.n_vars())
190  // continue;
191 
193 
194  // Recreate any user or internal constraints
195  sys.reinit_constraints();
196 
197  sys.prolong_vectors();
198  }
199  mesh_changed = true;
200  dof_constraints_created = true;
201  }
202 
203  if (this->_refine_in_reinit)
204  {
205  // Don't override any user refinement settings
206  MeshRefinement mesh_refine(_mesh);
207  mesh_refine.face_level_mismatch_limit() = 0; // unlimited
208  mesh_refine.overrefined_boundary_limit() = -1; // unlimited
209  mesh_refine.underrefined_boundary_limit() = -1; // unlimited
210 
211  // Try to coarsen the mesh, then restrict each system's vectors
212  // if necessary
213  if (mesh_refine.coarsen_elements())
214  {
215  for (unsigned int i=0; i != this->n_systems(); ++i)
216  {
217  System & sys = this->get_system(i);
218  if (!dof_constraints_created)
219  {
221  sys.reinit_constraints();
222  }
223  sys.restrict_vectors();
224  }
225  mesh_changed = true;
226  dof_constraints_created = true;
227  }
228 
229  // Once vectors are all restricted, we can delete
230  // children of coarsened elements
231  if (mesh_changed)
232  this->get_mesh().contract();
233 
234  // Try to refine the mesh, then prolong each system's vectors
235  // if necessary
236  if (mesh_refine.refine_elements())
237  {
238  for (unsigned int i=0; i != this->n_systems(); ++i)
239  {
240  System & sys = this->get_system(i);
241  if (!dof_constraints_created)
242  {
244  sys.reinit_constraints();
245  }
246  sys.prolong_vectors();
247  }
248  mesh_changed = true;
249  // dof_constraints_created = true;
250  }
251  }
252 
253  // If the mesh has changed, systems will need to create new dof
254  // constraints and update their global solution vectors
255  if (mesh_changed)
256  {
257  for (unsigned int i=0; i != this->n_systems(); ++i)
258  this->get_system(i).reinit();
259  }
260 #endif // #ifdef LIBMESH_ENABLE_AMR
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 (unsigned int i=0; i != this->n_systems(); ++i)
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  dof_map.prepare_send_list();
298  }
299 }
300 
301 
302 
303 
305 {
306  LOG_SCOPE("update()", "EquationSystems");
307 
308  // Localize each system's vectors
309  for (unsigned int i=0; i != this->n_systems(); ++i)
310  this->get_system(i).update();
311 }
312 
313 
314 
315 System & EquationSystems::add_system (const std::string & sys_type,
316  const std::string & name)
317 {
318  // If the user already built a system with this name, we'll
319  // trust them and we'll use it. That way they can pre-add
320  // non-standard derived system classes, and if their restart file
321  // has some non-standard sys_type we won't throw an error.
322  if (_systems.count(name))
323  {
324  return this->get_system(name);
325  }
326  // Build a basic System
327  else if (sys_type == "Basic")
328  this->add_system<System> (name);
329 
330  // Build a Newmark system
331  else if (sys_type == "Newmark")
332  this->add_system<NewmarkSystem> (name);
333 
334  // Build an Explicit system
335  else if ((sys_type == "Explicit"))
336  this->add_system<ExplicitSystem> (name);
337 
338  // Build an Implicit system
339  else if ((sys_type == "Implicit") ||
340  (sys_type == "Steady" ))
341  this->add_system<ImplicitSystem> (name);
342 
343  // build a transient implicit linear system
344  else if ((sys_type == "Transient") ||
345  (sys_type == "TransientImplicit") ||
346  (sys_type == "TransientLinearImplicit"))
347  this->add_system<TransientLinearImplicitSystem> (name);
348 
349  // build a transient implicit nonlinear system
350  else if (sys_type == "TransientNonlinearImplicit")
351  this->add_system<TransientNonlinearImplicitSystem> (name);
352 
353  // build a transient explicit system
354  else if (sys_type == "TransientExplicit")
355  this->add_system<TransientExplicitSystem> (name);
356 
357  // build a linear implicit system
358  else if (sys_type == "LinearImplicit")
359  this->add_system<LinearImplicitSystem> (name);
360 
361  // build a nonlinear implicit system
362  else if (sys_type == "NonlinearImplicit")
363  this->add_system<NonlinearImplicitSystem> (name);
364 
365  // build a Reduced Basis Construction system
366  else if (sys_type == "RBConstruction")
367  this->add_system<RBConstruction> (name);
368 
369  // build a transient Reduced Basis Construction system
370  else if (sys_type == "TransientRBConstruction")
371  this->add_system<TransientRBConstruction> (name);
372 
373 #ifdef LIBMESH_HAVE_SLEPC
374  // build an eigen system
375  else if (sys_type == "Eigen")
376  this->add_system<EigenSystem> (name);
377  else if (sys_type == "TransientEigenSystem")
378  this->add_system<TransientEigenSystem> (name);
379 #endif
380 
381 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
382  // build a frequency system
383  else if (sys_type == "Frequency")
384  this->add_system<FrequencySystem> (name);
385 #endif
386 
387  else
388  libmesh_error_msg("ERROR: Unknown system type: " << sys_type);
389 
390  // Return a reference to the new system
391  //return (*this)(name);
392  return this->get_system(name);
393 }
394 
395 
396 
397 
398 
399 
400 #ifdef LIBMESH_ENABLE_DEPRECATED
401 void EquationSystems::delete_system (const std::string & name)
402 {
403  libmesh_deprecated();
404 
405  if (!_systems.count(name))
406  libmesh_error_msg("ERROR: no system named " << name);
407 
408  delete _systems[name];
409 
410  _systems.erase (name);
411 }
412 #endif
413 
414 
415 
417 {
418  libmesh_assert (this->n_systems());
419 
420  for (unsigned int i=0; i != this->n_systems(); ++i)
421  this->get_system(i).solve();
422 }
423 
424 
425 
427 {
428  libmesh_assert (this->n_systems());
429 
430  for (unsigned int i=0; i != this->n_systems(); ++i)
431  this->get_system(i).sensitivity_solve(parameters_in);
432 }
433 
434 
435 
436 void EquationSystems::adjoint_solve (const QoISet & qoi_indices)
437 {
438  libmesh_assert (this->n_systems());
439 
440  for (unsigned int i=this->n_systems(); i != 0; --i)
441  this->get_system(i-1).adjoint_solve(qoi_indices);
442 }
443 
444 
445 
446 void EquationSystems::build_variable_names (std::vector<std::string> & var_names,
447  const FEType * type,
448  const std::set<std::string> * system_names) const
449 {
450  unsigned int var_num=0;
451 
452  const_system_iterator pos = _systems.begin();
453  const const_system_iterator end = _systems.end();
454 
455  // Need to size var_names by scalar variables plus all the
456  // vector components for all the vector variables
457  //Could this be replaced by a/some convenience methods?[PB]
458  {
459  unsigned int n_scalar_vars = 0;
460  unsigned int n_vector_vars = 0;
461 
462  for (; pos != end; ++pos)
463  {
464  // Check current system is listed in system_names, and skip pos if not
465  bool use_current_system = (system_names == libmesh_nullptr);
466  if (!use_current_system)
467  use_current_system = system_names->count(pos->first);
468  if (!use_current_system)
469  continue;
470 
471  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
472  {
473  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
474  n_vector_vars++;
475  else
476  n_scalar_vars++;
477  }
478  }
479 
480  // Here, we're assuming the number of vector components is the same
481  // as the mesh dimension. Will break for mixed dimension meshes.
482  unsigned int dim = this->get_mesh().mesh_dimension();
483  unsigned int nv = n_scalar_vars + dim*n_vector_vars;
484 
485  // We'd better not have more than dim*his->n_vars() (all vector variables)
486  libmesh_assert_less_equal ( nv, dim*this->n_vars() );
487 
488  // Here, we're assuming the number of vector components is the same
489  // as the mesh dimension. Will break for mixed dimension meshes.
490 
491  var_names.resize( nv );
492  }
493 
494  // reset
495  pos = _systems.begin();
496 
497  for (; pos != end; ++pos)
498  {
499  // Check current system is listed in system_names, and skip pos if not
500  bool use_current_system = (system_names == libmesh_nullptr);
501  if (!use_current_system)
502  use_current_system = system_names->count(pos->first);
503  if (!use_current_system)
504  continue;
505 
506  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
507  {
508  std::string var_name = pos->second->variable_name(vn);
509  FEType fe_type = pos->second->variable_type(vn);
510 
511  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type);
512 
513  // Filter on the type if requested
514  if (type == libmesh_nullptr || (type && *type == fe_type))
515  {
516  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
517  {
518  switch(n_vec_dim)
519  {
520  case 0:
521  case 1:
522  var_names[var_num++] = var_name;
523  break;
524  case 2:
525  var_names[var_num++] = var_name+"_x";
526  var_names[var_num++] = var_name+"_y";
527  break;
528  case 3:
529  var_names[var_num++] = var_name+"_x";
530  var_names[var_num++] = var_name+"_y";
531  var_names[var_num++] = var_name+"_z";
532  break;
533  default:
534  libmesh_error_msg("Invalid dim in build_variable_names");
535  }
536  }
537  else
538  var_names[var_num++] = var_name;
539  }
540  }
541  }
542  // Now resize again in case we filtered any names
543  var_names.resize(var_num);
544 }
545 
546 
547 
548 void EquationSystems::build_solution_vector (std::vector<Number> &,
549  const std::string &,
550  const std::string &) const
551 {
552  // TODO:[BSK] re-implement this from the method below
553  libmesh_not_implemented();
554 }
555 
556 
557 
558 
560 EquationSystems::build_parallel_solution_vector(const std::set<std::string> * system_names) const
561 {
562  LOG_SCOPE("build_parallel_solution_vector()", "EquationSystems");
563 
564  // This function must be run on all processors at once
565  parallel_object_only();
566 
567  const unsigned int dim = _mesh.mesh_dimension();
568  const dof_id_type nn = _mesh.n_nodes();
569 
570  // We'd better have a contiguous node numbering
571  libmesh_assert_equal_to (nn, _mesh.max_node_id());
572 
573  // allocate storage to hold
574  // (number_of_nodes)*(number_of_variables) entries.
575  // We have to differentiate between between scalar and vector
576  // variables. We intercept vector variables and treat each
577  // component as a scalar variable (consistently with build_solution_names).
578 
579  unsigned int nv = 0;
580 
581  //Could this be replaced by a/some convenience methods?[PB]
582  {
583  unsigned int n_scalar_vars = 0;
584  unsigned int n_vector_vars = 0;
585  const_system_iterator pos = _systems.begin();
586  const const_system_iterator end = _systems.end();
587 
588  for (; pos != end; ++pos)
589  {
590  // Check current system is listed in system_names, and skip pos if not
591  bool use_current_system = (system_names == libmesh_nullptr);
592  if (!use_current_system)
593  use_current_system = system_names->count(pos->first);
594  if (!use_current_system)
595  continue;
596 
597  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
598  {
599  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
600  n_vector_vars++;
601  else
602  n_scalar_vars++;
603  }
604  }
605  // Here, we're assuming the number of vector components is the same
606  // as the mesh dimension. Will break for mixed dimension meshes.
607  nv = n_scalar_vars + dim*n_vector_vars;
608  }
609 
610  // Get the number of local nodes
611  dof_id_type n_local_nodes = cast_int<dof_id_type>
614 
615  // Create a NumericVector to hold the parallel solution
617  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
618  parallel_soln.init(nn*nv, n_local_nodes*nv, false, PARALLEL);
619 
620  // Create a NumericVector to hold the "repeat_count" for each node - this is essentially
621  // the number of elements contributing to that node's value
623  NumericVector<Number> & repeat_count = *repeat_count_ptr;
624  repeat_count.init(nn*nv, n_local_nodes*nv, false, PARALLEL);
625 
626  repeat_count.close();
627 
628  unsigned int var_num=0;
629 
630  // For each system in this EquationSystems object,
631  // update the global solution and if we are on processor 0,
632  // loop over the elements and build the nodal solution
633  // from the element solution. Then insert this nodal solution
634  // into the vector passed to build_solution_vector.
635  const_system_iterator pos = _systems.begin();
636  const const_system_iterator end = _systems.end();
637 
638  for (; pos != end; ++pos)
639  {
640  // Check current system is listed in system_names, and skip pos if not
641  bool use_current_system = (system_names == libmesh_nullptr);
642  if (!use_current_system)
643  use_current_system = system_names->count(pos->first);
644  if (!use_current_system)
645  continue;
646 
647  const System & system = *(pos->second);
648  const unsigned int nv_sys = system.n_vars();
649  const unsigned int sys_num = system.number();
650 
651  //Could this be replaced by a/some convenience methods?[PB]
652  unsigned int n_scalar_vars = 0;
653  unsigned int n_vector_vars = 0;
654  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
655  {
656  if (FEInterface::field_type(pos->second->variable_type(vn)) == TYPE_VECTOR)
657  n_vector_vars++;
658  else
659  n_scalar_vars++;
660  }
661 
662  // Here, we're assuming the number of vector components is the same
663  // as the mesh dimension. Will break for mixed dimension meshes.
664  unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
665 
666  // Update the current_local_solution
667  {
668  System & non_const_sys = const_cast<System &>(system);
669  // We used to simply call non_const_sys.solution->close()
670  // here, but that is not allowed when the solution vector is
671  // locked read-only, for example when printing the solution
672  // during the middle of a solve... So try to be a bit
673  // more careful about calling close() unnecessarily.
674  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
675  if (!non_const_sys.solution->closed())
676  non_const_sys.solution->close();
677  non_const_sys.update();
678  }
679 
680  NumericVector<Number> & sys_soln(*system.current_local_solution);
681 
682  std::vector<Number> elem_soln; // The finite element solution
683  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
684  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
685 
686  unsigned var_inc = 0;
687  for (unsigned int var=0; var<nv_sys; var++)
688  {
689  const FEType & fe_type = system.variable_type(var);
690  const Variable & var_description = system.variable(var);
691  const DofMap & dof_map = system.get_dof_map();
692 
693  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type );
694 
695  for (const auto & elem : _mesh.active_local_element_ptr_range())
696  {
697  if (var_description.active_on_subdomain(elem->subdomain_id()))
698  {
699  dof_map.dof_indices (elem, dof_indices, var);
700 
701  elem_soln.resize(dof_indices.size());
702 
703  for (std::size_t i=0; i<dof_indices.size(); i++)
704  elem_soln[i] = sys_soln(dof_indices[i]);
705 
707  fe_type,
708  elem,
709  elem_soln,
710  nodal_soln);
711 
712 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
713  // infinite elements should be skipped...
714  if (!elem->infinite())
715 #endif
716  {
717  libmesh_assert_equal_to (nodal_soln.size(), n_vec_dim*elem->n_nodes());
718 
719  for (unsigned int n=0; n<elem->n_nodes(); n++)
720  {
721  for (unsigned int d=0; d < n_vec_dim; d++)
722  {
723  // For vector-valued elements, all components are in nodal_soln. For each
724  // node, the components are stored in order, i.e. node_0 -> s0_x, s0_y, s0_z
725  parallel_soln.add(nv*(elem->node_id(n)) + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
726 
727  // Increment the repeat count for this position
728  repeat_count.add(nv*(elem->node_id(n)) + (var_inc+d + var_num), 1);
729  }
730  }
731  }
732  }
733  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:
734  for (unsigned int n=0; n<elem->n_nodes(); n++)
735  // Only do this if this variable has NO DoFs at this node... it might have some from an adjoining element...
736  if (!elem->node_ptr(n)->n_dofs(sys_num, var))
737  for (unsigned int d=0; d < n_vec_dim; d++)
738  repeat_count.add(nv*(elem->node_id(n)) + (var+d + var_num), 1);
739 
740  } // end loop over elements
741  var_inc += n_vec_dim;
742  } // end loop on variables in this system
743 
744  var_num += nv_sys_split;
745  } // end loop over systems
746 
747  parallel_soln.close();
748  repeat_count.close();
749 
750  // Divide to get the average value at the nodes
751  parallel_soln /= repeat_count;
752 
753  return UniquePtr<NumericVector<Number>>(parallel_soln_ptr.release());
754 }
755 
756 
757 
758 void EquationSystems::build_solution_vector (std::vector<Number> & soln,
759  const std::set<std::string> * system_names) const
760 {
761  LOG_SCOPE("build_solution_vector()", "EquationSystems");
762 
763  // Call the parallel implementation
764  UniquePtr<NumericVector<Number>> parallel_soln =
765  this->build_parallel_solution_vector(system_names);
766 
767  // Localize the NumericVector into the provided std::vector.
768  parallel_soln->localize_to_one(soln);
769 }
770 
771 
772 
773 void EquationSystems::get_solution (std::vector<Number> & soln,
774  std::vector<std::string> & names) const
775 {
776  // This function must be run on all processors at once
777  parallel_object_only();
778 
779  libmesh_assert (this->n_systems());
780 
781  const dof_id_type ne = _mesh.n_elem();
782 
783  libmesh_assert_equal_to (ne, _mesh.max_elem_id());
784 
785  // If the names vector has entries, we will only populate the soln vector
786  // with names included in that list. Note: The names vector may be
787  // reordered upon exiting this function
788  std::vector<std::string> filter_names = names;
789  bool is_filter_names = !filter_names.empty();
790 
791  soln.clear();
792  names.clear();
793 
794  const FEType type(CONSTANT, MONOMIAL);
795 
796  dof_id_type nv = 0;
797 
798  // Find the total number of variables to output
799  std::vector<std::vector<unsigned>> do_output(_systems.size());
800  {
801  const_system_iterator pos = _systems.begin();
802  const const_system_iterator end = _systems.end();
803  unsigned sys_ctr = 0;
804 
805  for (; pos != end; ++pos, ++sys_ctr)
806  {
807  const System & system = *(pos->second);
808  const unsigned int nv_sys = system.n_vars();
809 
810  do_output[sys_ctr].resize(nv_sys);
811 
812  for (unsigned int var=0; var < nv_sys; ++var)
813  {
814  if (system.variable_type(var) != type ||
815  (is_filter_names && std::find(filter_names.begin(), filter_names.end(), system.variable_name(var)) == filter_names.end()))
816  continue;
817 
818  // Otherwise, this variable should be output
819  nv++;
820  do_output[sys_ctr][var] = 1;
821  }
822  }
823  }
824 
825  // If there are no variables to write out don't do anything...
826  if (!nv)
827  return;
828 
829  // We can handle the case where there are NULLs in the Elem vector
830  // by just having extra zeros in the solution vector.
831  numeric_index_type parallel_soln_global_size = ne*nv;
832 
833  numeric_index_type div = parallel_soln_global_size / this->n_processors();
834  numeric_index_type mod = parallel_soln_global_size % this->n_processors();
835 
836  // Initialize all processors to the average size.
837  numeric_index_type parallel_soln_local_size = div;
838 
839  // The first "mod" processors get an extra entry.
840  if (this->processor_id() < mod)
841  parallel_soln_local_size = div+1;
842 
843  // Create a NumericVector to hold the parallel solution
845  NumericVector<Number> & parallel_soln = *parallel_soln_ptr;
846  parallel_soln.init(parallel_soln_global_size,
847  parallel_soln_local_size,
848  /*fast=*/false,
849  /*ParallelType=*/PARALLEL);
850 
851  dof_id_type var_num = 0;
852 
853  // For each system in this EquationSystems object,
854  // update the global solution and collect the
855  // CONSTANT MONOMIALs. The entries are in variable-major
856  // format.
857  const_system_iterator pos = _systems.begin();
858  const const_system_iterator end = _systems.end();
859  unsigned sys_ctr = 0;
860 
861  for (; pos != end; ++pos, ++sys_ctr)
862  {
863  const System & system = *(pos->second);
864  const unsigned int nv_sys = system.n_vars();
865 
866  // Update the current_local_solution
867  {
868  System & non_const_sys = const_cast<System &>(system);
869  // We used to simply call non_const_sys.solution->close()
870  // here, but that is not allowed when the solution vector is
871  // locked read-only, for example when printing the solution
872  // during during the middle of a solve... So try to be a bit
873  // more careful about calling close() unnecessarily.
874  libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
875  if (!non_const_sys.solution->closed())
876  non_const_sys.solution->close();
877  non_const_sys.update();
878  }
879 
880  NumericVector<Number> & sys_soln(*system.current_local_solution);
881 
882  // The DOF indices for the finite element
883  std::vector<dof_id_type> dof_indices;
884 
885  // Loop over the variable names and load them in order
886  for (unsigned int var=0; var < nv_sys; ++var)
887  {
888  // Skip this variable if we are not outputting it.
889  if (!do_output[sys_ctr][var])
890  continue;
891 
892  names.push_back(system.variable_name(var));
893 
894  const Variable & variable = system.variable(var);
895  const DofMap & dof_map = system.get_dof_map();
896 
897  for (const auto & elem : _mesh.active_local_element_ptr_range())
898  if (variable.active_on_subdomain(elem->subdomain_id()))
899  {
900  dof_map.dof_indices (elem, dof_indices, var);
901 
902  libmesh_assert_equal_to (1, dof_indices.size());
903 
904  parallel_soln.set((ne*var_num)+elem->id(), sys_soln(dof_indices[0]));
905  }
906 
907  var_num++;
908  } // end loop on variables in this system
909  } // end loop over systems
910 
911  parallel_soln.close();
912 
913  parallel_soln.localize_to_one(soln);
914 }
915 
916 
917 
919  const std::set<std::string> * system_names) const
920 {
921  LOG_SCOPE("build_discontinuous_solution_vector()", "EquationSystems");
922 
923  libmesh_assert (this->n_systems());
924 
925  const unsigned int dim = _mesh.mesh_dimension();
926 
927  // Get the number of variables (nv) by counting the number of variables
928  // in each system listed in system_names
929  unsigned int nv = 0;
930 
931  {
932  const_system_iterator pos = _systems.begin();
933  const const_system_iterator end = _systems.end();
934 
935  for (; pos != end; ++pos)
936  {
937  // Check current system is listed in system_names, and skip pos if not
938  bool use_current_system = (system_names == libmesh_nullptr);
939  if (!use_current_system)
940  use_current_system = system_names->count(pos->first);
941  if (!use_current_system)
942  continue;
943 
944  const System & system = *(pos->second);
945  nv += system.n_vars();
946  }
947  }
948 
949  unsigned int tw=0;
950 
951  // get the total weight
952  for (const auto & elem : _mesh.active_element_ptr_range())
953  tw += elem->n_nodes();
954 
955  // Only if we are on processor zero, allocate the storage
956  // to hold (number_of_nodes)*(number_of_variables) entries.
957  if (_mesh.processor_id() == 0)
958  soln.resize(tw*nv);
959 
960  std::vector<Number> sys_soln;
961 
962 
963  unsigned int var_num=0;
964 
965  // For each system in this EquationSystems object,
966  // update the global solution and if we are on processor 0,
967  // loop over the elements and build the nodal solution
968  // from the element solution. Then insert this nodal solution
969  // into the vector passed to build_solution_vector.
970  {
971  const_system_iterator pos = _systems.begin();
972  const const_system_iterator end = _systems.end();
973 
974  for (; pos != end; ++pos)
975  {
976  // Check current system is listed in system_names, and skip pos if not
977  bool use_current_system = (system_names == libmesh_nullptr);
978  if (!use_current_system)
979  use_current_system = system_names->count(pos->first);
980  if (!use_current_system)
981  continue;
982 
983  const System & system = *(pos->second);
984  const unsigned int nv_sys = system.n_vars();
985 
986  system.update_global_solution (sys_soln, 0);
987 
988  if (_mesh.processor_id() == 0)
989  {
990  std::vector<Number> elem_soln; // The finite element solution
991  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
992  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
993 
994  for (unsigned int var=0; var<nv_sys; var++)
995  {
996  const FEType & fe_type = system.variable_type(var);
997 
998  unsigned int nn=0;
999 
1000  for (auto & elem : _mesh.active_element_ptr_range())
1001  {
1002  system.get_dof_map().dof_indices (elem, dof_indices, var);
1003 
1004  elem_soln.resize(dof_indices.size());
1005 
1006  for (std::size_t i=0; i<dof_indices.size(); i++)
1007  elem_soln[i] = sys_soln[dof_indices[i]];
1008 
1010  fe_type,
1011  elem,
1012  elem_soln,
1013  nodal_soln);
1014 
1015 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1016  // infinite elements should be skipped...
1017  if (!elem->infinite())
1018 #endif
1019  {
1020  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
1021 
1022  for (unsigned int n=0; n<elem->n_nodes(); n++)
1023  {
1024  soln[nv*(nn++) + (var + var_num)] +=
1025  nodal_soln[n];
1026  }
1027  }
1028  }
1029  }
1030  }
1031 
1032  var_num += nv_sys;
1033  }
1034  }
1035 }
1036 
1037 
1038 
1040  const Real threshold,
1041  const bool verbose) const
1042 {
1043  // safety check, whether we handle at least the same number
1044  // of systems
1045  std::vector<bool> os_result;
1046 
1047  if (this->n_systems() != other_es.n_systems())
1048  {
1049  if (verbose)
1050  {
1051  libMesh::out << " Fatal difference. This system handles "
1052  << this->n_systems() << " systems," << std::endl
1053  << " while the other system handles "
1054  << other_es.n_systems()
1055  << " systems." << std::endl
1056  << " Aborting comparison." << std::endl;
1057  }
1058  return false;
1059  }
1060  else
1061  {
1062  // start comparing each system
1063  const_system_iterator pos = _systems.begin();
1064  const const_system_iterator end = _systems.end();
1065 
1066  for (; pos != end; ++pos)
1067  {
1068  const std::string & sys_name = pos->first;
1069  const System & system = *(pos->second);
1070 
1071  // get the other system
1072  const System & other_system = other_es.get_system (sys_name);
1073 
1074  os_result.push_back (system.compare (other_system, threshold, verbose));
1075 
1076  }
1077 
1078  }
1079 
1080 
1081  // sum up the results
1082  if (os_result.size()==0)
1083  return true;
1084  else
1085  {
1086  bool os_identical;
1087  unsigned int n = 0;
1088  do
1089  {
1090  os_identical = os_result[n];
1091  n++;
1092  }
1093  while (os_identical && n<os_result.size());
1094  return os_identical;
1095  }
1096 }
1097 
1098 
1099 
1100 std::string EquationSystems::get_info () const
1101 {
1102  std::ostringstream oss;
1103 
1104  oss << " EquationSystems\n"
1105  << " n_systems()=" << this->n_systems() << '\n';
1106 
1107  // Print the info for the individual systems
1108  const_system_iterator pos = _systems.begin();
1109  const const_system_iterator end = _systems.end();
1110 
1111  for (; pos != end; ++pos)
1112  oss << pos->second->get_info();
1113 
1114 
1115  // // Possibly print the parameters
1116  // if (!this->parameters.empty())
1117  // {
1118  // oss << " n_parameters()=" << this->n_parameters() << '\n';
1119  // oss << " Parameters:\n";
1120 
1121  // for (std::map<std::string, Real>::const_iterator
1122  // param = _parameters.begin(); param != _parameters.end();
1123  // ++param)
1124  // oss << " "
1125  // << "\""
1126  // << param->first
1127  // << "\""
1128  // << "="
1129  // << param->second
1130  // << '\n';
1131  // }
1132 
1133  return oss.str();
1134 }
1135 
1136 
1137 
1138 void EquationSystems::print_info (std::ostream & os) const
1139 {
1140  os << this->get_info()
1141  << std::endl;
1142 }
1143 
1144 
1145 
1146 std::ostream & operator << (std::ostream & os,
1147  const EquationSystems & es)
1148 {
1149  es.print_info(os);
1150  return os;
1151 }
1152 
1153 
1154 
1155 unsigned int EquationSystems::n_vars () const
1156 {
1157  unsigned int tot=0;
1158 
1159  const_system_iterator pos = _systems.begin();
1160  const const_system_iterator end = _systems.end();
1161 
1162  for (; pos != end; ++pos)
1163  tot += pos->second->n_vars();
1164 
1165  return tot;
1166 }
1167 
1168 
1169 
1170 std::size_t EquationSystems::n_dofs () const
1171 {
1172  std::size_t tot=0;
1173 
1174  const_system_iterator pos = _systems.begin();
1175  const const_system_iterator end = _systems.end();
1176 
1177  for (; pos != end; ++pos)
1178  tot += pos->second->n_dofs();
1179 
1180  return tot;
1181 }
1182 
1183 
1184 
1185 
1187 {
1188  std::size_t tot=0;
1189 
1190  const_system_iterator pos = _systems.begin();
1191  const const_system_iterator end = _systems.end();
1192 
1193  for (; pos != end; ++pos)
1194  tot += pos->second->n_active_dofs();
1195 
1196  return tot;
1197 }
1198 
1199 
1201 {
1202  // All the nodes
1203  for (auto & node : _mesh.node_ptr_range())
1204  node->add_system();
1205 
1206  // All the elements
1207  for (auto & elem : _mesh.element_ptr_range())
1208  elem->add_system();
1209 }
1210 
1211 } // 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:39
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const =0
Creates a local copy of the global vector in v_local only on processor proc_id.
virtual bool is_serial() const
Definition: mesh_base.h:140
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
Data structure for specifying which Parameters should be independent variables in a parameter sensiti...
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
ImplicitSystem & sys
bool refine_elements()
Only refines the user-requested elements.
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:169
virtual ~EquationSystems()
Destructor.
processor_id_type n_processors() const
virtual dof_id_type max_node_id() const =0
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
void clean_refinement_flags()
Sets the refinement flag to Elem::DO_NOTHING for each element in the mesh.
static const Real TOLERANCE
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
Call adjoint_solve on all the individual equation systems.
const Parallel::Communicator & _communicator
friend std::ostream & operator<<(std::ostream &os, const EquationSystems &es)
Same as above, but allows you to also use stream syntax.
UniquePtr< NumericVector< Number > > build_parallel_solution_vector(const std::set< std::string > *system_names=libmesh_nullptr) const
A version of build_solution_vector which is appropriate for "parallel" output formats like Nemesis...
The libMesh namespace provides an interface to certain functionality in the library.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
Real distance(const Point &p)
void build_solution_vector(std::vector< Number > &soln, const std::string &system_name, const std::string &variable_name="all_vars") const
Fill the input vector soln with the solution values for the system named name.
void build_discontinuous_solution_vector(std::vector< Number > &soln, const std::set< std::string > *system_names=libmesh_nullptr) const
Fill the input vector soln with solution values.
virtual std::string get_info() const
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void prolong_vectors()
Prolong vectors after the mesh has refined.
Definition: system.C:375
MeshBase & _mesh
The mesh data structure.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
bool coarsen_elements()
Only coarsens the user-requested elements.
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
signed char & underrefined_boundary_limit()
If underrefined_boundary_limit is set to a nonnegative value, then refinement and coarsening will pro...
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
void allgather()
Serializes a distributed mesh and its associated degree of freedom numbering for all systems...
virtual void sensitivity_solve(const ParameterVector &parameters)
Call sensitivity_solve on all the individual equation systems.
std::map< std::string, System * >::const_iterator const_system_iterator
Typedef for constant system iterators.
virtual dof_id_type max_elem_id() const =0
void update_global_solution(std::vector< Number > &global_soln) const
Fill the input vector global_soln so that it contains the global solution on all processors.
Definition: system.C:662
This class defines the notion of a variable in the system.
Definition: variable.h:49
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
virtual SimpleRange< element_iterator > element_ptr_range()=0
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual node_iterator local_nodes_end()=0
std::map< std::string, System * >::iterator system_iterator
Typedef for system iterators.
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2114
virtual bool compare(const EquationSystems &other_es, const Real threshold, const bool verbose) const
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual void reinit()
Reinitialize all the systems.
void delete_system(const std::string &name)
Remove the system named name from the systems array.
virtual SimpleRange< node_iterator > node_ptr_range()=0
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
void distribute_dofs(MeshBase &)
Distribute dofs on the current mesh.
Definition: dof_map.C:882
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void restrict_vectors()
Restrict vectors after the mesh has coarsened.
Definition: system.C:322
bool _refine_in_reinit
Flag for whether to call coarsen/refine in reinit().
unsigned char & face_level_mismatch_limit()
If face_level_mismatch_limit is set to a nonzero value, then refinement and coarsening will produce m...
This class forms the base class for all other classes that are expected to be implemented in parallel...
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:414
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
unsigned int n_systems() const
virtual void solve()
Call solve on all the individual equation systems.
virtual bool compare(const System &other_system, const Real threshold, const bool verbose) const
Definition: system.C:531
std::size_t n_dofs() const
virtual node_iterator local_nodes_begin()=0
Iterate over local nodes (nodes whose processor_id() matches the current processor).
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool contract()=0
Delete subactive (i.e.
bool verify(const T &r, const Communicator &comm=Communicator_World)
unsigned int number() const
Definition: system.h:2006
T & set(const std::string &)
Definition: parameters.h:469
const Parallel::Communicator & comm() const
OStreamProxy out
virtual void clear()
Clears internal data structures & frees any allocated memory.
Definition: parameters.h:321
void build_variable_names(std::vector< std::string > &var_names, const FEType *type=libmesh_nullptr, const std::set< std::string > *system_names=libmesh_nullptr) const
Fill the input vector var_names with the names of the variables for each system.
signed char & overrefined_boundary_limit()
If overrefined_boundary_limit is set to a nonnegative value, then refinement and coarsening will prod...
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
Parameters parameters
Data structure holding arbitrary parameters.
unsigned int n_vars() const
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
void prepare_send_list()
Takes the _send_list vector (which may have duplicate entries) and sorts it.
Definition: dof_map.C:1635
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)
Build the nodal soln from the element soln.
Definition: fe_interface.C:526
virtual void init()
Initialize all the systems.
std::size_t n_active_dofs() const
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
virtual dof_id_type n_nodes() const =0
virtual dof_id_type n_elem() const =0
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...
void update()
Updates local values for all the systems.
void get_solution(std::vector< Number > &soln, std::vector< std::string > &names) const
Retrieve the solution data for CONSTANT MONOMIALs.
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
std::map< std::string, System * > _systems
Data structure holding the systems.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917