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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo