libMesh
system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 #include <sstream> // for std::ostringstream
22 
23 
24 // Local includes
25 #include "libmesh/dof_map.h"
26 #include "libmesh/equation_systems.h"
27 #include "libmesh/libmesh_logging.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/parameter_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"
41 #include "libmesh/fe_compute_data.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"
47 #include "libmesh/parallel_algebra.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 
359 #ifdef LIBMESH_ENABLE_GHOSTED
360  current_local_solution->init(this->n_dofs(),
361  this->n_local_dofs(), send_list,
362  false, GHOSTED);
363 #else
364  current_local_solution->init(this->n_dofs());
365 #endif
366 
368  solution->localize (*current_local_solution, send_list);
369 
370 #endif // LIBMESH_ENABLE_AMR
371 }
372 
373 
374 
376 {
377 #ifdef LIBMESH_ENABLE_AMR
378  // Currently project_vector handles both restriction and prolongation
379  this->restrict_vectors();
380 #endif
381 }
382 
383 
384 
386 {
387  //If no variables have been added to this system
388  //don't do anything
389  if (!this->n_vars())
390  return;
391 
392  // Constraints get handled in EquationSystems::reinit now
393  // _dof_map->create_dof_constraints(this->get_mesh());
394 
395  // Update the solution based on the projected
396  // current_local_solution.
397  solution->init (this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
398 
399  libmesh_assert_equal_to (solution->size(), current_local_solution->size());
400  // Not true with ghosted vectors
401  // libmesh_assert_equal_to (solution->size(), current_local_solution->local_size());
402 
403  const dof_id_type first_local_dof = solution->first_local_index();
404  const dof_id_type local_size = solution->local_size();
405 
406  for (dof_id_type i=0; i<local_size; i++)
407  solution->set(i+first_local_dof,
408  (*current_local_solution)(i+first_local_dof));
409 
410  solution->close();
411 }
412 
413 
415 {
416 #ifdef LIBMESH_ENABLE_CONSTRAINTS
418  user_constrain();
420 #endif
422 }
423 
424 
426 {
427  libmesh_assert(solution->closed());
428 
429  const std::vector<dof_id_type> & send_list = _dof_map->get_send_list ();
430 
431  // Check sizes
432  libmesh_assert_equal_to (current_local_solution->size(), solution->size());
433  // More processors than elements => empty send_list
434  // libmesh_assert (!send_list.empty());
435  libmesh_assert_less_equal (send_list.size(), solution->size());
436 
437  // Create current_local_solution from solution. This will
438  // put a local copy of solution into current_local_solution.
439  // Only the necessary values (specified by the send_list)
440  // are copied to minimize communication
441  solution->localize (*current_local_solution, send_list);
442 }
443 
444 
445 
447 {
448  parallel_object_only();
449 
450  // If this system is empty... don't do anything!
451  if (!this->n_vars())
452  return;
453 
454  const std::vector<dof_id_type> & send_list = this->get_dof_map().get_send_list ();
455 
456  // Check sizes
457  libmesh_assert_equal_to (current_local_solution->size(), solution->size());
458  // Not true with ghosted vectors
459  // libmesh_assert_equal_to (current_local_solution->local_size(), solution->size());
460  // libmesh_assert (!send_list.empty());
461  libmesh_assert_less_equal (send_list.size(), solution->size());
462 
463  // Create current_local_solution from solution. This will
464  // put a local copy of solution into current_local_solution.
465  solution->localize (*current_local_solution, send_list);
466 }
467 
468 
469 
471  const SubsetSolveMode /*subset_solve_mode*/)
472 {
473  if (subset != libmesh_nullptr)
474  libmesh_not_implemented();
475 }
476 
477 
478 
480 {
481  // Log how long the user's assembly code takes
482  LOG_SCOPE("assemble()", "System");
483 
484  // Call the user-specified assembly function
485  this->user_assembly();
486 }
487 
488 
489 
490 void System::assemble_qoi (const QoISet & qoi_indices)
491 {
492  // Log how long the user's assembly code takes
493  LOG_SCOPE("assemble_qoi()", "System");
494 
495  // Call the user-specified quantity of interest function
496  this->user_QOI(qoi_indices);
497 }
498 
499 
500 
501 void System::assemble_qoi_derivative(const QoISet & qoi_indices,
502  bool include_liftfunc,
503  bool apply_constraints)
504 {
505  // Log how long the user's assembly code takes
506  LOG_SCOPE("assemble_qoi_derivative()", "System");
507 
508  // Call the user-specified quantity of interest function
509  this->user_QOI_derivative(qoi_indices, include_liftfunc,
510  apply_constraints);
511 }
512 
513 
514 
515 void System::qoi_parameter_sensitivity (const QoISet & qoi_indices,
516  const ParameterVector & parameters,
517  SensitivityData & sensitivities)
518 {
519  // Forward sensitivities are more efficient for Nq > Np
520  if (qoi_indices.size(*this) > parameters.size())
521  forward_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
522  // Adjoint sensitivities are more efficient for Np > Nq,
523  // and an adjoint may be more reusable than a forward
524  // solution sensitivity in the Np == Nq case.
525  else
526  adjoint_qoi_parameter_sensitivity(qoi_indices, parameters, sensitivities);
527 }
528 
529 
530 
531 bool System::compare (const System & other_system,
532  const Real threshold,
533  const bool verbose) const
534 {
535  // we do not care for matrices, but for vectors
537  libmesh_assert (other_system._is_initialized);
538 
539  if (verbose)
540  {
541  libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
542  libMesh::out << " comparing matrices not supported." << std::endl;
543  libMesh::out << " comparing names...";
544  }
545 
546  // compare the name: 0 means identical
547  const int name_result = _sys_name.compare(other_system.name());
548  if (verbose)
549  {
550  if (name_result == 0)
551  libMesh::out << " identical." << std::endl;
552  else
553  libMesh::out << " names not identical." << std::endl;
554  libMesh::out << " comparing solution vector...";
555  }
556 
557 
558  // compare the solution: -1 means identical
559  const int solu_result = solution->compare (*other_system.solution.get(),
560  threshold);
561 
562  if (verbose)
563  {
564  if (solu_result == -1)
565  libMesh::out << " identical up to threshold." << std::endl;
566  else
567  libMesh::out << " first difference occurred at index = "
568  << solu_result << "." << std::endl;
569  }
570 
571 
572  // safety check, whether we handle at least the same number
573  // of vectors
574  std::vector<int> ov_result;
575 
576  if (this->n_vectors() != other_system.n_vectors())
577  {
578  if (verbose)
579  {
580  libMesh::out << " Fatal difference. This system handles "
581  << this->n_vectors() << " add'l vectors," << std::endl
582  << " while the other system handles "
583  << other_system.n_vectors()
584  << " add'l vectors." << std::endl
585  << " Aborting comparison." << std::endl;
586  }
587  return false;
588  }
589  else if (this->n_vectors() == 0)
590  {
591  // there are no additional vectors...
592  ov_result.clear ();
593  }
594  else
595  {
596  // compare other vectors
597  for (const_vectors_iterator pos = _vectors.begin();
598  pos != _vectors.end(); ++pos)
599  {
600  if (verbose)
601  libMesh::out << " comparing vector \""
602  << pos->first << "\" ...";
603 
604  // assume they have the same name
605  const NumericVector<Number> & other_system_vector =
606  other_system.get_vector(pos->first);
607 
608  ov_result.push_back(pos->second->compare (other_system_vector,
609  threshold));
610 
611  if (verbose)
612  {
613  if (ov_result[ov_result.size()-1] == -1)
614  libMesh::out << " identical up to threshold." << std::endl;
615  else
616  libMesh::out << " first difference occurred at" << std::endl
617  << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
618  }
619 
620  }
621 
622  } // finished comparing additional vectors
623 
624 
625  bool overall_result;
626 
627  // sum up the results
628  if ((name_result==0) && (solu_result==-1))
629  {
630  if (ov_result.size()==0)
631  overall_result = true;
632  else
633  {
634  bool ov_identical;
635  unsigned int n = 0;
636  do
637  {
638  ov_identical = (ov_result[n]==-1);
639  n++;
640  }
641  while (ov_identical && n<ov_result.size());
642  overall_result = ov_identical;
643  }
644  }
645  else
646  overall_result = false;
647 
648  if (verbose)
649  {
650  libMesh::out << " finished comparisons, ";
651  if (overall_result)
652  libMesh::out << "found no differences." << std::endl << std::endl;
653  else
654  libMesh::out << "found differences." << std::endl << std::endl;
655  }
656 
657  return overall_result;
658 }
659 
660 
661 
662 void System::update_global_solution (std::vector<Number> & global_soln) const
663 {
664  global_soln.resize (solution->size());
665 
666  solution->localize (global_soln);
667 }
668 
669 
670 
671 void System::update_global_solution (std::vector<Number> & global_soln,
672  const processor_id_type dest_proc) const
673 {
674  global_soln.resize (solution->size());
675 
676  solution->localize_to_one (global_soln, dest_proc);
677 }
678 
679 
680 
681 NumericVector<Number> & System::add_vector (const std::string & vec_name,
682  const bool projections,
683  const ParallelType type)
684 {
685  // Return the vector if it is already there.
686  if (this->have_vector(vec_name))
687  return *(_vectors[vec_name]);
688 
689  // Otherwise build the vector
690  NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
691  _vectors.insert (std::make_pair (vec_name, buf));
692  _vector_projections.insert (std::make_pair (vec_name, projections));
693 
694  _vector_types.insert (std::make_pair (vec_name, type));
695 
696  // Vectors are primal by default
697  _vector_is_adjoint.insert (std::make_pair (vec_name, -1));
698 
699  // Initialize it if necessary
700  if (_is_initialized)
701  {
702  if (type == GHOSTED)
703  {
704 #ifdef LIBMESH_ENABLE_GHOSTED
705  buf->init (this->n_dofs(), this->n_local_dofs(),
706  _dof_map->get_send_list(), false,
707  GHOSTED);
708 #else
709  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
710 #endif
711  }
712  else
713  buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
714  }
715 
716  return *buf;
717 }
718 
719 void System::remove_vector (const std::string & vec_name)
720 {
721  vectors_iterator pos = _vectors.find(vec_name);
722 
723  //Return if the vector does not exist
724  if (pos == _vectors.end())
725  return;
726 
727  delete pos->second;
728 
729  _vectors.erase(pos);
730 
731  _vector_projections.erase(vec_name);
732  _vector_is_adjoint.erase(vec_name);
733  _vector_types.erase(vec_name);
734 }
735 
736 const NumericVector<Number> * System::request_vector (const std::string & vec_name) const
737 {
738  const_vectors_iterator pos = _vectors.find(vec_name);
739 
740  if (pos == _vectors.end())
741  return libmesh_nullptr;
742 
743  return pos->second;
744 }
745 
746 
747 
748 NumericVector<Number> * System::request_vector (const std::string & vec_name)
749 {
750  vectors_iterator pos = _vectors.find(vec_name);
751 
752  if (pos == _vectors.end())
753  return libmesh_nullptr;
754 
755  return pos->second;
756 }
757 
758 
759 
760 const NumericVector<Number> * System::request_vector (const unsigned int vec_num) const
761 {
764  unsigned int num = 0;
765  while ((num<vec_num) && (v!=v_end))
766  {
767  num++;
768  ++v;
769  }
770  if (v==v_end)
771  return libmesh_nullptr;
772  return v->second;
773 }
774 
775 
776 
777 NumericVector<Number> * System::request_vector (const unsigned int vec_num)
778 {
780  vectors_iterator v_end = vectors_end();
781  unsigned int num = 0;
782  while ((num<vec_num) && (v!=v_end))
783  {
784  num++;
785  ++v;
786  }
787  if (v==v_end)
788  return libmesh_nullptr;
789  return v->second;
790 }
791 
792 
793 
794 const NumericVector<Number> & System::get_vector (const std::string & vec_name) const
795 {
796  // Make sure the vector exists
797  const_vectors_iterator pos = _vectors.find(vec_name);
798 
799  if (pos == _vectors.end())
800  libmesh_error_msg("ERROR: vector " << vec_name << " does not exist in this system!");
801 
802  return *(pos->second);
803 }
804 
805 
806 
807 NumericVector<Number> & System::get_vector (const std::string & vec_name)
808 {
809  // Make sure the vector exists
810  vectors_iterator pos = _vectors.find(vec_name);
811 
812  if (pos == _vectors.end())
813  libmesh_error_msg("ERROR: vector " << vec_name << " does not exist in this system!");
814 
815  return *(pos->second);
816 }
817 
818 
819 
820 const NumericVector<Number> & System::get_vector (const unsigned int vec_num) const
821 {
824  unsigned int num = 0;
825  while ((num<vec_num) && (v!=v_end))
826  {
827  num++;
828  ++v;
829  }
830  libmesh_assert (v != v_end);
831  return *(v->second);
832 }
833 
834 
835 
836 NumericVector<Number> & System::get_vector (const unsigned int vec_num)
837 {
839  vectors_iterator v_end = vectors_end();
840  unsigned int num = 0;
841  while ((num<vec_num) && (v!=v_end))
842  {
843  num++;
844  ++v;
845  }
846  libmesh_assert (v != v_end);
847  return *(v->second);
848 }
849 
850 
851 
852 const std::string & System::vector_name (const unsigned int vec_num) const
853 {
856  unsigned int num = 0;
857  while ((num<vec_num) && (v!=v_end))
858  {
859  num++;
860  ++v;
861  }
862  libmesh_assert (v != v_end);
863  return v->first;
864 }
865 
866 const std::string & System::vector_name (const NumericVector<Number> & vec_reference) const
867 {
870 
871  for (; v != v_end; ++v)
872  {
873  // Check if the current vector is the one whose name we want
874  if (&vec_reference == v->second)
875  break; // exit loop if it is
876  }
877 
878  // Before returning, make sure we didnt loop till the end and not find any match
879  libmesh_assert (v != v_end);
880 
881  // Return the string associated with the current vector
882  return v->first;
883 }
884 
885 
886 
887 void System::set_vector_preservation (const std::string & vec_name,
888  bool preserve)
889 {
890  _vector_projections[vec_name] = preserve;
891 }
892 
893 
894 
895 bool System::vector_preservation (const std::string & vec_name) const
896 {
897  if (_vector_projections.find(vec_name) == _vector_projections.end())
898  return false;
899 
900  return _vector_projections.find(vec_name)->second;
901 }
902 
903 
904 
905 void System::set_vector_as_adjoint (const std::string & vec_name,
906  int qoi_num)
907 {
908  // We reserve -1 for vectors which get primal constraints, -2 for
909  // vectors which get no constraints
910  libmesh_assert_greater_equal(qoi_num, -2);
911  _vector_is_adjoint[vec_name] = qoi_num;
912 }
913 
914 
915 
916 int System::vector_is_adjoint (const std::string & vec_name) const
917 {
918  libmesh_assert(_vector_is_adjoint.find(vec_name) !=
919  _vector_is_adjoint.end());
920 
921  return _vector_is_adjoint.find(vec_name)->second;
922 }
923 
924 
925 
927 {
928  std::ostringstream sensitivity_name;
929  sensitivity_name << "sensitivity_solution" << i;
930 
931  return this->add_vector(sensitivity_name.str());
932 }
933 
934 
935 
937 {
938  std::ostringstream sensitivity_name;
939  sensitivity_name << "sensitivity_solution" << i;
940 
941  return this->get_vector(sensitivity_name.str());
942 }
943 
944 
945 
947 {
948  std::ostringstream sensitivity_name;
949  sensitivity_name << "sensitivity_solution" << i;
950 
951  return this->get_vector(sensitivity_name.str());
952 }
953 
954 
955 
957 {
958  return this->add_vector("weighted_sensitivity_solution");
959 }
960 
961 
962 
964 {
965  return this->get_vector("weighted_sensitivity_solution");
966 }
967 
968 
969 
971 {
972  return this->get_vector("weighted_sensitivity_solution");
973 }
974 
975 
976 
978 {
979  std::ostringstream adjoint_name;
980  adjoint_name << "adjoint_solution" << i;
981 
982  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
983  this->set_vector_as_adjoint(adjoint_name.str(), i);
984  return returnval;
985 }
986 
987 
988 
990 {
991  std::ostringstream adjoint_name;
992  adjoint_name << "adjoint_solution" << i;
993 
994  return this->get_vector(adjoint_name.str());
995 }
996 
997 
998 
1000 {
1001  std::ostringstream adjoint_name;
1002  adjoint_name << "adjoint_solution" << i;
1003 
1004  return this->get_vector(adjoint_name.str());
1005 }
1006 
1007 
1008 
1010 {
1011  std::ostringstream adjoint_name;
1012  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1013 
1014  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1015  this->set_vector_as_adjoint(adjoint_name.str(), i);
1016  return returnval;
1017 }
1018 
1019 
1020 
1022 {
1023  std::ostringstream adjoint_name;
1024  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1025 
1026  return this->get_vector(adjoint_name.str());
1027 }
1028 
1029 
1030 
1032 {
1033  std::ostringstream adjoint_name;
1034  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1035 
1036  return this->get_vector(adjoint_name.str());
1037 }
1038 
1039 
1040 
1042 {
1043  std::ostringstream adjoint_rhs_name;
1044  adjoint_rhs_name << "adjoint_rhs" << i;
1045 
1046  return this->add_vector(adjoint_rhs_name.str(), false);
1047 }
1048 
1049 
1050 
1052 {
1053  std::ostringstream adjoint_rhs_name;
1054  adjoint_rhs_name << "adjoint_rhs" << i;
1055 
1056  return this->get_vector(adjoint_rhs_name.str());
1057 }
1058 
1059 
1060 
1061 const NumericVector<Number> & System::get_adjoint_rhs (unsigned int i) const
1062 {
1063  std::ostringstream adjoint_rhs_name;
1064  adjoint_rhs_name << "adjoint_rhs" << i;
1065 
1066  return this->get_vector(adjoint_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->add_vector(sensitivity_rhs_name.str(), false);
1077 }
1078 
1079 
1080 
1082 {
1083  std::ostringstream sensitivity_rhs_name;
1084  sensitivity_rhs_name << "sensitivity_rhs" << i;
1085 
1086  return this->get_vector(sensitivity_rhs_name.str());
1087 }
1088 
1089 
1090 
1092 {
1093  std::ostringstream sensitivity_rhs_name;
1094  sensitivity_rhs_name << "sensitivity_rhs" << i;
1095 
1096  return this->get_vector(sensitivity_rhs_name.str());
1097 }
1098 
1099 
1100 
1101 unsigned int System::add_variable (const std::string & var,
1102  const FEType & type,
1103  const std::set<subdomain_id_type> * const active_subdomains)
1104 {
1105  libmesh_assert(!this->is_initialized());
1106 
1107  // Make sure the variable isn't there already
1108  // or if it is, that it's the type we want
1109  for (unsigned int v=0; v<this->n_vars(); v++)
1110  if (this->variable_name(v) == var)
1111  {
1112  if (this->variable_type(v) == type)
1113  return _variables[v].number();
1114 
1115  libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1116  }
1117 
1118  // Optimize for VariableGroups here - if the user is adding multiple
1119  // variables of the same FEType and subdomain restriction, catch
1120  // that here and add them as members of the same VariableGroup.
1121  //
1122  // start by setting this flag to whatever the user has requested
1123  // and then consider the conditions which should negate it.
1124  bool should_be_in_vg = this->identify_variable_groups();
1125 
1126  // No variable groups, nothing to add to
1127  if (!this->n_variable_groups())
1128  should_be_in_vg = false;
1129 
1130  else
1131  {
1132  VariableGroup & vg(_variable_groups.back());
1133 
1134  // get a pointer to their subdomain restriction, if any.
1135  const std::set<subdomain_id_type> * const
1136  their_active_subdomains (vg.implicitly_active() ?
1137  libmesh_nullptr : &vg.active_subdomains());
1138 
1139  // Different types?
1140  if (vg.type() != type)
1141  should_be_in_vg = false;
1142 
1143  // they are restricted, we aren't?
1144  if (their_active_subdomains && !active_subdomains)
1145  should_be_in_vg = false;
1146 
1147  // they aren't restricted, we are?
1148  if (!their_active_subdomains && active_subdomains)
1149  should_be_in_vg = false;
1150 
1151  if (their_active_subdomains && active_subdomains)
1152  // restricted to different sets?
1153  if (*their_active_subdomains != *active_subdomains)
1154  should_be_in_vg = false;
1155 
1156  // OK, after all that, append the variable to the vg if none of the conditions
1157  // were violated
1158  if (should_be_in_vg)
1159  {
1160  const unsigned short curr_n_vars = cast_int<unsigned short>
1161  (this->n_vars());
1162 
1163  vg.append (var);
1164 
1165  _variables.push_back(vg(vg.n_variables()-1));
1166  _variable_numbers[var] = curr_n_vars;
1167  return curr_n_vars;
1168  }
1169  }
1170 
1171  // otherwise, fall back to adding a single variable group
1172  return this->add_variables (std::vector<std::string>(1, var),
1173  type,
1174  active_subdomains);
1175 }
1176 
1177 
1178 
1179 unsigned int System::add_variable (const std::string & var,
1180  const Order order,
1181  const FEFamily family,
1182  const std::set<subdomain_id_type> * const active_subdomains)
1183 {
1184  return this->add_variable(var,
1185  FEType(order, family),
1186  active_subdomains);
1187 }
1188 
1189 
1190 
1191 unsigned int System::add_variables (const std::vector<std::string> & vars,
1192  const FEType & type,
1193  const std::set<subdomain_id_type> * const active_subdomains)
1194 {
1195  libmesh_assert(!this->is_initialized());
1196 
1197  // Make sure the variable isn't there already
1198  // or if it is, that it's the type we want
1199  for (std::size_t ov=0; ov<vars.size(); ov++)
1200  for (unsigned int v=0; v<this->n_vars(); v++)
1201  if (this->variable_name(v) == vars[ov])
1202  {
1203  if (this->variable_type(v) == type)
1204  return _variables[v].number();
1205 
1206  libmesh_error_msg("ERROR: incompatible variable " << vars[ov] << " has already been added for this system!");
1207  }
1208 
1209  const unsigned short curr_n_vars = cast_int<unsigned short>
1210  (this->n_vars());
1211 
1212  const unsigned int next_first_component = this->n_components();
1213 
1214  // Add the variable group to the list
1215  _variable_groups.push_back((active_subdomains == libmesh_nullptr) ?
1216  VariableGroup(this, vars, curr_n_vars,
1217  next_first_component, type) :
1218  VariableGroup(this, vars, curr_n_vars,
1219  next_first_component, type, *active_subdomains));
1220 
1221  const VariableGroup & vg (_variable_groups.back());
1222 
1223  // Add each component of the group individually
1224  for (std::size_t v=0; v<vars.size(); v++)
1225  {
1226  _variables.push_back (vg(v));
1227  _variable_numbers[vars[v]] = cast_int<unsigned short>
1228  (curr_n_vars+v);
1229  }
1230 
1231  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1232 
1233  // BSK - Defer this now to System::init_data() so we can detect
1234  // VariableGroups 12/28/2012
1235  // // Add the variable group to the _dof_map
1236  // _dof_map->add_variable_group (vg);
1237 
1238  // Return the number of the new variable
1239  return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1240 }
1241 
1242 
1243 
1244 unsigned int System::add_variables (const std::vector<std::string> & vars,
1245  const Order order,
1246  const FEFamily family,
1247  const std::set<subdomain_id_type> * const active_subdomains)
1248 {
1249  return this->add_variables(vars,
1250  FEType(order, family),
1251  active_subdomains);
1252 }
1253 
1254 
1255 
1256 bool System::has_variable (const std::string & var) const
1257 {
1258  return _variable_numbers.count(var);
1259 }
1260 
1261 
1262 
1263 unsigned short int System::variable_number (const std::string & var) const
1264 {
1265  // Make sure the variable exists
1266  std::map<std::string, unsigned short int>::const_iterator
1267  pos = _variable_numbers.find(var);
1268 
1269  if (pos == _variable_numbers.end())
1270  libmesh_error_msg("ERROR: variable " << var << " does not exist in this system!");
1271 
1272  libmesh_assert_equal_to (_variables[pos->second].name(), var);
1273 
1274  return pos->second;
1275 }
1276 
1277 
1278 void System::get_all_variable_numbers(std::vector<unsigned int> & all_variable_numbers) const
1279 {
1280  all_variable_numbers.resize(n_vars());
1281 
1282  // Make sure the variable exists
1283  std::map<std::string, unsigned short int>::const_iterator
1284  it = _variable_numbers.begin();
1285  std::map<std::string, unsigned short int>::const_iterator
1286  it_end = _variable_numbers.end();
1287 
1288  unsigned int count = 0;
1289  for ( ; it != it_end; ++it)
1290  {
1291  all_variable_numbers[count] = it->second;
1292  count++;
1293  }
1294 }
1295 
1296 
1297 void System::local_dof_indices(const unsigned int var,
1298  std::set<dof_id_type> & var_indices) const
1299 {
1300  // Make sure the set is clear
1301  var_indices.clear();
1302 
1303  std::vector<dof_id_type> dof_indices;
1304 
1305  const dof_id_type
1306  first_local = this->get_dof_map().first_dof(),
1307  end_local = this->get_dof_map().end_dof();
1308 
1309  // Begin the loop over the elements
1310  for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1311  {
1312  this->get_dof_map().dof_indices (elem, dof_indices, var);
1313 
1314  for (std::size_t i=0; i<dof_indices.size(); i++)
1315  {
1316  dof_id_type dof = dof_indices[i];
1317 
1318  //If the dof is owned by the local processor
1319  if (first_local <= dof && dof < end_local)
1320  var_indices.insert(dof_indices[i]);
1321  }
1322  }
1323 }
1324 
1325 
1326 
1328  unsigned int var_num) const
1329 {
1330  /* Make sure the call makes sense. */
1331  libmesh_assert_less (var_num, this->n_vars());
1332 
1333  /* Get a reference to the mesh. */
1334  const MeshBase & mesh = this->get_mesh();
1335 
1336  /* Check which system we are. */
1337  const unsigned int sys_num = this->number();
1338 
1339  // Loop over nodes.
1340  for (const auto & node : mesh.local_node_ptr_range())
1341  {
1342  unsigned int n_comp = node->n_comp(sys_num,var_num);
1343  for (unsigned int i=0; i<n_comp; i++)
1344  {
1345  const dof_id_type index = node->dof_number(sys_num,var_num,i);
1346  v.set(index,0.0);
1347  }
1348  }
1349 
1350  // Loop over elements.
1351  for (const auto & elem : mesh.active_local_element_ptr_range())
1352  {
1353  unsigned int n_comp = elem->n_comp(sys_num,var_num);
1354  for (unsigned int i=0; i<n_comp; i++)
1355  {
1356  const dof_id_type index = elem->dof_number(sys_num,var_num,i);
1357  v.set(index,0.0);
1358  }
1359  }
1360 }
1361 
1362 
1363 
1365  unsigned int var,
1366  FEMNormType norm_type) const
1367 {
1368  std::set<dof_id_type> var_indices;
1369  local_dof_indices(var, var_indices);
1370 
1371  if (norm_type == DISCRETE_L1)
1372  return v.subset_l1_norm(var_indices);
1373  if (norm_type == DISCRETE_L2)
1374  return v.subset_l2_norm(var_indices);
1375  if (norm_type == DISCRETE_L_INF)
1376  return v.subset_linfty_norm(var_indices);
1377  else
1378  libmesh_error_msg("Invalid norm_type = " << norm_type);
1379 }
1380 
1381 
1382 
1384  unsigned int var,
1385  FEMNormType norm_type,
1386  std::set<unsigned int> * skip_dimensions) const
1387 {
1388  //short circuit to save time
1389  if (norm_type == DISCRETE_L1 ||
1390  norm_type == DISCRETE_L2 ||
1391  norm_type == DISCRETE_L_INF)
1392  return discrete_var_norm(v,var,norm_type);
1393 
1394  // Not a discrete norm
1395  std::vector<FEMNormType> norms(this->n_vars(), L2);
1396  std::vector<Real> weights(this->n_vars(), 0.0);
1397  norms[var] = norm_type;
1398  weights[var] = 1.0;
1399  Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1400  return val;
1401 }
1402 
1403 
1404 
1406  const SystemNorm & norm,
1407  std::set<unsigned int> * skip_dimensions) const
1408 {
1409  // This function must be run on all processors at once
1410  parallel_object_only();
1411 
1412  LOG_SCOPE ("calculate_norm()", "System");
1413 
1414  // Zero the norm before summation
1415  Real v_norm = 0.;
1416 
1417  if (norm.is_discrete())
1418  {
1419  //Check to see if all weights are 1.0 and all types are equal
1420  FEMNormType norm_type0 = norm.type(0);
1421  unsigned int check_var = 0;
1422  for (; check_var != this->n_vars(); ++check_var)
1423  if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1424  break;
1425 
1426  //All weights were 1.0 so just do the full vector discrete norm
1427  if (check_var == this->n_vars())
1428  {
1429  if (norm_type0 == DISCRETE_L1)
1430  return v.l1_norm();
1431  if (norm_type0 == DISCRETE_L2)
1432  return v.l2_norm();
1433  if (norm_type0 == DISCRETE_L_INF)
1434  return v.linfty_norm();
1435  else
1436  libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
1437  }
1438 
1439  for (unsigned int var=0; var != this->n_vars(); ++var)
1440  {
1441  // Skip any variables we don't need to integrate
1442  if (norm.weight(var) == 0.0)
1443  continue;
1444 
1445  v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1446  }
1447 
1448  return v_norm;
1449  }
1450 
1451  // Localize the potentially parallel vector
1453  local_v->init(v.size(), true, SERIAL);
1454  v.localize (*local_v, _dof_map->get_send_list());
1455 
1456  // I'm not sure how best to mix Hilbert norms on some variables (for
1457  // which we'll want to square then sum then square root) with norms
1458  // like L_inf (for which we'll just want to take an absolute value
1459  // and then sum).
1460  bool using_hilbert_norm = true,
1461  using_nonhilbert_norm = true;
1462 
1463  // Loop over all variables
1464  for (unsigned int var=0; var != this->n_vars(); ++var)
1465  {
1466  // Skip any variables we don't need to integrate
1467  Real norm_weight_sq = norm.weight_sq(var);
1468  if (norm_weight_sq == 0.0)
1469  continue;
1470  Real norm_weight = norm.weight(var);
1471 
1472  // Check for unimplemented norms (rather than just returning 0).
1473  FEMNormType norm_type = norm.type(var);
1474  if ((norm_type==H1) ||
1475  (norm_type==H2) ||
1476  (norm_type==L2) ||
1477  (norm_type==H1_SEMINORM) ||
1478  (norm_type==H2_SEMINORM))
1479  {
1480  if (!using_hilbert_norm)
1481  libmesh_not_implemented();
1482  using_nonhilbert_norm = false;
1483  }
1484  else if ((norm_type==L1) ||
1485  (norm_type==L_INF) ||
1486  (norm_type==W1_INF_SEMINORM) ||
1487  (norm_type==W2_INF_SEMINORM))
1488  {
1489  if (!using_nonhilbert_norm)
1490  libmesh_not_implemented();
1491  using_hilbert_norm = false;
1492  }
1493  else
1494  libmesh_not_implemented();
1495 
1496  const FEType & fe_type = this->get_dof_map().variable_type(var);
1497 
1498  // Allow space for dims 0-3, even if we don't use them all
1499  std::vector<FEBase *> fe_ptrs(4,libmesh_nullptr);
1500  std::vector<QBase *> q_rules(4,libmesh_nullptr);
1501 
1502  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1503 
1504  // Prepare finite elements for each dimension present in the mesh
1505  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
1506  d_it != elem_dims.end(); ++d_it)
1507  {
1508  if (skip_dimensions && skip_dimensions->find(*d_it) != skip_dimensions->end())
1509  continue;
1510 
1511  q_rules[*d_it] =
1512  fe_type.default_quadrature_rule (*d_it).release();
1513 
1514  // Construct finite element object
1515 
1516  fe_ptrs[*d_it] = FEBase::build(*d_it, fe_type).release();
1517 
1518  // Attach quadrature rule to FE object
1519  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it]);
1520  }
1521 
1522  std::vector<dof_id_type> dof_indices;
1523 
1524  // Begin the loop over the elements
1525  for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1526  {
1527  const unsigned int dim = elem->dim();
1528 
1529 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1530 
1531  // One way for implementing this would be to exchange the fe with the FEInterface- class.
1532  // However, it needs to be discussed whether integral-norms make sense for infinite elements.
1533  // or in which sense they could make sense.
1534  if (elem->infinite() )
1535  libmesh_not_implemented();
1536 
1537 #endif
1538 
1539  if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1540  continue;
1541 
1542  FEBase * fe = fe_ptrs[dim];
1543  QBase * qrule = q_rules[dim];
1544  libmesh_assert(fe);
1545  libmesh_assert(qrule);
1546 
1547  const std::vector<Real> & JxW = fe->get_JxW();
1548  const std::vector<std::vector<Real>> * phi = libmesh_nullptr;
1549  if (norm_type == H1 ||
1550  norm_type == H2 ||
1551  norm_type == L2 ||
1552  norm_type == L1 ||
1553  norm_type == L_INF)
1554  phi = &(fe->get_phi());
1555 
1556  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1557  if (norm_type == H1 ||
1558  norm_type == H2 ||
1559  norm_type == H1_SEMINORM ||
1560  norm_type == W1_INF_SEMINORM)
1561  dphi = &(fe->get_dphi());
1562 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1563  const std::vector<std::vector<RealTensor>> * d2phi = libmesh_nullptr;
1564  if (norm_type == H2 ||
1565  norm_type == H2_SEMINORM ||
1566  norm_type == W2_INF_SEMINORM)
1567  d2phi = &(fe->get_d2phi());
1568 #endif
1569 
1570  fe->reinit (elem);
1571 
1572  this->get_dof_map().dof_indices (elem, dof_indices, var);
1573 
1574  const unsigned int n_qp = qrule->n_points();
1575 
1576  const unsigned int n_sf = cast_int<unsigned int>
1577  (dof_indices.size());
1578 
1579  // Begin the loop over the Quadrature points.
1580  for (unsigned int qp=0; qp<n_qp; qp++)
1581  {
1582  if (norm_type == L1)
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 += norm_weight *
1588  JxW[qp] * std::abs(u_h);
1589  }
1590 
1591  if (norm_type == L_INF)
1592  {
1593  Number u_h = 0.;
1594  for (unsigned int i=0; i != n_sf; ++i)
1595  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1596  v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1597  }
1598 
1599  if (norm_type == H1 ||
1600  norm_type == H2 ||
1601  norm_type == L2)
1602  {
1603  Number u_h = 0.;
1604  for (unsigned int i=0; i != n_sf; ++i)
1605  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1606  v_norm += norm_weight_sq *
1607  JxW[qp] * TensorTools::norm_sq(u_h);
1608  }
1609 
1610  if (norm_type == H1 ||
1611  norm_type == H2 ||
1612  norm_type == H1_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 += norm_weight_sq *
1618  JxW[qp] * grad_u_h.norm_sq();
1619  }
1620 
1621  if (norm_type == W1_INF_SEMINORM)
1622  {
1623  Gradient grad_u_h;
1624  for (unsigned int i=0; i != n_sf; ++i)
1625  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1626  v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1627  }
1628 
1629 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1630  if (norm_type == H2 ||
1631  norm_type == H2_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 += norm_weight_sq *
1637  JxW[qp] * hess_u_h.norm_sq();
1638  }
1639 
1640  if (norm_type == W2_INF_SEMINORM)
1641  {
1642  Tensor hess_u_h;
1643  for (unsigned int i=0; i != n_sf; ++i)
1644  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1645  v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
1646  }
1647 #endif
1648  }
1649  }
1650 
1651  // Need to delete the FE and quadrature objects to prevent a memory leak
1652  for (std::size_t i=0; i<fe_ptrs.size(); i++)
1653  if (fe_ptrs[i])
1654  delete fe_ptrs[i];
1655 
1656  for (std::size_t i=0; i<q_rules.size(); i++)
1657  if (q_rules[i])
1658  delete q_rules[i];
1659  }
1660 
1661  if (using_hilbert_norm)
1662  {
1663  this->comm().sum(v_norm);
1664  v_norm = std::sqrt(v_norm);
1665  }
1666  else
1667  {
1668  this->comm().max(v_norm);
1669  }
1670 
1671  return v_norm;
1672 }
1673 
1674 
1675 
1676 std::string System::get_info() const
1677 {
1678  std::ostringstream oss;
1679 
1680 
1681  const std::string & sys_name = this->name();
1682 
1683  oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1684  << " Type \"" << this->system_type() << "\"\n"
1685  << " Variables=";
1686 
1687  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1688  {
1689  const VariableGroup & vg_description (this->variable_group(vg));
1690 
1691  if (vg_description.n_variables() > 1) oss << "{ ";
1692  for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
1693  oss << "\"" << vg_description.name(vn) << "\" ";
1694  if (vg_description.n_variables() > 1) oss << "} ";
1695  }
1696 
1697  oss << '\n';
1698 
1699  oss << " Finite Element Types=";
1700 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1701  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1702  oss << "\""
1703  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1704  << "\" ";
1705 #else
1706  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1707  {
1708  oss << "\""
1709  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1710  << "\", \""
1711  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1712  << "\" ";
1713  }
1714 
1715  oss << '\n' << " Infinite Element Mapping=";
1716  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1717  oss << "\""
1718  << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1719  << "\" ";
1720 #endif
1721 
1722  oss << '\n';
1723 
1724  oss << " Approximation Orders=";
1725  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1726  {
1727 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1728  oss << "\""
1729  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1730  << "\" ";
1731 #else
1732  oss << "\""
1733  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1734  << "\", \""
1735  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1736  << "\" ";
1737 #endif
1738  }
1739 
1740  oss << '\n';
1741 
1742  oss << " n_dofs()=" << this->n_dofs() << '\n';
1743  oss << " n_local_dofs()=" << this->n_local_dofs() << '\n';
1744 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1745  oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1746  oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1747 #endif
1748 
1749  oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1750  oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1751  // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1752 
1753  oss << this->get_dof_map().get_info();
1754 
1755  return oss.str();
1756 }
1757 
1758 
1759 
1761  const std::string & name))
1762 {
1764 
1766  {
1767  libmesh_here();
1768  libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1769  << std::endl;
1770 
1772  }
1773 
1775 }
1776 
1777 
1778 
1780 {
1782  {
1783  libmesh_here();
1784  libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1785  << std::endl;
1786 
1788  }
1789 
1790  _init_system_object = &init_in;
1791 }
1792 
1793 
1794 
1796  const std::string & name))
1797 {
1799 
1801  {
1802  libmesh_here();
1803  libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1804  << std::endl;
1805 
1807  }
1808 
1810 }
1811 
1812 
1813 
1815 {
1817  {
1818  libmesh_here();
1819  libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1820  << std::endl;
1821 
1823  }
1824 
1825  _assemble_system_object = &assemble_in;
1826 }
1827 
1828 
1829 
1831  const std::string & name))
1832 {
1834 
1836  {
1837  libmesh_here();
1838  libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1839  << std::endl;
1840 
1842  }
1843 
1845 }
1846 
1847 
1848 
1850 {
1852  {
1853  libmesh_here();
1854  libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1855  << std::endl;
1856 
1858  }
1859 
1860  _constrain_system_object = &constrain;
1861 }
1862 
1863 
1864 
1866  const std::string &,
1867  const QoISet &))
1868 {
1870 
1872  {
1873  libmesh_here();
1874  libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1875  << std::endl;
1876 
1878  }
1879 
1881 }
1882 
1883 
1884 
1886 {
1888  {
1889  libmesh_here();
1890  libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1891  << std::endl;
1892 
1894  }
1895 
1896  _qoi_evaluate_object = &qoi_in;
1897 }
1898 
1899 
1900 
1901 void System::attach_QOI_derivative(void fptr(EquationSystems &, const std::string &,
1902  const QoISet &, bool, bool))
1903 {
1905 
1907  {
1908  libmesh_here();
1909  libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1910  << std::endl;
1911 
1913  }
1914 
1916 }
1917 
1918 
1919 
1921 {
1923  {
1924  libmesh_here();
1925  libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1926  << std::endl;
1927 
1929  }
1930 
1931  _qoi_evaluate_derivative_object = &qoi_derivative;
1932 }
1933 
1934 
1935 
1937 {
1938  // Call the user-provided initialization function,
1939  // if it was provided
1941  this->_init_system_function (_equation_systems, this->name());
1942 
1943  // ...or the user-provided initialization object.
1946 }
1947 
1948 
1949 
1951 {
1952  // Call the user-provided assembly function,
1953  // if it was provided
1956 
1957  // ...or the user-provided assembly object.
1960 }
1961 
1962 
1963 
1965 {
1966  // Call the user-provided constraint function,
1967  // if it was provided
1970 
1971  // ...or the user-provided constraint object.
1974 }
1975 
1976 
1977 
1978 void System::user_QOI (const QoISet & qoi_indices)
1979 {
1980  // Call the user-provided quantity of interest function,
1981  // if it was provided
1983  this->_qoi_evaluate_function(_equation_systems, this->name(), qoi_indices);
1984 
1985  // ...or the user-provided QOI function object.
1987  this->_qoi_evaluate_object->qoi(qoi_indices);
1988 }
1989 
1990 
1991 
1992 void System::user_QOI_derivative(const QoISet & qoi_indices,
1993  bool include_liftfunc,
1994  bool apply_constraints)
1995 {
1996  // Call the user-provided quantity of interest derivative,
1997  // if it was provided
2000  (_equation_systems, this->name(), qoi_indices, include_liftfunc,
2001  apply_constraints);
2002 
2003  // ...or the user-provided QOI derivative function object.
2006  (qoi_indices, include_liftfunc, apply_constraints);
2007 }
2008 
2009 
2010 
2011 Number System::point_value(unsigned int var, const Point & p, const bool insist_on_success) const
2012 {
2013  // This function must be called on every processor; there's no
2014  // telling where in the partition p falls.
2015  parallel_object_only();
2016 
2017  // And every processor had better agree about which point we're
2018  // looking for
2019 #ifndef NDEBUG
2020  libmesh_assert(this->comm().verify(p(0)));
2021 #if LIBMESH_DIM > 1
2022  libmesh_assert(this->comm().verify(p(1)));
2023 #endif
2024 #if LIBMESH_DIM > 2
2025  libmesh_assert(this->comm().verify(p(2)));
2026 #endif
2027 #endif // NDEBUG
2028 
2029  // Get a reference to the mesh object associated with the system object that calls this function
2030  const MeshBase & mesh = this->get_mesh();
2031 
2032  // Use an existing PointLocator or create a new one
2033  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2034  PointLocatorBase & locator = *locator_ptr;
2035 
2036  if (!insist_on_success || !mesh.is_serial())
2037  locator.enable_out_of_mesh_mode();
2038 
2039  // Get a pointer to the element that contains P
2040  const Elem * e = locator(p);
2041 
2042  Number u = 0;
2043 
2044  if (e && this->get_dof_map().is_evaluable(*e, var))
2045  u = point_value(var, p, *e);
2046 
2047  // If I have an element containing p, then let's let everyone know
2048  processor_id_type lowest_owner =
2049  (e && (e->processor_id() == this->processor_id())) ?
2050  this->processor_id() : this->n_processors();
2051  this->comm().min(lowest_owner);
2052 
2053  // Everybody should get their value from a processor that was able
2054  // to compute it.
2055  // If nobody admits owning the point, we have a problem.
2056  if (lowest_owner != this->n_processors())
2057  this->comm().broadcast(u, lowest_owner);
2058  else
2059  libmesh_assert(!insist_on_success);
2060 
2061  return u;
2062 }
2063 
2064 Number System::point_value(unsigned int var, const Point & p, const Elem & e) const
2065 {
2066  // Ensuring that the given point is really in the element is an
2067  // expensive assert, but as long as debugging is turned on we might
2068  // as well try to catch a particularly nasty potential error
2070 
2071  // Get the dof map to get the proper indices for our computation
2072  const DofMap & dof_map = this->get_dof_map();
2073 
2074  // Make sure we can evaluate on this element.
2075  libmesh_assert (dof_map.is_evaluable(e, var));
2076 
2077  // Need dof_indices for phi[i][j]
2078  std::vector<dof_id_type> dof_indices;
2079 
2080  // Fill in the dof_indices for our element
2081  dof_map.dof_indices (&e, dof_indices, var);
2082 
2083  // Get the no of dofs associated with this point
2084  const unsigned int num_dofs = cast_int<unsigned int>
2085  (dof_indices.size());
2086 
2087  FEType fe_type = dof_map.variable_type(var);
2088 
2089  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h.
2090  Point coor = FEInterface::inverse_map(e.dim(), fe_type, &e, p);
2091 
2092  // get the shape function value via the FEInterface to also handle the case
2093  // of infinite elements correcly, the shape function is not fe->phi().
2094  FEComputeData fe_data(this->get_equation_systems(), coor);
2095  FEInterface::compute_data(e.dim(), fe_type, &e, fe_data);
2096 
2097  // Get ready to accumulate a value
2098  Number u = 0;
2099 
2100  for (unsigned int l=0; l<num_dofs; l++)
2101  {
2102  u += fe_data.shape[l]*this->current_solution (dof_indices[l]);
2103  }
2104 
2105  return u;
2106 }
2107 
2108 
2109 
2110 Number System::point_value(unsigned int var, const Point & p, const Elem * e) const
2111 {
2112  libmesh_assert(e);
2113  return this->point_value(var, p, *e);
2114 }
2115 
2116 
2117 
2118 Gradient System::point_gradient(unsigned int var, const Point & p, const bool insist_on_success) const
2119 {
2120  // This function must be called on every processor; there's no
2121  // telling where in the partition p falls.
2122  parallel_object_only();
2123 
2124  // And every processor had better agree about which point we're
2125  // looking for
2126 #ifndef NDEBUG
2127  libmesh_assert(this->comm().verify(p(0)));
2128 #if LIBMESH_DIM > 1
2129  libmesh_assert(this->comm().verify(p(1)));
2130 #endif
2131 #if LIBMESH_DIM > 2
2132  libmesh_assert(this->comm().verify(p(2)));
2133 #endif
2134 #endif // NDEBUG
2135 
2136  // Get a reference to the mesh object associated with the system object that calls this function
2137  const MeshBase & mesh = this->get_mesh();
2138 
2139  // Use an existing PointLocator or create a new one
2140  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2141  PointLocatorBase & locator = *locator_ptr;
2142 
2143  if (!insist_on_success || !mesh.is_serial())
2144  locator.enable_out_of_mesh_mode();
2145 
2146  // Get a pointer to the element that contains P
2147  const Elem * e = locator(p);
2148 
2149  Gradient grad_u;
2150 
2151  if (e && this->get_dof_map().is_evaluable(*e, var))
2152  grad_u = point_gradient(var, p, *e);
2153 
2154  // If I have an element containing p, then let's let everyone know
2155  processor_id_type lowest_owner =
2156  (e && (e->processor_id() == this->processor_id())) ?
2157  this->processor_id() : this->n_processors();
2158  this->comm().min(lowest_owner);
2159 
2160  // Everybody should get their value from a processor that was able
2161  // to compute it.
2162  // If nobody admits owning the point, we may have a problem.
2163  if (lowest_owner != this->n_processors())
2164  this->comm().broadcast(grad_u, lowest_owner);
2165  else
2166  libmesh_assert(!insist_on_success);
2167 
2168  return grad_u;
2169 }
2170 
2171 
2172 Gradient System::point_gradient(unsigned int var, const Point & p, const Elem & e) const
2173 {
2174  // Ensuring that the given point is really in the element is an
2175  // expensive assert, but as long as debugging is turned on we might
2176  // as well try to catch a particularly nasty potential error
2178 
2179 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2180  if (e.infinite())
2181  libmesh_not_implemented();
2182 #endif
2183 
2184  // Get the dof map to get the proper indices for our computation
2185  const DofMap & dof_map = this->get_dof_map();
2186 
2187  // Make sure we can evaluate on this element.
2188  libmesh_assert (dof_map.is_evaluable(e, var));
2189 
2190  // Need dof_indices for phi[i][j]
2191  std::vector<dof_id_type> dof_indices;
2192 
2193  // Fill in the dof_indices for our element
2194  dof_map.dof_indices (&e, dof_indices, var);
2195 
2196  // Get the no of dofs associated with this point
2197  const unsigned int num_dofs = cast_int<unsigned int>
2198  (dof_indices.size());
2199 
2200  FEType fe_type = dof_map.variable_type(var);
2201 
2202  // Build a FE again so we can calculate u(p)
2203  UniquePtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2204 
2205  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2206  // Build a vector of point co-ordinates to send to reinit
2207  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2208 
2209  // Get the values of the shape function derivatives
2210  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
2211 
2212  // Reinitialize the element and compute the shape function values at coor
2213  fe->reinit (&e, &coor);
2214 
2215  // Get ready to accumulate a gradient
2216  Gradient grad_u;
2217 
2218  for (unsigned int l=0; l<num_dofs; l++)
2219  {
2220  grad_u.add_scaled (dphi[l][0], this->current_solution (dof_indices[l]));
2221  }
2222 
2223  return grad_u;
2224 }
2225 
2226 
2227 
2228 Gradient System::point_gradient(unsigned int var, const Point & p, const Elem * e) const
2229 {
2230  libmesh_assert(e);
2231  return this->point_gradient(var, p, *e);
2232 }
2233 
2234 
2235 
2236 // We can only accumulate a hessian with --enable-second
2237 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2238 Tensor System::point_hessian(unsigned int var, const Point & p, const bool insist_on_success) const
2239 {
2240  // This function must be called on every processor; there's no
2241  // telling where in the partition p falls.
2242  parallel_object_only();
2243 
2244  // And every processor had better agree about which point we're
2245  // looking for
2246 #ifndef NDEBUG
2247  libmesh_assert(this->comm().verify(p(0)));
2248 #if LIBMESH_DIM > 1
2249  libmesh_assert(this->comm().verify(p(1)));
2250 #endif
2251 #if LIBMESH_DIM > 2
2252  libmesh_assert(this->comm().verify(p(2)));
2253 #endif
2254 #endif // NDEBUG
2255 
2256  // Get a reference to the mesh object associated with the system object that calls this function
2257  const MeshBase & mesh = this->get_mesh();
2258 
2259  // Use an existing PointLocator or create a new one
2260  UniquePtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2261  PointLocatorBase & locator = *locator_ptr;
2262 
2263  if (!insist_on_success || !mesh.is_serial())
2264  locator.enable_out_of_mesh_mode();
2265 
2266  // Get a pointer to the element that contains P
2267  const Elem * e = locator(p);
2268 
2269  Tensor hess_u;
2270 
2271  if (e && this->get_dof_map().is_evaluable(*e, var))
2272  hess_u = point_hessian(var, p, *e);
2273 
2274  // If I have an element containing p, then let's let everyone know
2275  processor_id_type lowest_owner =
2276  (e && (e->processor_id() == this->processor_id())) ?
2277  this->processor_id() : this->n_processors();
2278  this->comm().min(lowest_owner);
2279 
2280  // Everybody should get their value from a processor that was able
2281  // to compute it.
2282  // If nobody admits owning the point, we may have a problem.
2283  if (lowest_owner != this->n_processors())
2284  this->comm().broadcast(hess_u, lowest_owner);
2285  else
2286  libmesh_assert(!insist_on_success);
2287 
2288  return hess_u;
2289 }
2290 
2291 Tensor System::point_hessian(unsigned int var, const Point & p, const Elem & e) const
2292 {
2293  // Ensuring that the given point is really in the element is an
2294  // expensive assert, but as long as debugging is turned on we might
2295  // as well try to catch a particularly nasty potential error
2297 
2298 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2299  if (e.infinite())
2300  libmesh_not_implemented();
2301 #endif
2302 
2303  // Get the dof map to get the proper indices for our computation
2304  const DofMap & dof_map = this->get_dof_map();
2305 
2306  // Make sure we can evaluate on this element.
2307  libmesh_assert (dof_map.is_evaluable(e, var));
2308 
2309  // Need dof_indices for phi[i][j]
2310  std::vector<dof_id_type> dof_indices;
2311 
2312  // Fill in the dof_indices for our element
2313  dof_map.dof_indices (&e, dof_indices, var);
2314 
2315  // Get the no of dofs associated with this point
2316  const unsigned int num_dofs = cast_int<unsigned int>
2317  (dof_indices.size());
2318 
2319  FEType fe_type = dof_map.variable_type(var);
2320 
2321  // Build a FE again so we can calculate u(p)
2322  UniquePtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2323 
2324  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2325  // Build a vector of point co-ordinates to send to reinit
2326  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2327 
2328  // Get the values of the shape function derivatives
2329  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
2330 
2331  // Reinitialize the element and compute the shape function values at coor
2332  fe->reinit (&e, &coor);
2333 
2334  // Get ready to accumulate a hessian
2335  Tensor hess_u;
2336 
2337  for (unsigned int l=0; l<num_dofs; l++)
2338  {
2339  hess_u.add_scaled (d2phi[l][0], this->current_solution (dof_indices[l]));
2340  }
2341 
2342  return hess_u;
2343 }
2344 
2345 
2346 
2347 Tensor System::point_hessian(unsigned int var, const Point & p, const Elem * e) const
2348 {
2349  libmesh_assert(e);
2350  return this->point_hessian(var, p, *e);
2351 }
2352 
2353 
2354 
2355 #else
2356 Tensor System::point_hessian(unsigned int, const Point &, const bool) const
2357 {
2358  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2359 
2360  // Avoid compiler warnings
2361  return Tensor();
2362 }
2363 
2364 Tensor System::point_hessian(unsigned int, const Point &, const Elem &) const
2365 {
2366  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2367 
2368  // Avoid compiler warnings
2369  return Tensor();
2370 }
2371 
2372 Tensor System::point_hessian(unsigned int, const Point &, const Elem *) const
2373 {
2374  libmesh_error_msg("We can only accumulate a hessian with --enable-second");
2375 
2376  // Avoid compiler warnings
2377  return Tensor();
2378 }
2379 
2380 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
2381 
2382 } // namespace libMesh
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Allows one to set the QoI index controlling whether the vector identified by vec_name represents a so...
Definition: system.C:905
Number point_value(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2011
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
vectors_iterator vectors_end()
End of vectors container.
Definition: system.h:2238
FEFamily family
The type of finite element.
Definition: fe_type.h:203
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1545
dof_id_type n_constrained_dofs() const
Definition: system.C:155
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:281
std::map< std::string, unsigned short int > _variable_numbers
The variable numbers corresponding to user-specified names, useful for name-based lookups...
Definition: system.h:1903
double abs(double a)
const FEType & type() const
Definition: variable.h:119
std::string get_info() const
Gets summary info about the sparsity bandwidth and constraints.
Definition: dof_map.C:2635
This is the EquationSystems class.
virtual void clear()
Clear all the data structures associated with the system.
Definition: system.C:203
Abstract base class to be used for system initialization.
Definition: system.h:95
virtual bool is_serial() const
Definition: mesh_base.h:140
Assembly * _assemble_system_object
Object that assembles the system.
Definition: system.h:1822
virtual numeric_index_type size() const =0
bool _basic_system_only
Holds true if the components of more advanced system types (e.g.
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
Holds the type of a vector.
Definition: system.h:1933
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
bool _is_initialized
true when additional vectors and variables do not require immediate initialization, false otherwise.
Definition: system.h:1952
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Function to impose constraints.
Definition: system.h:1827
NumericVector< Number > & add_adjoint_solution(unsigned int i=0)
Definition: system.C:977
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
virtual void forward_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Solves for parameter sensitivities using the forward method.
Definition: system.h:2300
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
Data structure for specifying which Parameters should be independent variables in a parameter sensiti...
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
bool is_evaluable(const DofObjectSubclass &obj, unsigned int var_num=libMesh::invalid_uint) const
Definition: dof_map.C:2357
Real norm_sq() const
Definition: type_tensor.h:1270
void add_scaled(const TypeTensor< T2 > &, const T)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:777
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
int vector_is_adjoint(const std::string &vec_name) const
Definition: system.C:916
virtual Real l2_norm() const =0
void min(T &r) const
Take a local variable and replace it with the minimum of it&#39;s values on all processors.
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
virtual void reinit()
Reinitializes degrees of freedom and other required data on the current mesh.
Definition: system.C:385
unsigned int dim
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
virtual void initialize()=0
Initialization function.
NumericVector< Number > & get_sensitivity_solution(unsigned int i=0)
Definition: system.C:936
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
Real weight(unsigned int var) const
Definition: system_norm.h:292
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
bool has_variable(const std::string &var) const
Definition: system.C:1256
const unsigned int _sys_number
The number associated with this system.
Definition: system.h:1887
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
Calls user qoi derivative function.
Definition: system.C:501
NumericVector< Number > & get_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1081
virtual void assemble()=0
Assembly function.
Constraint * _constrain_system_object
Object that constrains the system.
Definition: system.h:1833
virtual void assemble()
Prepares matrix and _dof_map for matrix assembly.
Definition: system.C:479
virtual void init_data()
Initializes the data for the system.
Definition: system.C:260
processor_id_type n_processors() const
void zero_variable(NumericVector< Number > &v, unsigned int var_num) const
Zeroes all dofs in v that correspond to variable number var_num.
Definition: system.C:1327
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
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
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
bool vector_preservation(const std::string &vec_name) const
Definition: system.C:895
std::map< std::string, NumericVector< Number > * >::const_iterator const_vectors_iterator
Definition: system.h:749
virtual void user_initialization()
Calls user&#39;s attached initialization function, or is overridden by the user in derived classes...
Definition: system.C:1936
Numeric vector.
Definition: dof_map.h:66
std::vector< Variable > _variables
The Variable in this System.
Definition: system.h:1892
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:43
NumericVector< Number > & add_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1071
void attach_QOI_derivative(void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints))
Register a user function for evaluating derivatives of a quantity of interest with respect to test fu...
Definition: system.C:1901
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
Definition: system.h:2226
void attach_QOI_derivative_object(QOIDerivative &qoi_derivative)
Register a user object for evaluating derivatives of a quantity of interest with respect to test func...
Definition: system.C:1920
virtual void adjoint_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Solves for parameter sensitivities using the adjoint method.
Definition: system.h:2291
const VariableGroup & variable_group(unsigned int vg) const
Return a constant reference to VariableGroup vg.
Definition: system.h:2124
Abstract base class to be used for system assembly.
Definition: system.h:119
NumericVector< Number > & add_weighted_sensitivity_solution()
Definition: system.C:956
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)
Function to evaluate quantity of interest.
Definition: system.h:1838
virtual Real l1_norm() const =0
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:1657
void init()
Initializes degrees of freedom on the current mesh.
Definition: system.C:235
long double max(long double a, double b)
UniquePtr< DofMap > _dof_map
Data structure describing the relationship between nodes, variables, etc...
Definition: system.h:1865
const std::string & name() const
Definition: system.h:1998
NumericVector< Number > & get_weighted_sensitivity_solution()
Definition: system.C:963
std::string get_info() const
Definition: system.C:1676
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void prolong_vectors()
Prolong vectors after the mesh has refined.
Definition: system.C:375
virtual void user_QOI(const QoISet &qoi_indices)
Calls user&#39;s attached quantity of interest function, or is overridden by the user in derived classes...
Definition: system.C:1978
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void attach_QOI_function(void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
Register a user function for evaluating the quantities of interest, whose values should be placed in ...
Definition: system.C:1865
virtual bool infinite() const =0
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
QOIDerivative * _qoi_evaluate_derivative_object
Object to compute derivatives of quantities of interest.
Definition: system.h:1859
void remove_vector(const std::string &vec_name)
Removes the additional vector vec_name from this system.
Definition: system.C:719
NumericVector< Number > & add_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1009
virtual void qoi_derivative(const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)=0
Quantity of interest derivative function.
virtual void qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
Solves for the derivative of each of the system&#39;s quantities of interest q in qoi[qoi_indices] with r...
Definition: system.C:515
void attach_constraint_object(Constraint &constrain)
Register a user object for imposing constraints.
Definition: system.C:1849
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const NumericVector< Number > * request_vector(const std::string &vec_name) const
Definition: system.C:736
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
void attach_init_object(Initialization &init)
Register a user class to use to initialize the system.
Definition: system.C:1779
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
FEFamily
defines an enum for finite element families.
void update_global_solution(std::vector< Number > &global_soln) const
Fill the input vector global_soln so that it contains the global solution on all processors.
Definition: system.C:662
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)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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
ParallelType
Defines an enum for parallel data structure types.
Data structure for holding completed parameter sensitivity calculations.
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
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
The VariableGroup in this System.
Definition: system.h:1897
bool is_initialized()
Definition: system.h:2070
This is a base class for classes which represent subsets of the dofs of a System. ...
Definition: system_subset.h:42
unsigned int n_variables() const
Definition: variable.h:217
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
void attach_constraint_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function for imposing constraints.
Definition: system.C:1830
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)
Lets the appropriate child of FEBase compute the requested data for the input specified in data...
Definition: fe_interface.C:827
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
virtual void qoi(const QoISet &qoi_indices)=0
Quantity of interest function.
QOI * _qoi_evaluate_object
Object to compute quantities of interest.
Definition: system.h:1845
Initialization * _init_system_object
Object that initializes the system.
Definition: system.h:1811
System & operator=(const System &)
This isn&#39;t a copyable object, so let&#39;s make sure nobody tries.
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:569
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:257
Tensor point_hessian(unsigned int var, const Point &p, const bool insist_on_success=true) const
Definition: system.C:2238
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Function that initializes the system.
Definition: system.h:1805
This is the base class for point locators.
virtual void restrict_vectors()
Restrict vectors after the mesh has coarsened.
Definition: system.C:322
void attach_init_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in initializing the system.
Definition: system.C:1760
void broadcast(T &data, const unsigned int root_id=0) const
Take a local value and broadcast it to all processors.
Real norm_sq() const
Definition: type_vector.h:940
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:206
void append(const std::string &var_name)
Appends a variable to the group.
Definition: variable.h:276
std::map< std::string, int > _vector_is_adjoint
Holds non-negative if a vector by that name should be projected using adjoint constraints/BCs, -1 if primal.
Definition: system.h:1928
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
This class implements reference counting.
This class forms the base class for all other classes that are expected to be implemented in parallel...
This class defines a logically grouped set of variables in the system.
Definition: variable.h:172
NumericVector< Number > & add_adjoint_rhs(unsigned int i=0)
Definition: system.C:1041
virtual Real linfty_norm() const =0
virtual void restrict_solve_to(const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)
After calling this method, any solve will be restricted to the given subdomain.
Definition: system.C:470
virtual void user_QOI_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
Calls user&#39;s attached quantity of interest derivative function, or is overridden by the user in deriv...
Definition: system.C:1992
virtual void reinit_constraints()
Reinitializes the constraints for this system.
Definition: system.C:414
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:2118
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:1383
void attach_assemble_object(Assembly &assemble)
Register a user object to use in assembling the system matrix and RHS.
Definition: system.C:1814
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Vector iterator typedefs.
Definition: system.h:748
FEFamily radial_family
For InfFE, family contains the radial shape family, while base_family contains the approximation type...
Definition: fe_type.h:249
std::size_t size() const
bool _solution_projection
Holds true if the solution vector should be projected onto a changed grid, false if it should be zero...
Definition: system.h:1940
EquationSystems & _equation_systems
Constant reference to the EquationSystems object used for the simulation.
Definition: system.h:1871
void local_dof_indices(const unsigned int var, std::set< dof_id_type > &var_indices) const
Fills the std::set with the degrees of freedom on the local processor corresponding the the variable ...
Definition: system.C:1297
virtual void constrain()=0
Constraint function.
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:531
Real weight_sq(unsigned int var) const
Definition: system_norm.h:338
void attach_QOI_object(QOI &qoi)
Register a user object for evaluating the quantities of interest, whose values should be placed in Sy...
Definition: system.C:1885
virtual std::string system_type() const
Definition: system.h:471
virtual void user_assembly()
Calls user&#39;s attached assembly function, or is overridden by the user in derived classes.
Definition: system.C:1950
dof_id_type n_local_constrained_dofs() const
Definition: system.C:170
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
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)
Allows one to set the boolean controlling whether the vector identified by vec_name should be "preser...
Definition: system.C:887
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
void process_constraints(MeshBase &)
Postprocesses any constrained degrees of freedom to be constrained only in terms of unconstrained dof...
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
Some systems need an arbitrary number of vectors.
Definition: system.h:1916
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:290
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
OStreamProxy out
const std::string _sys_name
A name associated with this system.
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)
Adds the variable var to the list of variables for this system.
Definition: system.C:1191
virtual unsigned int dim() const =0
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
virtual void re_update()
Re-update the local values when the mesh has changed.
Definition: system.C:446
std::map< std::string, bool > _vector_projections
Holds true if a vector by that name should be projected onto a changed grid, false if it should be ze...
Definition: system.h:1922
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Abstract base class to be used for system constraints.
Definition: system.h:143
dof_id_type n_dofs() const
Definition: system.C:148
void create_dof_constraints(const MeshBase &, Real time=0)
Rebuilds the raw degree of freedom and DofObject constraints.
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
Sets v(i) = value.
void prepare_send_list()
Takes the _send_list vector (which may have duplicate entries) and sorts it.
Definition: dof_map.C:1635
SubsetSolveMode
defines an enum for the question what happens to the dofs outside the given subset when a system is s...
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
virtual void user_constrain()
Calls user&#39;s attached constraint function, or is overridden by the user in derived classes...
Definition: system.C:1964
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet())
Calls user qoi function.
Definition: system.C:490
NumericVector< Number > & add_sensitivity_solution(unsigned int i=0)
Definition: system.C:926
const std::string & name(unsigned int v) const
Definition: variable.h:246
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Function that assembles the system.
Definition: system.h:1816
virtual ~System()
Destructor.
Definition: system.C:118
virtual void enable_out_of_mesh_mode()=0
Enables out-of-mesh mode.
System(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
Definition: system.C:59
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
void get_all_variable_numbers(std::vector< unsigned int > &all_variable_numbers) const
Fills all_variable_numbers with all the variable numbers for the variables that have been added to th...
Definition: system.C:1278
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:394
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
MeshBase & _mesh
Constant reference to the mesh data structure used for the simulation.
Definition: system.h:1877
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
unsigned int n_points() const
Definition: quadrature.h:113
UniquePtr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:534
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
Real discrete_var_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type) const
Finds the discrete norm for the entries in the vector corresponding to Dofs associated with var...
Definition: system.C:1364
const std::string & vector_name(const unsigned int vec_num) const
Definition: system.C:852
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const
This class forms the foundation from which generic finite elements may be derived.
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition: elem.C:2448
Abstract base class to be used for quantities of interest.
Definition: system.h:167
NumericVector< Number > & get_weighted_sensitivity_adjoint_solution(unsigned int i=0)
Definition: system.C:1021
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
FEMNormType
defines an enum for norms defined on vectors of finite element coefficients
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)
Function to evaluate quantity of interest derivative.
Definition: system.h:1850
Abstract base class to be used for derivatives of quantities of interest.
Definition: system.h:191
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1051