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