DOpE
pdeproblemcontainer.h
Go to the documentation of this file.
1 
24 #ifndef _PDEProblemContainer_H_
25 #define _PDEProblemContainer_H_
26 
27 #include "dopeexceptionhandler.h"
28 #include "outputhandler.h"
29 #include "functionalinterface.h"
30 #include "dofhandler_wrapper.h"
31 #include "fevalues_wrapper.h"
32 #include "function_wrapper.h"
33 #include "statespacetimehandler.h"
34 #include "primaldirichletdata.h"
35 #include "elementdatacontainer.h"
36 #include "facedatacontainer.h"
37 #include "stateproblem.h"
39 #include <deal.II/multigrid/mg_dof_handler.h>
40 #include "dopetypes.h"
41 #include "dwrdatacontainer.h"
42 
43 #include <lac/vector.h>
44 #include <lac/full_matrix.h>
45 #include <grid/tria_iterator.h>
46 #include <dofs/dof_handler.h>
47 #include <dofs/dof_accessor.h>
48 #include <dofs/dof_tools.h>
49 #include <fe/fe_system.h>
50 #include <fe/fe_values.h>
51 #include <fe/mapping.h>
52 #include <base/quadrature_lib.h>
53 #include <lac/block_sparsity_pattern.h>
54 #include <lac/compressed_simple_sparsity_pattern.h>
55 
56 // Multi-level routines (step-16 in deal.II)
57 #include <deal.II/multigrid/mg_dof_handler.h>
58 #include <deal.II/multigrid/mg_constrained_dofs.h>
59 #include <deal.II/multigrid/multigrid.h>
60 #include <deal.II/multigrid/mg_transfer.h>
61 #include <deal.II/multigrid/mg_tools.h>
62 #include <deal.II/multigrid/mg_coarse.h>
63 #include <deal.II/multigrid/mg_smoother.h>
64 #include <deal.II/multigrid/mg_matrix.h>
65 
66 
67 
68 #include <assert.h>
69 #include <string>
70 #include <vector>
71 
72 namespace DOpE
73 {
74  //Predeclaration necessary
75  template<typename VECTOR>
76  class DOpEOutputHandler;
77  template<typename VECTOR>
78  class DOpEExceptionHandler;
80 
97  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
98  int dealdim, template<int, int> class FE = dealii::FESystem,
99  template<int, int> class DH = dealii::DoFHandler>
101  {
102  public:
103  PDEProblemContainer(PDE& pde,
105 
106  /******************************************************/
107 
108  virtual
110 
111  /******************************************************/
112 
113  virtual std::string
114  GetName() const
115  {
116  return "PDEProblemContainer";
117  }
118 
119  /******************************************************/
123  StateProblem<
124  PDEProblemContainer<PDE, DD, SPARSITYPATTERN, VECTOR, dealdim, FE,
125  DH>, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim>&
127  {
128  if (_state_problem == NULL)
129  {
130  _state_problem = new StateProblem<
131  PDEProblemContainer<PDE, DD, SPARSITYPATTERN, VECTOR, dealdim,
132  FE, DH>, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim>(*this,
133  this->GetPDE());
134  }
135  return *_state_problem;
136  }
137 
138  //TODO This is Pfush needed to split into different subproblems and allow optproblem to
139  //be substituted as any of these problems. Can be removed once the splitting is complete.
142  {
143  return *this;
144  }
145  /******************************************************/
146 
154  void
155  ReInit(std::string algo_type);
156 
157  /******************************************************/
158 
159  void
161  {
162  _OutputHandler = OH;
163  }
164 
165  /******************************************************/
166 
167  void
169  {
170  _ExceptionHandler = OH;
171  }
172 
173  /******************************************************/
174 
182  void
183  SetType(std::string type, unsigned int num = 0);
184 
185  /******************************************************/
186 
206  template<typename DATACONTAINER>
207  double
208  ElementFunctional(const DATACONTAINER& edc);
209 
210  /******************************************************/
211 
224  double
226  const std::map<std::string, const dealii::Vector<double>*> &param_values,
227  const std::map<std::string, const VECTOR*> &domain_values);
228 
229  /******************************************************/
230 
250  template<typename FACEDATACONTAINER>
251  double
252  BoundaryFunctional(const FACEDATACONTAINER& fdc);
253 
254  /******************************************************/
255 
262  template<typename FACEDATACONTAINER>
263  double
264  FaceFunctional(const FACEDATACONTAINER& fdc);
265 
266  /******************************************************/
273  double
275  const std::map<std::string, const dealii::Vector<double>*> &values,
276  const std::map<std::string, const VECTOR*> &block_values);
277 
278  /******************************************************/
279 
303  template<typename DATACONTAINER>
304  void
305  ElementEquation(const DATACONTAINER& edc,
306  dealii::Vector<double> &local_vector, double scale,
307  double scale_ico);
308 
309  /******************************************************/
310 
311  /******************************************************/
312 
318  template<typename DATACONTAINER>
319  void
320  ElementTimeEquation(const DATACONTAINER& dc,
321  dealii::Vector<double> &local_vector, double scale = 1.);
322 
323  /******************************************************/
324 
331  template<typename DATACONTAINER>
332  void
333  ElementTimeEquationExplicit(const DATACONTAINER& dc,
334  dealii::Vector<double> &local_vector, double scale = 1.);
335 
336  /******************************************************/
337 
349  template<typename DATACONTAINER>
350  void
351  ElementRhs(const DATACONTAINER& dc,
352  dealii::Vector<double> &local_vector, double scale = 1.);
353 
354  /******************************************************/
355 
368  void
369  PointRhs(
370  const std::map<std::string, const dealii::Vector<double>*> &param_values,
371  const std::map<std::string, const VECTOR*> &domain_values,
372  VECTOR& rhs_vector, double scale = 1.);
373 
374  /******************************************************/
375 
400  template<typename DATACONTAINER>
401  void
402  ElementMatrix(const DATACONTAINER& dc,
403  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
404  double scale_ico = 1.);
405 
406  /******************************************************/
407 
419  template<typename DATACONTAINER>
420  void
421  ElementTimeMatrix(const DATACONTAINER& dc,
422  dealii::FullMatrix<double> &local_entry_matrix);
423 
424  /******************************************************/
425 
435  template<typename DATACONTAINER>
436  void
437  ElementTimeMatrixExplicit(const DATACONTAINER& dc,
438  dealii::FullMatrix<double> &local_entry_matrix);
439 
440  /******************************************************/
441 
448  template<typename FACEDATACONTAINER>
449  void
450  FaceEquation(const FACEDATACONTAINER& dc,
451  dealii::Vector<double> &local_vector, double scale,
452  double scale_ico);
453 
454  /******************************************************/
461  template<typename FACEDATACONTAINER>
462  void
463  InterfaceEquation(const FACEDATACONTAINER& dc,
464  dealii::Vector<double> &local_vector, double scale,
465  double scale_ico);
466 
467  /******************************************************/
474  template<typename FACEDATACONTAINER>
475  void
476  FaceRhs(const FACEDATACONTAINER& dc,
477  dealii::Vector<double> &local_vector, double scale = 1.);
478 
479  /******************************************************/
480 
487  template<typename FACEDATACONTAINER>
488  void
489  FaceMatrix(const FACEDATACONTAINER& dc,
490  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
491  double scale_ico = 1.);
492 
493  /******************************************************/
500  template<typename FACEDATACONTAINER>
501  void
502  InterfaceMatrix(const FACEDATACONTAINER& dc,
503  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
504  double scale_ico = 1.);
505 
506  /******************************************************/
507 
513  template<typename FACEDATACONTAINER>
514  void
515  BoundaryEquation(const FACEDATACONTAINER& dc,
516  dealii::Vector<double> &local_vector, double scale,
517  double scale_ico);
518 
519  /******************************************************/
526  template<typename FACEDATACONTAINER>
527  void
528  BoundaryRhs(const FACEDATACONTAINER& dc,
529  dealii::Vector<double> &local_vector, double scale = 1.);
530 
531  /******************************************************/
532 
539  template<typename FACEDATACONTAINER>
540  void
541  BoundaryMatrix(const FACEDATACONTAINER& dc,
542  dealii::FullMatrix<double> &local_matrix, double scale = 1.,
543  double scale_ico = 1.);
544 
545  /******************************************************/
546 
547  const FE<dealdim, dealdim>&
548  GetFESystem() const;
549 
550  /******************************************************/
557  bool
558  HasFaces() const;
559  /******************************************************/
563  bool
564  HasPoints() const;
565 
566  /******************************************************/
573  bool
574  HasInterfaces() const;
575 
576  /******************************************************/
577 
578  dealii::UpdateFlags
579  GetUpdateFlags() const;
580 
581  /******************************************************/
582 
583  dealii::UpdateFlags
584  GetFaceUpdateFlags() const;
585 
586  /******************************************************/
587 
588  void
589  SetDirichletBoundaryColors(unsigned int color,
590  const std::vector<bool>& comp_mask, const DD* values);
591 
592  /******************************************************/
593 
594  const std::vector<unsigned int>&
595  GetDirichletColors() const;
596  const std::vector<bool>&
597  GetDirichletCompMask(unsigned int color) const;
598 
599  /******************************************************/
600 
601  const dealii::Function<dealdim> &
602  GetDirichletValues(unsigned int color,
603  const std::map<std::string, const dealii::Vector<double>*> &param_values,
604  const std::map<std::string, const VECTOR*> &domain_values) const;
605 
606  /******************************************************/
607 
608  void
609  SetInitialValues(const dealii::Function<dealdim>* values)
610  {
611  _initial_values = values;
612  }
613  const dealii::Function<dealdim>&
615  {
616  return *_initial_values;
617  }
618 
619  /******************************************************/
620 
621  void
622  SetBoundaryEquationColors(unsigned int color);
623  const std::vector<unsigned int>&
625  void
626  SetBoundaryFunctionalColors(unsigned int color);
627  const std::vector<unsigned int>&
629 
630  /******************************************************/
631 
640  void
643  VECTOR, dealdim>* F)
644  {
645  _aux_functionals.push_back(F);
646  if (_functional_position.find(F->GetName())
647  != _functional_position.end())
648  {
649  throw DOpEException(
650  "You cant add two functionals with the same name.",
651  "PDEProblemContainer::AddFunctional");
652  }
653  _functional_position[F->GetName()] = _aux_functionals.size() - 1;
654  //-1 because we are in the pde case and have therefore no cost functional
655  }
656 
657  /******************************************************/
658 
669  void
670  SetFunctionalForErrorEstimation(std::string functional_name)
671  {
672  bool found = false;
673  //we go through all aux functionals.
674  for (unsigned int i = 0; i < this->GetNFunctionals(); i++)
675  {
676  if (_aux_functionals[i]->GetName() == functional_name)
677  {
678  //if the names match, we have found our functional.
679  found = true;
680  _functional_for_ee_num = i;
681  break;
682  }
683  }
684  //If we have not found a functional with the given name,
685  //we throw an error.
686  if (!found)
687  {
688  throw DOpEException(
689  "Can't find functional " + functional_name
690  + " in _aux_functionals",
691  "PDEProblemContainer::SetFunctionalForErrorEstimation");
692  }
693  }
694 
695  /******************************************************/
696 
697  unsigned int
699  {
700  return _aux_functionals.size();
701  }
702 
703  /******************************************************/
704 
705  unsigned int
706  GetStateNBlocks() const;
707 
708  /******************************************************/
709 
710  unsigned int
711  GetNBlocks() const;
712 
713  /******************************************************/
714 
715  unsigned int
716  GetDoFsPerBlock(unsigned int b) const;
717 
718  /******************************************************/
719 
720  const std::vector<unsigned int>&
721  GetDoFsPerBlock() const;
722 
723  /******************************************************/
724 
725  const dealii::ConstraintMatrix&
726  GetDoFConstraints() const;
727 
728  /******************************************************/
729 
730  std::string
731  GetDoFType() const;
732  std::string
733  GetFunctionalType() const;
734  std::string
735  GetFunctionalName() const;
736 
737  /******************************************************/
738 
739  bool
740  NeedTimeFunctional() const;
741 
742  /******************************************************/
743 
746  {
747  assert(_ExceptionHandler);
748  return _ExceptionHandler;
749  }
750 
751  /******************************************************/
752 
755  {
756  assert(_OutputHandler);
757  return _OutputHandler;
758  }
759 
760  /******************************************************/
761 
762  void
763  SetTime(double time, const TimeIterator& interval);
764 
765  /******************************************************/
766 
769  {
770  return _STH;
771  }
772 
773  /******************************************************/
774 
777  {
778  return _STH;
779  }
780 
781  /******************************************************/
782 
783  void
784  ComputeSparsityPattern(SPARSITYPATTERN & sparsity) const;
785 
786  /******************************************************/
787  /*
788  * Experimental status:
789  * Needed for MG prec
790  */
791  void
792  ComputeMGSparsityPattern(dealii::MGLevelObject<dealii::BlockSparsityPattern> & mg_sparsity_patterns,
793  unsigned int n_levels) const;
794 
795  /******************************************************/
796  /*
797  * Experimental status:
798  * Needed for MG prec
799  */
800  void
801  ComputeMGSparsityPattern(dealii::MGLevelObject<dealii::SparsityPattern> & mg_sparsity_patterns,
802  unsigned int n_levels) const;
803 
804  /******************************************************/
805 
806  template<typename INTEGRATOR>
807  void
808  AddAuxiliaryToIntegrator(INTEGRATOR& /*integrator*/)
809  {
810 
811  }
812 
813  /******************************************************/
814 
815  template<typename INTEGRATOR>
816  void
817  DeleteAuxiliaryFromIntegrator(INTEGRATOR& /*integrator*/)
818  {
819 
820  }
821 
822  /******************************************************/
823 
824  unsigned int
826  {
827  return this->GetPDE().GetStateNBlocks();
828  }
829 
830  /******************************************************/
831 
832  std::vector<unsigned int>&
834  {
835  return this->GetPDE().GetStateBlockComponent();
836  }
837 
838  /******************************************************/
839 
846  const std::map<std::string, unsigned int>&
848  {
849  return _functional_position;
850  }
851 
852  /******************************************************/
853  private:
854  DOpEExceptionHandler<VECTOR>* _ExceptionHandler;
855  DOpEOutputHandler<VECTOR>* _OutputHandler;
856 
857  std::string _algo_type;
858 
859  std::vector<
861  VECTOR, dealdim>*> _aux_functionals;
862  std::map<std::string, unsigned int> _functional_position;
863 
864  unsigned int _functional_for_ee_num;
866 
867  std::vector<unsigned int> _dirichlet_colors;
868  std::vector<std::vector<bool> > _dirichlet_comps;
869  std::vector<PrimalDirichletData<DD, VECTOR, dealdim>*> _primal_dirichlet_values;
870  const dealii::Function<dealdim>* _zero_dirichlet_values;
871 
872  const dealii::Function<dealdim>* _initial_values;
873 
874  std::vector<unsigned int> _state_boundary_equation_colors;
875  std::vector<unsigned int> _adjoint_boundary_equation_colors;
876 
877  std::vector<unsigned int> _boundary_functional_colors;
878 
879  StateProblem<
880  PDEProblemContainer<PDE, DD, SPARSITYPATTERN, VECTOR, dealdim, FE,
881  DH>, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim>* _state_problem;
882 
883  friend class StateProblem<
884  PDEProblemContainer<PDE, DD, SPARSITYPATTERN, VECTOR, dealdim, FE,
885  DH>, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim> ;
886  };
887  /******************************************************/
888 
889  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
890  int dealdim, template<int, int> class FE, template<int, int> class DH>
892  PDE& pde,
894  ProblemContainerInternal<PDE>(pde), _STH(&STH), _state_problem(NULL)
895  {
896  _ExceptionHandler = NULL;
897  _OutputHandler = NULL;
898  _zero_dirichlet_values = new ZeroFunction<dealdim>(
899  this->GetPDE().GetStateNComponents());
900  _algo_type = "";
901  _functional_for_ee_num = dealii::numbers::invalid_unsigned_int;
902  }
903 
904  /******************************************************/
905 
906  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
907  int dealdim, template<int, int> class FE, template<int, int> class DH>
909  {
910  if (_zero_dirichlet_values != NULL)
911  {
912  delete _zero_dirichlet_values;
913  }
914  for (unsigned int i = 0; i < _primal_dirichlet_values.size(); i++)
915  {
916  if (_primal_dirichlet_values[i] != NULL)
917  delete _primal_dirichlet_values[i];
918  }
919  if (_state_problem != NULL)
920  {
921  delete _state_problem;
922  }
923  }
924 
925  /******************************************************/
926 
927  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
928  int dealdim, template<int, int> class FE, template<int, int> class DH>
929  void
931  std::string algo_type)
932  {
933  if (_state_problem != NULL)
934  {
935  delete _state_problem;
936  _state_problem = NULL;
937  }
938 
939  if (_algo_type != algo_type && _algo_type != "")
940  {
941  throw DOpEException("Conflicting Algorithms!",
942  "PDEProblemContainer::ReInit");
943  }
944  else
945  {
946  _algo_type = algo_type;
947  this->SetTypeInternal("");
948 
949  if (_algo_type == "reduced")
950  {
951  GetSpaceTimeHandler()->ReInit(this->GetPDE().GetStateNBlocks(),
952  this->GetPDE().GetStateBlockComponent());
953  }
954  else
955  {
956  throw DOpEException("Unknown Algorithm " + _algo_type,
957  "PDEProblemContainer::ReInit");
958  }
959  }
960  }
961 
962  /******************************************************/
963 
964  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
965  int dealdim, template<int, int> class FE, template<int, int> class DH>
966  void
968  std::string type, unsigned int num)
969  {
970  if (this->GetType() != type || this->GetTypeNum() != num)
971  {
972  this->SetTypeInternal(type);
973  this->SetTypeNumInternal(num);
974  this->GetPDE().SetProblemType(type);
975  if (_functional_for_ee_num != dealii::numbers::invalid_unsigned_int)
976  _aux_functionals[_functional_for_ee_num]->SetProblemType(type);
977  }
978  //Nothing to do.
979  }
980 
981  /******************************************************/
982 
983  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
984  int dealdim, template<int, int> class FE, template<int, int> class DH>
985  template<typename DATACONTAINER>
986  double
988  const DATACONTAINER& edc)
989  {
990 
991  if (this->GetType() == "cost_functional")
992  {
993  return 0;
994  }
995  else if (this->GetType() == "aux_functional")
996  {
997  // state values in quadrature points
998  return _aux_functionals[this->GetTypeNum()]->ElementValue(edc);
999  }
1000  else if (this->GetType() == "error_evaluation")
1001  {
1002  return _aux_functionals[_functional_for_ee_num]->ElementValue(edc);
1003  }
1004  else
1005  {
1006  throw DOpEException("Not implemented",
1007  "PDEProblemContainer::ElementFunctional");
1008  }
1009  }
1010 
1011  /******************************************************/
1012  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1013  int dealdim, template<int, int> class FE, template<int, int> class DH>
1014  double
1016  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1017  const std::map<std::string, const VECTOR*> &domain_values)
1018  {
1019  if (this->GetType() == "cost_functional")
1020  {
1021  return 0.;
1022  } //endif cost_functional
1023  else if (this->GetType() == "aux_functional")
1024  {
1025  // state values in quadrature points
1026  return _aux_functionals[this->GetTypeNum()]->PointValue(
1027  this->GetSpaceTimeHandler()->GetStateDoFHandler(),
1028  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1029  domain_values);
1030 
1031  } //endif aux_functional
1032  else if (this->GetType() == "error_evaluation")
1033  {
1034  return _aux_functionals[_functional_for_ee_num]->PointValue(
1035  this->GetSpaceTimeHandler()->GetStateDoFHandler(),
1036  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1037  domain_values);
1038  }
1039  else
1040  {
1041  throw DOpEException("Not implemented",
1042  "PDEProblemContainer::PointFunctional");
1043  }
1044  }
1045 
1046  /******************************************************/
1047 
1048  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1049  int dealdim, template<int, int> class FE, template<int, int> class DH>
1050  template<typename FACEDATACONTAINER>
1051  double
1053  const FACEDATACONTAINER& fdc)
1054  {
1055  if (this->GetType() == "cost_functional")
1056  {
1057  // state values in quadrature points
1058  return 0.;
1059  }
1060  else if (this->GetType() == "aux_functional")
1061  {
1062  // state values in quadrature points
1063  return _aux_functionals[this->GetTypeNum()]->BoundaryValue(fdc);
1064  }
1065  else if (this->GetType() == "error_evaluation")
1066  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1067  {
1068  return _aux_functionals[_functional_for_ee_num]->BoundaryValue(fdc);
1069  }
1070  else
1071  {
1072  throw DOpEException("Not implemented",
1073  "PDEProblemContainer::BoundaryFunctional");
1074  }
1075  }
1076 
1077  /******************************************************/
1078 
1079  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1080  int dealdim, template<int, int> class FE, template<int, int> class DH>
1081  template<typename FACEDATACONTAINER>
1082  double
1084  const FACEDATACONTAINER& fdc)
1085  {
1086  if (this->GetType() == "cost_functional")
1087  {
1088  // state values in quadrature points
1089  return 0.;
1090  }
1091  else if (this->GetType() == "aux_functional")
1092  {
1093  // state values in quadrature points
1094  return _aux_functionals[this->GetTypeNum()]->FaceValue(fdc);
1095  }
1096  else if (this->GetType() == "error_evaluation")
1097  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1098  {
1099  return _aux_functionals[_functional_for_ee_num]->FaceValue(fdc);
1100  }
1101  else
1102  {
1103  throw DOpEException("Not implemented",
1104  "PDEProblemContainer::FaceFunctional");
1105  }
1106  }
1107 
1108  /******************************************************/
1109 
1110  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1111  int dealdim, template<int, int> class FE, template<int, int> class DH>
1112  double
1114  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1115  const std::map<std::string, const VECTOR*> &domain_values)
1116  {
1117  if (this->GetType() == "cost_functional")
1118  {
1119  // state values in quadrature points
1120  return 0.;
1121  }
1122  else if (this->GetType() == "aux_functional")
1123  {
1124  // state values in quadrature points
1125  return _aux_functionals[this->GetTypeNum()]->AlgebraicValue(
1126  param_values, domain_values);
1127  }
1128  else if (this->GetType() == "error_evaluation")
1129  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1130  {
1131  return _aux_functionals[_functional_for_ee_num]->AlgebraicValue(
1132  param_values, domain_values);
1133  }
1134  else
1135  {
1136  throw DOpEException("Not implemented",
1137  "PDEProblemContainer::AlgebraicFunctional");
1138  }
1139  }
1140 
1141  /******************************************************/
1142 
1143  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1144  int dealdim, template<int, int> class FE, template<int, int> class DH>
1145  template<typename DATACONTAINER>
1146  void
1148  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1149  double scale, double scale_ico)
1150  {
1151  if (this->GetType() == "state")
1152  {
1153  // state values in quadrature points
1154  this->GetPDE().ElementEquation(edc, local_vector, scale, scale_ico);
1155  }
1156  else if (this->GetType() == "adjoint_for_ee")
1157  {
1158  // state values in quadrature points
1159  this->GetPDE().ElementEquation_U(edc, local_vector, scale,
1160  scale_ico);
1161  }
1162  else
1163  {
1164  throw DOpEException("Not implemented",
1165  "PDEProblemContainer::ElementEquation");
1166  }
1167  }
1168 
1169  /******************************************************/
1170 
1171  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1172  int dealdim, template<int, int> class FE, template<int, int> class DH>
1173  template<typename DATACONTAINER>
1174  void
1176  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1177  double scale)
1178  {
1179 
1180  if (this->GetType() == "state")
1181  {
1182  // state values in quadrature points
1183  this->GetPDE().ElementTimeEquation(edc, local_vector, scale);
1184  }
1185  else if (this->GetType() == "adjoint_for_ee")
1186  {
1187  throw DOpEException("Not implemented",
1188  "OptProblem::ElementTimeEquation");
1189  }
1190  else
1191  {
1192  throw DOpEException("Not implemented",
1193  "PDEProblemContainer::ElementTimeEquation");
1194  }
1195  }
1196 
1197  /******************************************************/
1198 
1199  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1200  int dealdim, template<int, int> class FE, template<int, int> class DH>
1201  template<typename DATACONTAINER>
1202  void
1204  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1205  double scale)
1206  {
1207 
1208  if (this->GetType() == "state")
1209  {
1210  // state values in quadrature points
1211  this->GetPDE().ElementTimeEquationExplicit(edc, local_vector,
1212  scale);
1213  }
1214  else if (this->GetType() == "adjoint_for_ee")
1215  {
1216  throw DOpEException("Not implemented",
1217  "OptProblem::ElementTimeEquation");
1218  }
1219  else
1220  {
1221  throw DOpEException("Not implemented",
1222  "PDEProblemContainer::ElementTimeEquation");
1223  }
1224  }
1225 
1226  /******************************************************/
1227 
1228  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1229  int dealdim, template<int, int> class FE, template<int, int> class DH>
1230  template<typename FACEDATACONTAINER>
1231  void
1233  const FACEDATACONTAINER& fdc,
1234  dealii::Vector<double> &local_vector, double scale,
1235  double scale_ico)
1236  {
1237  if (this->GetType() == "state")
1238  {
1239  // state values in quadrature points
1240  this->GetPDE().FaceEquation(fdc, local_vector, scale, scale_ico);
1241  }
1242  else if (this->GetType() == "adjoint_for_ee")
1243  {
1244  // state values in quadrature points
1245  this->GetPDE().FaceEquation_U(fdc, local_vector, scale,
1246  scale_ico);
1247  }
1248  else
1249  {
1250  throw DOpEException("Not implemented",
1251  "PDEProblemContainer::FaceEquation");
1252  }
1253  }
1254 
1255  /******************************************************/
1256 
1257  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1258  int dealdim, template<int, int> class FE, template<int, int> class DH>
1259  template<typename FACEDATACONTAINER>
1260  void
1262  const FACEDATACONTAINER& fdc,
1263  dealii::Vector<double> &local_vector, double scale,
1264  double scale_ico)
1265  {
1266  if (this->GetType() == "state")
1267  {
1268  // state values in quadrature points
1269  this->GetPDE().InterfaceEquation(fdc, local_vector, scale,
1270  scale_ico);
1271  }
1272  else if (this->GetType() == "adjoint_for_ee")
1273  {
1274  // state values in quadrature points
1275  this->GetPDE().InterfaceEquation_U(fdc, local_vector, scale,
1276  scale_ico);
1277  }
1278  else
1279  {
1280  throw DOpEException("Not implemented",
1281  "PDEProblemContainer::InterfaceEquation");
1282  }
1283  }
1284  /******************************************************/
1285 
1286  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1287  int dealdim, template<int, int> class FE, template<int, int> class DH>
1288  template<typename FACEDATACONTAINER>
1289  void
1291  const FACEDATACONTAINER& fdc,
1292  dealii::Vector<double> &local_vector, double scale,
1293  double scale_ico)
1294  {
1295  if (this->GetType() == "state")
1296  {
1297  // state values in quadrature points
1298  this->GetPDE().BoundaryEquation(fdc, local_vector, scale,
1299  scale_ico);
1300  }
1301  else if (this->GetType() == "adjoint_for_ee")
1302  {
1303  // state values in quadrature points
1304  this->GetPDE().BoundaryEquation_U(fdc, local_vector, scale,
1305  scale_ico);
1306  }
1307  else
1308  {
1309  throw DOpEException("Not implemented",
1310  "PDEProblemContainer::ElementBoundaryEquation");
1311  }
1312  }
1313 
1314  /******************************************************/
1315 
1316  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1317  int dealdim, template<int, int> class FE, template<int, int> class DH>
1318  template<typename DATACONTAINER>
1319  void
1321  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1322  double scale)
1323  {
1324 
1325  if (this->GetType() == "state")
1326  {
1327  // state values in quadrature points
1328  this->GetPDE().ElementRightHandSide(edc, local_vector, scale);
1329  }
1330  else if (this->GetType() == "adjoint_for_ee")
1331  {
1332  //TODO currently, pointvalue is not working for dual error estimation!
1333  //values of the derivative of the functional for error estimation.
1334  //Check, if we have to evaluate an integral over a domain.
1335  if (_aux_functionals[_functional_for_ee_num]->GetType().find("domain")
1336  != std::string::npos)
1337  _aux_functionals[_functional_for_ee_num]->ElementValue_U(edc,
1338  local_vector, scale);
1339  }
1340  else
1341  {
1342  throw DOpEException("Not implemented",
1343  "PDEProblemContainer::ElementRhs");
1344  }
1345  }
1346 
1347  /******************************************************/
1348 
1349  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1350  int dealdim, template<int, int> class FE, template<int, int> class DH>
1351  void
1353  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1354  const std::map<std::string, const VECTOR*> &domain_values,
1355  VECTOR& rhs_vector, double scale)
1356  {
1357  if (this->GetType() == "adjoint_for_ee")
1358  {
1359  //values of the derivative of the functional for error estimation
1360  if (_aux_functionals[_functional_for_ee_num]->GetType().find("point")
1361  != std::string::npos)
1362  _aux_functionals[_functional_for_ee_num]->PointValue_U(
1363  this->GetSpaceTimeHandler()->GetStateDoFHandler(),
1364  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1365  domain_values, rhs_vector, scale);
1366  }
1367  else
1368  {
1369  throw DOpEException("Not implemented", "OptProblem::ElementRhs");
1370  }
1371  }
1372 
1373  /******************************************************/
1374 
1375  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1376  int dealdim, template<int, int> class FE, template<int, int> class DH>
1377  template<typename FACEDATACONTAINER>
1378  void
1380  const FACEDATACONTAINER& fdc,
1381  dealii::Vector<double> &local_vector, double scale)
1382  {
1383  if (this->GetType() == "state")
1384  {
1385  // state values in face quadrature points
1386  this->GetPDE().FaceRightHandSide(fdc, local_vector, scale);
1387  }
1388  else if (this->GetType() == "adjoint_for_ee")
1389  {
1390  //values of the derivative of the functional for error estimation
1391  if (_aux_functionals[_functional_for_ee_num]->GetType().find("face")
1392  != std::string::npos)
1393  _aux_functionals[_functional_for_ee_num]->FaceValue_U(fdc,
1394  local_vector, scale);
1395  }
1396  else
1397  {
1398  throw DOpEException("Not implemented",
1399  "PDEProblemContainer::ElementFaceRhs");
1400  }
1401  }
1402 
1403  /******************************************************/
1404 
1405  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1406  int dealdim, template<int, int> class FE, template<int, int> class DH>
1407  template<typename FACEDATACONTAINER>
1408  void
1410  const FACEDATACONTAINER& fdc,
1411  dealii::Vector<double> &local_vector, double scale)
1412  {
1413  if (this->GetType() == "state")
1414  {
1415  // state values in face quadrature points
1416  this->GetPDE().BoundaryRightHandSide(fdc, local_vector, scale);
1417  }
1418  else if (this->GetType() == "adjoint_for_ee")
1419  {
1420  //values of the derivative of the functional for error estimation
1421  if (_aux_functionals[_functional_for_ee_num]->GetType().find(
1422  "boundary") != std::string::npos)
1423  _aux_functionals[_functional_for_ee_num]->BoundaryValue_U(fdc,
1424  local_vector, scale);
1425  }
1426  else
1427  {
1428  throw DOpEException("Not implemented",
1429  "PDEProblemContainer::ElementBoundaryRhs");
1430  }
1431  }
1432 
1433  /******************************************************/
1434 
1435  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1436  int dealdim, template<int, int> class FE, template<int, int> class DH>
1437  template<typename DATACONTAINER>
1438  void
1440  const DATACONTAINER& edc,
1441  dealii::FullMatrix<double> &local_entry_matrix, double scale,
1442  double scale_ico)
1443  {
1444 
1445  if (this->GetType() == "state")
1446  {
1447  // state values in quadrature points
1448  this->GetPDE().ElementMatrix(edc, local_entry_matrix, scale, scale_ico);
1449  }
1450  else if (this->GetType() == "adjoint_for_ee")
1451  {
1452  // state values in quadrature points
1453  this->GetPDE().ElementMatrix_T(edc, local_entry_matrix, scale,
1454  scale_ico);
1455  }
1456  else
1457  {
1458  throw DOpEException("Not implemented",
1459  "PDEProblemContainer::ElementMatrix");
1460  }
1461 
1462  }
1463 
1464  /******************************************************/
1465 
1466  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1467  int dealdim, template<int, int> class FE, template<int, int> class DH>
1468  template<typename DATACONTAINER>
1469  void
1471  const DATACONTAINER& edc, FullMatrix<double> &local_entry_matrix)
1472  {
1473 
1474  if (this->GetType() == "state")
1475  {
1476  // state values in quadrature points
1477  this->GetPDE().ElementTimeMatrix(edc, local_entry_matrix);
1478  }
1479  else if (this->GetType() == "dual_for_ee")
1480  {
1481  throw DOpEException("Not implemented",
1482  "PDEProblemContainer::ElementTimeMatrix");
1483  }
1484  else
1485  {
1486  throw DOpEException("Not implemented",
1487  "PDEProblemContainer::ElementTimeMatrix");
1488  }
1489 
1490  }
1491 
1492  /******************************************************/
1493 
1494  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1495  int dealdim, template<int, int> class FE, template<int, int> class DH>
1496  template<typename DATACONTAINER>
1497  void
1499  const DATACONTAINER& edc,
1500  dealii::FullMatrix<double> &local_entry_matrix)
1501  {
1502 
1503  if (this->GetType() == "state")
1504  {
1505  // state values in quadrature points
1506  this->GetPDE().ElementTimeMatrixExplicit(edc, local_entry_matrix);
1507  }
1508  else if (this->GetType() == "dual_for_ee")
1509  {
1510  throw DOpEException("Not implemented",
1511  "PDEProblemContainer::ElementTimeMatrix");
1512  }
1513  else
1514  {
1515  throw DOpEException("Not implemented",
1516  "PDEProblemContainer::ElementTimeMatrix");
1517  }
1518 
1519  }
1520 
1521  /******************************************************/
1522 
1523  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1524  int dealdim, template<int, int> class FE, template<int, int> class DH>
1525  template<typename FACEDATACONTAINER>
1526  void
1528  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
1529  double scale, double scale_ico)
1530  {
1531  if (this->GetType() == "state")
1532  {
1533  // state values in face quadrature points
1534  this->GetPDE().FaceMatrix(fdc, local_entry_matrix, scale, scale_ico);
1535  }
1536  else if (this->GetType() == "adjoint_for_ee")
1537  {
1538  this->GetPDE().FaceMatrix_T(fdc, local_entry_matrix, scale,
1539  scale_ico);
1540  }
1541  else
1542  {
1543  throw DOpEException("Not implemented",
1544  "PDEProblemContainer::NewtonFaceMatrix");
1545  }
1546 
1547  }
1548 
1549  /******************************************************/
1550 
1551  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1552  int dealdim, template<int, int> class FE, template<int, int> class DH>
1553  template<typename FACEDATACONTAINER>
1554  void
1556  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
1557  double scale, double scale_ico)
1558  {
1559  if (this->GetType() == "state")
1560  {
1561  // state values in face quadrature points
1562  this->GetPDE().InterfaceMatrix(fdc, local_entry_matrix, scale,
1563  scale_ico);
1564  }
1565  else if (this->GetType() == "adjoint_for_ee")
1566  {
1567  this->GetPDE().InterfaceMatrix_T(fdc, local_entry_matrix, scale,
1568  scale_ico);
1569  }
1570  else
1571  {
1572  throw DOpEException("Not implemented",
1573  "PDEProblemContainer::NewtonInterfaceMatrix");
1574  }
1575  }
1576 
1577  /******************************************************/
1578 
1579  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1580  int dealdim, template<int, int> class FE, template<int, int> class DH>
1581  template<typename FACEDATACONTAINER>
1582  void
1584  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_matrix,
1585  double scale, double scale_ico)
1586  {
1587  if (this->GetType() == "state")
1588  {
1589  // state values in face quadrature points
1590  this->GetPDE().BoundaryMatrix(fdc, local_matrix, scale,
1591  scale_ico);
1592  }
1593  else if (this->GetType() == "adjoint_for_ee")
1594  {
1595  // state values in quadrature points
1596  this->GetPDE().BoundaryMatrix_T(fdc, local_matrix, scale,
1597  scale_ico);
1598  }
1599  else
1600  {
1601  throw DOpEException("Not implemented",
1602  "PDEProblemContainer::ElementBoundaryMatrix");
1603  }
1604 
1605  }
1606 
1607  /******************************************************/
1608 
1609  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1610  int dealdim, template<int, int> class FE, template<int, int> class DH>
1611  std::string
1613  {
1614  if (this->GetType() == "state" || this->GetType() == "adjoint_for_ee"
1615  || this->GetType() == "error_evaluation")
1616  {
1617  return "state";
1618  }
1619  else
1620  {
1621  throw DOpEException("Unknown Type:" + this->GetType(),
1622  "PDEProblemContainer::GetDoFType");
1623  }
1624  }
1625 
1626  /******************************************************/
1627 
1628  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1629  int dealdim, template<int, int> class FE, template<int, int> class DH>
1630  const FE<dealdim, dealdim>&
1632  {
1633  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee")
1634  {
1635  return this->GetSpaceTimeHandler()->GetFESystem("state");
1636  }
1637  else
1638  {
1639  throw DOpEException("Unknown Type:" + this->GetType(),
1640  "PDEProblemContainer::GetFESystem");
1641  }
1642  }
1643 
1644  /******************************************************/
1645 
1646  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1647  int dealdim, template<int, int> class FE, template<int, int> class DH>
1648  UpdateFlags
1650  {
1651 
1652  UpdateFlags r;
1653  if (this->GetType().find("aux_functional") != std::string::npos)
1654  {
1655  r = _aux_functionals[this->GetTypeNum()]->GetUpdateFlags();
1656  }
1657  else
1658  {
1659  r = this->GetPDE().GetUpdateFlags();
1660  if (this->GetType() == "adjoint_for_ee"
1661  || this->GetType() == "error_evaluation")
1662  {
1663  if (_functional_for_ee_num != dealii::numbers::invalid_unsigned_int)
1664  r = r | _aux_functionals[_functional_for_ee_num]->GetUpdateFlags();
1665  }
1666  }
1667  return r | update_JxW_values;
1668  }
1669 
1670  /******************************************************/
1671 
1672  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1673  int dealdim, template<int, int> class FE, template<int, int> class DH>
1674  UpdateFlags
1676  {
1677  UpdateFlags r;
1678  if (this->GetType().find("aux_functional") != std::string::npos)
1679  {
1680  r = _aux_functionals[this->GetTypeNum()]->GetFaceUpdateFlags();
1681  }
1682  else
1683  {
1684  r = this->GetPDE().GetFaceUpdateFlags();
1685  if (this->GetType() == "adjoint_for_ee"
1686  || this->GetType() == "error_evaluation")
1687  {
1688  if (_functional_for_ee_num != dealii::numbers::invalid_unsigned_int)
1689  r =
1690  r
1691  | _aux_functionals[_functional_for_ee_num]->GetFaceUpdateFlags();
1692  }
1693  }
1694  return r | update_JxW_values;
1695  }
1696 
1697  /******************************************************/
1698 
1699  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1700  int dealdim, template<int, int> class FE, template<int, int> class DH>
1701  std::string
1703  {
1704  if (this->GetType() == "aux_functional")
1705  {
1706  return _aux_functionals[this->GetTypeNum()]->GetType();
1707  }
1708  else if (this->GetType() == "error_evaluation")
1709  {
1710  return _aux_functionals[_functional_for_ee_num]->GetType();
1711  }
1712  return "none";
1713  }
1714 
1715  /******************************************************/
1716 
1717  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1718  int dealdim, template<int, int> class FE, template<int, int> class DH>
1719  std::string
1721  {
1722  if (this->GetType() == "aux_functional")
1723  {
1724  return _aux_functionals[this->GetTypeNum()]->GetName();
1725  }
1726  else if (this->GetType() == "error_evaluation")
1727  {
1728  return _aux_functionals[_functional_for_ee_num]->GetName();
1729  }
1730  return "";
1731  }
1732 
1733  /******************************************************/
1734 
1735  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1736  int dealdim, template<int, int> class FE, template<int, int> class DH>
1737  void
1739  double time, const TimeIterator& interval)
1740  {
1741  GetSpaceTimeHandler()->SetInterval(interval);
1742 
1743  { //Zeit an Dirichlet Werte uebermitteln
1744  for (unsigned int i = 0; i < _primal_dirichlet_values.size(); i++)
1745  _primal_dirichlet_values[i]->SetTime(time);
1746  for (unsigned int i = 0; i < _aux_functionals.size(); i++)
1747  _aux_functionals[i]->SetTime(time);
1748  //PDE
1749  this->GetPDE().SetTime(time);
1750  }
1751  }
1752 
1753  /******************************************************/
1754 
1755  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1756  int dealdim, template<int, int> class FE, template<int, int> class DH>
1757  void
1759  SPARSITYPATTERN & sparsity) const
1760  {
1761  if (this->GetType() == "state" || this->GetType() == "adjoint_for_ee")
1762  {
1763  this->GetSpaceTimeHandler()->ComputeStateSparsityPattern(sparsity);
1764  }
1765  else
1766  {
1767  throw DOpEException("Unknown type " + this->GetType(),
1768  "PDEProblemContainer::ComputeSparsityPattern");
1769  }
1770  }
1771 
1772  /******************************************************/
1773 
1774  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1775  int dealdim, template<int, int> class FE, template<int, int> class DH>
1776  void
1778  dealii::MGLevelObject<dealii::BlockSparsityPattern> & mg_sparsity_patterns,
1779  unsigned int n_levels) const
1780  {
1781  if (this->GetType() == "state")
1782  {
1783  this->GetSpaceTimeHandler()->ComputeMGStateSparsityPattern(mg_sparsity_patterns, n_levels);
1784  }
1785  else
1786  {
1787  throw DOpEException("Unknown type " + this->GetType(),
1788  "PDEProblemContainer::ComputeMGSparsityPattern");
1789  }
1790  }
1791 
1792  /******************************************************/
1793 
1794  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1795  int dealdim, template<int, int> class FE, template<int, int> class DH>
1796  void
1798  dealii::MGLevelObject<dealii::SparsityPattern> & mg_sparsity_patterns,
1799  unsigned int n_levels) const
1800  {
1801  if (this->GetType() == "state")
1802  {
1803  this->GetSpaceTimeHandler()->ComputeMGStateSparsityPattern(mg_sparsity_patterns, n_levels);
1804  }
1805  else
1806  {
1807  throw DOpEException("Unknown type " + this->GetType(),
1808  "PDEProblemContainer::ComputeMGSparsityPattern");
1809  }
1810  }
1811 
1812 
1813 
1814 
1815  /******************************************************/
1816 
1817  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1818  int dealdim, template<int, int> class FE, template<int, int> class DH>
1819  bool
1821  {
1822  if (this->GetType().find("aux_functional") != std::string::npos)
1823  {
1824  return _aux_functionals[this->GetTypeNum()]->HasFaces();
1825  }
1826  else
1827  {
1828  if ((this->GetType() == "state"))
1829  {
1830  return this->GetPDE().HasFaces();
1831  }
1832  else if (this->GetType() == "adjoint_for_ee")
1833  {
1834  return this->GetPDE().HasFaces()
1835  || _aux_functionals[_functional_for_ee_num]->HasFaces();
1836  }
1837  else
1838  {
1839  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
1840  "PDEProblemContainer::HasFaces");
1841  }
1842  }
1843  }
1844 
1845  /******************************************************/
1846 
1847  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1848  int dealdim, template<int, int> class FE, template<int, int> class DH>
1849  bool
1851  {
1852  if ((this->GetType() == "state") || this->GetType() == "aux_functional")
1853  {
1854  //We dont need PointRhs in these cases
1855  return false;
1856  }
1857  else if (this->GetType() == "adjoint_for_ee")
1858  {
1859  return _aux_functionals[_functional_for_ee_num]->HasPoints();
1860  }
1861  else
1862  {
1863  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
1864  "PDEProblemContainer::HasFaces");
1865  }
1866 
1867  }
1868 
1869  /******************************************************/
1870 
1871  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1872  int dealdim, template<int, int> class FE, template<int, int> class DH>
1873  bool
1875  {
1876  if (this->GetType().find("aux_functional") != std::string::npos)
1877  {
1878  return false;
1879  }
1880  else
1881  {
1882  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee")
1883  {
1884  return this->GetPDE().HasInterfaces();
1885  }
1886  else
1887  {
1888  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
1889  "PDEProblemContainer::HasFaces");
1890  }
1891  }
1892  }
1893 
1894  /******************************************************/
1895 
1896  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1897  int dealdim, template<int, int> class FE, template<int, int> class DH>
1898  void
1900  unsigned int color, const std::vector<bool>& comp_mask,
1901  const DD* values)
1902  {
1903  assert(values);
1904  assert(values->n_components() == comp_mask.size());
1905 
1906  unsigned int comp = _dirichlet_colors.size();
1907  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
1908  {
1909  if (_dirichlet_colors[i] == color)
1910  {
1911  comp = i;
1912  break;
1913  }
1914  }
1915  if (comp != _dirichlet_colors.size())
1916  {
1917  std::stringstream s;
1918  s << "DirichletColor" << color << " has multiple occurrences !";
1919  throw DOpEException(s.str(),
1920  "PDEProblemContainer::SetDirichletBoundary");
1921  }
1922  _dirichlet_colors.push_back(color);
1923  _dirichlet_comps.push_back(comp_mask);
1925  DD, VECTOR, dealdim>(*values);
1926  _primal_dirichlet_values.push_back(data);
1927  }
1928 
1929  /******************************************************/
1930 
1931  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1932  int dealdim, template<int, int> class FE, template<int, int> class DH>
1933  const std::vector<unsigned int>&
1935  {
1936  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee")
1937  {
1938  return _dirichlet_colors;
1939  }
1940  else
1941  {
1942  throw DOpEException("Unknown Type:" + this->GetType(),
1943  "PDEProblemContainer::GetDirichletColors");
1944  }
1945  }
1946 
1947  /******************************************************/
1948 
1949  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1950  int dealdim, template<int, int> class FE, template<int, int> class DH>
1951  const std::vector<bool>&
1953  unsigned int color) const
1954  {
1955  if ((this->GetType() == "state" || this->GetType() == "adjoint_for_ee"))
1956  {
1957  unsigned int comp = _dirichlet_colors.size();
1958  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
1959  {
1960  if (_dirichlet_colors[i] == color)
1961  {
1962  comp = i;
1963  break;
1964  }
1965  }
1966  if (comp == _dirichlet_colors.size())
1967  {
1968  std::stringstream s;
1969  s << "DirichletColor" << color << " has not been found !";
1970  throw DOpEException(s.str(),
1971  "PDEProblemContainer::GetDirichletCompMask");
1972  }
1973  return _dirichlet_comps[comp];
1974  }
1975  else
1976  {
1977  throw DOpEException("Unknown Type:" + this->GetType(),
1978  "PDEProblemContainer::GetDirichletCompMask");
1979  }
1980  }
1981 
1982  /******************************************************/
1983 
1984  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
1985  int dealdim, template<int, int> class FE, template<int, int> class DH>
1986  const Function<dealdim>&
1988  unsigned int color,
1989  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1990  const std::map<std::string, const VECTOR*> &domain_values) const
1991  {
1992 
1993  unsigned int col = _dirichlet_colors.size();
1994  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee")
1995  {
1996  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
1997  {
1998  if (_dirichlet_colors[i] == color)
1999  {
2000  col = i;
2001  break;
2002  }
2003  }
2004  if (col == _dirichlet_colors.size())
2005  {
2006  std::stringstream s;
2007  s << "DirichletColor" << color << " has not been found !";
2008  throw DOpEException(s.str(),
2009  "PDEProblemContainer::GetDirichletValues");
2010  }
2011  }
2012  else
2013  {
2014  throw DOpEException("Unknown Type:" + this->GetType(),
2015  "PDEProblemContainer::GetDirichletValues");
2016  }
2017 
2018  if (this->GetType() == "state")
2019  {
2020  _primal_dirichlet_values[col]->ReInit(param_values, domain_values,
2021  color);
2022  return *(_primal_dirichlet_values[col]);
2023  }
2024  else if (this->GetType() == "adjoint_for_ee"
2025  || (this->GetType() == "adjoint_hessian"))
2026  {
2027  return *(_zero_dirichlet_values);
2028  }
2029  else
2030  {
2031  throw DOpEException("Unknown Type:" + this->GetType(),
2032  "PDEProblemContainer::GetDirichletValues");
2033  }
2034  }
2035 
2036  /******************************************************/
2037 
2038  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2039  int dealdim, template<int, int> class FE, template<int, int> class DH>
2040  const std::vector<unsigned int>&
2042  {
2043  if (this->GetType() == "state")
2044  {
2045  return _state_boundary_equation_colors;
2046  }
2047  else if (this->GetType() == "adjoint_for_ee")
2048  {
2049  return _adjoint_boundary_equation_colors;
2050  }
2051  else
2052  {
2053  throw DOpEException("Unknown Type:" + this->GetType(),
2054  "PDEProblemContainer::GetBoundaryEquationColors");
2055  }
2056  }
2057 
2058  /******************************************************/
2059 
2060  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2061  int dealdim, template<int, int> class FE, template<int, int> class DH>
2062  void
2064  unsigned int color)
2065  {
2066  { //State Boundary Equation colors are simply inserted
2067  unsigned int comp = _state_boundary_equation_colors.size();
2068  for (unsigned int i = 0; i < _state_boundary_equation_colors.size();
2069  ++i)
2070  {
2071  if (_state_boundary_equation_colors[i] == color)
2072  {
2073  comp = i;
2074  break;
2075  }
2076  }
2077  if (comp != _state_boundary_equation_colors.size())
2078  {
2079  std::stringstream s;
2080  s << "Boundary Equation Color" << color
2081  << " has multiple occurences !";
2082  throw DOpEException(s.str(),
2083  "PDEProblemContainer::SetBoundaryEquationColors");
2084  }
2085  _state_boundary_equation_colors.push_back(color);
2086  }
2087  { //For the adjoint they are added with the boundary functional colors
2088  unsigned int comp = _adjoint_boundary_equation_colors.size();
2089  for (unsigned int i = 0; i < _adjoint_boundary_equation_colors.size();
2090  ++i)
2091  {
2092  if (_adjoint_boundary_equation_colors[i] == color)
2093  {
2094  comp = i;
2095  break;
2096  }
2097  }
2098  if (comp != _adjoint_boundary_equation_colors.size())
2099  {
2100  //Seems this color is already added, however it might have been a functional color
2101  //so we don't do anything.
2102  }
2103  else
2104  {
2105  _adjoint_boundary_equation_colors.push_back(color);
2106  }
2107 
2108  }
2109  }
2110 
2111  /******************************************************/
2112 
2113  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2114  int dealdim, template<int, int> class FE, template<int, int> class DH>
2115  const std::vector<unsigned int>&
2117  {
2118  //FIXME cost_functional?? This is pdeproblemcontainer, we should not have a cost functional! ~cg
2119  if (this->GetType() == "cost_functional"
2120  || this->GetType() == "aux_functional"
2121  || this->GetType() == "error_evaluation")
2122  {
2123  return _boundary_functional_colors;
2124  }
2125  else
2126  {
2127  throw DOpEException("Unknown Type:" + this->GetType(),
2128  "PDEProblemContainer::GetBoundaryFunctionalColors");
2129  }
2130  }
2131 
2132  /******************************************************/
2133 
2134  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2135  int dealdim, template<int, int> class FE, template<int, int> class DH>
2136  void
2138  unsigned int color)
2139  {
2140  { //Boundary Functional colors are simply inserted
2141  unsigned int comp = _boundary_functional_colors.size();
2142  for (unsigned int i = 0; i < _boundary_functional_colors.size(); ++i)
2143  {
2144  if (_boundary_functional_colors[i] == color)
2145  {
2146  comp = i;
2147  break;
2148  }
2149  }
2150  if (comp != _boundary_functional_colors.size())
2151  {
2152  std::stringstream s;
2153  s << "Boundary Functional Color" << color
2154  << " has multiple occurences !";
2155  throw DOpEException(s.str(),
2156  "PDEProblemContainer::SetBoundaryFunctionalColors");
2157  }
2158  _boundary_functional_colors.push_back(color);
2159  }
2160  { //For the adjoint they are added to the boundary equation colors
2161  unsigned int comp = _adjoint_boundary_equation_colors.size();
2162  for (unsigned int i = 0; i < _adjoint_boundary_equation_colors.size();
2163  ++i)
2164  {
2165  if (_adjoint_boundary_equation_colors[i] == color)
2166  {
2167  comp = i;
2168  break;
2169  }
2170  }
2171  if (comp != _adjoint_boundary_equation_colors.size())
2172  {
2173  //Seems this color is already added, however it might have been a equation color
2174  //so we don't do anything.
2175  }
2176  else
2177  {
2178  _adjoint_boundary_equation_colors.push_back(color);
2179  }
2180  }
2181  }
2182 
2183  /******************************************************/
2184 
2185  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2186  int dealdim, template<int, int> class FE, template<int, int> class DH>
2187  unsigned int
2189  {
2190  return this->GetPDE().GetStateNBlocks();
2191  }
2192 
2193  /******************************************************/
2194 
2195  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2196  int dealdim, template<int, int> class FE, template<int, int> class DH>
2197  unsigned int
2199  {
2200  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee"))
2201  {
2202  return this->GetStateNBlocks();
2203  }
2204  else
2205  {
2206  throw DOpEException("Unknown Type:" + this->GetType(),
2207  "PDEProblemContainer::GetNBlocks");
2208  }
2209  }
2210 
2211  /******************************************************/
2212 
2213  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2214  int dealdim, template<int, int> class FE, template<int, int> class DH>
2215  unsigned int
2217  unsigned int b) const
2218  {
2219  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee"))
2220  {
2221  return GetSpaceTimeHandler()->GetStateDoFsPerBlock(b);
2222  }
2223  else
2224  {
2225  throw DOpEException("Unknown Type:" + this->GetType(),
2226  "PDEProblemContainer::GetDoFsPerBlock");
2227  }
2228  }
2229 
2230  /******************************************************/
2231 
2232  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2233  int dealdim, template<int, int> class FE, template<int, int> class DH>
2234  const std::vector<unsigned int>&
2236  {
2237  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee"))
2238  {
2239  return GetSpaceTimeHandler()->GetStateDoFsPerBlock();
2240  }
2241  else
2242  {
2243  throw DOpEException("Unknown Type:" + this->GetType(),
2244  "PDEProblemContainer::GetDoFsPerBlock");
2245  }
2246  }
2247 
2248  /******************************************************/
2249 
2250  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2251  int dealdim, template<int, int> class FE, template<int, int> class DH>
2252  const dealii::ConstraintMatrix&
2254  {
2255 // std::cout << "Constraints:"
2256 // << GetSpaceTimeHandler()->GetStateDoFConstraints().n_constraints()
2257 // << " and type is " << this->GetType() << std::endl;
2258  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee"))
2259  {
2260  return GetSpaceTimeHandler()->GetStateDoFConstraints();
2261  }
2262  else
2263  {
2264  throw DOpEException("Unknown Type:" + this->GetType(),
2265  "PDEProblemContainer::GetDoFConstraints");
2266  }
2267  }
2268 
2269  /******************************************************/
2270 
2271  template<typename PDE, typename DD, typename SPARSITYPATTERN, typename VECTOR,
2272  int dealdim, template<int, int> class FE, template<int, int> class DH>
2273  bool
2275  {
2276  if (this->GetType() == "cost_functional")
2277  return false;
2278  else if (this->GetType() == "aux_functional")
2279  return _aux_functionals[this->GetTypeNum()]->NeedTime();
2280  else if (this->GetType() == "error_evaluation")
2281  return _aux_functionals[_functional_for_ee_num]->NeedTime();
2282  else
2283  throw DOpEException("Not implemented",
2284  "PDEProblemContainer::NeedTimeFunctional");
2285  }
2286 
2287 }
2288 #endif
Definition: stateproblem.h:48
StateProblem< PDEProblemContainer< PDE, DD, SPARSITYPATTERN, VECTOR, dealdim, FE, DH >, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim > & GetStateProblem()
Definition: pdeproblemcontainer.h:126
double BoundaryFunctional(const FACEDATACONTAINER &fdc)
Definition: pdeproblemcontainer.h:1052
void ElementMatrix(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)
Definition: pdeproblemcontainer.h:1439
void RegisterExceptionHandler(DOpEExceptionHandler< VECTOR > *OH)
Definition: pdeproblemcontainer.h:168
void SetInitialValues(const dealii::Function< dealdim > *values)
Definition: pdeproblemcontainer.h:609
bool HasInterfaces() const
Definition: pdeproblemcontainer.h:1874
PDEProblemContainer< PDE, DD, SPARSITYPATTERN, VECTOR, dealdim, FE, DH > & GetBaseProblem()
Definition: pdeproblemcontainer.h:141
void ReInit(std::string algo_type)
Definition: pdeproblemcontainer.h:930
unsigned int GetStateNBlocks()
Definition: pdeproblemcontainer.h:825
void RegisterOutputHandler(DOpEOutputHandler< VECTOR > *OH)
Definition: pdeproblemcontainer.h:160
const dealii::Function< dealdim > & GetDirichletValues(unsigned int color, const std::map< std::string, const dealii::Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values) const
Definition: pdeproblemcontainer.h:1987
void FaceRhs(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1379
dealii::UpdateFlags GetFaceUpdateFlags() const
Definition: pdeproblemcontainer.h:1675
virtual std::string GetName() const
Definition: pdeproblemcontainer.h:114
unsigned int GetStateNBlocks() const
Definition: pdeproblemcontainer.h:2188
std::string GetDoFType() const
Definition: pdeproblemcontainer.h:1612
void SetBoundaryFunctionalColors(unsigned int color)
Definition: pdeproblemcontainer.h:2137
Definition: elementdatacontainer.h:56
unsigned int GetNFunctionals() const
Definition: pdeproblemcontainer.h:698
void ComputeSparsityPattern(SPARSITYPATTERN &sparsity) const
Definition: pdeproblemcontainer.h:1758
void FaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: pdeproblemcontainer.h:1232
Definition: problemcontainer_internal.h:38
Definition: optproblemcontainer.h:70
const std::vector< unsigned int > & GetDirichletColors() const
Definition: pdeproblemcontainer.h:1934
std::vector< unsigned int > & GetStateBlockComponent()
Definition: pdeproblemcontainer.h:833
void ComputeMGSparsityPattern(dealii::MGLevelObject< dealii::BlockSparsityPattern > &mg_sparsity_patterns, unsigned int n_levels) const
Definition: pdeproblemcontainer.h:1777
const dealii::ConstraintMatrix & GetDoFConstraints() const
Definition: pdeproblemcontainer.h:2253
void AddFunctional(FunctionalInterface< ElementDataContainer, FaceDataContainer, DH, VECTOR, dealdim > *F)
Definition: pdeproblemcontainer.h:641
std::string GetFunctionalType() const
Definition: pdeproblemcontainer.h:1702
void SetTime(double time, const TimeIterator &interval)
Definition: pdeproblemcontainer.h:1738
const StateSpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dealdim > * GetSpaceTimeHandler() const
Definition: pdeproblemcontainer.h:768
Definition: timeiterator.h:63
bool NeedTimeFunctional() const
Definition: pdeproblemcontainer.h:2274
const PDE & GetPDE() const
Definition: problemcontainer_internal.h:103
unsigned int GetNBlocks() const
Definition: pdeproblemcontainer.h:2198
const dealii::Function< dealdim > & GetInitialValues() const
Definition: pdeproblemcontainer.h:614
const std::vector< unsigned int > & GetDoFsPerBlock() const
Definition: pdeproblemcontainer.h:2235
StateSpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dealdim > * GetSpaceTimeHandler()
Definition: pdeproblemcontainer.h:776
Definition: facedatacontainer.h:48
void ElementTimeEquationExplicit(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1203
void SetBoundaryEquationColors(unsigned int color)
Definition: pdeproblemcontainer.h:2063
dealii::UpdateFlags GetUpdateFlags() const
Definition: pdeproblemcontainer.h:1649
bool HasPoints() const
Definition: pdeproblemcontainer.h:1850
const std::map< std::string, unsigned int > & GetFunctionalPosition() const
Definition: pdeproblemcontainer.h:847
std::string GetFunctionalName() const
Definition: pdeproblemcontainer.h:1720
Definition: statespacetimehandler.h:59
double ElementFunctional(const DATACONTAINER &edc)
Definition: pdeproblemcontainer.h:987
const std::vector< unsigned int > & GetBoundaryEquationColors() const
Definition: pdeproblemcontainer.h:2041
Definition: primaldirichletdata.h:41
void ElementRhs(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1320
void ElementEquation(const DATACONTAINER &edc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: pdeproblemcontainer.h:1147
Definition: functionalinterface.h:53
void SetDirichletBoundaryColors(unsigned int color, const std::vector< bool > &comp_mask, const DD *values)
Definition: pdeproblemcontainer.h:1899
void SetFunctionalForErrorEstimation(std::string functional_name)
Definition: pdeproblemcontainer.h:670
void BoundaryRhs(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1409
bool HasFaces() const
Definition: pdeproblemcontainer.h:1820
void AddAuxiliaryToIntegrator(INTEGRATOR &)
Definition: pdeproblemcontainer.h:808
void ElementTimeMatrixExplicit(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix)
Definition: pdeproblemcontainer.h:1498
DOpEExceptionHandler< VECTOR > * GetExceptionHandler()
Definition: pdeproblemcontainer.h:745
PDEProblemContainer(PDE &pde, StateSpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dealdim > &STH)
Definition: pdeproblemcontainer.h:891
double FaceFunctional(const FACEDATACONTAINER &fdc)
Definition: pdeproblemcontainer.h:1083
void InterfaceMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)
void ElementTimeMatrix(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix)
virtual ~PDEProblemContainer()
Definition: pdeproblemcontainer.h:908
void BoundaryMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_matrix, double scale=1., double scale_ico=1.)
const FE< dealdim, dealdim > & GetFESystem() const
Definition: pdeproblemcontainer.h:1631
void BoundaryEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: pdeproblemcontainer.h:1290
void ElementTimeEquation(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1175
void PointRhs(const std::map< std::string, const dealii::Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values, VECTOR &rhs_vector, double scale=1.)
Definition: pdeproblemcontainer.h:1352
Definition: pdeproblemcontainer.h:100
const std::vector< bool > & GetDirichletCompMask(unsigned int color) const
Definition: pdeproblemcontainer.h:1952
double AlgebraicFunctional(const std::map< std::string, const dealii::Vector< double > * > &values, const std::map< std::string, const VECTOR * > &block_values)
Definition: pdeproblemcontainer.h:1113
Definition: dopeexception.h:35
void InterfaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: pdeproblemcontainer.h:1261
void DeleteAuxiliaryFromIntegrator(INTEGRATOR &)
Definition: pdeproblemcontainer.h:817
void SetType(std::string type, unsigned int num=0)
Definition: pdeproblemcontainer.h:967
const std::vector< unsigned int > & GetBoundaryFunctionalColors() const
Definition: pdeproblemcontainer.h:2116
DOpEOutputHandler< VECTOR > * GetOutputHandler()
Definition: pdeproblemcontainer.h:754
double PointFunctional(const std::map< std::string, const dealii::Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values)
Definition: pdeproblemcontainer.h:1015
Definition: optproblemcontainer.h:72
void FaceMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)