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