system.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 
20 // C++ includes
21 #include <sstream> // for std::ostringstream
22 
23 
24 // Local includes
25 #include "libmesh/dof_map.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/numeric_vector.h"
31 #include "libmesh/point.h" // For point_value
32 #include "libmesh/point_locator_base.h" // For point_value
33 #include "libmesh/qoi_set.h"
34 #include "libmesh/string_to_enum.h"
35 #include "libmesh/system.h"
36 #include "libmesh/system_norm.h"
37 #include "libmesh/utility.h"
38 #include "libmesh/elem.h"
39 #include "libmesh/fe_type.h"
40 
41 // includes for calculate_norm, point_*
42 #include "libmesh/fe_base.h"
43 #include "libmesh/fe_interface.h"
44 #include "libmesh/parallel.h"
46 #include "libmesh/quadrature.h"
47 #include "libmesh/tensor_value.h"
48 #include "libmesh/vector_value.h"
49 #include "libmesh/tensor_tools.h"
50 
51 namespace libMesh
52 {
53 
54 
55 // ------------------------------------------------------------
56 // System implementation
58  const std::string & name_in,
59  const unsigned int number_in) :
60 
61  ParallelObject (es),
62  assemble_before_solve (true),
63  use_fixed_solution (false),
64  extra_quadrature_order (0),
65  solution (NumericVector<Number>::build(this->comm())),
66  current_local_solution (NumericVector<Number>::build(this->comm())),
67  time (0.),
68  qoi (0),
69  _init_system_function (libmesh_nullptr),
70  _init_system_object (libmesh_nullptr),
71  _assemble_system_function (libmesh_nullptr),
72  _assemble_system_object (libmesh_nullptr),
73  _constrain_system_function (libmesh_nullptr),
74  _constrain_system_object (libmesh_nullptr),
75  _qoi_evaluate_function (libmesh_nullptr),
76  _qoi_evaluate_object (libmesh_nullptr),
77  _qoi_evaluate_derivative_function (libmesh_nullptr),
78  _qoi_evaluate_derivative_object (libmesh_nullptr),
79  _dof_map (new DofMap(number_in, es.get_mesh())),
80  _equation_systems (es),
81  _mesh (es.get_mesh()),
82  _sys_name (name_in),
83  _sys_number (number_in),
84  _active (true),
85  _solution_projection (true),
86  _basic_system_only (false),
87  _is_initialized (false),
88  _identify_variable_groups (true),
89  _additional_data_written (false),
90  adjoint_already_solved (false),
91  _hide_output (false)
92 {
93 }
94 
95 
96 
97 // No copy construction of System objects!
98 System::System (const System & other) :
100  ParallelObject(other),
102  _mesh(other._mesh),
103  _sys_number(other._sys_number)
104 {
105  libmesh_not_implemented();
106 }
107 
108 
109 
111 {
112  libmesh_not_implemented();
113 }
114 
115 
117 {
118  // Null-out the function pointers. Since this
119  // class is getting destructed it is pointless,
120  // but a good habit.
124 
127 
128  // libmesh_nullptr-out user-provided objects.
134 
135  // Clear data
136  // Note: although clear() is virtual, C++ only calls
137  // the clear() of the base class in the destructor.
138  // Thus we add System namespace to make it clear.
139  System::clear ();
140 
141  libmesh_exceptionless_assert (!libMesh::closed());
142 }
143 
144 
145 
147 {
148  return _dof_map->n_dofs();
149 }
150 
151 
152 
154 {
155 #ifdef LIBMESH_ENABLE_CONSTRAINTS
156 
157  return _dof_map->n_constrained_dofs();
158 
159 #else
160 
161  return 0;
162 
163 #endif
164 }
165 
166 
167 
169 {
170 #ifdef LIBMESH_ENABLE_CONSTRAINTS
171 
172  return _dof_map->n_local_constrained_dofs();
173 
174 #else
175 
176  return 0;
177 
178 #endif
179 }
180 
181 
182 
184 {
185  return _dof_map->n_dofs_on_processor (this->processor_id());
186 }
187 
188 
189 
190 Number System::current_solution (const dof_id_type global_dof_number) const
191 {
192  // Check the sizes
193  libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
194  libmesh_assert_less (global_dof_number, current_local_solution->size());
195 
196  return (*current_local_solution)(global_dof_number);
197 }
198 
199 
200 
202 {
203  _variables.clear();
204 
205  _variable_numbers.clear();
206 
207  _dof_map->clear ();
208 
209  solution->clear ();
210 
211  current_local_solution->clear ();
212 
213  // clear any user-added vectors
214  {
215  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
216  {
217  pos->second->clear ();
218  delete pos->second;
219  pos->second = libmesh_nullptr;
220  }
221 
222  _vectors.clear();
223  _vector_projections.clear();
224  _vector_is_adjoint.clear();
225  _vector_types.clear();
226  _is_initialized = false;
227  }
228 
229 }
230 
231 
232 
234 {
235  // Calling init() twice on the same system currently works evil
236  // magic, whether done directly or via EquationSystems::read()
237  libmesh_assert(!this->is_initialized());
238 
239  // First initialize any required data:
240  // either only the basic System data
241  if (_basic_system_only)
243  // or all the derived class' data too
244  else
245  this->init_data();
246 
247  // If no variables have been added to this system
248  // don't do anything
249  if (!this->n_vars())
250  return;
251 
252  // Then call the user-provided initialization function
253  this->user_initialization();
254 }
255 
256 
257 
259 {
260  MeshBase & mesh = this->get_mesh();
261 
262  // Add all variable groups to our underlying DofMap
263  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
264  _dof_map->add_variable_group(this->variable_group(vg));
265 
266  // Distribute the degrees of freedom on the mesh
267  _dof_map->distribute_dofs (mesh);
268 
269  // Recreate any user or internal constraints
270  this->reinit_constraints();
271 
272  // And clean up the send_list before we first use it
273  _dof_map->prepare_send_list();
274 
275  // Resize the solution conformal to the current mesh
276  solution->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
277 
278  // Resize the current_local_solution for the current mesh
279 #ifdef LIBMESH_ENABLE_GHOSTED
280  current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
281  _dof_map->get_send_list(), false,
282  GHOSTED);
283 #else
284  current_local_solution->init (this->n_dofs(), false, SERIAL);
285 #endif
286 
287  // from now on, adding additional vectors or variables can't be done
288  // without immediately initializing them
289  _is_initialized = true;
290 
291  // initialize & zero other vectors, if necessary
292  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
293  {
294  ParallelType type = _vector_types[pos->first];
295 
296  if (type == GHOSTED)
297  {
298 #ifdef LIBMESH_ENABLE_GHOSTED
299  pos->second->init (this->n_dofs(), this->n_local_dofs(),
300  _dof_map->get_send_list(), false,
301  GHOSTED);
302 #else
303  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
304 #endif
305  }
306  else if (type == SERIAL)
307  {
308  pos->second->init (this->n_dofs(), false, type);
309  }
310  else
311  {
312  libmesh_assert_equal_to(type, PARALLEL);
313  pos->second->init (this->n_dofs(), this->n_local_dofs(), false, type);
314  }
315  }
316 }
317 
318 
319 
321 {
322 #ifdef LIBMESH_ENABLE_AMR
323  // Restrict the _vectors on the coarsened cells
324  for (vectors_iterator pos = _vectors.begin(); pos != _vectors.end(); ++pos)
325  {
326  NumericVector<Number> * v = pos->second;
327 
328  if (_vector_projections[pos->first])
329  {
330  this->project_vector (*v, this->vector_is_adjoint(pos->first));
331  }
332  else
333  {
334  ParallelType type = _vector_types[pos->first];
335 
336  if (type == GHOSTED)
337  {
338 #ifdef LIBMESH_ENABLE_GHOSTED
339  pos->second->init (this->n_dofs(), this->n_local_dofs(),
340  _dof_map->get_send_list(), false,
341  GHOSTED);
342 #else
343  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
344 #endif
345  }
346  else
347  pos->second->init (this->n_dofs(), this->n_local_dofs(), false, type);
348  }
349  }
350 
351  const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
352 
353  // Restrict the solution on the coarsened cells
355  this->project_vector (*solution);
356 
357 #ifdef LIBMESH_ENABLE_GHOSTED
358  current_local_solution->init(this->n_dofs(),
359  this->n_local_dofs(), send_list,
360  false, GHOSTED);
361 #else
362  current_local_solution->init(this->n_dofs());
363 #endif
364 
366  solution->localize (*current_local_solution, send_list);
367 
368 #endif // LIBMESH_ENABLE_AMR
369 }
370 
371 
372 
374 {
375 #ifdef LIBMESH_ENABLE_AMR
376  // Currently project_vector handles both restriction and prolongation
377  this->restrict_vectors();
378 #endif
379 }
380 
381 
382 
384 {
385  //If no variables have been added to this system
386  //don't do anything
387  if (!this->n_vars())
388  return;
389 
390  // Constraints get handled in EquationSystems::reinit now
391  // _dof_map->create_dof_constraints(this->get_mesh());
392 
393  // Update the solution based on the projected
394  // current_local_solution.
395  solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
396 
397  libmesh_assert_equal_to (solution->size(), current_local_solution->size());
398  // Not true with ghosted vectors
399  // libmesh_assert_equal_to (solution->size(), current_local_solution->local_size());
400 
401  const dof_id_type first_local_dof = solution->first_local_index();
402  const dof_id_type local_size = solution->local_size();
403 
404  for (dof_id_type i=0; i<local_size; i++)
405  solution->set(i+first_local_dof,
406  (*current_local_solution)(i+first_local_dof));
407 
408  solution->close();
409 }
410 
411 
413 {
414 #ifdef LIBMESH_ENABLE_CONSTRAINTS
416  user_constrain();
418 #endif
420 }
421 
422 
424 {
425  libmesh_assert(solution->closed());
426 
427  const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
428 
429  // Check sizes
430  libmesh_assert_equal_to (current_local_solution->size(), solution->size());
431  // More processors than elements => empty send_list
432  // libmesh_assert (!send_list.empty());
433  libmesh_assert_less_equal (send_list.size(), solution->size());
434 
435  // Create current_local_solution from solution. This will
436  // put a local copy of solution into current_local_solution.
437  // Only the necessary values (specified by the send_list)
438  // are copied to minimize communication
439  solution->localize (*current_local_solution, send_list);
440 }
441 
442 
443 
445 {
446  parallel_object_only();
447 
448  // If this system is empty... don't do anything!
449  if (!this->n_vars())
450  return;
451 
452  const std::vector<dof_id_type> & send_list = this->get_dof_map().get_send_list ();
453 
454  // Check sizes
455  libmesh_assert_equal_to (current_local_solution->size(), solution->size());
456  // Not true with ghosted vectors
457  // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
458  // libmesh_assert (!send_list.empty());
459  libmesh_assert_less_equal (send_list.size(), solution->size());
460 
461  // Create current_local_solution from solution. This will
462  // put a local copy of solution into current_local_solution.
463  solution->localize (*current_local_solution, send_list);
464 }
465 
466 
467 
469  const SubsetSolveMode /*subset_solve_mode*/)
470 {
471  if (subset != libmesh_nullptr)
472  libmesh_not_implemented();
473 }
474 
475 
476 
478 {
479  // Log how long the user's assembly code takes
480  LOG_SCOPE("assemble()", "System");
481 
482  // Call the user-specified assembly function
483  this->user_assembly();
484 }
485 
486 
487 
488 void System::assemble_qoi (const QoISet & qoi_indices)
489 {
490  // Log how long the user's assembly code takes
491  LOG_SCOPE("assemble_qoi()", "System");
492 
493  // Call the user-specified quantity of interest function
494  this->user_QOI(qoi_indices);
495 }
496 
497 
498 
499 void System::assemble_qoi_derivative(const QoISet & qoi_indices,
500  bool include_liftfunc,
501  bool apply_constraints)
502 {
503  // Log how long the user's assembly code takes
504  LOG_SCOPE("assemble_qoi_derivative()", "System");
505 
506  // Call the user-specified quantity of interest function
507  this->user_QOI_derivative(qoi_indices, include_liftfunc,
508  apply_constraints);
509 }
510 
511 
512 
513 void System::qoi_parameter_sensitivity (const QoISet & qoi_indices,
514  const ParameterVector & parameters,
515  SensitivityData & sensitivities)
516 {
517  // Forward sensitivities are more efficient for Nq > Np
518  if (qoi_indices.size(*this) > parameters.size())
519  forward_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
520  // Adjoint sensitivities are more efficient for Np > Nq,
521  // and an adjoint may be more reusable than a forward
522  // solution sensitivity in the Np == Nq case.
523  else
524  adjoint_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
525 }
526 
527 
528 
529 bool System::compare (const System & other_system,
530  const Real threshold,
531  const bool verbose) const
532 {
533  // we do not care for matrices, but for vectors
535  libmesh_assert (other_system._is_initialized);
536 
537  if (verbose)
538  {
539  libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
540  libMesh::out << " comparing matrices not supported." << std::endl;
541  libMesh::out << " comparing names...";
542  }
543 
544  // compare the name: 0 means identical
545  const int name_result = _sys_name.compare(other_system.name());
546  if (verbose)
547  {
548  if (name_result == 0)
549  libMesh::out << " identical." << std::endl;
550  else
551  libMesh::out << " names not identical." << std::endl;
552  libMesh::out << " comparing solution vector...";
553  }
554 
555 
556  // compare the solution: -1 means identical
557  const int solu_result = solution->compare (*other_system.solution.get(),
558  threshold);
559 
560  if (verbose)
561  {
562  if (solu_result == -1)
563  libMesh::out << " identical up to threshold." << std::endl;
564  else
565  libMesh::out << " first difference occurred at index = "
566  << solu_result << "." << std::endl;
567  }
568 
569 
570  // safety check, whether we handle at least the same number
571  // of vectors
572  std::vector<int> ov_result;
573 
574  if (this->n_vectors() != other_system.n_vectors())
575  {
576  if (verbose)
577  {
578  libMesh::out << " Fatal difference. This system handles "
579  << this->n_vectors() << " add'l vectors," << std::endl
580  << " while the other system handles "
581  << other_system.n_vectors()
582  << " add'l vectors." << std::endl
583  << " Aborting comparison." << std::endl;
584  }
585  return false;
586  }
587  else if (this->n_vectors() == 0)
588  {
589  // there are no additional vectors...
590  ov_result.clear ();
591  }
592  else
593  {
594  // compare other vectors
595  for (const_vectors_iterator pos = _vectors.begin();
596  pos != _vectors.end(); ++pos)
597  {
598  if (verbose)
599  libMesh::out << " comparing vector \""
600  << pos->first << "\" ...";
601 
602  // assume they have the same name
603  const NumericVector<Number> & other_system_vector =
604  other_system.get_vector(pos->first);
605 
606  ov_result.push_back(pos->second->compare (other_system_vector,
607  threshold));
608 
609  if (verbose)
610  {
611  if (ov_result[ov_result.size()-1] == -1)
612  libMesh::out << " identical up to threshold." << std::endl;
613  else
614  libMesh::out << " first difference occurred at" << std::endl
615  << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
616  }
617 
618  }
619 
620  } // finished comparing additional vectors
621 
622 
623  bool overall_result;
624 
625  // sum up the results
626  if ((name_result==0) && (solu_result==-1))
627  {
628  if (ov_result.size()==0)
629  overall_result = true;
630  else
631  {
632  bool ov_identical;
633  unsigned int n = 0;
634  do
635  {
636  ov_identical = (ov_result[n]==-1);
637  n++;
638  }
639  while (ov_identical && n<ov_result.size());
640  overall_result = ov_identical;
641  }
642  }
643  else
644  overall_result = false;
645 
646  if (verbose)
647  {
648  libMesh::out << " finished comparisons, ";
649  if (overall_result)
650  libMesh::out << "found no differences." << std::endl << std::endl;
651  else
652  libMesh::out << "found differences." << std::endl << std::endl;
653  }
654 
655  return overall_result;
656 }
657 
658 
659 
660 void System::update_global_solution (std::vector<Number> & global_soln) const
661 {
662  global_soln.resize (solution->size());
663 
664  solution->localize (global_soln);
665 }
666 
667 
668 
669 void System::update_global_solution (std::vector<Number> & global_soln,
670  const processor_id_type dest_proc) const
671 {
672  global_soln.resize (solution->size());
673 
674  solution->localize_to_one (global_soln, dest_proc);
675 }
676 
677 
678 
679 NumericVector<Number> & System::add_vector (const std::string & vec_name,
680  const bool projections,
681  const ParallelType type)
682 {
683  // Return the vector if it is already there.
684  if (this->have_vector(vec_name))
685  return *(_vectors[vec_name]);
686 
687  // Otherwise build the vector
688  NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
689  _vectors.insert (std::make_pair (vec_name, buf));
690  _vector_projections.insert (std::make_pair (vec_name, projections));
691 
692  _vector_types.insert (std::make_pair (vec_name, type));
693 
694  // Vectors are primal by default
695  _vector_is_adjoint.insert (std::make_pair (vec_name, -1));
696 
697  // Initialize it if necessary
698  if (_is_initialized)
699  {
700  if (type == GHOSTED)
701  {
702 #ifdef LIBMESH_ENABLE_GHOSTED
703  buf->init (this->n_dofs(), this->n_local_dofs(),
704  _dof_map->get_send_list(), false,
705  GHOSTED);
706 #else
707  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
708 #endif
709  }
710  else
711  buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
712  }
713 
714  return *buf;
715 }
716 
717 void System::remove_vector (const std::string & vec_name)
718 {
719  vectors_iterator pos = _vectors.find(vec_name);
720 
721  //Return if the vector does not exist
722  if (pos == _vectors.end())
723  return;
724 
725  delete pos->second;
726 
727  _vectors.erase(pos);
728 
729  _vector_projections.erase(vec_name);
730  _vector_is_adjoint.erase(vec_name);
731  _vector_types.erase(vec_name);
732 }
733 
734 const NumericVector<Number> * System::request_vector (const std::string & vec_name) const
735 {
736  const_vectors_iterator pos = _vectors.find(vec_name);
737 
738  if (pos == _vectors.end())
739  return libmesh_nullptr;
740 
741  return pos->second;
742 }
743 
744 
745 
746 NumericVector<Number> * System::request_vector (const std::string & vec_name)
747 {
748  vectors_iterator pos = _vectors.find(vec_name);
749 
750  if (pos == _vectors.end())
751  return libmesh_nullptr;
752 
753  return pos->second;
754 }
755 
756 
757 
758 const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
759 {
762  unsigned int num = 0;
763  while ((num<vec_num) && (v!=v_end))
764  {
765  num++;
766  ++v;
767  }
768  if (v==v_end)
769  return libmesh_nullptr;
770  return v->second;
771 }
772 
773 
774 
775 NumericVector<Number> * System::request_vector (const unsigned int vec_num)
776 {
778  vectors_iterator v_end = vectors_end();
779  unsigned int num = 0;
780  while ((num<vec_num) && (v!=v_end))
781  {
782  num++;
783  ++v;
784  }
785  if (v==v_end)
786  return libmesh_nullptr;
787  return v->second;
788 }
789 
790 
791 
792 const NumericVector<Number> & System::get_vector (const std::string & vec_name) const
793 {
794  // Make sure the vector exists
795  const_vectors_iterator pos = _vectors.find(vec_name);
796 
797  if (pos == _vectors.end())
798  libmesh_error_msg("ERROR: vector " << vec_name << " does not exist in this system!");
799 
800  return *(pos->second);
801 }
802 
803 
804 
805 NumericVector<Number> & System::get_vector (const std::string & vec_name)
806 {
807  // Make sure the vector exists
808  vectors_iterator pos = _vectors.find(vec_name);
809 
810  if (pos == _vectors.end())
811  libmesh_error_msg("ERROR: vector " << vec_name << " does not exist in this system!");
812 
813  return *(pos->second);
814 }
815 
816 
817 
818 const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
819 {
822  unsigned int num = 0;
823  while ((num<vec_num) && (v!=v_end))
824  {
825  num++;
826  ++v;
827  }
828  libmesh_assert (v != v_end);
829  return *(v->second);
830 }
831 
832 
833 
834 NumericVector<Number> & System::get_vector (const unsigned int vec_num)
835 {
837  vectors_iterator v_end = vectors_end();
838  unsigned int num = 0;
839  while ((num<vec_num) && (v!=v_end))
840  {
841  num++;
842  ++v;
843  }
844  libmesh_assert (v != v_end);
845  return *(v->second);
846 }
847 
848 
849 
850 const std::string & System::vector_name (const unsigned int vec_num) const
851 {
854  unsigned int num = 0;
855  while ((num<vec_num) && (v!=v_end))
856  {
857  num++;
858  ++v;
859  }
860  libmesh_assert (v != v_end);
861  return v->first;
862 }
863 
864 const std::string & System::vector_name (const NumericVector<Number> & vec_reference) const
865 {
868 
869  for (; v != v_end; ++v)
870  {
871  // Check if the current vector is the one whose name we want
872  if (&vec_reference == v->second)
873  break; // exit loop if it is
874  }
875 
876  // Before returning, make sure we didnt loop till the end and not find any match
877  libmesh_assert (v != v_end);
878 
879  // Return the string associated with the current vector
880  return v->first;
881 }
882 
883 
884 
885 void System::set_vector_preservation (const std::string & vec_name,
886  bool preserve)
887 {
888  _vector_projections[vec_name] = preserve;
889 }
890 
891 
892 
893 bool System::vector_preservation (const std::string & vec_name) const
894 {
895  if (_vector_projections.find(vec_name) == _vector_projections.end())
896  return false;
897 
898  return _vector_projections.find(vec_name)->second;
899 }
900 
901 
902 
903 void System::set_vector_as_adjoint (const std::string & vec_name,
904  int qoi_num)
905 {
906  // We reserve -1 for vectors which get primal constraints, -2 for
907  // vectors which get no constraints
908  libmesh_assert_greater_equal(qoi_num, -2);
909  _vector_is_adjoint[vec_name] = qoi_num;
910 }
911 
912 
913 
914 int System::vector_is_adjoint (const std::string & vec_name) const
915 {
916  libmesh_assert(_vector_is_adjoint.find(vec_name) !=
917  _vector_is_adjoint.end());
918 
919  return _vector_is_adjoint.find(vec_name)->second;
920 }
921 
922 
923 
925 {
926  std::ostringstream sensitivity_name;
927  sensitivity_name << "sensitivity_solution" << i;
928 
929  return this->add_vector(sensitivity_name.str());
930 }
931 
932 
933 
935 {
936  std::ostringstream sensitivity_name;
937  sensitivity_name << "sensitivity_solution" << i;
938 
939  return this->get_vector(sensitivity_name.str());
940 }
941 
942 
943 
945 {
946  std::ostringstream sensitivity_name;
947  sensitivity_name << "sensitivity_solution" << i;
948 
949  return this->get_vector(sensitivity_name.str());
950 }
951 
952 
953 
955 {
956  return this->add_vector("weighted_sensitivity_solution");
957 }
958 
959 
960 
962 {
963  return this->get_vector("weighted_sensitivity_solution");
964 }
965 
966 
967 
969 {
970  return this->get_vector("weighted_sensitivity_solution");
971 }
972 
973 
974 
976 {
977  std::ostringstream adjoint_name;
978  adjoint_name << "adjoint_solution" << i;
979 
980  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
981  this->set_vector_as_adjoint(adjoint_name.str(), i);
982  return returnval;
983 }
984 
985 
986 
988 {
989  std::ostringstream adjoint_name;
990  adjoint_name << "adjoint_solution" << i;
991 
992  return this->get_vector(adjoint_name.str());
993 }
994 
995 
996 
998 {
999  std::ostringstream adjoint_name;
1000  adjoint_name << "adjoint_solution" << i;
1001 
1002  return this->get_vector(adjoint_name.str());
1003 }
1004 
1005 
1006 
1008 {
1009  std::ostringstream adjoint_name;
1010  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1011 
1012  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1013  this->set_vector_as_adjoint(adjoint_name.str(), i);
1014  return returnval;
1015 }
1016 
1017 
1018 
1020 {
1021  std::ostringstream adjoint_name;
1022  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1023 
1024  return this->get_vector(adjoint_name.str());
1025 }
1026 
1027 
1028 
1030 {
1031  std::ostringstream adjoint_name;
1032  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1033 
1034  return this->get_vector(adjoint_name.str());
1035 }
1036 
1037 
1038 
1040 {
1041  std::ostringstream adjoint_rhs_name;
1042  adjoint_rhs_name << "adjoint_rhs" << i;
1043 
1044  return this->add_vector(adjoint_rhs_name.str(), false);
1045 }
1046 
1047 
1048 
1050 {
1051  std::ostringstream adjoint_rhs_name;
1052  adjoint_rhs_name << "adjoint_rhs" << i;
1053 
1054  return this->get_vector(adjoint_rhs_name.str());
1055 }
1056 
1057 
1058 
1059 const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
1060 {
1061  std::ostringstream adjoint_rhs_name;
1062  adjoint_rhs_name << "adjoint_rhs" << i;
1063 
1064  return this->get_vector(adjoint_rhs_name.str());
1065 }
1066 
1067 
1068 
1070 {
1071  std::ostringstream sensitivity_rhs_name;
1072  sensitivity_rhs_name << "sensitivity_rhs" << i;
1073 
1074  return this->add_vector(sensitivity_rhs_name.str(), false);
1075 }
1076 
1077 
1078 
1080 {
1081  std::ostringstream sensitivity_rhs_name;
1082  sensitivity_rhs_name << "sensitivity_rhs" << i;
1083 
1084  return this->get_vector(sensitivity_rhs_name.str());
1085 }
1086 
1087 
1088 
1090 {
1091  std::ostringstream sensitivity_rhs_name;
1092  sensitivity_rhs_name << "sensitivity_rhs" << i;
1093 
1094  return this->get_vector(sensitivity_rhs_name.str());
1095 }
1096 
1097 
1098 
1099 unsigned int System::add_variable (const std::string & var,
1100  const FEType & type,
1101  const std::set<subdomain_id_type> * const active_subdomains)
1102 {
1103  libmesh_assert(!this->is_initialized());
1104 
1105  // Make sure the variable isn't there already
1106  // or if it is, that it's the type we want
1107  for (unsigned int v=0; v<this->n_vars(); v++)
1108  if (this->variable_name(v) == var)
1109  {
1110  if (this->variable_type(v) == type)
1111  return _variables[v].number();
1112 
1113  libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1114  }
1115 
1116  // Optimize for VariableGroups here - if the user is adding multiple
1117  // variables of the same FEType and subdomain restriction, catch
1118  // that here and add them as members of the same VariableGroup.
1119  //
1120  // start by setting this flag to whatever the user has requested
1121  // and then consider the conditions which should negate it.
1122  bool should_be_in_vg = this->identify_variable_groups();
1123 
1124  // No variable groups, nothing to add to
1125  if (!this->n_variable_groups())
1126  should_be_in_vg = false;
1127 
1128  else
1129  {
1130  VariableGroup & vg(_variable_groups.back());
1131 
1132  // get a pointer to their subdomain restriction, if any.
1133  const std::set<subdomain_id_type> * const
1134  their_active_subdomains (vg.implicitly_active() ?
1135  libmesh_nullptr : &vg.active_subdomains());
1136 
1137  // Different types?
1138  if (vg.type() != type)
1139  should_be_in_vg = false;
1140 
1141  // they are restricted, we aren't?
1142  if (their_active_subdomains && !active_subdomains)
1143  should_be_in_vg = false;
1144 
1145  // they aren't restricted, we are?
1146  if (!their_active_subdomains && active_subdomains)
1147  should_be_in_vg = false;
1148 
1149  if (their_active_subdomains && active_subdomains)
1150  // restricted to different sets?
1151  if (*their_active_subdomains != *active_subdomains)
1152  should_be_in_vg = false;
1153 
1154  // OK, after all that, append the variable to the vg if none of the conditions
1155  // were violated
1156  if (should_be_in_vg)
1157  {
1158  const unsigned short curr_n_vars = cast_int<unsigned short>
1159  (this->n_vars());
1160 
1161  vg.append (var);
1162 
1163  _variables.push_back(vg(vg.n_variables()-1));
1164  _variable_numbers[var] = curr_n_vars;
1165  return curr_n_vars;
1166  }
1167  }
1168 
1169  // otherwise, fall back to adding a single variable group
1170  return this->add_variables (std::vector<std::string>(1, var),
1171  type,
1172  active_subdomains);
1173 }
1174 
1175 
1176 
1177 unsigned int System::add_variable (const std::string & var,
1178  const Order order,
1179  const FEFamily family,
1180  const std::set<subdomain_id_type> * const active_subdomains)
1181 {
1182  return this->add_variable(var,
1183  FEType(order, family),
1184  active_subdomains);
1185 }
1186 
1187 
1188 
1189 unsigned int System::add_variables (const std::vector<std::string> & vars,
1190  const FEType & type,
1191  const std::set<subdomain_id_type> * const active_subdomains)
1192 {
1193  libmesh_assert(!this->is_initialized());
1194 
1195  // Make sure the variable isn't there already
1196  // or if it is, that it's the type we want
1197  for (std::size_t ov=0; ov<vars.size(); ov++)
1198  for (unsigned int v=0; v<this->n_vars(); v++)
1199  if (this->variable_name(v) == vars[ov])
1200  {
1201  if (this->variable_type(v) == type)
1202  return _variables[v].number();
1203 
1204  libmesh_error_msg("ERROR: incompatible variable " << vars[ov] << " has already been added for this system!");
1205  }
1206 
1207  const unsigned short curr_n_vars = cast_int<unsigned short>
1208  (this->n_vars());
1209 
1210  const unsigned int next_first_component = this->n_components();
1211 
1212  // Add the variable group to the list
1213  _variable_groups.push_back((active_subdomains == libmesh_nullptr) ?
1214  VariableGroup(this, vars, curr_n_vars,
1215  next_first_component, type) :
1216  VariableGroup(this, vars, curr_n_vars,
1217  next_first_component, type, *active_subdomains));
1218 
1219  const VariableGroup & vg (_variable_groups.back());
1220 
1221  // Add each component of the group individually
1222  for (std::size_t v=0; v<vars.size(); v++)
1223  {
1224  _variables.push_back (vg(v));
1225  _variable_numbers[vars[v]] = cast_int<unsigned short>
1226  (curr_n_vars+v);
1227  }
1228 
1229  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1230 
1231  // BSK - Defer this now to System::init_data() so we can detect
1232  // VariableGroups 12/28/2012
1233  // // Add the variable group to the _dof_map
1234  // _dof_map->add_variable_group (vg);
1235 
1236  // Return the number of the new variable
1237  return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1238 }
1239 
1240 
1241 
1242 unsigned int System::add_variables (const std::vector<std::string> & vars,
1243  const Order order,
1244  const FEFamily family,
1245  const std::set<subdomain_id_type> * const active_subdomains)
1246 {
1247  return this->add_variables(vars,
1248  FEType(order, family),
1249  active_subdomains);
1250 }
1251 
1252 
1253 
1254 bool System::has_variable (const std::string & var) const
1255 {
1256  return _variable_numbers.count(var);
1257 }
1258 
1259 
1260 
1261 unsigned short int System::variable_number (const std::string & var) const
1262 {
1263  // Make sure the variable exists
1264  std::map<std::string, unsigned short int>::const_iterator
1265  pos = _variable_numbers.find(var);
1266 
1267  if (pos == _variable_numbers.end())
1268  libmesh_error_msg("ERROR: variable " << var << " does not exist in this system!");
1269 
1270  libmesh_assert_equal_to (_variables[pos->second].name(), var);
1271 
1272  return pos->second;
1273 }
1274 
1275 
1276 void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
1277 {
1278  all_variable_numbers.resize(n_vars());
1279 
1280  // Make sure the variable exists
1281  std::map<std::string, unsigned short int>::const_iterator
1282  it = _variable_numbers.begin();
1283  std::map<std::string, unsigned short int>::const_iterator
1284  it_end = _variable_numbers.end();
1285 
1286  unsigned int count = 0;
1287  for ( ; it != it_end; ++it)
1288  {
1289  all_variable_numbers[count] = it->second;
1290  count++;
1291  }
1292 }
1293 
1294 
1295 void System::local_dof_indices(const unsigned int var,
1296  std::set<dof_id_type> & var_indices) const
1297 {
1298  // Make sure the set is clear
1299  var_indices.clear();
1300 
1301  std::vector<dof_id_type> dof_indices;
1302 
1303  // Begin the loop over the elements
1306  const MeshBase::const_element_iterator end_el =
1308 
1309  const dof_id_type
1310  first_local = this->get_dof_map().first_dof(),
1311  end_local = this->get_dof_map().end_dof();
1312 
1313  for ( ; el != end_el; ++el)
1314  {
1315  const Elem * elem = *el;
1316  this->get_dof_map().dof_indices (elem, dof_indices, var);
1317 
1318  for (std::size_t i=0; i<dof_indices.size(); i++)
1319  {
1320  dof_id_type dof = dof_indices[i];
1321 
1322  //If the dof is owned by the local processor
1323  if (first_local <= dof && dof < end_local)
1324  var_indices.insert(dof_indices[i]);
1325  }
1326  }
1327 }
1328 
1329 
1330 
1332  unsigned int var_num) const
1333 {
1334  /* Make sure the call makes sense. */
1335  libmesh_assert_less (var_num, this->n_vars());
1336 
1337  /* Get a reference to the mesh. */
1338  const MeshBase & mesh = this->get_mesh();
1339 
1340  /* Check which system we are. */
1341  const unsigned int sys_num = this->number();
1342 
1343  /* Loop over nodes. */
1344  {
1346  const MeshBase::const_node_iterator end_it = mesh.local_nodes_end();
1347  for ( ; it != end_it; ++it)
1348  {
1349  const Node * node = *it;
1350  unsigned int n_comp = node->n_comp(sys_num,var_num);
1351  for (unsigned int i=0; i<n_comp; i++)
1352  {
1353  const dof_id_type index = node->dof_number(sys_num,var_num,i);
1354  v.set(index,0.0);
1355  }
1356  }
1357  }
1358 
1359  /* Loop over elements. */
1360  {
1363  for ( ; it != end_it; ++it)
1364  {
1365  const Elem * elem = *it;
1366  unsigned int n_comp = elem->n_comp(sys_num,var_num);
1367  for (unsigned int i=0; i<n_comp; i++)
1368  {
1369  const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1370  v.set(index,0.0);
1371  }
1372  }
1373  }
1374 }
1375 
1376 
1377 
1379  unsigned int var,
1380  FEMNormType norm_type) const
1381 {
1382  std::set<dof_id_type> var_indices;
1383  local_dof_indices(var, var_indices);
1384 
1385  if (norm_type == DISCRETE_L1)
1386  return v.subset_l1_norm(var_indices);
1387  if (norm_type == DISCRETE_L2)
1388  return v.subset_l2_norm(var_indices);
1389  if (norm_type == DISCRETE_L_INF)
1390  return v.subset_linfty_norm(var_indices);
1391  else
1392  libmesh_error_msg("Invalid norm_type = " << norm_type);
1393 }
1394 
1395 
1396 
1398  unsigned int var,
1399  FEMNormType norm_type,
1400  std::set<unsigned int> * skip_dimensions) const
1401 {
1402  //short circuit to save time
1403  if (norm_type == DISCRETE_L1 ||
1404  norm_type == DISCRETE_L2 ||
1405  norm_type == DISCRETE_L_INF)
1406  return discrete_var_norm(v,var,norm_type);
1407 
1408  // Not a discrete norm
1409  std::vector<FEMNormType> norms(this->n_vars(), L2);
1410  std::vector<Real> weights(this->n_vars(), 0.0);
1411  norms[var] = norm_type;
1412  weights[var] = 1.0;
1413  Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1414  return val;
1415 }
1416 
1417 
1418 
1420  const SystemNorm & norm,
1421  std::set<unsigned int> * skip_dimensions) const
1422 {
1423  // This function must be run on all processors at once
1424  parallel_object_only();
1425 
1426  LOG_SCOPE ("calculate_norm()", "System");
1427 
1428  // Zero the norm before summation
1429  Real v_norm = 0.;
1430 
1431  if (norm.is_discrete())
1432  {
1433  //Check to see if all weights are 1.0 and all types are equal
1434  FEMNormType norm_type0 = norm.type(0);
1435  unsigned int check_var = 0;
1436  for (; check_var != this->n_vars(); ++check_var)
1437  if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1438  break;
1439 
1440  //All weights were 1.0 so just do the full vector discrete norm
1441  if (check_var == this->n_vars())
1442  {
1443  if (norm_type0 == DISCRETE_L1)
1444  return v.l1_norm();
1445  if (norm_type0 == DISCRETE_L2)
1446  return v.l2_norm();
1447  if (norm_type0 == DISCRETE_L_INF)
1448  return v.linfty_norm();
1449  else
1450  libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
1451  }
1452 
1453  for (unsigned int var=0; var != this->n_vars(); ++var)
1454  {
1455  // Skip any variables we don't need to integrate
1456  if (norm.weight(var) == 0.0)
1457  continue;
1458 
1459  v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1460  }
1461 
1462  return v_norm;
1463  }
1464 
1465  // Localize the potentially parallel vector
1467  local_v->init(v.size(), true, SERIAL);
1468  v.localize (*local_v, _dof_map->get_send_list());
1469 
1470  // I'm not sure how best to mix Hilbert norms on some variables (for
1471  // which we'll want to square then sum then square root) with norms
1472  // like L_inf (for which we'll just want to take an absolute value
1473  // and then sum).
1474  bool using_hilbert_norm = true,
1475  using_nonhilbert_norm = true;
1476 
1477  // Loop over all variables
1478  for (unsigned int var=0; var != this->n_vars(); ++var)
1479  {
1480  // Skip any variables we don't need to integrate
1481  Real norm_weight_sq = norm.weight_sq(var);
1482  if (norm_weight_sq == 0.0)
1483  continue;
1484  Real norm_weight = norm.weight(var);
1485 
1486  // Check for unimplemented norms (rather than just returning 0).
1487  FEMNormType norm_type = norm.type(var);
1488  if ((norm_type==H1) ||
1489  (norm_type==H2) ||
1490  (norm_type==L2) ||
1491  (norm_type==H1_SEMINORM) ||
1492  (norm_type==H2_SEMINORM))
1493  {
1494  if (!using_hilbert_norm)
1495  libmesh_not_implemented();
1496  using_nonhilbert_norm = false;
1497  }
1498  else if ((norm_type==L1) ||
1499  (norm_type==L_INF) ||
1500  (norm_type==W1_INF_SEMINORM) ||
1501  (norm_type==W2_INF_SEMINORM))
1502  {
1503  if (!using_nonhilbert_norm)
1504  libmesh_not_implemented();
1505  using_hilbert_norm = false;
1506  }
1507  else
1508  libmesh_not_implemented();
1509 
1510  const FEType & fe_type = this->get_dof_map().variable_type(var);
1511 
1512  // Allow space for dims 0-3, even if we don't use them all
1513  std::vector<FEBase *> fe_ptrs(4,libmesh_nullptr);
1514  std::vector<QBase *> q_rules(4,libmesh_nullptr);
1515 
1516  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1517 
1518  // Prepare finite elements for each dimension present in the mesh
1519  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
1520  d_it != elem_dims.end(); ++d_it)
1521  {
1522  if (skip_dimensions && skip_dimensions->find(*d_it) != skip_dimensions->end())
1523  continue;
1524 
1525  q_rules[*d_it] =
1526  fe_type.default_quadrature_rule (*d_it).release();
1527 
1528  // Construct finite element object
1529 
1530  fe_ptrs[*d_it] = FEBase::build(*d_it, fe_type).release();
1531 
1532  // Attach quadrature rule to FE object
1533  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it]);
1534  }
1535 
1536  std::vector<dof_id_type> dof_indices;
1537 
1538  // Begin the loop over the elements
1541  const MeshBase::const_element_iterator end_el =
1543 
1544  for ( ; el != end_el; ++el)
1545  {
1546  const Elem * elem = *el;
1547  const unsigned int dim = elem->dim();
1548 
1549  if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1550  continue;
1551 
1552  FEBase * fe = fe_ptrs[dim];
1553  QBase * qrule = q_rules[dim];
1554  libmesh_assert(fe);
1555  libmesh_assert(qrule);
1556 
1557  const std::vector<Real> & JxW = fe->get_JxW();
1558  const std::vector<std::vector<Real> > * phi = libmesh_nullptr;
1559  if (norm_type == H1 ||
1560  norm_type == H2 ||
1561  norm_type == L2 ||
1562  norm_type == L1 ||
1563  norm_type == L_INF)
1564  phi = &(fe->get_phi());
1565 
1566  const std::vector<std::vector<RealGradient> > * dphi = libmesh_nullptr;
1567  if (norm_type == H1 ||
1568  norm_type == H2 ||
1569  norm_type == H1_SEMINORM ||
1570  norm_type == W1_INF_SEMINORM)
1571  dphi = &(fe->get_dphi());
1572 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1573  const std::vector<std::vector<RealTensor> > * d2phi = libmesh_nullptr;
1574  if (norm_type == H2 ||
1575  norm_type == H2_SEMINORM ||
1576  norm_type == W2_INF_SEMINORM)
1577  d2phi = &(fe->get_d2phi());
1578 #endif
1579 
1580  fe->reinit (elem);
1581 
1582  this->get_dof_map().dof_indices (elem, dof_indices, var);
1583 
1584  const unsigned int n_qp = qrule->n_points();
1585 
1586  const unsigned int n_sf = cast_int<unsigned int>
1587  (dof_indices.size());
1588 
1589  // Begin the loop over the Quadrature points.
1590  for (unsigned int qp=0; qp<n_qp; qp++)
1591  {
1592  if (norm_type == L1)
1593  {
1594  Number u_h = 0.;
1595  for (unsigned int i=0; i != n_sf; ++i)
1596  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1597  v_norm += norm_weight *
1598  JxW[qp] * std::abs(u_h);
1599  }
1600 
1601  if (norm_type == L_INF)
1602  {
1603  Number u_h = 0.;
1604  for (unsigned int i=0; i != n_sf; ++i)
1605  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1606  v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1607  }
1608 
1609  if (norm_type == H1 ||
1610  norm_type == H2 ||
1611  norm_type == L2)
1612  {
1613  Number u_h = 0.;
1614  for (unsigned int i=0; i != n_sf; ++i)
1615  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1616  v_norm += norm_weight_sq *
1617  JxW[qp] * TensorTools::norm_sq(u_h);
1618  }
1619 
1620  if (norm_type == H1 ||
1621  norm_type == H2 ||
1622  norm_type == H1_SEMINORM)
1623  {
1624  Gradient grad_u_h;
1625  for (unsigned int i=0; i != n_sf; ++i)
1626  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1627  v_norm += norm_weight_sq *
1628  JxW[qp] * grad_u_h.norm_sq();
1629  }
1630 
1631  if (norm_type == W1_INF_SEMINORM)
1632  {
1633  Gradient grad_u_h;
1634  for (unsigned int i=0; i != n_sf; ++i)
1635  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1636  v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1637  }
1638 
1639 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1640  if (norm_type == H2 ||
1641  norm_type == H2_SEMINORM)
1642  {
1643  Tensor hess_u_h;
1644  for (unsigned int i=0; i != n_sf; ++i)
1645  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1646  v_norm += norm_weight_sq *
1647  JxW[qp] * hess_u_h.norm_sq();
1648  }
1649 
1650  if (norm_type == W2_INF_SEMINORM)
1651  {
1652  Tensor hess_u_h;
1653  for (unsigned int i=0; i != n_sf; ++i)
1654  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1655  v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
1656  }
1657 #endif
1658  }
1659  }
1660 
1661  // Need to delete the FE and quadrature objects to prevent a memory leak
1662  for (std::size_t i=0; i<fe_ptrs.size(); i++)
1663  if (fe_ptrs[i])
1664  delete fe_ptrs[i];
1665 
1666  for (std::size_t i=0; i<q_rules.size(); i++)
1667  if (q_rules[i])
1668  delete q_rules[i];
1669  }
1670 
1671  if (using_hilbert_norm)
1672  {
1673  this->comm().sum(v_norm);
1674  v_norm = std::sqrt(v_norm);
1675  }
1676  else
1677  {
1678  this->comm().max(v_norm);
1679  }
1680 
1681  return v_norm;
1682 }
1683 
1684 
1685 
1686 std::string System::get_info() const
1687 {
1688  std::ostringstream oss;
1689 
1690 
1691  const std::string & sys_name = this->name();
1692 
1693  oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1694  << " Type \"" << this->system_type() << "\"\n"
1695  << " Variables=";
1696 
1697  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1698  {
1699  const VariableGroup & vg_description (this->variable_group(vg));
1700 
1701  if (vg_description.n_variables() > 1) oss << "{ ";
1702  for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
1703  oss << "\"" << vg_description.name(vn) << "\" ";
1704  if (vg_description.n_variables() > 1) oss << "} ";
1705  }
1706 
1707  oss << '\n';
1708 
1709  oss << " Finite Element Types=";
1710 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1711  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1712  oss << "\""
1713  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1714  << "\" ";
1715 #else
1716  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1717  {
1718  oss << "\""
1719  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1720  << "\", \""
1721  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1722  << "\" ";
1723  }
1724 
1725  oss << '\n' << " Infinite Element Mapping=";
1726  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1727  oss << "\""
1728  << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1729  << "\" ";
1730 #endif
1731 
1732  oss << '\n';
1733 
1734  oss << " Approximation Orders=";
1735  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1736  {
1737 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1738  oss << "\""
1739  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1740  << "\" ";
1741 #else
1742  oss << "\""
1743  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1744  << "\", \""
1745  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1746  << "\" ";
1747 #endif
1748  }
1749 
1750  oss << '\n';
1751 
1752  oss << " n_dofs()=" << this->n_dofs() << '\n';
1753  oss << " n_local_dofs()=" << this->n_local_dofs() << '\n';
1754 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1755  oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1756  oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1757 #endif
1758 
1759  oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1760  oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1761  // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1762 
1763  oss << this->get_dof_map().get_info();
1764 
1765  return oss.str();
1766 }
1767 
1768 
1769 
1771  const std::string & name))
1772 {
1773  libmesh_assert(fptr);
1774 
1776  {
1777  libmesh_here();
1778  libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1779  << std::endl;
1780 
1782  }
1783 
1784  _init_system_function = fptr;
1785 }
1786 
1787 
1788 
1790 {
1792  {
1793  libmesh_here();
1794  libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1795  << std::endl;
1796 
1798  }
1799 
1800  _init_system_object = &init_in;
1801 }
1802 
1803 
1804 
1806  const std::string & name))
1807 {
1808  libmesh_assert(fptr);
1809 
1811  {
1812  libmesh_here();
1813  libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1814  << std::endl;
1815 
1817  }
1818 
1820 }
1821 
1822 
1823 
1825 {
1827  {
1828  libmesh_here();
1829  libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1830  << std::endl;
1831 
1833  }
1834 
1835  _assemble_system_object = &assemble_in;
1836 }
1837 
1838 
1839 
1841  const std::string & name))
1842 {
1843  libmesh_assert(fptr);
1844 
1846  {
1847  libmesh_here();
1848  libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1849  << std::endl;
1850 
1852  }
1853 
1855 }
1856 
1857 
1858 
1860 {
1862  {
1863  libmesh_here();
1864  libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1865  << std::endl;
1866 
1868  }
1869 
1870  _constrain_system_object = &constrain;
1871 }
1872 
1873 
1874 
1876  const std::string &,
1877  const QoISet &))
1878 {
1879  libmesh_assert(fptr);
1880 
1882  {
1883  libmesh_here();
1884  libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1885  << std::endl;
1886 
1888  }
1889 
1890  _qoi_evaluate_function = fptr;
1891 }
1892 
1893 
1894 
1896 {
1898  {
1899  libmesh_here();
1900  libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1901  << std::endl;
1902 
1904  }
1905 
1906  _qoi_evaluate_object = &qoi_in;
1907 }
1908 
1909 
1910 
1911 void System::attach_QOI_derivative(void fptr(EquationSystems &, const std::string &,
1912  const QoISet &, bool, bool))
1913 {
1914  libmesh_assert(fptr);
1915 
1917  {
1918  libmesh_here();
1919  libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1920  << std::endl;
1921 
1923  }
1924 
1926 }
1927 
1928 
1929 
1931 {
1933  {
1934  libmesh_here();
1935  libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1936  << std::endl;
1937 
1939  }
1940 
1941  _qoi_evaluate_derivative_object = &qoi_derivative;
1942 }
1943 
1944 
1945 
1947 {
1948  // Call the user-provided initialization function,
1949  // if it was provided
1951  this->_init_system_function (_equation_systems, this->name());
1952 
1953  // ...or the user-provided initialization object.
1956 }
1957 
1958 
1959 
1961 {
1962  // Call the user-provided assembly function,
1963  // if it was provided
1966 
1967  // ...or the user-provided assembly object.
1970 }
1971 
1972 
1973 
1975 {
1976  // Call the user-provided constraint function,
1977  // if it was provided
1980 
1981  // ...or the user-provided constraint object.
1984 }
1985 
1986 
1987 
1988 void System::user_QOI (const QoISet & qoi_indices)
1989 {
1990  // Call the user-provided quantity of interest function,
1991  // if it was provided
1993  this->_qoi_evaluate_function(_equation_systems, this->name(), qoi_indices);
1994 
1995  // ...or the user-provided QOI function object.
1997  this->_qoi_evaluate_object->qoi(qoi_indices);
1998 }
1999 
2000 
2001 
2002 void System::user_QOI_derivative(const QoISet & qoi_indices,
2003  bool include_liftfunc,
2004  bool apply_constraints)
2005 {
2006  // Call the user-provided quantity of interest derivative,
2007  // if it was provided
2010  (_equation_systems, this->name(), qoi_indices, include_liftfunc,
2011  apply_constraints);
2012 
2013  // ...or the user-provided QOI derivative function object.
2016  (qoi_indices, include_liftfunc, apply_constraints);
2017 }
2018 
2019 
2020 
2021 Number System::point_value(unsigned int var, const Point & p, const bool insist_on_success) const
2022 {
2023  // This function must be called on every processor; there's no
2024  // telling where in the partition p falls.
2025  parallel_object_only();
2026 
2027  // And every processor had better agree about which point we're
2028  // looking for
2029 #ifndef NDEBUG
2030  libmesh_assert(this->comm().verify(p(0)));
2031 #if LIBMESH_DIM > 1
2032  libmesh_assert(this->comm().verify(p(1)));
2033 #endif
2034 #if LIBMESH_DIM > 2
2035  libmesh_assert(this->comm().verify(p(2)));
2036 #endif
2037 #endif // NDEBUG
2038 
2039  // Get a reference to the mesh object associated with the system object that calls this function
2040  const MeshBase & mesh = this->get_mesh();
2041 
2042  // Use an existing PointLocator or create a new one
2043  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2044  PointLocatorBase & locator = *locator_ptr;
2045 
2046  if (!insist_on_success || !mesh.is_serial())
2047  locator.enable_out_of_mesh_mode();
2048 
2049  // Get a pointer to the element that contains P
2050  const Elem * e = locator(p);
2051 
2052  Number u = 0;
2053 
2054  if (e && this->get_dof_map().is_evaluable(*e, var))
2055  u = point_value(var, p, *e);
2056 
2057  // If I have an element containing p, then let's let everyone know
2058  processor_id_type lowest_owner =
2059  (e && (e->processor_id() == this->processor_id())) ?
2060  this->processor_id() : this->n_processors();
2061  this->comm().min(lowest_owner);
2062 
2063  // Everybody should get their value from a processor that was able
2064  // to compute it.
2065  // If nobody admits owning the point, we have a problem.
2066  if (lowest_owner != this->n_processors())
2067  this->comm().broadcast(u, lowest_owner);
2068  else
2069  libmesh_assert(!insist_on_success);
2070 
2071  return u;
2072 }
2073 
2074 Number System::point_value(unsigned int var, const Point & p, const Elem & e) const
2075 {
2076  // Ensuring that the given point is really in the element is an
2077  // expensive assert, but as long as debugging is turned on we might
2078  // as well try to catch a particularly nasty potential error
2080 
2081  // Get the dof map to get the proper indices for our computation
2082  const DofMap & dof_map = this->get_dof_map();
2083 
2084  // Make sure we can evaluate on this element.
2085  libmesh_assert (dof_map.is_evaluable(e, var));
2086 
2087  // Need dof_indices for phi[i][j]
2088  std::vector<dof_id_type> dof_indices;
2089 
2090  // Fill in the dof_indices for our element
2091  dof_map.dof_indices (&e, dof_indices, var);
2092 
2093  // Get the no of dofs associated with this point
2094  const unsigned int num_dofs = cast_int<unsigned int>
2095  (dof_indices.size());
2096 
2097  FEType fe_type = dof_map.variable_type(var);
2098 
2099  // Build a FE so we can calculate u(p)
2100  UniquePtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2101 
2102  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2103  // Build a vector of point co-ordinates to send to reinit
2104  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2105 
2106  // Get the shape function values
2107  const std::vector<std::vector<Real> > & phi = fe->get_phi();
2108 
2109  // Reinitialize the element and compute the shape function values at coor
2110  fe->reinit (&e, &coor);
2111 
2112  // Get ready to accumulate a value
2113  Number u = 0;
2114 
2115  for (unsigned int l=0; l<num_dofs; l++)
2116  {
2117  u += phi[l][0]*this->current_solution (dof_indices[l]);
2118  }
2119 
2120  return u;
2121 }
2122 
2123 
2124 
2125 Number System::point_value(unsigned int var, const Point & p, const Elem * e) const
2126 {
2127  libmesh_assert(e);
2128  return this->point_value(var, p, *e);
2129 }
2130 
2131 
2132 
2133 Gradient System::point_gradient(unsigned int var, const Point & p, const bool insist_on_success) const
2134 {
2135  // This function must be called on every processor; there's no
2136  // telling where in the partition p falls.
2137  parallel_object_only();
2138 
2139  // And every processor had better agree about which point we're
2140  // looking for
2141 #ifndef NDEBUG
2142  libmesh_assert(this->comm().verify(p(0)));
2143 #if LIBMESH_DIM > 1
2144  libmesh_assert(this->comm().verify(p(1)));
2145 #endif
2146 #if LIBMESH_DIM > 2
2147  libmesh_assert(this->comm().verify(p(2)));
2148 #endif
2149 #endif // NDEBUG
2150 
2151  // Get a reference to the mesh object associated with the system object that calls this function
2152  const MeshBase & mesh = this->get_mesh();
2153 
2154  // Use an existing PointLocator or create a new one
2155  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2156  PointLocatorBase & locator = *locator_ptr;
2157 
2158  if (!insist_on_success || !mesh.is_serial())
2159  locator.enable_out_of_mesh_mode();
2160 
2161  // Get a pointer to the element that contains P
2162  const Elem * e = locator(p);
2163 
2164  Gradient grad_u;
2165 
2166  if (e && this->get_dof_map().is_evaluable(*e, var))
2167  grad_u = point_gradient(var, p, *e);
2168 
2169  // If I have an element containing p, then let's let everyone know
2170  processor_id_type lowest_owner =
2171  (e && (e->processor_id() == this->processor_id())) ?
2172  this->processor_id() : this->n_processors();
2173  this->comm().min(lowest_owner);
2174 
2175  // Everybody should get their value from a processor that was able
2176  // to compute it.
2177  // If nobody admits owning the point, we may have a problem.
2178  if (lowest_owner != this->n_processors())
2179  this->comm().broadcast(grad_u, lowest_owner);
2180  else
2181  libmesh_assert(!insist_on_success);
2182 
2183  return grad_u;
2184 }
2185 
2186 
2187 Gradient System::point_gradient(unsigned int var, const Point & p, const Elem & e) const
2188 {
2189  // Ensuring that the given point is really in the element is an
2190  // expensive assert, but as long as debugging is turned on we might
2191  // as well try to catch a particularly nasty potential error
2193 
2194  // Get the dof map to get the proper indices for our computation
2195  const DofMap & dof_map = this->get_dof_map();
2196 
2197  // Make sure we can evaluate on this element.
2198  libmesh_assert (dof_map.is_evaluable(e, var));
2199 
2200  // Need dof_indices for phi[i][j]
2201  std::vector<dof_id_type> dof_indices;
2202 
2203  // Fill in the dof_indices for our element
2204  dof_map.dof_indices (&e, dof_indices, var);
2205 
2206  // Get the no of dofs associated with this point
2207  const unsigned int num_dofs = cast_int<unsigned int>
2208  (dof_indices.size());
2209 
2210  FEType fe_type = dof_map.variable_type(var);
2211 
2212  // Build a FE again so we can calculate u(p)
2213  UniquePtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2214 
2215  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2216  // Build a vector of point co-ordinates to send to reinit
2217  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2218 
2219  // Get the values of the shape function derivatives
2220  const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();
2221 
2222  // Reinitialize the element and compute the shape function values at coor
2223  fe->reinit (&e, &coor);
2224 
2225  // Get ready to accumulate a gradient
2226  Gradient grad_u;
2227 
2228  for (unsigned int l=0; l<num_dofs; l++)
2229  {
2230  grad_u.add_scaled (dphi[l][0], this->current_solution (dof_indices[l]));
2231  }
2232 
2233  return grad_u;
2234 }
2235 
2236 
2237 
2238 Gradient System::point_gradient(unsigned int var, const Point & p, const Elem * e) const
2239 {
2240  libmesh_assert(e);
2241  return this->point_gradient(var, p, *e);
2242 }
2243 
2244 
2245 
2246 // We can only accumulate a hessian with --enable-second
2247 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2248 Tensor System::point_hessian(unsigned int var, const Point & p, const bool insist_on_success) const
2249 {
2250  // This function must be called on every processor; there's no
2251  // telling where in the partition p falls.
2252  parallel_object_only();
2253 
2254  // And every processor had better agree about which point we're
2255  // looking for
2256 #ifndef NDEBUG
2257  libmesh_assert(this->comm().verify(p(0)));
2258 #if LIBMESH_DIM > 1
2259  libmesh_assert(this->comm().verify(p(1)));
2260 #endif
2261 #if LIBMESH_DIM > 2
2262  libmesh_assert(this->comm().verify(p(2)));
2263 #endif
2264 #endif // NDEBUG
2265 
2266  // Get a reference to the mesh object associated with the system object that calls this function
2267  const MeshBase & mesh = this->get_mesh();
2268 
2269  // Use an existing PointLocator or create a new one
2270  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2271  PointLocatorBase & locator = *locator_ptr;
2272 
2273  if (!insist_on_success || !mesh.is_serial())
2274  locator.enable_out_of_mesh_mode();
2275 
2276  // Get a pointer to the element that contains P
2277  const Elem * e = locator(p);
2278 
2279  Tensor hess_u;
2280 
2281  if (e && this->get_dof_map().is_evaluable(*e, var))
2282  hess_u = point_hessian(var, p, *e);
2283 
2284  // If I have an element containing p, then let's let everyone know
2285  processor_id_type lowest_owner =
2286  (e && (e->processor_id() == this->processor_id())) ?
2287  this->processor_id() : this->n_processors();
2288  this->comm().min(lowest_owner);
2289 
2290  // Everybody should get their value from a processor that was able
2291  // to compute it.
2292  // If nobody admits owning the point, we may have a problem.
2293  if (lowest_owner != this->n_processors())
2294  this->comm().broadcast(hess_u, lowest_owner);
2295  else
2296  libmesh_assert(!insist_on_success);
2297 
2298  return hess_u;
2299 }
2300 
2301 Tensor System::point_hessian(unsigned int var, const Point & p, const Elem & e) const
2302 {
2303  // Ensuring that the given point is really in the element is an
2304  // expensive assert, but as long as debugging is turned on we might
2305  // as well try to catch a particularly nasty potential error
2307 
2308  // Get the dof map to get the proper indices for our computation
2309  const DofMap & dof_map = this->get_dof_map();
2310 
2311  // Make sure we can evaluate on this element.
2312  libmesh_assert (dof_map.is_evaluable(e, var));
2313 
2314  // Need dof_indices for phi[i][j]
2315  std::vector<dof_id_type> dof_indices;
2316 
2317  // Fill in the dof_indices for our element
2318  dof_map.dof_indices (&e, dof_indices, var);
2319 
2320  // Get the no of dofs associated with this point
2321  const unsigned int num_dofs = cast_int<unsigned int>
2322  (dof_indices.size());
2323 
2324  FEType fe_type = dof_map.variable_type(var);
2325 
2326  // Build a FE again so we can calculate u(p)
2327  UniquePtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2328 
2329  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2330  // Build a vector of point co-ordinates to send to reinit
2331  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2332 
2333  // Get the values of the shape function derivatives
2334  const std::vector<std::vector<RealTensor> > & d2phi = fe->get_d2phi();
2335 
2336  // Reinitialize the element and compute the shape function values at coor
2337  fe->reinit (&e, &coor);
2338 
2339  // Get ready to accumulate a hessian
2340  Tensor hess_u;
2341 
2342  for (unsigned int l=0; l<num_dofs; l++)
2343  {
2344  hess_u.add_scaled (d2phi[l][0], this->current_solution (dof_indices[l]));
2345  }
2346 
2347  return hess_u;
2348 }
2349 
2350 
2351 
2352 Tensor System::point_hessian(unsigned int var, const Point & p, const Elem * e) const
2353 {
2354  libmesh_assert(e);
2355  return this->point_hessian(var, p, *e);
2356 }
2357 
2358 
2359 
2360 #else
2361 Tensor System::point_hessian(unsigned int, const Point &, const bool) const
2362 {
2363  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2364 
2365  // Avoid compiler warnings
2366  return Tensor();
2367 }
2368 
2369 Tensor System::point_hessian(unsigned int, const Point &, const Elem &) const
2370 {
2371  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2372 
2373  // Avoid compiler warnings
2374  return Tensor();
2375 }
2376 
2377 Tensor System::point_hessian(unsigned int, const Point &, const Elem *) const
2378 {
2379  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2380 
2381  // Avoid compiler warnings
2382  return Tensor();
2383 }
2384 
2385 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2386 
2387 } // namespace libMesh
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Definition: system.C:903
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2021
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:178
vectors_iterator vectors_end()
Definition: system.h:2236
FEFamily family
Definition: fe_type.h:203
dof_id_type n_constrained_dofs() const
Definition: system.C:153
bool closed()
Definition: libmesh.C:281
std::map< std::string, unsigned short int > _variable_numbers
Definition: system.h:1901
double abs(double a)
const FEType & type() const
Definition: variable.h:119
std::string get_info() const
Definition: dof_map.C:2764
Manages multiples systems of equations.
virtual void clear()
Definition: system.C:201
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
virtual bool is_serial() const
Definition: mesh_base.h:137
Assembly * _assemble_system_object
Definition: system.h:1820
virtual numeric_index_type size() const =0
bool _basic_system_only
Definition: system.h:1944
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2132
std::map< std::string, ParallelType > _vector_types
Definition: system.h:1931
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
bool _is_initialized
Definition: system.h:1950
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1825
NumericVector< Number > & add_adjoint_solution(unsigned int i=0)
Definition: system.C:975
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
virtual void forward_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Definition: system.h:2298
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1686
Specifies parameters for parameter sensitivity calculations.
Real norm_sq() const
Definition: type_tensor.h:1262
void add_scaled(const TypeTensor< T2 > &, const T)
Definition: type_tensor.h:773
void add_scaled(const TypeVector< T2 > &, const T)
Definition: type_vector.h:620
int vector_is_adjoint(const std::string &vec_name) const
Definition: system.C:914
virtual Real l2_norm() const =0
Used to specify quantities of interest in a simulation.
Definition: qoi_set.h:45
virtual void reinit()
Definition: system.C:383
NumericVector< Number > & get_sensitivity_solution(unsigned int i=0)
Definition: system.C:934
OrderWrapper radial_order
Definition: fe_type.h:236
Real weight(unsigned int var) const
Definition: system_norm.h:292
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1099
bool has_variable(const std::string &var) const
Definition: system.C:1254
const unsigned int _sys_number
Definition: system.h:1885
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
Definition: system.C:499
NumericVector< Number > & get_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1079
virtual void assemble()=0
Constraint * _constrain_system_object
Definition: system.h:1831
virtual void assemble()
Definition: system.C:477
virtual void init_data()
Definition: system.C:258
processor_id_type n_processors() const
void zero_variable(NumericVector< Number > &v, unsigned int var_num) const
Definition: system.C:1331
The base class for all geometric element types.
Definition: elem.h:86
unsigned int size(const System &sys) const
Definition: qoi_set.C:35
Real norm() const
Definition: type_vector.h:903
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
UniquePtr< NumericVector< Number > > current_local_solution
Definition: system.h:1533
const class libmesh_nullptr_t libmesh_nullptr
bool vector_preservation(const std::string &vec_name) const
Definition: system.C:893
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:802
std::map< std::string, NumericVector< Number > * >::const_iterator const_vectors_iterator
Definition: system.h:749
virtual void user_initialization()
Definition: system.C:1946
std::vector< Variable > _variables
Definition: system.h:1890
OrderWrapper order
Definition: fe_type.h:197
NumericVector< Number > & add_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1069
void attach_QOI_derivative(void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints))
Definition: system.C:1911
vectors_iterator vectors_begin()
Definition: system.h:2224
void attach_QOI_derivative_object(QOIDerivative &qoi_derivative)
Definition: system.C:1930
virtual void adjoint_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Definition: system.h:2289
const VariableGroup & variable_group(unsigned int vg) const
Definition: system.h:2122
NumericVector< Number > & add_weighted_sensitivity_solution()
Definition: system.C:954
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2162
void(* _qoi_evaluate_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices)
Definition: system.h:1836
virtual Real l1_norm() const =0
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:1646
void init()
Definition: system.C:233
long double max(long double a, double b)
UniquePtr< DofMap > _dof_map
Definition: system.h:1863
const std::string & name() const
Definition: system.h:1996
NumericVector< Number > & get_weighted_sensitivity_solution()
Definition: system.C:961
std::string get_info() const
Definition: system.C:1686
Base class for Mesh.
Definition: mesh_base.h:67
virtual void prolong_vectors()
Definition: system.C:373
virtual void user_QOI(const QoISet &qoi_indices)
Definition: system.C:1988
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void attach_QOI_function(void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
Definition: system.C:1875
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
QOIDerivative * _qoi_evaluate_derivative_object
Definition: system.h:1857
void remove_vector(const std::string &vec_name)
Definition: system.C:717
NumericVector< Number > & add_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1007
virtual void qoi_derivative(const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)=0
virtual void qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Definition: system.C:513
void attach_constraint_object(Constraint &constrain)
Definition: system.C:1859
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:167
const NumericVector< Number > * request_vector(const std::string &vec_name) const
Definition: system.C:734
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1261
void attach_init_object(Initialization &init)
Definition: system.C:1789
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:575
void update_global_solution(std::vector< Number > &global_soln) const
Definition: system.C:660
Real norm() const
Definition: type_tensor.h:1193
unsigned int n_variable_groups() const
Definition: system.h:2092
const MeshBase & get_mesh() const
Definition: system.h:2012
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:679
const DofMap & get_dof_map() const
Definition: system.h:2028
unsigned int n_components() const
Definition: system.h:2100
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:231
Holds completed parameter sensitivity calculations.
bool is_discrete() const
Definition: system_norm.h:254
bool have_vector(const std::string &vec_name) const
Definition: system.h:2204
virtual node_iterator local_nodes_end()=0
virtual element_iterator active_local_elements_begin()=0
std::vector< VariableGroup > _variable_groups
Definition: system.h:1895
bool is_initialized()
Definition: system.h:2068
unsigned int n_variables() const
Definition: variable.h:217
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Definition: system.C:1805
void attach_constraint_function(void fptr(EquationSystems &es, const std::string &name))
Definition: system.C:1840
virtual Real subset_l1_norm(const std::set< numeric_index_type > &indices) const
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
virtual void qoi(const QoISet &qoi_indices)=0
QOI * _qoi_evaluate_object
Definition: system.h:1843
Initialization * _init_system_object
Definition: system.h:1809
System & operator=(const System &)
Definition: system.C:110
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:556
InfMapType inf_map
Definition: fe_type.h:257
Tensor point_hessian(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2248
UniquePtr< NumericVector< Number > > solution
Definition: system.h:1521
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1803
virtual void restrict_vectors()
Definition: system.C:320
void attach_init_function(void fptr(EquationSystems &es, const std::string &name))
Definition: system.C:1770
void broadcast(T &data, const unsigned int root_id=0) const
Real norm_sq() const
Definition: type_vector.h:932
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:203
void append(const std::string &var_name)
Definition: variable.h:276
std::map< std::string, int > _vector_is_adjoint
Definition: system.h:1926
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:792
NumericVector< Number > & add_adjoint_rhs(unsigned int i=0)
Definition: system.C:1039
virtual Real linfty_norm() const =0
virtual void restrict_solve_to(const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)
Definition: system.C:468
virtual void user_QOI_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
Definition: system.C:2002
virtual void reinit_constraints()
Definition: system.C:412
virtual unsigned int n_matrices() const
Definition: system.h:2218
Gradient point_gradient(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2133
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=libmesh_nullptr) const
Definition: system.C:1397
void attach_assemble_object(Assembly &assemble)
Definition: system.C:1824
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Definition: system.h:748
FEFamily radial_family
Definition: fe_type.h:249
std::size_t size() const
bool _solution_projection
Definition: system.h:1938
EquationSystems & _equation_systems
Definition: system.h:1869
void local_dof_indices(const unsigned int var, std::set< dof_id_type > &var_indices) const
Definition: system.C:1295
virtual void constrain()=0
dof_id_type n_local_dofs() const
Definition: system.C:183
virtual bool compare(const System &other_system, const Real threshold, const bool verbose) const
Definition: system.C:529
Real weight_sq(unsigned int var) const
Definition: system_norm.h:338
void attach_QOI_object(QOI &qoi)
Definition: system.C:1895
virtual std::string system_type() const
Definition: system.h:471
virtual void user_assembly()
Definition: system.C:1960
dof_id_type n_local_constrained_dofs() const
Definition: system.C:168
virtual node_iterator local_nodes_begin()=0
virtual void update()
Definition: system.C:423
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_vector_preservation(const std::string &vec_name, bool preserve)
Definition: system.C:885
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:190
void process_constraints(MeshBase &)
bool verify(const T &r, const Communicator &comm=Communicator_World)
bool identify_variable_groups() const
Definition: system.h:2180
NumberTensorValue Tensor
unsigned int number() const
Definition: system.h:2004
std::map< std::string, NumericVector< Number > * > _vectors
Definition: system.h:1914
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:290
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
const std::string _sys_name
Definition: system.h:1880
unsigned int add_variables(const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1189
virtual unsigned int dim() const =0
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:987
virtual void re_update()
Definition: system.C:444
std::map< std::string, bool > _vector_projections
Definition: system.h:1920
dof_id_type n_dofs() const
Definition: system.C:146
void create_dof_constraints(const MeshBase &, Real time=0)
unsigned int n_vars() const
Definition: system.h:2084
FEMNormType type(unsigned int var) const
Definition: system_norm.h:268
virtual void set(const numeric_index_type i, const T value)=0
void prepare_send_list()
Definition: dof_map.C:1723
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
virtual void user_constrain()
Definition: system.C:1974
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet())
Definition: system.C:488
NumericVector< Number > & add_sensitivity_solution(unsigned int i=0)
Definition: system.C:924
const std::string & name(unsigned int v) const
Definition: variable.h:246
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1814
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:772
virtual ~System()
Definition: system.C:116
virtual void enable_out_of_mesh_mode()=0
System(EquationSystems &es, const std::string &name, const unsigned int number)
Definition: system.C:57
OStreamProxy out(std::cout)
void get_all_variable_numbers(std::vector< unsigned int > &all_variable_numbers) const
Definition: system.C:1276
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:394
virtual element_iterator active_local_elements_end()=0
A geometric point in (x,y,z) space.
Definition: point.h:38
MeshBase & _mesh
Definition: system.h:1875
virtual void localize(std::vector< T > &v_local) const =0
unsigned int n_points() const
Definition: quadrature.h:113
UniquePtr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:531
Base class for all quadrature families and orders.
Definition: quadrature.h:53
Real discrete_var_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type) const
Definition: system.C:1378
const std::string & vector_name(const unsigned int vec_num) const
Definition: system.C:850
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition: elem.C:2469
NumericVector< Number > & get_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1019
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:686
unsigned int n_vectors() const
Definition: system.h:2212
void(* _qoi_evaluate_derivative_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)
Definition: system.h:1848
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2005
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1049