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