equation_systems.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
24 #include "libmesh/fe_interface.h"
28 #include "libmesh/newmark_system.h"
32 #include "libmesh/eigen_system.h"
33 #include "libmesh/parallel.h"
35 #include "libmesh/dof_map.h"
36 #include "libmesh/mesh_base.h"
37 #include "libmesh/elem.h"
39 
40 // Include the systems before this one to avoid
41 // overlapping forward declarations.
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  _mesh_data (mesh_data)
58 {
59  // Set default parameters
60  this->parameters.set<Real> ("linear solver tolerance") = TOLERANCE * TOLERANCE;
61  this->parameters.set<unsigned int>("linear solver maximum iterations") = 5000;
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 = NULL;
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  // Distribute the mesh if possible
101  if (this->n_processors() > 1)
103 
104  // Tell all the \p DofObject entities how many systems
105  // there are.
106  {
108  const MeshBase::node_iterator node_end = _mesh.nodes_end();
109 
110  for ( ; node_it != node_end; ++node_it)
111  (*node_it)->set_n_systems(n_sys);
112 
114  const MeshBase::element_iterator elem_end = _mesh.elements_end();
115 
116  for ( ; elem_it != elem_end; ++elem_it)
117  (*elem_it)->set_n_systems(n_sys);
118  }
119 
120  for (unsigned int i=0; i != this->n_systems(); ++i)
121  this->get_system(i).init();
122 
123 #ifdef LIBMESH_ENABLE_AMR
124  MeshRefinement mesh_refine(_mesh);
125  mesh_refine.clean_refinement_flags();
126 #endif
127 }
128 
129 
130 
132 {
133  parallel_object_only();
134 
135  libmesh_assert_not_equal_to (this->n_systems(), 0);
136 
137 #ifdef DEBUG
138  // Make sure all the \p DofObject entities know how many systems
139  // there are.
140  {
141  // All the nodes
143  const MeshBase::node_iterator node_end = _mesh.nodes_end();
144 
145  for ( ; node_it != node_end; ++node_it)
146  {
147  Node *node = *node_it;
148  libmesh_assert_equal_to (node->n_systems(), this->n_systems());
149  }
150 
151  // All the elements
153  const MeshBase::element_iterator elem_end = _mesh.elements_end();
154 
155  for ( ; elem_it != elem_end; ++elem_it)
156  {
157  Elem *elem = *elem_it;
158  libmesh_assert_equal_to (elem->n_systems(), this->n_systems());
159  }
160  }
161 #endif
162 
163  // Localize each system's vectors
164  for (unsigned int i=0; i != this->n_systems(); ++i)
165  this->get_system(i).re_update();
166 
167 #ifdef LIBMESH_ENABLE_AMR
168 
169  bool dof_constraints_created = false;
170  bool mesh_changed = false;
171 
172  // FIXME: For backwards compatibility, assume
173  // refine_and_coarsen_elements or refine_uniformly have already
174  // been called
175  {
176  for (unsigned int i=0; i != this->n_systems(); ++i)
177  {
178  System &sys = this->get_system(i);
179 
180  // Don't do anything if the system doesn't have any variables in it
181  if(!sys.n_vars())
182  continue;
183 
185 
186  // Recreate any hanging node constraints
188 
189  // Apply any user-defined constraints
190  sys.user_constrain();
191 
192  // Expand any recursive constraints
194 
195  // And clean up the send_list before we use it again
197 
198  sys.prolong_vectors();
199  }
200  mesh_changed = true;
201  dof_constraints_created = true;
202  }
203 
204  // FIXME: Where should the user set maintain_level_one now??
205  // Don't override previous settings, for now
206 
207  MeshRefinement mesh_refine(_mesh);
208 
209  mesh_refine.face_level_mismatch_limit() = false;
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  {
222  sys.user_constrain();
225 
226  }
227  sys.restrict_vectors();
228  }
229  mesh_changed = true;
230  dof_constraints_created = true;
231  }
232 
233  // Once vectors are all restricted, we can delete
234  // children of coarsened elements
235  if (mesh_changed)
236  this->get_mesh().contract();
237 
238  // Try to refine the mesh, then prolong each system's vectors
239  // if necessary
240  if (mesh_refine.refine_elements())
241  {
242  for (unsigned int i=0; i != this->n_systems(); ++i)
243  {
244  System &sys = this->get_system(i);
245  if (!dof_constraints_created)
246  {
249  sys.user_constrain();
252 
253  }
254  sys.prolong_vectors();
255  }
256  mesh_changed = true;
257  dof_constraints_created = true;
258  }
259 
260  // If the mesh has changed, systems will need to create new dof
261  // constraints and update their global solution vectors
262  if (mesh_changed)
263  {
264  for (unsigned int i=0; i != this->n_systems(); ++i)
265  this->get_system(i).reinit();
266  }
267 #endif // #ifdef LIBMESH_ENABLE_AMR
268 }
269 
270 
271 
273 {
274  // A serial mesh means nothing needs to be done
275  if (_mesh.is_serial())
276  return;
277 
278  const unsigned int n_sys = this->n_systems();
279 
280  libmesh_assert_not_equal_to (n_sys, 0);
281 
282  // Gather the mesh
283  _mesh.allgather();
284 
285  // Tell all the \p DofObject entities how many systems
286  // there are.
287  {
289  const MeshBase::node_iterator node_end = _mesh.nodes_end();
290 
291  for ( ; node_it != node_end; ++node_it)
292  (*node_it)->set_n_systems(n_sys);
293 
295  const MeshBase::element_iterator elem_end = _mesh.elements_end();
296 
297  for ( ; elem_it != elem_end; ++elem_it)
298  (*elem_it)->set_n_systems(n_sys);
299  }
300 
301  // And distribute each system's dofs
302  for (unsigned int i=0; i != this->n_systems(); ++i)
303  {
304  System &sys = this->get_system(i);
305  DofMap &dof_map = sys.get_dof_map();
306  dof_map.distribute_dofs(_mesh);
307 
308 #ifdef LIBMESH_ENABLE_CONSTRAINTS
309  // The user probably won't need constraint equations or the
310  // send_list after an allgather, but let's keep it in consistent
311  // shape just in case.
312  dof_map.create_dof_constraints(_mesh, sys.time);
313  sys.user_constrain();
314  dof_map.process_constraints(_mesh);
315 #endif
316  dof_map.prepare_send_list();
317  }
318 }
319 
320 
321 
322 
324 {
325  START_LOG("update()","EquationSystems");
326 
327  // Localize each system's vectors
328  for (unsigned int i=0; i != this->n_systems(); ++i)
329  this->get_system(i).update();
330 
331  STOP_LOG("update()","EquationSystems");
332 }
333 
334 
335 
336 System & EquationSystems::add_system (const std::string& sys_type,
337  const std::string& name)
338 {
339  // If the user already built a system with this name, we'll
340  // trust them and we'll use it. That way they can pre-add
341  // non-standard derived system classes, and if their restart file
342  // has some non-standard sys_type we won't throw an error.
343  if (_systems.count(name))
344  {
345  return this->get_system(name);
346  }
347  // Build a basic System
348  else if (sys_type == "Basic")
349  this->add_system<System> (name);
350 
351  // Build a Newmark system
352  else if (sys_type == "Newmark")
353  this->add_system<NewmarkSystem> (name);
354 
355  // Build an Explicit system
356  else if ((sys_type == "Explicit"))
357  this->add_system<ExplicitSystem> (name);
358 
359  // Build an Implicit system
360  else if ((sys_type == "Implicit") ||
361  (sys_type == "Steady" ))
362  this->add_system<ImplicitSystem> (name);
363 
364  // build a transient implicit linear system
365  else if ((sys_type == "Transient") ||
366  (sys_type == "TransientImplicit") ||
367  (sys_type == "TransientLinearImplicit"))
368  this->add_system<TransientLinearImplicitSystem> (name);
369 
370  // build a transient implicit nonlinear system
371  else if (sys_type == "TransientNonlinearImplicit")
372  this->add_system<TransientNonlinearImplicitSystem> (name);
373 
374  // build a transient explicit system
375  else if (sys_type == "TransientExplicit")
376  this->add_system<TransientExplicitSystem> (name);
377 
378  // build a linear implicit system
379  else if (sys_type == "LinearImplicit")
380  this->add_system<LinearImplicitSystem> (name);
381 
382  // build a nonlinear implicit system
383  else if (sys_type == "NonlinearImplicit")
384  this->add_system<NonlinearImplicitSystem> (name);
385 
386  // build a Reduced Basis Construction system
387  else if (sys_type == "RBConstruction")
388  this->add_system<RBConstruction> (name);
389 
390  // build a transient Reduced Basis Construction system
391  else if (sys_type == "TransientRBConstruction")
392  this->add_system<TransientRBConstruction> (name);
393 
394 #ifdef LIBMESH_HAVE_SLEPC
395  // build an eigen system
396  else if (sys_type == "Eigen")
397  this->add_system<EigenSystem> (name);
398 #endif
399 
400 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
401  // build a frequency system
402  else if (sys_type == "Frequency")
403  this->add_system<FrequencySystem> (name);
404 #endif
405 
406  else
407  {
408  libMesh::err << "ERROR: Unknown system type: "
409  << sys_type
410  << std::endl;
411  libmesh_error();
412  }
413 
414  // Return a reference to the new system
415  //return (*this)(name);
416  return this->get_system(name);
417 }
418 
419 
420 
421 
422 
423 
424 void EquationSystems::delete_system (const std::string& name)
425 {
426  libmesh_deprecated();
427 
428  if (!_systems.count(name))
429  {
430  libMesh::err << "ERROR: no system named "
431  << name << std::endl;
432 
433  libmesh_error();
434  }
435 
436  delete _systems[name];
437 
438  _systems.erase (name);
439 }
440 
441 
442 
444 {
445  libmesh_assert (this->n_systems());
446 
447  for (unsigned int i=0; i != this->n_systems(); ++i)
448  this->get_system(i).solve();
449 }
450 
451 
452 
454 {
455  libmesh_assert (this->n_systems());
456 
457  for (unsigned int i=0; i != this->n_systems(); ++i)
458  this->get_system(i).sensitivity_solve(parameters_in);
459 }
460 
461 
462 
463 void EquationSystems::adjoint_solve (const QoISet& qoi_indices)
464 {
465  libmesh_assert (this->n_systems());
466 
467  for (unsigned int i=this->n_systems(); i != 0; --i)
468  this->get_system(i-1).adjoint_solve(qoi_indices);
469 }
470 
471 
472 
473 void EquationSystems::build_variable_names (std::vector<std::string>& var_names,
474  const FEType *type,
475  const std::set<std::string>* system_names) const
476 {
477  libmesh_assert (this->n_systems());
478 
479  unsigned int var_num=0;
480 
481  const_system_iterator pos = _systems.begin();
482  const const_system_iterator end = _systems.end();
483 
484  // Need to size var_names by scalar variables plus all the
485  // vector components for all the vector variables
486  //Could this be replaced by a/some convenience methods?[PB]
487  {
488  unsigned int n_scalar_vars = 0;
489  unsigned int n_vector_vars = 0;
490 
491  for (; pos != end; ++pos)
492  {
493  // Check current system is listed in system_names, and skip pos if not
494  bool use_current_system = (system_names == NULL);
495  if(!use_current_system)
496  {
497  use_current_system = std::find( system_names->begin(), system_names->end(), pos->first ) != system_names->end();
498  }
499  if(!use_current_system)
500  {
501  continue;
502  }
503 
504  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
505  {
506  if( FEInterface::field_type(pos->second->variable_type(vn)) ==
507  TYPE_VECTOR )
508  n_vector_vars++;
509  else
510  n_scalar_vars++;
511  }
512  }
513 
514  // Here, we're assuming the number of vector components is the same
515  // as the mesh dimension. Will break for mixed dimension meshes.
516  unsigned int dim = this->get_mesh().mesh_dimension();
517  unsigned int nv = n_scalar_vars + dim*n_vector_vars;
518 
519  // We'd better not have more than dim*his->n_vars() (all vector variables)
520  libmesh_assert_less_equal ( nv, dim*this->n_vars() );
521 
522  // Here, we're assuming the number of vector components is the same
523  // as the mesh dimension. Will break for mixed dimension meshes.
524 
525  var_names.resize( nv );
526  }
527 
528  // reset
529  pos = _systems.begin();
530 
531  for (; pos != end; ++pos)
532  {
533  // Check current system is listed in system_names, and skip pos if not
534  bool use_current_system = (system_names == NULL);
535  if(!use_current_system)
536  {
537  use_current_system = std::find( system_names->begin(), system_names->end(), pos->first ) != system_names->end();
538  }
539  if(!use_current_system)
540  {
541  continue;
542  }
543 
544  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
545  {
546  std::string var_name = pos->second->variable_name(vn);
547  FEType fe_type = pos->second->variable_type(vn);
548 
549  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type);
550 
551  // Filter on the type if requested
552  if (type == NULL || (type && *type == fe_type))
553  {
554  if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
555  {
556  switch(n_vec_dim)
557  {
558  case 0:
559  case 1:
560  var_names[var_num++] = var_name;
561  break;
562  case 2:
563  var_names[var_num++] = var_name+"_x";
564  var_names[var_num++] = var_name+"_y";
565  break;
566  case 3:
567  var_names[var_num++] = var_name+"_x";
568  var_names[var_num++] = var_name+"_y";
569  var_names[var_num++] = var_name+"_z";
570  break;
571  default:
572  libMesh::err << "Invalid dim in build_variable_names" << std::endl;
573  libmesh_error();
574  }
575  }
576  else
577  var_names[var_num++] = var_name;
578  }
579  }
580  }
581  // Now resize again in case we filtered any names
582  var_names.resize(var_num);
583 }
584 
585 
586 
587 void EquationSystems::build_solution_vector (std::vector<Number>&,
588  const std::string&,
589  const std::string&) const
590 {
591  //TODO:[BSK] re-implement this from the method below
592  libmesh_error();
593 
594 // // Get a reference to the named system
595 // const System& system = this->get_system(system_name);
596 
597 // // Get the number associated with the variable_name we are passed
598 // const unsigned short int variable_num = system.variable_number(variable_name);
599 
600 // // Get the dimension of the current mesh
601 // const unsigned int dim = _mesh.mesh_dimension();
602 
603 // // If we're on processor 0, allocate enough memory to hold the solution.
604 // // Since we're only looking at one variable, there will be one solution value
605 // // for each node in the mesh.
606 // if (_mesh.processor_id() == 0)
607 // soln.resize(_mesh.n_nodes());
608 
609 // // Vector to hold the global solution from all processors
610 // std::vector<Number> sys_soln;
611 
612 // // Update the global solution from all processors
613 // system.update_global_solution (sys_soln, 0);
614 
615 // // Temporary vector to store the solution on an individual element.
616 // std::vector<Number> elem_soln;
617 
618 // // The FE solution interpolated to the nodes
619 // std::vector<Number> nodal_soln;
620 
621 // // The DOF indices for the element
622 // std::vector<dof_id_type> dof_indices;
623 
624 // // Determine the finite/infinite element type used in this system
625 // const FEType& fe_type = system.variable_type(variable_num);
626 
627 // // Define iterators to iterate over all the elements of the mesh
628 // const_active_elem_iterator it (_mesh.elements_begin());
629 // const const_active_elem_iterator end(_mesh.elements_end());
630 
631 // // Loop over elements
632 // for ( ; it != end; ++it)
633 // {
634 // // Convenient shortcut to the element pointer
635 // const Elem* elem = *it;
636 
637 // // Fill the dof_indices vector for this variable
638 // system.get_dof_map().dof_indices(elem,
639 // dof_indices,
640 // variable_num);
641 
642 // // Resize the element solution vector to fit the
643 // // dof_indices for this element.
644 // elem_soln.resize(dof_indices.size());
645 
646 // // Transfer the system solution to the element
647 // // solution by mapping it through the dof_indices vector.
648 // for (unsigned int i=0; i<dof_indices.size(); i++)
649 // elem_soln[i] = sys_soln[dof_indices[i]];
650 
651 // // Using the FE interface, compute the nodal_soln
652 // // for the current elemnt type given the elem_soln
653 // FEInterface::nodal_soln (dim,
654 // fe_type,
655 // elem,
656 // elem_soln,
657 // nodal_soln);
658 
659 // // Sanity check -- make sure that there are the same number
660 // // of entries in the nodal_soln as there are nodes in the
661 // // element!
662 // libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
663 
664 // // Copy the nodal solution over into the correct place in
665 // // the global soln vector which will be returned to the user.
666 // for (unsigned int n=0; n<elem->n_nodes(); n++)
667 // soln[elem->node(n)] = nodal_soln[n];
668 // }
669 }
670 
671 
672 
673 
674 void EquationSystems::build_solution_vector (std::vector<Number>& soln,
675  const std::set<std::string>* system_names) const
676 {
677  START_LOG("build_solution_vector()", "EquationSystems");
678 
679  // This function must be run on all processors at once
680  parallel_object_only();
681 
682  libmesh_assert (this->n_systems());
683 
684  const unsigned int dim = _mesh.mesh_dimension();
685  const dof_id_type nn = _mesh.n_nodes();
686 
687  // We'd better have a contiguous node numbering
688  libmesh_assert_equal_to (nn, _mesh.max_node_id());
689 
690  // allocate storage to hold
691  // (number_of_nodes)*(number_of_variables) entries.
692  // We have to differentiate between between scalar and vector
693  // variables. We intercept vector variables and treat each
694  // component as a scalar variable (consistently with build_solution_names).
695 
696  unsigned int nv = 0;
697 
698  //Could this be replaced by a/some convenience methods?[PB]
699  {
700  unsigned int n_scalar_vars = 0;
701  unsigned int n_vector_vars = 0;
702  const_system_iterator pos = _systems.begin();
703  const const_system_iterator end = _systems.end();
704 
705  for (; pos != end; ++pos)
706  {
707  // Check current system is listed in system_names, and skip pos if not
708  bool use_current_system = (system_names == NULL);
709  if(!use_current_system)
710  {
711  use_current_system = std::find( system_names->begin(), system_names->end(), pos->first ) != system_names->end();
712  }
713  if(!use_current_system)
714  {
715  continue;
716  }
717 
718  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
719  {
720  if( FEInterface::field_type(pos->second->variable_type(vn)) ==
721  TYPE_VECTOR )
722  n_vector_vars++;
723  else
724  n_scalar_vars++;
725  }
726  }
727  // Here, we're assuming the number of vector components is the same
728  // as the mesh dimension. Will break for mixed dimension meshes.
729  nv = n_scalar_vars + dim*n_vector_vars;
730  }
731 
732  // Get the number of elements that share each node. We will
733  // compute the average value at each node. This is particularly
734  // useful for plotting discontinuous data.
737 
738  // Get the number of local nodes
739  dof_id_type n_local_nodes = std::distance(_mesh.local_nodes_begin(), _mesh.local_nodes_end());
740 
741  // Create a NumericVector to hold the parallel solution
743  NumericVector<Number> &parallel_soln = *parallel_soln_ptr;
744  parallel_soln.init(nn*nv, n_local_nodes*nv, false, PARALLEL);
745 
746  // Create a NumericVector to hold the "repeat_count" for each node - this is essentially
747  // the number of elements contributing to that node's value
749  NumericVector<Number> &repeat_count = *repeat_count_ptr;
750  repeat_count.init(nn*nv, n_local_nodes*nv, false, PARALLEL);
751 
752  repeat_count.close();
753 
754  unsigned int var_num=0;
755 
756  // For each system in this EquationSystems object,
757  // update the global solution and if we are on processor 0,
758  // loop over the elements and build the nodal solution
759  // from the element solution. Then insert this nodal solution
760  // into the vector passed to build_solution_vector.
761  const_system_iterator pos = _systems.begin();
762  const const_system_iterator end = _systems.end();
763 
764  for (; pos != end; ++pos)
765  {
766  // Check current system is listed in system_names, and skip pos if not
767  bool use_current_system = (system_names == NULL);
768  if(!use_current_system)
769  {
770  use_current_system = std::find( system_names->begin(), system_names->end(), pos->first ) != system_names->end();
771  }
772  if(!use_current_system)
773  {
774  continue;
775  }
776 
777  const System& system = *(pos->second);
778  const unsigned int nv_sys = system.n_vars();
779  const unsigned int sys_num = system.number();
780 
781  //Could this be replaced by a/some convenience methods?[PB]
782  unsigned int n_scalar_vars = 0;
783  unsigned int n_vector_vars = 0;
784  for (unsigned int vn=0; vn<pos->second->n_vars(); vn++)
785  {
786  if( FEInterface::field_type(pos->second->variable_type(vn)) ==
787  TYPE_VECTOR )
788  n_vector_vars++;
789  else
790  n_scalar_vars++;
791  }
792 
793  // Here, we're assuming the number of vector components is the same
794  // as the mesh dimension. Will break for mixed dimension meshes.
795  unsigned int nv_sys_split = n_scalar_vars + dim*n_vector_vars;
796 
797  // Update the current_local_solution
798  {
799  System & non_const_sys = const_cast<System &>(system);
800  non_const_sys.solution->close();
801  non_const_sys.update();
802  }
803 
804  NumericVector<Number> & sys_soln(*system.current_local_solution);
805 
806  std::vector<Number> elem_soln; // The finite element solution
807  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
808  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
809 
810  for (unsigned int var=0; var<nv_sys; var++)
811  {
812  const FEType& fe_type = system.variable_type(var);
813  const Variable &var_description = system.variable(var);
814  const DofMap &dof_map = system.get_dof_map();
815 
816  unsigned int n_vec_dim = FEInterface::n_vec_dim( pos->second->get_mesh(), fe_type );
817 
820 
821  for ( ; it != end_elem; ++it)
822  {
823  const Elem* elem = *it;
824 
825  if (var_description.active_on_subdomain((*it)->subdomain_id()))
826  {
827  dof_map.dof_indices (elem, dof_indices, var);
828 
829  elem_soln.resize(dof_indices.size());
830 
831  for (unsigned int i=0; i<dof_indices.size(); i++)
832  elem_soln[i] = sys_soln(dof_indices[i]);
833 
835  fe_type,
836  elem,
837  elem_soln,
838  nodal_soln);
839 
840 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
841  // infinite elements should be skipped...
842  if (!elem->infinite())
843 #endif
844  {
845  libmesh_assert_equal_to (nodal_soln.size(), n_vec_dim*elem->n_nodes());
846 
847  for (unsigned int n=0; n<elem->n_nodes(); n++)
848  {
849  for( unsigned int d=0; d < n_vec_dim; d++ )
850  {
851  // For vector-valued elements, all components are in nodal_soln. For each
852  // node, the components are stored in order, i.e. node_0 -> s0_x, s0_y, s0_z
853  parallel_soln.add(nv*(elem->node(n)) + (var+d + var_num), nodal_soln[n_vec_dim*n+d]);
854 
855  // Increment the repeat count for this position
856  repeat_count.add(nv*(elem->node(n)) + (var+d + var_num), 1);
857  }
858  }
859  }
860  }
861  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:
862  for (unsigned int n=0; n<elem->n_nodes(); n++)
863  // Only do this if this variable has NO DoFs at this node... it might have some from an ajoining element...
864  if(!elem->get_node(n)->n_dofs(sys_num, var))
865  for( unsigned int d=0; d < n_vec_dim; d++ )
866  repeat_count.add(nv*(elem->node(n)) + (var+d + var_num), 1);
867 
868  } // end loop over elements
869  } // end loop on variables in this system
870 
871  var_num += nv_sys_split;
872  } // end loop over systems
873 
874  parallel_soln.close();
875  repeat_count.close();
876 
877  // Divide to get the average value at the nodes
878  parallel_soln /= repeat_count;
879 
880  parallel_soln.localize_to_one(soln);
881 
882  STOP_LOG("build_solution_vector()", "EquationSystems");
883 }
884 
885 
886 void EquationSystems::get_solution (std::vector<Number>& soln,
887  std::vector<std::string> & names ) const
888 {
889  // This function must be run on all processors at once
890  parallel_object_only();
891 
892  libmesh_assert (this->n_systems());
893 
894  const dof_id_type ne = _mesh.n_elem();
895 
896  libmesh_assert_equal_to (ne, _mesh.max_elem_id());
897 
898  // Get the number of local elements
899  dof_id_type n_local_elems = std::distance(_mesh.local_elements_begin(), _mesh.local_elements_end());
900 
901  // If the names vector has entries, we will only populate the soln vector
902  // with names included in that list. Note: The names vector may be
903  // reordered upon exiting this function
904  std::vector<std::string> filter_names = names;
905  bool is_filter_names = ! filter_names.empty();
906 
907  soln.clear();
908  names.clear();
909 
910  const FEType type(CONSTANT, MONOMIAL);
911 
912  dof_id_type nv = 0;
913 
914  // Find the total number of variables to output
915  {
916  const_system_iterator pos = _systems.begin();
917  const const_system_iterator end = _systems.end();
918 
919  for (; pos != end; ++pos)
920  {
921  const System& system = *(pos->second);
922  const unsigned int nv_sys = system.n_vars();
923 
924  for (unsigned int var=0; var < nv_sys; ++var)
925  {
926  if ( system.variable_type( var ) != type ||
927  ( is_filter_names && std::find(filter_names.begin(), filter_names.end(), system.variable_name( var )) == filter_names.end()) )
928  continue;
929 
930  nv++;
931  }
932  }
933  }
934 
935  if(!nv) // If there are no variables to write out don't do anything...
936  return;
937 
938  // Create a NumericVector to hold the parallel solution
940  NumericVector<Number> &parallel_soln = *parallel_soln_ptr;
941  parallel_soln.init(ne*nv, n_local_elems*nv, false, PARALLEL);
942 
943  dof_id_type var_num = 0;
944 
945  // For each system in this EquationSystems object,
946  // update the global solution and collect the
947  // CONSTANT MONOMIALs. The entries are in variable-major
948  // format.
949  const_system_iterator pos = _systems.begin();
950  const const_system_iterator end = _systems.end();
951 
952  for (; pos != end; ++pos)
953  {
954  const System& system = *(pos->second);
955  const unsigned int nv_sys = system.n_vars();
956 
957  // Update the current_local_solution
958  {
959  System & non_const_sys = const_cast<System &>(system);
960  non_const_sys.solution->close();
961  non_const_sys.update();
962  }
963 
964  NumericVector<Number> & sys_soln(*system.current_local_solution);
965 
966  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
967 
968  // Loop over the variable names and load them in order
969  for (unsigned int var=0; var < nv_sys; ++var)
970  {
971  if ( system.variable_type( var ) != type ||
972  ( is_filter_names && std::find(filter_names.begin(), filter_names.end(), system.variable_name( var )) == filter_names.end()) )
973  continue;
974 
975  names.push_back( system.variable_name( var ) );
976 
977  const Variable & variable = system.variable(var);
978  const DofMap & dof_map = system.get_dof_map();
979 
982 
983  for ( ; it != end_elem; ++it)
984  {
985  if (variable.active_on_subdomain((*it)->subdomain_id()))
986  {
987  const Elem* elem = *it;
988 
989  dof_map.dof_indices (elem, dof_indices, var);
990 
991  libmesh_assert_equal_to ( 1, dof_indices.size() );
992 
993  parallel_soln.set((ne*var_num)+elem->id(), sys_soln(dof_indices[0]));
994  }
995  }
996 
997  var_num++;
998  } // end loop on variables in this system
999  } // end loop over systems
1000 
1001  parallel_soln.close();
1002 
1003  parallel_soln.localize_to_one(soln);
1004 }
1005 
1006 
1007 
1008 void EquationSystems::build_discontinuous_solution_vector (std::vector<Number>& soln) const
1009 {
1010  START_LOG("build_discontinuous_solution_vector()", "EquationSystems");
1011 
1012  libmesh_assert (this->n_systems());
1013 
1014  const unsigned int dim = _mesh.mesh_dimension();
1015  const unsigned int nv = this->n_vars();
1016  unsigned int tw=0;
1017 
1018  // get the total weight
1019  {
1022 
1023  for ( ; it != end; ++it)
1024  tw += (*it)->n_nodes();
1025  }
1026 
1027 
1028  // Only if we are on processor zero, allocate the storage
1029  // to hold (number_of_nodes)*(number_of_variables) entries.
1030  if (_mesh.processor_id() == 0)
1031  soln.resize(tw*nv);
1032 
1033  std::vector<Number> sys_soln;
1034 
1035 
1036  unsigned int var_num=0;
1037 
1038  // For each system in this EquationSystems object,
1039  // update the global solution and if we are on processor 0,
1040  // loop over the elements and build the nodal solution
1041  // from the element solution. Then insert this nodal solution
1042  // into the vector passed to build_solution_vector.
1043  const_system_iterator pos = _systems.begin();
1044  const const_system_iterator end = _systems.end();
1045 
1046  for (; pos != end; ++pos)
1047  {
1048  const System& system = *(pos->second);
1049  const unsigned int nv_sys = system.n_vars();
1050 
1051  system.update_global_solution (sys_soln, 0);
1052 
1053  if (_mesh.processor_id() == 0)
1054  {
1055  std::vector<Number> elem_soln; // The finite element solution
1056  std::vector<Number> nodal_soln; // The FE solution interpolated to the nodes
1057  std::vector<dof_id_type> dof_indices; // The DOF indices for the finite element
1058 
1059  for (unsigned int var=0; var<nv_sys; var++)
1060  {
1061  const FEType& fe_type = system.variable_type(var);
1062 
1065 
1066  unsigned int nn=0;
1067 
1068  for ( ; it != end_elem; ++it)
1069  {
1070  const Elem* elem = *it;
1071  system.get_dof_map().dof_indices (elem, dof_indices, var);
1072 
1073  elem_soln.resize(dof_indices.size());
1074 
1075  for (unsigned int i=0; i<dof_indices.size(); i++)
1076  elem_soln[i] = sys_soln[dof_indices[i]];
1077 
1079  fe_type,
1080  elem,
1081  elem_soln,
1082  nodal_soln);
1083 
1084 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1085  // infinite elements should be skipped...
1086  if (!elem->infinite())
1087 #endif
1088  {
1089  libmesh_assert_equal_to (nodal_soln.size(), elem->n_nodes());
1090 
1091  for (unsigned int n=0; n<elem->n_nodes(); n++)
1092  {
1093  soln[nv*(nn++) + (var + var_num)] +=
1094  nodal_soln[n];
1095  }
1096  }
1097  }
1098  }
1099  }
1100 
1101  var_num += nv_sys;
1102  }
1103 
1104  STOP_LOG("build_discontinuous_solution_vector()", "EquationSystems");
1105 }
1106 
1107 
1108 
1110  const Real threshold,
1111  const bool verbose) const
1112 {
1113  // safety check, whether we handle at least the same number
1114  // of systems
1115  std::vector<bool> os_result;
1116 
1117  if (this->n_systems() != other_es.n_systems())
1118  {
1119  if (verbose)
1120  {
1121  libMesh::out << " Fatal difference. This system handles "
1122  << this->n_systems() << " systems," << std::endl
1123  << " while the other system handles "
1124  << other_es.n_systems()
1125  << " systems." << std::endl
1126  << " Aborting comparison." << std::endl;
1127  }
1128  return false;
1129  }
1130  else
1131  {
1132  // start comparing each system
1133  const_system_iterator pos = _systems.begin();
1134  const const_system_iterator end = _systems.end();
1135 
1136  for (; pos != end; ++pos)
1137  {
1138  const std::string& sys_name = pos->first;
1139  const System& system = *(pos->second);
1140 
1141  // get the other system
1142  const System& other_system = other_es.get_system (sys_name);
1143 
1144  os_result.push_back (system.compare (other_system, threshold, verbose));
1145 
1146  }
1147 
1148  }
1149 
1150 
1151  // sum up the results
1152  if (os_result.size()==0)
1153  return true;
1154  else
1155  {
1156  bool os_identical;
1157  unsigned int n = 0;
1158  do
1159  {
1160  os_identical = os_result[n];
1161  n++;
1162  }
1163  while (os_identical && n<os_result.size());
1164  return os_identical;
1165  }
1166 }
1167 
1168 
1169 
1170 std::string EquationSystems::get_info () const
1171 {
1172  std::ostringstream oss;
1173 
1174  oss << " EquationSystems\n"
1175  << " n_systems()=" << this->n_systems() << '\n';
1176 
1177  // Print the info for the individual systems
1178  const_system_iterator pos = _systems.begin();
1179  const const_system_iterator end = _systems.end();
1180 
1181  for (; pos != end; ++pos)
1182  oss << pos->second->get_info();
1183 
1184 
1185 // // Possibly print the parameters
1186 // if (!this->parameters.empty())
1187 // {
1188 // oss << " n_parameters()=" << this->n_parameters() << '\n';
1189 // oss << " Parameters:\n";
1190 
1191 // for (std::map<std::string, Real>::const_iterator
1192 // param = _parameters.begin(); param != _parameters.end();
1193 // ++param)
1194 // oss << " "
1195 // << "\""
1196 // << param->first
1197 // << "\""
1198 // << "="
1199 // << param->second
1200 // << '\n';
1201 // }
1202 
1203  return oss.str();
1204 }
1205 
1206 
1207 
1208 void EquationSystems::print_info (std::ostream& os) const
1209 {
1210  os << this->get_info()
1211  << std::endl;
1212 }
1213 
1214 
1215 
1216 std::ostream& operator << (std::ostream& os, const EquationSystems& es)
1217 {
1218  es.print_info(os);
1219  return os;
1220 }
1221 
1222 
1223 
1224 unsigned int EquationSystems::n_vars () const
1225 {
1226  unsigned int tot=0;
1227 
1228  const_system_iterator pos = _systems.begin();
1229  const const_system_iterator end = _systems.end();
1230 
1231  for (; pos != end; ++pos)
1232  tot += pos->second->n_vars();
1233 
1234  return tot;
1235 }
1236 
1237 
1238 
1239 std::size_t EquationSystems::n_dofs () const
1240 {
1241  std::size_t tot=0;
1242 
1243  const_system_iterator pos = _systems.begin();
1244  const const_system_iterator end = _systems.end();
1245 
1246  for (; pos != end; ++pos)
1247  tot += pos->second->n_dofs();
1248 
1249  return tot;
1250 }
1251 
1252 
1253 
1254 
1256 {
1257  std::size_t tot=0;
1258 
1259  const_system_iterator pos = _systems.begin();
1260  const const_system_iterator end = _systems.end();
1261 
1262  for (; pos != end; ++pos)
1263  tot += pos->second->n_active_dofs();
1264 
1265  return tot;
1266 }
1267 
1268 
1270 {
1271  // All the nodes
1273  const MeshBase::node_iterator node_end = _mesh.nodes_end();
1274 
1275  for ( ; node_it != node_end; ++node_it)
1276  (*node_it)->add_system();
1277 
1278  // All the elements
1280  const MeshBase::element_iterator elem_end = _mesh.elements_end();
1281 
1282  for ( ; elem_it != elem_end; ++elem_it)
1283  (*elem_it)->add_system();
1284 }
1285 
1286 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo