DOpE
optproblemcontainer.h
Go to the documentation of this file.
1 
24 #ifndef OptProblemContainer_H_
25 #define OptProblemContainer_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 "spacetimehandler.h"
34 #include "primaldirichletdata.h"
35 #include "tangentdirichletdata.h"
36 #include "transposeddirichletdatainterface.h"
37 #include "transposedgradientdirichletdata.h"
38 #include "transposedhessiandirichletdata.h"
39 #include "constraintvector.h"
40 #include "controlvector.h"
41 #include "statevector.h"
42 #include "elementdatacontainer.h"
43 #include "facedatacontainer.h"
44 #include "stateproblem.h"
45 #include "dopetypes.h"
46 #include "dwrdatacontainer.h"
47 #include "problemcontainer_internal.h"
48 
49 #include <deal.II/lac/vector.h>
50 #include <deal.II/lac/full_matrix.h>
51 #include <deal.II/grid/tria_iterator.h>
52 #include <deal.II/dofs/dof_handler.h>
53 #include <deal.II/dofs/dof_accessor.h>
54 #include <deal.II/dofs/dof_tools.h>
55 #include <deal.II/fe/fe_system.h>
56 #include <deal.II/fe/fe_values.h>
57 #include <deal.II/fe/mapping.h>
58 #include <deal.II/base/quadrature_lib.h>
59 #include <deal.II/lac/block_sparsity_pattern.h>
60 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
61 
62 #include <assert.h>
63 #include <string>
64 #include <vector>
65 
66 namespace DOpE
67 {
68  //Predeclaration necessary
69  template<typename VECTOR>
70  class DOpEOutputHandler;
71  template<typename VECTOR>
72  class DOpEExceptionHandler;
74 
96  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
97  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
98  typename VECTOR, int dopedim, int dealdim,
99  template<int, int> class FE = dealii::FESystem,
100  template<int, int> class DH = dealii::DoFHandler>
101  class OptProblemContainer : public ProblemContainerInternal<PDE>
102  {
103  public:
104  OptProblemContainer(FUNCTIONAL& functional, PDE& pde,
105  CONSTRAINTS& constraints,
106  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>& STH);
107 
108  /******************************************************/
109 
110  virtual ~OptProblemContainer();
111 
112  /******************************************************/
113 
114  virtual std::string
115  GetName() const
116  {
117  return "OptProblem";
118  }
119 
120  /******************************************************/
124  StateProblem<
125  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
126  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>,
127  PDE, DD, SPARSITYPATTERN, VECTOR, dealdim>&
128  GetStateProblem()
129  {
130  if (state_problem_ == NULL)
131  {
132  state_problem_ = new StateProblem<
133  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
134  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE,
135  DH>, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim>(
136  *this, this->GetPDE());
137  }
138  return *state_problem_;
139  }
140 
141  //TODO This is Pfush needed to split into different subproblems and allow optproblem to
142  //be substituted as any of these problems. Can be removed once the splitting is complete.
143  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
144  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>&
145  GetBaseProblem()
146  {
147  return *this;
148  }
149  /******************************************************/
150 
158  void
159  ReInit(std::string algo_type);
160 
161  /******************************************************/
162 
163  void
164  RegisterOutputHandler(DOpEOutputHandler<VECTOR>* OH)
165  {
166  OutputHandler_ = OH;
167  }
168 
169  /******************************************************/
170 
171  void
172  RegisterExceptionHandler(DOpEExceptionHandler<VECTOR>* OH)
173  {
174  ExceptionHandler_ = OH;
175  }
176 
177  /******************************************************/
178 
186  void
187  SetType(std::string type, unsigned int num = 0);
188 
189  /******************************************************/
190 
206  template<typename DATACONTAINER>
207  double
208  ElementFunctional(const DATACONTAINER& edc);
209 
210  /******************************************************/
211 
221  double
222  PointFunctional(
223  const std::map<std::string, const dealii::Vector<double>*> &param_values,
224  const std::map<std::string, const VECTOR*> &domain_values);
225 
226  /******************************************************/
227 
244  template<typename FACEDATACONTAINER>
245  double
246  BoundaryFunctional(const FACEDATACONTAINER& fdc);
247 
248  /******************************************************/
249 
256  template<typename FACEDATACONTAINER>
257  double
258  FaceFunctional(const FACEDATACONTAINER& fdc);
259 
260  /******************************************************/
267  double
268  AlgebraicFunctional(
269  const std::map<std::string, const dealii::Vector<double>*> &values,
270  const std::map<std::string, const VECTOR*> &block_values);
271 
278  void
279  AlgebraicResidual(VECTOR& residual,
280  const std::map<std::string, const dealii::Vector<double>*> &values,
281  const std::map<std::string, const VECTOR*> &block_values);
282 
283  /******************************************************/
284 
308  template<typename DATACONTAINER>
309  void
310  ElementEquation(const DATACONTAINER& edc,
311  dealii::Vector<double> &local_vector, double scale,
312  double scale_ico);
313 
314  /******************************************************/
315 
321  template<typename DATACONTAINER>
322  void
323  ElementTimeEquation(const DATACONTAINER& dc,
324  dealii::Vector<double> &local_vector, double scale = 1.);
325 
326  /******************************************************/
327 
337  template<typename DATACONTAINER>
338  void
339  ElementTimeEquationExplicit(const DATACONTAINER& dc,
340  dealii::Vector<double> &local_vector, double scale = 1.);
341 
342  /******************************************************/
343 
355  template<typename DATACONTAINER>
356  void
357  ElementRhs(const DATACONTAINER& dc,
358  dealii::Vector<double> &local_vector, double scale = 1.);
359 
360  /******************************************************/
361 
374  void
375  PointRhs(
376  const std::map<std::string, const dealii::Vector<double>*> &param_values,
377  const std::map<std::string, const VECTOR*> &domain_values,
378  VECTOR& rhs_vector, double scale = 1.);
379 
380  /******************************************************/
381 
405  template<typename DATACONTAINER>
406  void
407  ElementMatrix(const DATACONTAINER& dc,
408  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
409  double scale_ico = 1.);
410 
411  /******************************************************/
412 
424  template<typename DATACONTAINER>
425  void
426  ElementTimeMatrix(const DATACONTAINER& dc,
427  dealii::FullMatrix<double> &local_entry_matrix);
428 
429  /******************************************************/
430 
440  template<typename DATACONTAINER>
441  void
442  ElementTimeMatrixExplicit(const DATACONTAINER& dc,
443  dealii::FullMatrix<double> &local_entry_matrix);
444 
445  /******************************************************/
446 
453  template<typename FACEDATACONTAINER>
454  void
455  FaceEquation(const FACEDATACONTAINER& dc,
456  dealii::Vector<double> &local_vector, double scale,
457  double scale_ico);
458 
459  /******************************************************/
466  //FIXME maybe InterfaceEquation and InterfaceMatrix could get
467  //integrated into FaceEquation and FaceMatrix?
468  template<typename FACEDATACONTAINER>
469  void
470  InterfaceEquation(const FACEDATACONTAINER& dc,
471  dealii::Vector<double> &local_vector, double scale,
472  double scale_ico);
473 
474  /******************************************************/
481  template<typename FACEDATACONTAINER>
482  void
483  FaceRhs(const FACEDATACONTAINER& dc,
484  dealii::Vector<double> &local_vector, double scale = 1.);
485 
486  /******************************************************/
487 
494  template<typename FACEDATACONTAINER>
495  void
496  FaceMatrix(const FACEDATACONTAINER& dc,
497  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
498  double scale_ico = 1.);
499 
500  /******************************************************/
507  template<typename FACEDATACONTAINER>
508  void
509  InterfaceMatrix(const FACEDATACONTAINER& dc,
510  dealii::FullMatrix<double> &local_entry_matrix, double scale = 1.,
511  double scale_ico = 1.);
512 
513  /******************************************************/
514 
521  template<typename FACEDATACONTAINER>
522  void
523  BoundaryEquation(const FACEDATACONTAINER& dc,
524  dealii::Vector<double> &local_vector, double scale,
525  double scale_ico);
526 
527  /******************************************************/
528 
535  template<typename FACEDATACONTAINER>
536  void
537  BoundaryRhs(const FACEDATACONTAINER& dc,
538  dealii::Vector<double> &local_vector, double scale = 1.);
539 
540  /******************************************************/
541 
548  template<typename FACEDATACONTAINER>
549  void
550  BoundaryMatrix(const FACEDATACONTAINER& dc,
551  dealii::FullMatrix<double> &local_matrix, double scale = 1.,
552  double scale_ico = 1.);
553  /******************************************************/
554  void
555  ComputeLocalControlConstraints(VECTOR& constraints,
556  const std::map<std::string, const dealii::Vector<double>*> &values,
557  const std::map<std::string, const VECTOR*> &block_values);
558  /******************************************************/
559 
565  void
566  GetControlBoxConstraints(VECTOR& lb, VECTOR& ub) const
567  {
568  this->GetConstraints()->GetControlBoxConstraints(lb, ub);
569  }
570 
571  /******************************************************/
572 
573  const FE<dealdim, dealdim>&
574  GetFESystem() const;
575 
576  /******************************************************/
583  bool
584  HasFaces() const;
585 
586  /******************************************************/
593  bool
594  HasPoints() const;
595 
596  /******************************************************/
603  bool
604  HasInterfaces() const;
605 
606  /******************************************************/
607 
608  dealii::UpdateFlags
609  GetUpdateFlags() const;
610 
611  /******************************************************/
612 
613  dealii::UpdateFlags
614  GetFaceUpdateFlags() const;
615 
616  /******************************************************/
617 
618  void
619  SetControlDirichletBoundaryColors(unsigned int color,
620  const std::vector<bool>& comp_mask,
621  const DOpEWrapper::Function<dealdim>* values);
622 
623  /******************************************************/
624 
625  void
626  SetDirichletBoundaryColors(unsigned int color,
627  const std::vector<bool>& comp_mask, const DD* values);
628 
629  /******************************************************/
630 
631  const std::vector<unsigned int>&
632  GetDirichletColors() const;
633  const std::vector<unsigned int>&
634  GetTransposedDirichletColors() const;
635  const std::vector<bool>&
636  GetDirichletCompMask(unsigned int color) const;
637  const std::vector<bool>&
638  GetTransposedDirichletCompMask(unsigned int color) const;
639 
640  /******************************************************/
641 
642  const dealii::Function<dealdim> &
643  GetDirichletValues(unsigned int color,
644  const std::map<std::string, const dealii::Vector<double>*> &param_values,
645  const std::map<std::string, const VECTOR*> &domain_values) const;
646 
647  /******************************************************/
648 
649  const TransposedDirichletDataInterface<dealdim> &
650  GetTransposedDirichletValues(unsigned int color,
651  const std::map<std::string, const dealii::Vector<double>*> &param_values,
652  const std::map<std::string, const VECTOR*> &domain_values) const;
653 
654  /******************************************************/
655 
656  void
657  SetInitialValues(const dealii::Function<dealdim>* values)
658  {
659  initial_values_ = values;
660  }
661  const dealii::Function<dealdim>&
662  GetInitialValues() const
663  {
664  return *initial_values_;
665  }
666 
667  /******************************************************/
668 
669  void
670  SetControlBoundaryEquationColors(unsigned int color);
671  void
672  SetBoundaryEquationColors(unsigned int color);
673  const std::vector<unsigned int>&
674  GetBoundaryEquationColors() const;
675  void
676  SetBoundaryFunctionalColors(unsigned int color);
677  const std::vector<unsigned int>&
678  GetBoundaryFunctionalColors() const;
679 
680  /******************************************************/
681 
691  void
692  AddFunctional(FUNCTIONAL_INTERFACE* F)
693  {
694  aux_functionals_.push_back(F);
695  if (functional_position_.find(F->GetName())
696  != functional_position_.end())
697  {
698  throw DOpEException(
699  "You cant add two functionals with the same name.",
700  "OPTProblemContainer::AddFunctional");
701  }
702  functional_position_[F->GetName()] = aux_functionals_.size();
703  //remember! At functional_values_[0] we store always the cost functional!
704  }
705 
706  /******************************************************/
707 
717  void
718  SetFunctionalForErrorEstimation(std::string functional_name)
719  {
720  if (GetFunctional()->GetName() == functional_name)
721  {
722  functional_for_ee_is_cost_ = true;
723  functional_for_ee_num_ = dealii::numbers::invalid_unsigned_int;
724  }
725  else
726  {
727  functional_for_ee_is_cost_ = false;
728  bool found = false;
729  //we go through all aux functionals.
730  for (unsigned int i = 0; i < this->GetNFunctionals(); i++)
731  {
732  if (aux_functionals_[i]->GetName() == functional_name)
733  {
734  //if the names match, we have found our functional.
735  found = true;
736  functional_for_ee_num_ = i;
737  }
738  }
739  //If we have not found a functional with the given name,
740  //we throw an error.
741  if (!found)
742  {
743  throw DOpEException(
744  "Can't find functional " + functional_name
745  + " in aux_functionals_",
746  "Optproblem::SetFunctionalForErrorEstimation");
747  }
748  }
749  }
750 
751  /******************************************************/
752 
753  unsigned int
754  GetNFunctionals() const
755  {
756  return aux_functionals_.size();
757  }
758 
759  /******************************************************/
760 
761  unsigned int
762  GetControlNBlocks() const;
763 
764  /******************************************************/
765 
766  unsigned int
767  GetStateNBlocks() const;
768 
769  /******************************************************/
770 
771  unsigned int
772  GetNBlocks() const;
773 
774  /******************************************************/
775 
776  unsigned int
777  GetDoFsPerBlock(unsigned int b) const;
778 
779  /******************************************************/
780 
781  const std::vector<unsigned int>&
782  GetDoFsPerBlock() const;
783 
784  /******************************************************/
785 
786  const dealii::ConstraintMatrix&
787  GetDoFConstraints() const;
788 
789  /******************************************************/
790 
791  std::string
792  GetDoFType() const;
793  std::string
794  GetFunctionalType() const;
795  std::string
796  GetFunctionalName() const;
797  std::string
798  GetConstraintType() const;
799 
800  /******************************************************/
801 
802  bool
803  NeedTimeFunctional() const;
804 
805  /******************************************************/
806 
807  bool
808  HasControlInDirichletData() const;
809 
810  /******************************************************/
811 
812  DOpEExceptionHandler<VECTOR>*
813  GetExceptionHandler()
814  {
815  assert(ExceptionHandler_);
816  return ExceptionHandler_;
817  }
818 
819  /******************************************************/
820 
821  DOpEOutputHandler<VECTOR>*
822  GetOutputHandler()
823  {
824  assert(OutputHandler_);
825  return OutputHandler_;
826  }
827 
828  /******************************************************/
829 
839  void
840  SetTime(double time,
841  unsigned int time_dof_number,
842  const TimeIterator& interval, bool initial = false);
843 
844  /******************************************************/
845 
846  const SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>*
847  GetSpaceTimeHandler() const
848  {
849  return STH_;
850  }
851 
852  /******************************************************/
853 
854  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>*
855  GetSpaceTimeHandler()
856  {
857  return STH_;
858  }
859 
860  /******************************************************/
861 
862  void
863  ComputeSparsityPattern(SPARSITYPATTERN & sparsity) const;
864 
865  /******************************************************/
866 
867  void
868  PostProcessConstraints(ConstraintVector<VECTOR>& g) const;
869 
870  /******************************************************/
871 
872  void
873  AddAuxiliaryControl(const ControlVector<VECTOR>* c, std::string name);
874 
875  /******************************************************/
876 
877  void
878  AddAuxiliaryState(const StateVector<VECTOR>* c, std::string name);
879 
880  /******************************************************/
881 
882  void
883  AddAuxiliaryConstraint(const ConstraintVector<VECTOR>* c,
884  std::string name);
885 
886  /******************************************************/
887 
888  const ControlVector<VECTOR>*
889  GetAuxiliaryControl(std::string name) const;
890 
891  /******************************************************/
892 
893  const StateVector<VECTOR>*
894  GetAuxiliaryState(std::string name) const;
895 
896  /******************************************************/
897 
898  void
899  DeleteAuxiliaryControl(std::string name);
900 
901  /******************************************************/
902 
903  void
904  DeleteAuxiliaryState(std::string name);
905 
906  /******************************************************/
907 
908  void
909  DeleteAuxiliaryConstraint(std::string name);
910 
911  /*****************************************************************/
912 
913  const ConstraintVector<VECTOR>*
914  GetAuxiliaryConstraint(std::string name)
915  {
916  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
917  auxiliary_constraints_.find(name);
918  if (it == auxiliary_constraints_.end())
919  {
920  throw DOpEException("Could not find data" + name,
921  "OptProblemContainer::GetAuxiliaryConstraint");
922  }
923  return it->second;
924  }
925  /*****************************************************************/
934  template<typename INTEGRATOR>
935  void
936  AddAuxiliaryToIntegrator(INTEGRATOR& integrator)
937  {
938  {
939  //Only add control vector if the vecor is distributed in time, or
940  //we are calculating the initial value.
941  if((GetSpaceTimeHandler()->GetControlType()==DOpEtypes::ControlType::initial && initial_)
942  || (GetSpaceTimeHandler()->GetControlType()!=DOpEtypes::ControlType::initial))
943  {
944  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
945  auxiliary_controls_.begin();
946  for (; it != auxiliary_controls_.end(); it++)
947  {
948  if (dopedim == dealdim)
949  {
950  integrator.AddDomainData(it->first,
951  &(it->second->GetSpacialVector()));
952  }
953  else if (dopedim == 0)
954  {
955  integrator.AddParamData(it->first,
956  &(it->second->GetSpacialVectorCopy()));
957  }
958  else
959  {
960  throw DOpEException("dopedim not implemented",
961  "OptProblemContainer::AddAuxiliaryToIntegrator");
962  }
963  }
964  }
965  }
966  {
967  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
968  auxiliary_state_.begin();
969  for (; it != auxiliary_state_.end(); it++)
970  {
971  integrator.AddDomainData(it->first,
972  &(it->second->GetSpacialVector()));
973  }
974  }
975  {
976  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
977  auxiliary_constraints_.begin();
978  for (; it != auxiliary_constraints_.end(); it++)
979  {
980  integrator.AddDomainData(it->first + "_local",
981  &(it->second->GetSpacialVector("local")));
982  integrator.AddParamData(it->first + "_global",
983  &(it->second->GetGlobalConstraints()));
984  }
985  }
986  }
987 
988  /*****************************************************************/
1001  template<typename INTEGRATOR>
1002  void
1003  AddNextAuxiliaryToIntegrator(INTEGRATOR& integrator)
1004  {
1005  {
1006  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
1007  auxiliary_state_.find("state");
1008  if (it != auxiliary_state_.end())
1009  {
1010  integrator.AddDomainData("state_i+1",
1011  &(it->second->GetNextSpacialVector()));
1012  }
1013  }
1014  }
1015  /*****************************************************************/
1028  template<typename INTEGRATOR>
1029  void
1030  AddPreviousAuxiliaryToIntegrator(INTEGRATOR& integrator)
1031  {
1032  {
1033  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
1034  auxiliary_state_.find("state");
1035  if (it != auxiliary_state_.end())
1036  {
1037  integrator.AddDomainData("state_i-1",
1038  &(it->second->GetPreviousSpacialVector()));
1039  }
1040  }
1041  }
1042 
1043  /*****************************************************************/
1056  template<typename INTEGRATOR>
1057  void
1058  DeleteNextAuxiliaryFromIntegrator(INTEGRATOR& integrator)
1059  {
1060  {
1061  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
1062  auxiliary_state_.find("state");
1063  if (it != auxiliary_state_.end())
1064  {
1065  integrator.DeleteDomainData("state_i+1");
1066  }
1067  }
1068  }
1069  /*****************************************************************/
1082  template<typename INTEGRATOR>
1083  void
1084  DeletePreviousAuxiliaryFromIntegrator(INTEGRATOR& integrator)
1085  {
1086  {
1087  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
1088  auxiliary_state_.find("state");
1089  if (it != auxiliary_state_.end())
1090  {
1091  integrator.DeleteDomainData("state_i-1");
1092  }
1093  }
1094  }
1095 
1096 
1097  /******************************************************/
1106  template<typename INTEGRATOR>
1107  void
1108  DeleteAuxiliaryFromIntegrator(INTEGRATOR& integrator)
1109  {
1110  {
1111  //Only delete control vector if the vecor is distributed in time, or
1112  //we are calculating the initial value. Otherwise it would not have been added.
1113  if((GetSpaceTimeHandler()->GetControlType()==DOpEtypes::ControlType::initial && initial_)
1114  || (GetSpaceTimeHandler()->GetControlType()!=DOpEtypes::ControlType::initial))
1115  {
1116  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
1117  auxiliary_controls_.begin();
1118  for (; it != auxiliary_controls_.end(); it++)
1119  {
1120  if (dopedim == dealdim)
1121  {
1122  integrator.DeleteDomainData(it->first);
1123  }
1124  else if (dopedim == 0)
1125  {
1126  integrator.DeleteParamData(it->first);
1127  it->second->UnLockCopy();
1128  }
1129  else
1130  {
1131  throw DOpEException("dopedim not implemented",
1132  "OptProblemContainer::AddAuxiliaryToIntegrator");
1133  }
1134  }
1135  }
1136  }
1137  {
1138  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
1139  auxiliary_state_.begin();
1140  for (; it != auxiliary_state_.end(); it++)
1141  {
1142  integrator.DeleteDomainData(it->first);
1143  }
1144  }
1145  {
1146  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
1147  auxiliary_constraints_.begin();
1148  for (; it != auxiliary_constraints_.end(); it++)
1149  {
1150  integrator.DeleteDomainData(it->first + "_local");
1151  integrator.DeleteParamData(it->first + "_global");
1152  }
1153  }
1154  }
1155 
1156  /******************************************************/
1157  const std::map<std::string, unsigned int>&
1158  GetFunctionalPosition() const
1159  {
1160  return functional_position_;
1161  }
1162 
1163  unsigned int
1164  GetStateNBlocks()
1165  {
1166  return this->GetPDE().GetStateNBlocks();
1167  }
1168 
1169  /******************************************************/
1170 
1171  std::vector<unsigned int>&
1172  GetControlBlockComponent()
1173  {
1174  return this->GetPDE().GetControlBlockComponent();
1175  }
1176  /******************************************************/
1177 
1178  std::vector<unsigned int>&
1179  GetStateBlockComponent()
1180  {
1181  return this->GetPDE().GetStateBlockComponent();
1182  }
1183 
1184  /******************************************************/
1185  /******************************************************/
1186  /****For the initial values ***************/
1187  template<typename DATACONTAINER>
1188  void
1189  Init_ElementEquation(const DATACONTAINER& edc,
1190  dealii::Vector<double> &local_vector, double scale,
1191  double scale_ico)
1192  {
1193  if (this->GetType() == "adjoint" || this->GetType() == "tangent"
1194  || this->GetType() == "adjoint_hessian")
1195  {
1196  this->GetPDE().Init_ElementEquation(edc, local_vector, scale,
1197  scale_ico);
1198  }
1199  else
1200  {
1201  throw DOpEException("Not implemented",
1202  "OptProblemContainer::Init_ElementEquation");
1203  }
1204  }
1205 
1206  template<typename DATACONTAINER>
1207  void
1208  Init_ElementRhs(const DATACONTAINER& edc,
1209  dealii::Vector<double> &local_vector, double scale)
1210  {
1211  //FIXME: We should take care of the cases of boundary and face functionals...
1212  if (this->GetType() == "adjoint")
1213  {
1214  if (GetFunctional()->NeedTime())
1215  {
1216  if (GetFunctional()->GetType().find("timelocal")
1217  != std::string::npos)
1218  {
1219  if (GetFunctional()->GetType().find("domain")
1220  != std::string::npos)
1221  {
1222  GetFunctional()->ElementValue_U(edc, local_vector, scale);
1223  }
1224  }
1225  }
1226  }
1227  else if (this->GetType() == "tangent")
1228  {
1229  this->GetPDE().Init_ElementRhs_QT(edc, local_vector, scale);
1230  }
1231  else if (this->GetType() == "adjoint_hessian")
1232  {
1233  if (GetFunctional()->NeedTime())
1234  {
1235  if (GetFunctional()->GetType().find("timelocal")
1236  != std::string::npos)
1237  {
1238  if (GetFunctional()->GetType().find("domain")
1239  != std::string::npos)
1240  {
1241  GetFunctional()->ElementValue_UU(edc, local_vector, scale);
1242  }
1243  }
1244  }
1245  }
1246  else
1247  {
1248  throw DOpEException("Not implemented",
1249  "OptProblemContainer::Init_ElementRhs");
1250  }
1251  }
1252 
1253  void
1254  Init_PointRhs(
1255  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1256  const std::map<std::string, const VECTOR*> &domain_values,
1257  VECTOR& rhs_vector, double scale = 1.)
1258  {
1259  if (this->GetType() == "adjoint")
1260  {
1261  if (GetFunctional()->NeedTime())
1262  {
1263  if (GetFunctional()->GetType().find("timelocal")
1264  != std::string::npos)
1265  {
1266  if (GetFunctional()->GetType().find("point")
1267  != std::string::npos)
1268  {
1269  GetFunctional()->PointValue_U(param_values, domain_values,
1270  rhs_vector, scale);
1271  }
1272  }
1273  }
1274  }
1275  else if (this->GetType() == "tangent")
1276  {
1277  //Nothing to do for tangent, since no point sources are allowed in the initial condition.
1278  }
1279  else if (this->GetType() == "adjoint_hessian")
1280  {
1281  if (GetFunctional()->NeedTime())
1282  {
1283  if (GetFunctional()->GetType().find("timelocal")
1284  != std::string::npos)
1285  {
1286  if (GetFunctional()->GetType().find("point")
1287  != std::string::npos)
1288  {
1289  GetFunctional()->PointValue_UU(param_values, domain_values,
1290  rhs_vector, scale);
1291  }
1292  }
1293  }
1294  }
1295  else
1296  {
1297  throw DOpEException("Not implemented",
1298  "OptProblemContainer::Init_PointRhs");
1299  }
1300  }
1301 
1302  template<typename DATACONTAINER>
1303  void
1304  Init_ElementMatrix(const DATACONTAINER& edc,
1305  dealii::FullMatrix<double> &local_entry_matrix, double scale,
1306  double scale_ico)
1307  {
1308  if (this->GetType() == "adjoint" || this->GetType() == "tangent"
1309  || this->GetType() == "adjoint_hessian")
1310  {
1311  this->GetPDE().Init_ElementMatrix(edc, local_entry_matrix, scale,
1312  scale_ico);
1313  }
1314  else
1315  {
1316  throw DOpEException("Not implemented",
1317  "OptProblemContainer::Init_ElementMatrix");
1318  }
1319  }
1320 
1321  /******************************************************/
1322 
1327  bool
1328  EEFunctionalIsCost() const
1329  {
1330  return functional_for_ee_is_cost_;
1331  }
1332 
1333  template<typename ELEMENTITERATOR>
1334  bool AtInterface(ELEMENTITERATOR& element, unsigned int face) const;
1335 
1336  protected:
1337  FUNCTIONAL*
1338  GetFunctional();
1339  const FUNCTIONAL*
1340  GetFunctional() const;
1341  CONSTRAINTS *
1342  GetConstraints()
1343  {
1344  return constraints_;
1345  }
1346  const CONSTRAINTS *
1347  GetConstraints() const
1348  {
1349  return constraints_;
1350  }
1351 
1352  const VECTOR*
1353  GetBlockVector(const std::map<std::string, const VECTOR*>& values,
1354  std::string name)
1355  {
1356  typename std::map<std::string, const VECTOR*>::const_iterator it =
1357  values.find(name);
1358  if (it == values.end())
1359  {
1360  throw DOpEException("Did not find " + name,
1361  "OptProblemContainer::GetBlockVector");
1362  }
1363  return it->second;
1364  }
1365  const dealii::Vector<double>*
1366  GetVector(const std::map<std::string, const Vector<double>*>& values,
1367  std::string name)
1368  {
1369  typename std::map<std::string, const Vector<double>*>::const_iterator it =
1370  values.find(name);
1371  if (it == values.end())
1372  {
1373  throw DOpEException("Did not find " + name,
1374  "OptProblemContainer::GetVector");
1375  }
1376  return it->second;
1377  }
1378 
1379  /******************************************************/
1380  private:
1381  DOpEExceptionHandler<VECTOR>* ExceptionHandler_;
1382  DOpEOutputHandler<VECTOR>* OutputHandler_;
1383  std::string algo_type_;
1384 
1385  bool functional_for_ee_is_cost_;
1386  double c_interval_length_, interval_length_;
1387  unsigned int functional_for_ee_num_;
1388  std::vector<FUNCTIONAL_INTERFACE*> aux_functionals_;
1389  std::map<std::string, unsigned int> functional_position_;
1390  FUNCTIONAL* functional_;
1391  CONSTRAINTS* constraints_;
1392  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>* STH_;
1393 
1394  std::vector<unsigned int> control_dirichlet_colors_;
1395  std::vector<unsigned int> control_transposed_dirichlet_colors_;
1396  std::vector<std::vector<bool> > control_dirichlet_comps_;
1397  std::vector<const DOpEWrapper::Function<dealdim>*> control_dirichlet_values_;
1398  std::vector<
1399  TransposedGradientDirichletData<DD, VECTOR, dealdim>*> transposed_control_gradient_dirichlet_values_;
1400  std::vector<
1401  TransposedHessianDirichletData<DD, VECTOR, dealdim>*> transposed_control_hessian_dirichlet_values_;
1402 
1403  std::vector<unsigned int> dirichlet_colors_;
1404  std::vector<std::vector<bool> > dirichlet_comps_;
1405  std::vector<PrimalDirichletData<DD, VECTOR, dealdim>*> primal_dirichlet_values_;
1406  std::vector<TangentDirichletData<DD, VECTOR, dealdim>*> tangent_dirichlet_values_;
1407  const dealii::Function<dealdim>* zero_dirichlet_values_;
1408 
1409  const dealii::Function<dealdim>* initial_values_;
1410 
1411  std::vector<unsigned int> control_boundary_equation_colors_;
1412  std::vector<unsigned int> state_boundary_equation_colors_;
1413  std::vector<unsigned int> adjoint_boundary_equation_colors_;
1414 
1415  std::vector<unsigned int> boundary_functional_colors_;
1416 
1417  std::map<std::string, const ControlVector<VECTOR>*> auxiliary_controls_;
1418  std::map<std::string, const StateVector<VECTOR>*> auxiliary_state_;
1419  std::map<std::string, const ConstraintVector<VECTOR>*> auxiliary_constraints_;
1420 
1421  bool initial_; //Do we solve the problem at initial time?
1422 
1423  StateProblem<
1424  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1425  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>,
1426  PDE, DD, SPARSITYPATTERN, VECTOR, dealdim> * state_problem_;
1427 
1428  friend class StateProblem<
1429  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1430  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>,
1431  PDE, DD, SPARSITYPATTERN, VECTOR, dealdim> ;
1432  };
1433  /******************************************************/
1434 
1435  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1436  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1437  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1438  template<int, int> class DH>
1439  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1440  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::OptProblemContainer(
1441  FUNCTIONAL& functional, PDE& pde, CONSTRAINTS& constraints,
1442  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>& STH) :
1443  ProblemContainerInternal<PDE>(pde), functional_(&functional), constraints_(
1444  &constraints), STH_(&STH), state_problem_(NULL)
1445  {
1446  ExceptionHandler_ = NULL;
1447  OutputHandler_ = NULL;
1448  zero_dirichlet_values_ = new ZeroFunction<dealdim>(
1449  this->GetPDE().GetStateNComponents());
1450  algo_type_ = "";
1451  functional_position_[functional_->GetName()] = 0;
1452  //remember! At functional_values_[0] we store always the cost functional!
1453  functional_for_ee_num_ = dealii::numbers::invalid_unsigned_int;
1454  c_interval_length_ = 1.;
1455  interval_length_ = 1.;
1456  }
1457 
1458  /******************************************************/
1459 
1460  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1461  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1462  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1463  template<int, int> class DH>
1464  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1465  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::~OptProblemContainer()
1466  {
1467  if (zero_dirichlet_values_ != NULL)
1468  {
1469  delete zero_dirichlet_values_;
1470  }
1471 
1472  for (unsigned int i = 0;
1473  i < transposed_control_gradient_dirichlet_values_.size(); i++)
1474  {
1475  if (transposed_control_gradient_dirichlet_values_[i] != NULL)
1476  delete transposed_control_gradient_dirichlet_values_[i];
1477  }
1478  for (unsigned int i = 0;
1479  i < transposed_control_hessian_dirichlet_values_.size(); i++)
1480  {
1481  if (transposed_control_hessian_dirichlet_values_[i] != NULL)
1482  delete transposed_control_hessian_dirichlet_values_[i];
1483  }
1484  for (unsigned int i = 0; i < primal_dirichlet_values_.size(); i++)
1485  {
1486  if (primal_dirichlet_values_[i] != NULL)
1487  delete primal_dirichlet_values_[i];
1488  }
1489  for (unsigned int i = 0; i < tangent_dirichlet_values_.size(); i++)
1490  {
1491  if (tangent_dirichlet_values_[i] != NULL)
1492  delete tangent_dirichlet_values_[i];
1493  }
1494  if (state_problem_ != NULL)
1495  {
1496  delete state_problem_;
1497  }
1498  }
1499 
1500  /******************************************************/
1501 
1502  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1503  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1504  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1505  template<int, int> class DH>
1506  void
1507  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1508  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ReInit(
1509  std::string algo_type)
1510  {
1511  if (state_problem_ != NULL)
1512  {
1513  delete state_problem_;
1514  state_problem_ = NULL;
1515  }
1516 
1517  if (algo_type_ != algo_type && algo_type_ != "")
1518  {
1519  throw DOpEException("Conflicting Algorithms!",
1520  "OptProblemContainer::ReInit");
1521  }
1522  else
1523  {
1524  algo_type_ = algo_type;
1525  this->SetTypeInternal("");
1526 
1527  if (algo_type_ == "reduced")
1528  {
1529  GetSpaceTimeHandler()->ReInit(this->GetPDE().GetControlNBlocks(),
1530  this->GetPDE().GetControlBlockComponent(),
1531  this->GetPDE().GetStateNBlocks(),
1532  this->GetPDE().GetStateBlockComponent());
1533  }
1534  else
1535  {
1536  throw DOpEException("Unknown Algorithm " + algo_type_,
1537  "OptProblemContainer::ReInit");
1538  }
1539  }
1540  }
1541 
1542  /******************************************************/
1543 
1544  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1545  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1546  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1547  template<int, int> class DH>
1548  void
1549  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1550  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetType(
1551  std::string type, unsigned int num)
1552  {
1553  if (this->GetType() != type || this->GetTypeNum() != num)
1554  {
1555  this->SetTypeNumInternal(num);
1556  this->SetTypeInternal(type);
1557  this->GetPDE().SetProblemType(type);
1558  if (functional_for_ee_num_ != dealii::numbers::invalid_unsigned_int)
1559  aux_functionals_[functional_for_ee_num_]->SetProblemType(type);
1560  this->GetConstraints()->SetProblemType(type, num);
1561 
1562 #if dope_dimension > 0
1563  if(dealdim == dopedim)
1564  {
1565  //Prepare DoFHandlerPointer
1566 
1567  {
1568  if(this->GetType() == "state" ||this->GetType() == "adjoint"
1569  || this->GetType() == "adjoint_for_ee" || this->GetType() == "cost_functional"
1570  || this->GetType() == "aux_functional" || this->GetType() == "functional_for_ee"
1571  || this->GetType() == "tangent" || this->GetType() == "adjoint_hessian"
1572  || this->GetType() == "error_evaluation"
1573  || this->GetType().find("constraints") != std::string::npos)
1574  {
1575  GetSpaceTimeHandler()->SetDoFHandlerOrdering(1,0);
1576  }
1577  else if (this->GetType() == "gradient"||this->GetType() == "hessian"||this->GetType() == "hessian_inverse" || this->GetType() == "global_constraint_gradient"|| this->GetType() == "global_constraint_hessian")
1578  {
1579  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0,1);
1580  }
1581  else
1582  {
1583  throw DOpEException("problem_type_ : "+this->GetType()+" not implemented!", "OptProblemContainer::SetType");
1584  }
1585  }
1586  }
1587  else
1588  {
1589  throw DOpEException("dopedim not implemented", "OptProblemContainer::SetType");
1590  }
1591 #else
1592  //dopedim ==0
1593  {
1594  //Prepare DoFHandlerPointer
1595  {
1596 
1597  if (this->GetType() == "state" || this->GetType() == "adjoint"
1598  || this->GetType() == "adjoint_for_ee"
1599  || this->GetType() == "functional_for_ee"
1600  || this->GetType() == "cost_functional"
1601  || this->GetType() == "aux_functional"
1602  || this->GetType() == "tangent"
1603  || this->GetType() == "error_evaluation"
1604  || this->GetType() == "adjoint_hessian")
1605  {
1606  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0, 0);
1607  }
1608  else if (this->GetType() == "gradient"
1609  || this->GetType() == "hessian_inverse"
1610  || this->GetType() == "hessian")
1611  {
1612  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0, 0);
1613  }
1614  else
1615  {
1616  throw DOpEException(
1617  "problem_type_ : " + this->GetType() + " not implemented!",
1618  "OptProblemContainer::SetType");
1619  }
1620  }
1621  }
1622 #endif
1623  }
1624  }
1625 
1626  /******************************************************/
1627 
1628  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1629  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1630  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1631  template<int, int> class DH>
1632  template<typename DATACONTAINER>
1633  double
1634  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1635  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementFunctional(
1636  const DATACONTAINER& edc)
1637  {
1638 
1639  if (this->GetType() == "cost_functional")
1640  {
1641  // state values in quadrature points
1642  return GetFunctional()->ElementValue(edc);
1643  }
1644  else if (this->GetType() == "aux_functional")
1645  {
1646  // state values in quadrature points
1647  return aux_functionals_[this->GetTypeNum()]->ElementValue(edc);
1648  }
1649  else if (this->GetType() == "functional_for_ee")
1650  { // TODO is this correct? Should not be needed.
1651  return aux_functionals_[functional_for_ee_num_]->ElementValue(edc);
1652  }
1653  else if (this->GetType().find("constraints") != std::string::npos)
1654  {
1655  return GetConstraints()->ElementValue(edc);
1656  }
1657  else
1658  {
1659  throw DOpEException("Not implemented",
1660  "OptProblemContainer::ElementFunctional");
1661  }
1662  }
1663 
1664  /******************************************************/
1665  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1666  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1667  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1668  template<int, int> class DH>
1669  double
1670  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1671  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PointFunctional(
1672  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1673  const std::map<std::string, const VECTOR*> &domain_values)
1674  {
1675  if (this->GetType() == "cost_functional")
1676  {
1677  // state values in quadrature points
1678  return GetFunctional()->PointValue(
1679  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1680  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1681  domain_values);
1682 
1683  } //endif cost_functional
1684  else if (this->GetType() == "aux_functional")
1685  {
1686  // state values in quadrature points
1687  return aux_functionals_[this->GetTypeNum()]->PointValue(
1688  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1689  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1690  domain_values);
1691 
1692  } //endif aux_functional
1693  else if (this->GetType() == "functional_for_ee")
1694  {
1695  // TODO is this correct? Should not be needed.
1696  return aux_functionals_[functional_for_ee_num_]->PointValue(
1697  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1698  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1699  domain_values);
1700  } //endif functional_for_ee
1701  else if (this->GetType().find("constraints") != std::string::npos)
1702  {
1703  return GetConstraints()->PointValue(
1704  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1705  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1706  domain_values);
1707 
1708  } //endif constraints
1709  else
1710  {
1711  throw DOpEException("Not implemented",
1712  "OptProblemContainer::PointFunctional");
1713  }
1714  }
1715 
1716  /******************************************************/
1717 
1718  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1719  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1720  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1721  template<int, int> class DH>
1722  template<typename FACEDATACONTAINER>
1723  double
1724  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1725  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryFunctional(
1726  const FACEDATACONTAINER& fdc)
1727  {
1728  if (this->GetType() == "cost_functional")
1729  {
1730  // state values in quadrature points
1731  return GetFunctional()->BoundaryValue(fdc);
1732  }
1733  else if (this->GetType() == "aux_functional")
1734  {
1735  // state values in quadrature points
1736  return aux_functionals_[this->GetTypeNum()]->BoundaryValue(fdc);
1737  }
1738  else if (this->GetType() == "functional_for_ee")
1739  // TODO is this correct? Should not be needed.
1740  {
1741  return aux_functionals_[functional_for_ee_num_]->BoundaryValue(fdc);
1742  }
1743  else if (this->GetType().find("constraints") != std::string::npos)
1744  {
1745  return GetConstraints()->BoundaryValue(fdc);
1746  }
1747  else
1748  {
1749  throw DOpEException("Not implemented",
1750  "OptProblemContainer::BoundaryFunctional");
1751  }
1752  }
1753 
1754  /******************************************************/
1755 
1756  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1757  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1758  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1759  template<int, int> class DH>
1760  template<typename FACEDATACONTAINER>
1761  double
1762  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1763  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceFunctional(
1764  const FACEDATACONTAINER& fdc)
1765  {
1766  if (this->GetType() == "cost_functional")
1767  {
1768  // state values in quadrature points
1769  return GetFunctional()->FaceValue(fdc);
1770  }
1771  else if (this->GetType() == "aux_functional")
1772  {
1773  // state values in quadrature points
1774  return aux_functionals_[this->GetTypeNum()]->FaceValue(fdc);
1775  }
1776  else if (this->GetType() == "functional_for_ee")
1777  // TODO is this correct? Should not be needed.
1778  {
1779  return aux_functionals_[functional_for_ee_num_]->FaceValue(fdc);
1780  }
1781  else
1782  {
1783  throw DOpEException("Not implemented",
1784  "OptProblemContainer::FaceFunctional");
1785  }
1786  }
1787 
1788  /******************************************************/
1789 
1790  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1791  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1792  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1793  template<int, int> class DH>
1794  double
1795  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1796  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AlgebraicFunctional(
1797  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1798  const std::map<std::string, const VECTOR*> &domain_values)
1799  {
1800  if (this->GetType() == "cost_functional")
1801  {
1802  // state values in quadrature points
1803  return GetFunctional()->AlgebraicValue(param_values, domain_values);
1804  }
1805  else if (this->GetType() == "aux_functional")
1806  {
1807  // state values in quadrature points
1808  return aux_functionals_[this->GetTypeNum()]->AlgebraicValue(
1809  param_values, domain_values);
1810  }
1811  else if (this->GetType() == "functional_for_ee")
1812  // TODO is this correct? Should not be needed.
1813  {
1814  return aux_functionals_[functional_for_ee_num_]->AlgebraicValue(
1815  param_values, domain_values);
1816  }
1817  else
1818  {
1819  throw DOpEException("Not implemented",
1820  "OptProblemContainer::AlgebraicFunctional");
1821  }
1822  }
1823 
1824  /******************************************************/
1825 
1826  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1827  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1828  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1829  template<int, int> class DH>
1830  template<typename DATACONTAINER>
1831  void
1832  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1833  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementEquation(
1834  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1835  double scale, double scale_ico)
1836  {
1837 
1838  if (this->GetType() == "state")
1839  {
1840  // state values in quadrature points
1841  this->GetPDE().ElementEquation(edc, local_vector, scale*interval_length_, scale_ico*interval_length_);
1842  }
1843  else if ((this->GetType() == "adjoint")
1844  || (this->GetType() == "adjoint_for_ee"))
1845  {
1846  // state values in quadrature points
1847  this->GetPDE().ElementEquation_U(edc, local_vector, scale*interval_length_,
1848  scale_ico*interval_length_);
1849  }
1850  else if (this->GetType() == "adjoint_hessian")
1851  {
1852  // state values in quadrature points
1853  this->GetPDE().ElementEquation_UTT(edc, local_vector, scale*interval_length_,
1854  scale_ico*interval_length_);
1855  }
1856  else if (this->GetType() == "tangent")
1857  {
1858  // state values in quadrature points
1859  this->GetPDE().ElementEquation_UT(edc, local_vector, scale*interval_length_,
1860  scale_ico*interval_length_);
1861  }
1862  else if ((this->GetType() == "gradient")
1863  || (this->GetType() == "hessian"))
1864  {
1865  // control values in quadrature points
1866  this->GetPDE().ControlElementEquation(edc, local_vector, scale*c_interval_length_);
1867  }
1868  else
1869  {
1870  throw DOpEException("Not implemented",
1871  "OptProblemContainer::ElementEquation");
1872  }
1873  }
1874 
1875  /******************************************************/
1876 
1877  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1878  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1879  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1880  template<int, int> class DH>
1881  void
1882  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1883  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AlgebraicResidual(
1884  VECTOR& residual,
1885  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1886  const std::map<std::string, const VECTOR*> &domain_values)
1887  {
1888  if (this->GetType() == "gradient")
1889  {
1890  // state values in quadrature points
1891  return GetFunctional()->AlgebraicGradient_Q(residual, param_values,
1892  domain_values);
1893  }
1894  else
1895  {
1896  throw DOpEException("Not implemented",
1897  "OptProblemContainer::AlgebraicFunctional");
1898  }
1899  }
1900 
1901  /******************************************************/
1902 
1903  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1904  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1905  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1906  template<int, int> class DH>
1907  template<typename DATACONTAINER>
1908  void
1909  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1910  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeEquation(
1911  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1912  double scale)
1913  {
1914 
1915  if (this->GetType() == "state")
1916  {
1917  this->GetPDE().ElementTimeEquation(edc, local_vector, scale);
1918  }
1919  else if (this->GetType() == "adjoint"
1920  || this->GetType() == "adjoint_for_ee")
1921  {
1922  this->GetPDE().ElementTimeEquation_U(edc, local_vector, scale);
1923  }
1924  else if (this->GetType() == "adjoint_hessian")
1925  {
1926  this->GetPDE().ElementTimeEquation_UTT(edc, local_vector, scale);
1927  }
1928  else if (this->GetType() == "tangent")
1929  {
1930  this->GetPDE().ElementTimeEquation_UT(edc, local_vector, scale);
1931  }
1932  else if ((this->GetType() == "gradient")
1933  || (this->GetType() == "hessian"))
1934  {
1935  throw DOpEException("Not implemented",
1936  "OptProblemContainer::ElementTimeEquation");
1937  }
1938  else
1939  {
1940  throw DOpEException("Not implemented",
1941  "OptProblemContainer::ElementTimeEquation");
1942  }
1943  }
1944 
1945  /******************************************************/
1946 
1947  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1948  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1949  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1950  template<int, int> class DH>
1951  template<typename DATACONTAINER>
1952  void
1953  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1954  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeEquationExplicit(
1955  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1956  double scale)
1957  {
1958 
1959  if (this->GetType() == "state")
1960  {
1961  this->GetPDE().ElementTimeEquationExplicit(edc, local_vector,
1962  scale);
1963  }
1964  else if (this->GetType() == "adjoint"
1965  || this->GetType() == "adjoint_for_ee")
1966  {
1967  this->GetPDE().ElementTimeEquationExplicit_U(edc, local_vector,
1968  scale);
1969  }
1970  else if (this->GetType() == "adjoint_hessian")
1971  {
1972  this->GetPDE().ElementTimeEquationExplicit_UTT(edc, local_vector,
1973  scale);
1974  }
1975  else if (this->GetType() == "tangent")
1976  {
1977  this->GetPDE().ElementTimeEquationExplicit_UT(edc, local_vector,
1978  scale);
1979  }
1980  else if ((this->GetType() == "gradient")
1981  || (this->GetType() == "hessian"))
1982  {
1983  throw DOpEException("Not implemented",
1984  "OptProblemContainer::ElementTimeEquationExplicit");
1985  }
1986  else
1987  {
1988  throw DOpEException("Not implemented",
1989  "OptProblemContainer::ElementTimeEquationExplicit");
1990  }
1991  }
1992 
1993  /******************************************************/
1994 
1995  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1996  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1997  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1998  template<int, int> class DH>
1999  template<typename FACEDATACONTAINER>
2000  void
2001  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2002  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceEquation(
2003  const FACEDATACONTAINER& fdc,
2004  dealii::Vector<double> &local_vector, double scale,
2005  double scale_ico)
2006  {
2007  if (this->GetType() == "state")
2008  {
2009  // state values in quadrature points
2010  this->GetPDE().FaceEquation(fdc, local_vector, scale*interval_length_, scale_ico*interval_length_);
2011  }
2012  else if (this->GetType() == "adjoint"
2013  || this->GetType() == "adjoint_for_ee")
2014  {
2015  // state values in quadrature points
2016  this->GetPDE().FaceEquation_U(fdc, local_vector, scale*interval_length_,
2017  scale_ico*interval_length_);
2018  }
2019  else if (this->GetType() == "adjoint_hessian")
2020  {
2021  // state values in quadrature points
2022  this->GetPDE().FaceEquation_UTT(fdc, local_vector, scale*interval_length_,
2023  scale_ico*interval_length_);
2024  }
2025  else if (this->GetType() == "tangent")
2026  {
2027  // state values in quadrature points
2028  this->GetPDE().FaceEquation_UT(fdc, local_vector, scale*interval_length_,
2029  scale_ico*interval_length_);
2030  }
2031 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
2032 // {
2033 // // control values in quadrature points
2034 // this->GetPDE().ControlFaceEquation(fdc, local_vector, scale);
2035 // }
2036  else
2037  {
2038  throw DOpEException("Not implemented",
2039  "OptProblemContainer::FaceEquation");
2040  }
2041  }
2042 
2043  /******************************************************/
2044 
2045  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2046  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2047  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2048  template<int, int> class DH>
2049  template<typename FACEDATACONTAINER>
2050  void
2051  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2052  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::InterfaceEquation(
2053  const FACEDATACONTAINER& fdc,
2054  dealii::Vector<double> &local_vector, double scale,
2055  double scale_ico)
2056  {
2057  if (this->GetType() == "state")
2058  {
2059  this->GetPDE().InterfaceEquation(fdc, local_vector, scale*interval_length_,
2060  scale_ico*interval_length_);
2061  }
2062  else if (this->GetType() == "adjoint"
2063  || this->GetType() == "adjoint_for_ee")
2064  {
2065  // state values in quadrature points
2066  this->GetPDE().InterfaceEquation_U(fdc, local_vector, scale*interval_length_,
2067  scale_ico*interval_length_);
2068  }
2069  else
2070  {
2071  throw DOpEException("Not implemented",
2072  "OptProblemContainer::InterfaceEquation");
2073  }
2074  }
2075  /******************************************************/
2076 
2077  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2078  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2079  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2080  template<int, int> class DH>
2081  template<typename FACEDATACONTAINER>
2082  void
2083  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2084  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryEquation(
2085  const FACEDATACONTAINER& fdc,
2086  dealii::Vector<double> &local_vector, double scale,
2087  double scale_ico)
2088  {
2089  if (this->GetType() == "state")
2090  {
2091  // state values in quadrature points
2092  this->GetPDE().BoundaryEquation(fdc, local_vector, scale*interval_length_,
2093  scale_ico*interval_length_);
2094  }
2095  else if (this->GetType() == "adjoint"
2096  || this->GetType() == "adjoint_for_ee")
2097  {
2098  // state values in quadrature points
2099  this->GetPDE().BoundaryEquation_U(fdc, local_vector, scale*interval_length_,
2100  scale_ico*interval_length_);
2101  }
2102  else if (this->GetType() == "adjoint_hessian")
2103  {
2104  // state values in quadrature points
2105  this->GetPDE().BoundaryEquation_UTT(fdc, local_vector, scale*interval_length_,
2106  scale_ico*interval_length_);
2107  }
2108  else if (this->GetType() == "tangent")
2109  {
2110  // state values in quadrature points
2111  this->GetPDE().BoundaryEquation_UT(fdc, local_vector, scale*interval_length_,
2112  scale_ico*interval_length_);
2113  }
2114  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
2115  {
2116  // control values in quadrature points
2117  this->GetPDE().ControlBoundaryEquation(fdc, local_vector, scale*c_interval_length_);
2118  }
2119  else
2120  {
2121  throw DOpEException("Not implemented",
2122  "OptProblemContainer::ElementBoundaryEquation");
2123  }
2124  }
2125 
2126  /******************************************************/
2127 
2128  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2129  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2130  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2131  template<int, int> class DH>
2132  template<typename DATACONTAINER>
2133  void
2134  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2135  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementRhs(
2136  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
2137  double scale)
2138  {
2139  if (this->GetType() == "state")
2140  {
2141  // state values in quadrature points
2142  this->GetPDE().ElementRightHandSide(edc, local_vector, scale*interval_length_);
2143  }
2144  else if (this->GetType() == "adjoint")
2145  {
2146  // state values in quadrature points
2147  if(GetFunctional()->NeedTime())
2148  {
2149  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2150  {
2151  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2152  {
2153  GetFunctional()->ElementValue_U(edc, local_vector, scale*interval_length_);
2154  }
2155  else // Otherwise always local if(GetFunctional()->GetType().find("timelocal") != std::string::npos)
2156  {
2157  GetFunctional()->ElementValue_U(edc, local_vector, scale);
2158  }
2159 
2160  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2161  {
2162  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2163  "OptProblemContainer::ElementRhs");
2164  }
2165  }
2166  }
2167  }
2168  else if (this->GetType() == "adjoint_for_ee")
2169  {
2170  //values of the derivative of the functional for error estimation
2171  if (aux_functionals_[functional_for_ee_num_]->NeedTime())
2172  {
2173  if (aux_functionals_[functional_for_ee_num_]->GetType().find("domain") != std::string::npos)
2174  {
2175  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos)
2176  {
2177  aux_functionals_[functional_for_ee_num_]->ElementValue_U(edc,
2178  local_vector, scale*interval_length_);
2179  }
2180  else
2181  {
2182  aux_functionals_[functional_for_ee_num_]->ElementValue_U(edc, local_vector, scale);
2183  }
2184  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos && aux_functionals_[functional_for_ee_num_]->GetType().find("timelocal") != std::string::npos)
2185  {
2186  throw DOpEException("Conflicting functional types: "+ aux_functionals_[functional_for_ee_num_]->GetType(),
2187  "OptProblemContainer::ElementRhs");
2188  }
2189  }
2190  }
2191  }
2192  else if (this->GetType() == "tangent")
2193  {
2194  // state values in quadrature points
2195  scale *= -1;
2196  this->GetPDE().ElementEquation_QT(edc, local_vector, scale*interval_length_, scale*interval_length_);
2197  }
2198  else if (this->GetType() == "adjoint_hessian")
2199  {
2200  // state values in quadrature points
2201  if(GetFunctional()->NeedTime())
2202  {
2203  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2204  {
2205  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2206  {
2207  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2208  {
2209  GetFunctional()->ElementValue_UU(edc, local_vector, scale*interval_length_);
2210  GetFunctional()->ElementValue_QU(edc, local_vector, scale*interval_length_);
2211  }
2212  else
2213  {
2214  GetFunctional()->ElementValue_UU(edc, local_vector, scale);
2215  GetFunctional()->ElementValue_QU(edc, local_vector, scale);
2216  }
2217  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2218  {
2219  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2220  "OptProblemContainer::ElementRhs");
2221  }
2222  }
2223  }
2224  }
2225  scale *= -1;
2226  this->GetPDE().ElementEquation_UU(edc, local_vector, scale*interval_length_, scale*interval_length_);
2227  this->GetPDE().ElementEquation_QU(edc, local_vector, scale*interval_length_, scale*interval_length_);
2228  //TODO: make some example where this realy matters to check if this is right
2229  this->GetPDE().ElementTimeEquationExplicit_UU(edc, local_vector,
2230  scale);
2231  }
2232  else if (this->GetType() == "gradient")
2233  {
2234  if (GetSpaceTimeHandler()->GetControlType()
2235  == DOpEtypes::ControlType::initial && initial_)
2236  {
2237  this->GetPDE().Init_ElementRhs_Q(edc, local_vector, scale);
2238  }
2239  // state values in quadrature points
2240  if(GetFunctional()->NeedTime())
2241  {
2242  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2243  {
2244  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2245  {
2246  GetFunctional()->ElementValue_Q(edc, local_vector, scale*interval_length_);
2247  }
2248  else
2249  {
2250  GetFunctional()->ElementValue_Q(edc, local_vector, scale);
2251  }
2252  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2253  {
2254  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2255  "OptProblemContainer::ElementRhs");
2256  }
2257  }
2258  }
2259  scale *= -1;
2260  this->GetPDE().ElementEquation_Q(edc, local_vector, scale*interval_length_, scale*interval_length_);
2261  }
2262  else if (this->GetType() == "hessian")
2263  {
2264  if (GetSpaceTimeHandler()->GetControlType()
2265  == DOpEtypes::ControlType::initial && initial_)
2266  {
2267  this->GetPDE().Init_ElementRhs_QTT(edc, local_vector, scale);
2268  this->GetPDE().Init_ElementRhs_QQ(edc, local_vector, scale);
2269  }
2270  if(GetFunctional()->NeedTime())
2271  {
2272  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2273  {
2274  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2275  {
2276  GetFunctional()->ElementValue_QQ(edc, local_vector, scale*interval_length_);
2277  GetFunctional()->ElementValue_UQ(edc, local_vector, scale*interval_length_);
2278  }
2279  else
2280  {
2281  GetFunctional()->ElementValue_QQ(edc, local_vector, scale);
2282  GetFunctional()->ElementValue_UQ(edc, local_vector, scale);
2283  }
2284  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2285  {
2286  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2287  "OptProblemContainer::ElementRhs");
2288  }
2289  }
2290  }
2291 
2292  scale *= -1;
2293  this->GetPDE().ElementEquation_QTT(edc, local_vector, scale*interval_length_, scale*interval_length_);
2294  this->GetPDE().ElementEquation_UQ(edc, local_vector, scale*interval_length_, scale*interval_length_);
2295  this->GetPDE().ElementEquation_QQ(edc, local_vector, scale*interval_length_, scale*interval_length_);
2296  }
2297  else if (this->GetType() == "global_constraint_gradient")
2298  {
2299  assert(interval_length_==1.);
2300  GetConstraints()->ElementValue_Q(edc, local_vector, scale);
2301  }
2302  else if (this->GetType() == "global_constraint_hessian")
2303  {
2304  assert(interval_length_==1.);
2305  GetConstraints()->ElementValue_QQ(edc, local_vector, scale);
2306  }
2307  else
2308  {
2309  throw DOpEException("Not implemented",
2310  "OptProblemContainer::ElementRhs");
2311  }
2312 
2313  }
2314 
2315  /******************************************************/
2316 
2317  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2318  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2319  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2320  template<int, int> class DH>
2321  void
2322  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2323  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PointRhs(
2324  const std::map<std::string, const dealii::Vector<double>*> &param_values,
2325  const std::map<std::string, const VECTOR*> &domain_values,
2326  VECTOR& rhs_vector, double scale)
2327  {
2328 
2329  if (this->GetType() == "adjoint")
2330  {
2331  // state values in quadrature points
2332  if(GetFunctional()->NeedTime())
2333  {
2334  if (GetFunctional()->GetType().find("point") != std::string::npos)
2335  {
2336  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2337  {
2338  GetFunctional()->PointValue_U(
2339  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2340  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2341  domain_values, rhs_vector, scale*interval_length_);
2342  }
2343  else
2344  {
2345  GetFunctional()->PointValue_U(
2346  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2347  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2348  domain_values, rhs_vector, scale);
2349  }
2350  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2351  {
2352  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2353  "OptProblemContainer::PointRhs");
2354  }
2355  }
2356  }
2357  }
2358  else if (this->GetType() == "adjoint_for_ee")
2359  {
2360  //values of the derivative of the functional for error estimation
2361  if (aux_functionals_[functional_for_ee_num_]->NeedTime())
2362  {
2363  if (aux_functionals_[functional_for_ee_num_]->GetType().find("point") != std::string::npos)
2364  {
2365  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos)
2366  {
2367  aux_functionals_[functional_for_ee_num_]->PointValue_U(
2368  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2369  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2370  domain_values, rhs_vector, scale*interval_length_);
2371  }
2372  else
2373  {
2374  aux_functionals_[functional_for_ee_num_]->PointValue_U(
2375  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2376  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2377  domain_values, rhs_vector, scale);
2378  }
2379  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos && aux_functionals_[functional_for_ee_num_]->GetType().find("timelocal") != std::string::npos)
2380  {
2381  throw DOpEException("Conflicting functional types: "+ aux_functionals_[functional_for_ee_num_]->GetType(),
2382  "OptProblemContainer::PointRhs");
2383  }
2384  }
2385  }
2386  }
2387  else if (this->GetType() == "adjoint_hessian")
2388  {
2389  // state values in quadrature points
2390  if(GetFunctional()->NeedTime())
2391  {
2392  if (GetFunctional()->GetType().find("point") != std::string::npos)
2393  {
2394  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2395  {
2396  GetFunctional()->PointValue_UU(
2397  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2398  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2399  domain_values, rhs_vector, scale*interval_length_);
2400  GetFunctional()->PointValue_QU(
2401  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2402  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2403  domain_values, rhs_vector, scale*interval_length_);
2404  }
2405  else
2406  {
2407  GetFunctional()->PointValue_UU(
2408  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2409  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2410  domain_values, rhs_vector, scale);
2411  GetFunctional()->PointValue_QU(
2412  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2413  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2414  domain_values, rhs_vector, scale);
2415  }
2416  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2417  {
2418  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2419  "OptProblemContainer::PointRhs");
2420  }
2421  }
2422  }
2423  }
2424  else if (this->GetType() == "gradient")
2425  {
2426  // state values in quadrature points
2427  if(GetFunctional()->NeedTime())
2428  {
2429  if (GetFunctional()->GetType().find("point") != std::string::npos)
2430  {
2431  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2432  {
2433  GetFunctional()->PointValue_Q(
2434  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2435  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2436  domain_values, rhs_vector, scale*interval_length_);
2437  }
2438  else
2439  {
2440  GetFunctional()->PointValue_Q(
2441  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2442  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2443  domain_values, rhs_vector, scale);
2444  }
2445  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2446  {
2447  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2448  "OptProblemContainer::PointRhs");
2449  }
2450  }
2451  }
2452  }
2453  else if (this->GetType() == "hessian")
2454  {
2455  // state values in quadrature points
2456  if(GetFunctional()->NeedTime())
2457  {
2458  if (GetFunctional()->GetType().find("point") != std::string::npos)
2459  {
2460  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2461  {
2462  GetFunctional()->PointValue_QQ(
2463  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2464  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2465  domain_values, rhs_vector, scale*interval_length_);
2466  GetFunctional()->PointValue_UQ(
2467  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2468  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2469  domain_values, rhs_vector, scale*interval_length_);
2470  }
2471  else
2472  {
2473  GetFunctional()->PointValue_QQ(
2474  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2475  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2476  domain_values, rhs_vector, scale);
2477  GetFunctional()->PointValue_UQ(
2478  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2479  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2480  domain_values, rhs_vector, scale);
2481  }
2482  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2483  {
2484  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2485  "OptProblemContainer::PointRhs");
2486  }
2487  }
2488  }
2489  }
2490  else
2491  {
2492  throw DOpEException("Not implemented", "OptProblem::PointRhs");
2493  }
2494 
2495  }
2496 
2497  /******************************************************/
2498 
2499  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2500  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2501  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2502  template<int, int> class DH>
2503  template<typename FACEDATACONTAINER>
2504  void
2505  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2506  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceRhs(
2507  const FACEDATACONTAINER& fdc,
2508  dealii::Vector<double> &local_vector, double scale)
2509  {
2510 
2511  if (this->GetType() == "state")
2512  {
2513  // state values in face quadrature points
2514  this->GetPDE().FaceRightHandSide(fdc, local_vector, scale*interval_length_);
2515  }
2516  else if (this->GetType() == "adjoint")
2517  {
2518  // state values in quadrature points
2519  if(GetFunctional()->NeedTime())
2520  {
2521  if (GetFunctional()->GetType().find("face") != std::string::npos)
2522  {
2523  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2524  {
2525  GetFunctional()->FaceValue_U(fdc, local_vector, scale*interval_length_);
2526  }
2527  else
2528  {
2529  GetFunctional()->FaceValue_U(fdc, local_vector, scale);
2530  }
2531 
2532  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2533  {
2534  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2535  "OptProblemContainer::FaceRhs");
2536  }
2537  }
2538  }
2539  }
2540  else if (this->GetType() == "adjoint_for_ee")
2541  {
2542  //values of the derivative of the functional for error estimation
2543  if (aux_functionals_[functional_for_ee_num_]->NeedTime())
2544  {
2545  if (aux_functionals_[functional_for_ee_num_]->GetType().find("face")
2546  != std::string::npos)
2547  {
2548  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos)
2549  {
2550  aux_functionals_[functional_for_ee_num_]->FaceValue_U(fdc,
2551  local_vector, scale*interval_length_);
2552  }
2553  else
2554  {
2555  aux_functionals_[functional_for_ee_num_]->FaceValue_U(fdc,
2556  local_vector, scale);
2557  }
2558  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos && aux_functionals_[functional_for_ee_num_]->GetType().find("timelocal") != std::string::npos)
2559  {
2560  throw DOpEException("Conflicting functional types: "+ aux_functionals_[functional_for_ee_num_]->GetType(),
2561  "OptProblemContainer::FaceRhs");
2562  }
2563  }
2564  }
2565 
2566  }
2567  else if (this->GetType() == "tangent")
2568  {
2569  // state values in quadrature points
2570  scale *= -1;
2571  this->GetPDE().FaceEquation_QT(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2572  }
2573  else if (this->GetType() == "adjoint_hessian")
2574  {
2575  // state values in quadrature points
2576  if(GetFunctional()->NeedTime())
2577  {
2578  if (GetFunctional()->GetType().find("face") != std::string::npos)
2579  {
2580  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2581  {
2582  GetFunctional()->FaceValue_UU(fdc, local_vector, scale*interval_length_);
2583  GetFunctional()->FaceValue_QU(fdc, local_vector, scale*interval_length_);
2584  }
2585  else
2586  {
2587  GetFunctional()->FaceValue_UU(fdc, local_vector, scale);
2588  GetFunctional()->FaceValue_QU(fdc, local_vector, scale);
2589  }
2590  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2591  {
2592  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2593  "OptProblemContainer::FaceRhs");
2594  }
2595  }
2596  }
2597 
2598  scale *= -1;
2599  this->GetPDE().FaceEquation_UU(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2600  this->GetPDE().FaceEquation_QU(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2601  }
2602  else if (this->GetType() == "gradient")
2603  {
2604  // state values in quadrature points
2605  if(GetFunctional()->NeedTime())
2606  {
2607  if (GetFunctional()->GetType().find("face") != std::string::npos)
2608  {
2609  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2610  {
2611  GetFunctional()->FaceValue_Q(fdc, local_vector, scale*interval_length_);
2612  }
2613  else
2614  {
2615  GetFunctional()->FaceValue_Q(fdc, local_vector, scale);
2616  }
2617  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2618  {
2619  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2620  "OptProblemContainer::FaceRhs");
2621  }
2622  }
2623  }
2624 
2625  scale *= -1;
2626  this->GetPDE().FaceEquation_Q(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2627  }
2628  else if (this->GetType() == "hessian")
2629  {
2630  // state values in quadrature points
2631  if(GetFunctional()->NeedTime())
2632  {
2633  if (GetFunctional()->GetType().find("face") != std::string::npos)
2634  {
2635  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2636  {
2637  GetFunctional()->FaceValue_QQ(fdc, local_vector, scale*interval_length_);
2638  GetFunctional()->FaceValue_UQ(fdc, local_vector, scale*interval_length_);
2639  }
2640  else
2641  {
2642  GetFunctional()->FaceValue_QQ(fdc, local_vector, scale);
2643  GetFunctional()->FaceValue_UQ(fdc, local_vector, scale);
2644  }
2645  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2646  {
2647  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2648  "OptProblemContainer::FaceRhs");
2649  }
2650  }
2651  }
2652 
2653  scale *= -1;
2654  this->GetPDE().FaceEquation_QTT(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2655  this->GetPDE().FaceEquation_UQ(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2656  this->GetPDE().FaceEquation_QQ(fdc, local_vector, scale*interval_length_, scale*interval_length_);
2657  }
2658  else
2659  {
2660  throw DOpEException("Not implemented",
2661  "OptProblemContainer::FaceRhs");
2662  }
2663 
2664  }
2665 
2666  /******************************************************/
2667 
2668  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2669  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2670  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2671  template<int, int> class DH>
2672  template<typename FACEDATACONTAINER>
2673  void
2674  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2675  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryRhs(
2676  const FACEDATACONTAINER& fdc,
2677  dealii::Vector<double> &local_vector, double scale)
2678  {
2679 
2680  if (this->GetType() == "state")
2681  {
2682  // state values in face quadrature points
2683  this->GetPDE().BoundaryRightHandSide(fdc, local_vector, scale*interval_length_);
2684  }
2685  else if (this->GetType() == "adjoint")
2686  {
2687  // state values in quadrature points
2688  if(GetFunctional()->NeedTime())
2689  {
2690  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2691  {
2692  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2693  {
2694  GetFunctional()->BoundaryValue_U(fdc, local_vector, scale*interval_length_);
2695  }
2696  else
2697  {
2698  GetFunctional()->BoundaryValue_U(fdc, local_vector, scale);
2699  }
2700  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2701  {
2702  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2703  "OptProblemContainer::BoundaryRhs");
2704  }
2705  }
2706  }
2707  }
2708  else if (this->GetType() == "adjoint_for_ee")
2709  {
2710  //values of the derivative of the functional for error estimation
2711  if (aux_functionals_[functional_for_ee_num_]->NeedTime())
2712  {
2713  if (aux_functionals_[functional_for_ee_num_]->GetType().find(
2714  "boundary") != std::string::npos)
2715  {
2716  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos)
2717  {
2718  aux_functionals_[functional_for_ee_num_]->BoundaryValue_U(fdc,
2719  local_vector, scale*interval_length_);
2720  }
2721  else
2722  {
2723  aux_functionals_[functional_for_ee_num_]->BoundaryValue_U(fdc,
2724  local_vector, scale);
2725  }
2726  if(aux_functionals_[functional_for_ee_num_]->GetType().find("timedistributed") != std::string::npos && aux_functionals_[functional_for_ee_num_]->GetType().find("timelocal") != std::string::npos)
2727  {
2728  throw DOpEException("Conflicting functional types: "+ aux_functionals_[functional_for_ee_num_]->GetType(),
2729  "OptProblemContainer::BoundaryRhs");
2730  }
2731  }
2732  }
2733  }
2734  else if (this->GetType() == "tangent")
2735  {
2736  // state values in quadrature points
2737  scale *= -1;
2738  this->GetPDE().BoundaryEquation_QT(fdc, local_vector, scale*interval_length_,
2739  scale*interval_length_);
2740  }
2741  else if (this->GetType() == "adjoint_hessian")
2742  {
2743  // state values in quadrature points
2744  if(GetFunctional()->NeedTime())
2745  {
2746  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2747  {
2748  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2749  {
2750  GetFunctional()->BoundaryValue_UU(fdc, local_vector, scale*interval_length_);
2751  GetFunctional()->BoundaryValue_QU(fdc, local_vector, scale*interval_length_);
2752  }
2753  else
2754  {
2755  GetFunctional()->BoundaryValue_UU(fdc, local_vector, scale);
2756  GetFunctional()->BoundaryValue_QU(fdc, local_vector, scale);
2757  }
2758  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2759  {
2760  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2761  "OptProblemContainer::BoundaryRhs");
2762  }
2763  }
2764  }
2765 
2766  scale *= -1;
2767  this->GetPDE().BoundaryEquation_UU(fdc, local_vector, scale*interval_length_,
2768  scale*interval_length_);
2769  this->GetPDE().BoundaryEquation_QU(fdc, local_vector, scale*interval_length_,
2770  scale*interval_length_);
2771  }
2772  else if (this->GetType() == "gradient")
2773  {
2774  // state values in quadrature points
2775  if(GetFunctional()->NeedTime())
2776  {
2777  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2778  {
2779  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2780  {
2781  GetFunctional()->BoundaryValue_Q(fdc, local_vector, scale*interval_length_);
2782  }
2783  else
2784  {
2785  GetFunctional()->BoundaryValue_Q(fdc, local_vector, scale);
2786  }
2787  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2788  {
2789  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2790  "OptProblemContainer::BoundaryRhs");
2791  }
2792  }
2793  }
2794 
2795  scale *= -1;
2796  this->GetPDE().BoundaryEquation_Q(fdc, local_vector, scale*interval_length_,
2797  scale*interval_length_);
2798  }
2799  else if (this->GetType() == "hessian")
2800  {
2801  // state values in quadrature points
2802  if(GetFunctional()->NeedTime())
2803  {
2804  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2805  {
2806  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos)
2807  {
2808  GetFunctional()->BoundaryValue_QQ(fdc, local_vector, scale*interval_length_);
2809  GetFunctional()->BoundaryValue_UQ(fdc, local_vector, scale*interval_length_);
2810  }
2811  else
2812  {
2813  GetFunctional()->BoundaryValue_QQ(fdc, local_vector, scale);
2814  GetFunctional()->BoundaryValue_UQ(fdc, local_vector, scale);
2815  }
2816  if(GetFunctional()->GetType().find("timedistributed") != std::string::npos && GetFunctional()->GetType().find("timelocal") != std::string::npos)
2817  {
2818  throw DOpEException("Conflicting functional types: "+ GetFunctional()->GetType(),
2819  "OptProblemContainer::BoundaryRhs");
2820  }
2821  }
2822  }
2823 
2824  scale *= -1;
2825  this->GetPDE().BoundaryEquation_QTT(fdc, local_vector, scale*interval_length_,
2826  scale*interval_length_);
2827  this->GetPDE().BoundaryEquation_UQ(fdc, local_vector, scale*interval_length_,
2828  scale*interval_length_);
2829  this->GetPDE().BoundaryEquation_QQ(fdc, local_vector, scale*interval_length_,
2830  scale*interval_length_);
2831  }
2832  else
2833  {
2834  throw DOpEException("Not implemented",
2835  "OptProblemContainer::BoundaryRhs");
2836  }
2837 
2838  }
2839 
2840  /******************************************************/
2841 
2842  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2843  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2844  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2845  template<int, int> class DH>
2846  template<typename DATACONTAINER>
2847  void
2848  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2849  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementMatrix(
2850  const DATACONTAINER& edc,
2851  dealii::FullMatrix<double> &local_entry_matrix, double scale,
2852  double scale_ico)
2853  {
2854 
2855  if (this->GetType() == "state" || this->GetType() == "tangent")
2856  {
2857  // state values in quadrature points
2858  this->GetPDE().ElementMatrix(edc, local_entry_matrix, scale*interval_length_, scale_ico*interval_length_);
2859  }
2860  else if (this->GetType() == "adjoint"
2861  || this->GetType() == "adjoint_for_ee"
2862  || this->GetType() == "adjoint_hessian")
2863  {
2864  // state values in quadrature points
2865  this->GetPDE().ElementMatrix_T(edc, local_entry_matrix, scale*interval_length_,
2866  scale_ico*interval_length_);
2867  }
2868  else if ((this->GetType() == "gradient")
2869  || (this->GetType() == "hessian"))
2870  {
2871  // control values in quadrature points
2872  this->GetPDE().ControlElementMatrix(edc, local_entry_matrix, scale*c_interval_length_);
2873  }
2874  else
2875  {
2876  throw DOpEException("Not implemented",
2877  "OptProblemContainer::ElementMatrix");
2878  }
2879 
2880  }
2881 
2882  /******************************************************/
2883 
2884  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2885  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2886  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2887  template<int, int> class DH>
2888  template<typename DATACONTAINER>
2889  void
2890  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2891  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeMatrix(
2892  const DATACONTAINER& edc, FullMatrix<double> &local_entry_matrix)
2893  {
2894 
2895  if (this->GetType() == "state" || this->GetType() == "tangent")
2896  {
2897  // state values in quadrature points
2898  this->GetPDE().ElementTimeMatrix(edc, local_entry_matrix);
2899  }
2900  else if (this->GetType() == "adjoint"
2901  || this->GetType() == "adjoint_for_ee"
2902  || this->GetType() == "adjoint_hessian")
2903  {
2904  this->GetPDE().ElementTimeMatrix_T(edc, local_entry_matrix);
2905  }
2906  else if ((this->GetType() == "gradient")
2907  || (this->GetType() == "hessian"))
2908  {
2909  throw DOpEException("Not implemented",
2910  "OptProblemContainer::ElementTimeMatrix");
2911  }
2912  else
2913  {
2914  throw DOpEException("Not implemented",
2915  "OptProblemContainer::ElementTimeMatrix");
2916  }
2917 
2918  }
2919 
2920  /******************************************************/
2921 
2922  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2923  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2924  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2925  template<int, int> class DH>
2926  template<typename DATACONTAINER>
2927  void
2928  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2929  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeMatrixExplicit(
2930  const DATACONTAINER& edc,
2931  dealii::FullMatrix<double> &local_entry_matrix)
2932  {
2933 
2934  if (this->GetType() == "state" || this->GetType() == "tangent")
2935  {
2936  this->GetPDE().ElementTimeMatrixExplicit(edc, local_entry_matrix);
2937  }
2938  else if (this->GetType() == "adjoint"
2939  || this->GetType() == "adjoint_hessian")
2940  {
2941  this->GetPDE().ElementTimeMatrixExplicit_T(edc, local_entry_matrix);
2942  }
2943  else if ((this->GetType() == "gradient")
2944  || (this->GetType() == "hessian"))
2945  {
2946  throw DOpEException("Not implemented",
2947  "OptProblemContainer::ElementTimeMatrixExplicit");
2948  }
2949  else
2950  {
2951  throw DOpEException("Not implemented",
2952  "OptProblemContainer::ElementTimeMatrixExplicit");
2953  }
2954 
2955  }
2956 
2957  /******************************************************/
2958 
2959  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2960  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2961  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2962  template<int, int> class DH>
2963  template<typename FACEDATACONTAINER>
2964  void
2965  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2966  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceMatrix(
2967  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
2968  double scale, double scale_ico)
2969  {
2970  if (this->GetType() == "state" || this->GetType() == "tangent")
2971  {
2972  // state values in face quadrature points
2973  this->GetPDE().FaceMatrix(fdc, local_entry_matrix, scale*interval_length_, scale_ico*interval_length_);
2974  }
2975  else if (this->GetType() == "adjoint"
2976  || this->GetType() == "adjoint_for_ee"
2977  || this->GetType() == "adjoint_hessian")
2978  {
2979  // state values in quadrature points
2980  this->GetPDE().FaceMatrix_T(fdc, local_entry_matrix, scale*interval_length_,
2981  scale_ico*interval_length_);
2982  }
2983 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
2984 // {
2985 // // control values in quadrature points
2986 // this->GetPDE().ControlFaceMatrix(fdc, local_entry_matrix);
2987 // }
2988  else
2989  {
2990  throw DOpEException("Not implemented",
2991  "OptProblemContainer::NewtonFaceMatrix");
2992  }
2993 
2994  }
2995 
2996  /******************************************************/
2997 
2998  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2999  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3000  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3001  template<int, int> class DH>
3002  template<typename FACEDATACONTAINER>
3003  void
3004  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
3005  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::InterfaceMatrix(
3006  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
3007  double scale, double scale_ico)
3008  {
3009  if (this->GetType() == "state")
3010  {
3011  this->GetPDE().InterfaceMatrix(fdc, local_entry_matrix, scale*interval_length_,
3012  scale_ico*interval_length_);
3013  }
3014  else if (this->GetType() == "adjoint"
3015  || this->GetType() == "adjoint_for_ee"
3016  || this->GetType() == "adjoint_hessian")
3017  {
3018  this->GetPDE().InterfaceMatrix_T(fdc, local_entry_matrix, scale*interval_length_,
3019  scale_ico*interval_length_);
3020  }
3021  else
3022  {
3023  throw DOpEException("Not implemented",
3024  "OptProblemContainer::NewtonInterfaceMatrix");
3025  }
3026  }
3027 
3028  /******************************************************/
3029 
3030  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3031  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3032  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3033  template<int, int> class DH>
3034  template<typename FACEDATACONTAINER>
3035  void
3036  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
3037  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryMatrix(
3038  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_matrix,
3039  double scale, double scale_ico)
3040  {
3041  if (this->GetType() == "state" || this->GetType() == "tangent")
3042  {
3043  // state values in face quadrature points
3044  this->GetPDE().BoundaryMatrix(fdc, local_matrix, scale*interval_length_,
3045  scale_ico*interval_length_);
3046  }
3047  else if (this->GetType() == "adjoint"
3048  || this->GetType() == "adjoint_for_ee"
3049  || this->GetType() == "adjoint_hessian")
3050  {
3051  // state values in quadrature points
3052  this->GetPDE().BoundaryMatrix_T(fdc, local_matrix, scale*interval_length_,
3053  scale_ico*interval_length_);
3054  }
3055  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
3056  {
3057  // control values in quadrature points
3058  this->GetPDE().ControlBoundaryMatrix(fdc, local_matrix, scale*c_interval_length_);
3059  }
3060  else
3061  {
3062  throw DOpEException("Not implemented",
3063  "OptProblemContainer::ElementBoundaryMatrix");
3064  }
3065 
3066  }
3067 
3068  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3069  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3070  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3071  template<int, int> class DH>
3072  void
3073  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3074  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ComputeLocalControlConstraints(
3075  VECTOR& constraints,
3076  const std::map<std::string, const dealii::Vector<double>*> &/*values*/,
3077  const std::map<std::string, const VECTOR*> &block_values)
3078  {
3079  if (this->GetType() == "constraints")
3080  {
3081  if (this->GetSpaceTimeHandler()->GetNLocalConstraints() != 0)
3082  {
3083  const VECTOR& control = *GetBlockVector(block_values, "control");
3084  this->GetConstraints()->EvaluateLocalControlConstraints(control,
3085  constraints);
3086  }
3087  }
3088  else
3089  {
3090  throw DOpEException("Wrong problem type" + this->GetType(),
3091  "OptProblemContainer::ComputeLocalConstraints");
3092  }
3093  }
3094 
3095  /******************************************************/
3096 
3097  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3098  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3099  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3100  template<int, int> class DH>
3101  std::string
3102  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3103  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFType() const
3104  {
3105  if (this->GetType() == "state" || this->GetType() == "adjoint"
3106  || this->GetType() == "adjoint_for_ee" || this->GetType() == "tangent"
3107  || this->GetType() == "adjoint_hessian")
3108  {
3109  return "state";
3110  }
3111  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
3112  || (this->GetType() == "hessian_inverse"))
3113  {
3114  return "control";
3115  }
3116  else
3117  {
3118  throw DOpEException("Unknown Type:" + this->GetType(),
3119  "OptProblemContainer::GetDoFType");
3120  }
3121  }
3122 
3123  /******************************************************/
3124 
3125  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3126  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3127  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3128  template<int, int> class DH>
3129  const FE<dealdim, dealdim>&
3130  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3131  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFESystem() const
3132  {
3133  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3134  || this->GetType() == "adjoint_for_ee" || this->GetType() == "tangent"
3135  || this->GetType() == "adjoint_hessian")
3136  {
3137  return this->GetSpaceTimeHandler()->GetFESystem("state");
3138  }
3139  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
3140  || (this->GetType() == "global_constraint_gradient"))
3141  {
3142 #if dope_dimension > 0
3143  if(dopedim == dealdim)
3144  return this->GetSpaceTimeHandler()->GetFESystem("control");
3145  else
3146  throw DOpEException("Non matching dimensions!","OptProblemContainer::GetFESystem");
3147 #else
3148  return this->GetSpaceTimeHandler()->GetFESystem("state");
3149 #endif
3150  }
3151  else
3152  {
3153  throw DOpEException("Unknown Type:" + this->GetType(),
3154  "OptProblemContainer::GetFESystem");
3155  }
3156  }
3157 
3158  /******************************************************/
3159 
3160  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3161  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3162  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3163  template<int, int> class DH>
3164  UpdateFlags
3165  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3166  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetUpdateFlags() const
3167  {
3168 
3169  UpdateFlags r;
3170  if (this->GetType().find("aux_functional") != std::string::npos)
3171  {
3172  r = aux_functionals_[this->GetTypeNum()]->GetUpdateFlags();
3173  }
3174  else if (this->GetType().find("functional") != std::string::npos)
3175  {
3176  r = this->GetFunctional()->GetUpdateFlags();
3177  }
3178  else if (this->GetType().find("constraints") != std::string::npos)
3179  {
3180  r = this->GetConstraints()->GetUpdateFlags();
3181  }
3182  else if (this->GetType() == "functional_for_ee")
3183  {
3184  r = aux_functionals_[functional_for_ee_num_]->GetUpdateFlags();
3185  }
3186  else
3187  {
3188  r = this->GetPDE().GetUpdateFlags();
3189  if (this->GetType() == "adjoint_hessian" || this->GetType() == "adjoint"
3190  || (this->GetType() == "hessian")
3191  || (this->GetType() == "gradient"))
3192  {
3193  r = r | this->GetFunctional()->GetUpdateFlags();
3194  }
3195  else if (this->GetType() == "adjoint_for_ee")
3196  {
3197  r = r | aux_functionals_[functional_for_ee_num_]->GetUpdateFlags();
3198  }
3199  }
3200  return r | update_JxW_values;
3201  }
3202 
3203  /******************************************************/
3204 
3205  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3206  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3207  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3208  template<int, int> class DH>
3209  UpdateFlags
3210  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3211  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFaceUpdateFlags() const
3212  {
3213  UpdateFlags r;
3214  if (this->GetType().find("aux_functional") != std::string::npos)
3215  {
3216  r = aux_functionals_[this->GetTypeNum()]->GetFaceUpdateFlags();
3217  }
3218  else if (this->GetType().find("functional") != std::string::npos)
3219  {
3220  r = this->GetFunctional()->GetFaceUpdateFlags();
3221  }
3222  else if (this->GetType().find("constraints") != std::string::npos)
3223  {
3224  r = this->GetConstraints()->GetFaceUpdateFlags();
3225  }
3226  else if (this->GetType() == "functional_for_ee")
3227  {
3228  r = aux_functionals_[functional_for_ee_num_]->GetUpdateFlags();
3229  }
3230  else
3231  {
3232  r = this->GetPDE().GetFaceUpdateFlags();
3233  if (this->GetType() == "adjoint_hessian" || this->GetType() == "adjoint"
3234  || (this->GetType() == "hessian"))
3235  {
3236  r = r | this->GetFunctional()->GetFaceUpdateFlags();
3237  }
3238  else if (this->GetType() == "adjoint_for_ee")
3239  {
3240  r = r
3241  | aux_functionals_[functional_for_ee_num_]->GetFaceUpdateFlags();
3242  }
3243  }
3244  return r | update_JxW_values;
3245  }
3246 
3247  /******************************************************/
3248 
3249  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3250  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3251  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3252  template<int, int> class DH>
3253  std::string
3254  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3255  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctionalType() const
3256  {
3257  if (this->GetType() == "aux_functional")
3258  {
3259  return aux_functionals_[this->GetTypeNum()]->GetType();
3260  }
3261  else if (this->GetType() == "functional_for_ee")
3262  {
3263  return aux_functionals_[functional_for_ee_num_]->GetType();
3264  }
3265  return GetFunctional()->GetType();
3266  }
3267 
3268  /******************************************************/
3269 
3270  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3271  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3272  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3273  template<int, int> class DH>
3274  std::string
3275  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3276  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctionalName() const
3277  {
3278  if (this->GetType() == "aux_functional")
3279  {
3280  return aux_functionals_[this->GetTypeNum()]->GetName();
3281  }
3282  else if (this->GetType() == "functional_for_ee")
3283  {
3284  return aux_functionals_[functional_for_ee_num_]->GetName();
3285  }
3286  return GetFunctional()->GetName();
3287  }
3288 
3289  /******************************************************/
3290 
3291  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3292  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3293  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3294  template<int, int> class DH>
3295  std::string
3296  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3297  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetConstraintType() const
3298  {
3299  return GetConstraints()->GetType();
3300  }
3301 
3302  /******************************************************/
3303 
3304  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3305  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3306  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3307  template<int, int> class DH>
3308  void
3309  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3310  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetTime(double time,
3311  unsigned int time_dof_number,
3312  const TimeIterator& interval, bool initial)
3313  {
3314  GetSpaceTimeHandler()->SetInterval(interval);
3315  initial_ = initial;
3316  // GetSpaceTimeHandler()->SetTimeDoFNumber(time_point);
3317  interval_length_ = GetSpaceTimeHandler()->GetStepSize();
3318 
3319  if(GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::initial || GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::stationary )
3320  {
3321  c_interval_length_ = 1.;
3322  }
3323  else
3324  {
3325  c_interval_length_ = GetSpaceTimeHandler()->GetStepSize();
3326  }
3327 
3328  { //Zeit an Dirichlet Werte uebermitteln
3329  for (unsigned int i = 0;
3330  i < transposed_control_gradient_dirichlet_values_.size(); i++)
3331  transposed_control_gradient_dirichlet_values_[i]->SetTime(time);
3332  for (unsigned int i = 0;
3333  i < transposed_control_hessian_dirichlet_values_.size(); i++)
3334  transposed_control_hessian_dirichlet_values_[i]->SetTime(time);
3335  for (unsigned int i = 0; i < primal_dirichlet_values_.size(); i++)
3336  primal_dirichlet_values_[i]->SetTime(time);
3337  for (unsigned int i = 0; i < tangent_dirichlet_values_.size(); i++)
3338  tangent_dirichlet_values_[i]->SetTime(time);
3339  for (unsigned int i = 0; i < control_dirichlet_values_.size(); i++)
3340  control_dirichlet_values_[i]->SetTime(time);
3341  //Functionals
3342  GetFunctional()->SetTime(time);
3343  for (unsigned int i = 0; i < aux_functionals_.size(); i++)
3344  aux_functionals_[i]->SetTime(time);
3345  //PDE
3346  this->GetPDE().SetTime(time);
3347  }
3348  //Update Auxiliary Control, State and Constraint Vectors
3349  {
3350  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
3351  auxiliary_controls_.begin();
3352  for (; it != auxiliary_controls_.end(); it++)
3353  {
3354  it->second->SetTimeDoFNumber(time_dof_number);
3355  }
3356  }
3357  {
3358  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
3359  auxiliary_state_.begin();
3360  for (; it != auxiliary_state_.end(); it++)
3361  {
3362  it->second->SetTimeDoFNumber(time_dof_number, interval);
3363  }
3364  }
3365  {
3366  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
3367  auxiliary_constraints_.begin();
3368  for (; it != auxiliary_constraints_.end(); it++)
3369  {
3370  it->second->SetTimeDoFNumber(time_dof_number);
3371  }
3372  }
3373 
3374  }
3375 
3376  /******************************************************/
3377 
3378  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3379  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3380  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3381  template<int, int> class DH>
3382  void
3383  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3384  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ComputeSparsityPattern(
3385  SPARSITYPATTERN & sparsity) const
3386  {
3387  if (this->GetType() == "state" || this->GetType() == "tangent"
3388  || this->GetType() == "adjoint_for_ee" || this->GetType() == "adjoint"
3389  || this->GetType() == "adjoint_hessian")
3390  {
3391  this->GetSpaceTimeHandler()->ComputeStateSparsityPattern(sparsity);
3392  }
3393  else if ((this->GetType() == "gradient")
3394  || (this->GetType() == "hessian"))
3395  {
3396 #if dope_dimension > 0
3397  this->GetSpaceTimeHandler()->ComputeControlSparsityPattern(sparsity);
3398 #else
3399  throw DOpEException("Wrong dimension",
3400  "OptProblemContainer::ComputeSparsityPattern");
3401 #endif
3402  }
3403  else
3404  {
3405  throw DOpEException("Unknown type " + this->GetType(),
3406  "OptProblemContainer::ComputeSparsityPattern");
3407  }
3408  }
3409 
3410  /******************************************************/
3411 
3412  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3413  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3414  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3415  template<int, int> class DH>
3416  void
3417  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3418  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PostProcessConstraints(
3419  ConstraintVector<VECTOR>& g) const
3420  {
3421  return this->GetConstraints()->PostProcessConstraints(g);
3422  }
3423  /******************************************************/
3424 
3425  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3426  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3427  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3428  template<int, int> class DH>
3429  void
3430  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3431  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryControl(
3432  const ControlVector<VECTOR>* c, std::string name)
3433  {
3434  if (auxiliary_controls_.find(name) != auxiliary_controls_.end())
3435  {
3436  throw DOpEException(
3437  "Adding multiple Data with name " + name + " is prohibited!",
3438  "OptProblemContainer::AddAuxiliaryControl");
3439  }
3440  auxiliary_controls_.insert(
3441  std::pair<std::string, const ControlVector<VECTOR>*>(name, c));
3442  }
3443 
3444  /******************************************************/
3445 
3446  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3447  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3448  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3449  template<int, int> class DH>
3450  void
3451  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3452  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryState(
3453  const StateVector<VECTOR>* c, std::string name)
3454  {
3455  if (auxiliary_state_.find(name) != auxiliary_state_.end())
3456  {
3457  throw DOpEException(
3458  "Adding multiple Data with name " + name + " is prohibited!",
3459  "OptProblemContainer::AddAuxiliaryState");
3460  }
3461  auxiliary_state_.insert(
3462  std::pair<std::string, const StateVector<VECTOR>*>(name, c));
3463  }
3464  /******************************************************/
3465 
3466  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3467  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3468  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3469  template<int, int> class DH>
3470  void
3471  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3472  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryConstraint(
3473  const ConstraintVector<VECTOR>* c, std::string name)
3474  {
3475  if (auxiliary_constraints_.find(name) != auxiliary_constraints_.end())
3476  {
3477  throw DOpEException(
3478  "Adding multiple Data with name " + name + " is prohibited!",
3479  "OptProblemContainer::AddAuxiliaryConstraint");
3480  }
3481  auxiliary_constraints_.insert(
3482  std::pair<std::string, const ConstraintVector<VECTOR>*>(name, c));
3483  }
3484  /******************************************************/
3485 
3486  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3487  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3488  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3489  template<int, int> class DH>
3490  const ControlVector<VECTOR>*
3491  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3492  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetAuxiliaryControl(
3493  std::string name) const
3494  {
3495  typename std::map<std::string, const ControlVector<VECTOR> *>::const_iterator it =
3496  auxiliary_controls_.find(name);
3497  if (it == auxiliary_controls_.end())
3498  {
3499  throw DOpEException("Could not find Data with name " + name,
3500  "OptProblemContainer::GetAuxiliaryControl");
3501  }
3502  return it->second;
3503  }
3504 
3505  /******************************************************/
3506 
3507  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3508  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3509  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3510  template<int, int> class DH>
3511  const StateVector<VECTOR>*
3512  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3513  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetAuxiliaryState(
3514  std::string name) const
3515  {
3516  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
3517  auxiliary_state_.find(name);
3518  if (it == auxiliary_state_.end())
3519  {
3520  throw DOpEException("Could not find Data with name " + name,
3521  "OptProblemContainer::GetAuxiliaryState");
3522  }
3523  return it->second;
3524  }
3525 
3526  /******************************************************/
3527 
3528  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3529  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3530  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3531  template<int, int> class DH>
3532  void
3533  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3534  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryControl(
3535  std::string name)
3536  {
3537  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
3538  auxiliary_controls_.find(name);
3539  if (it == auxiliary_controls_.end())
3540  {
3541  throw DOpEException(
3542  "Deleting Data " + name + " is impossible! Data not found",
3543  "OptProblemContainer::DeleteAuxiliaryControl");
3544  }
3545  auxiliary_controls_.erase(it);
3546  }
3547 
3548  /******************************************************/
3549 
3550  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3551  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3552  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3553  template<int, int> class DH>
3554  void
3555  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3556  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryState(
3557  std::string name)
3558  {
3559  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
3560  auxiliary_state_.find(name);
3561  if (it == auxiliary_state_.end())
3562  {
3563  throw DOpEException(
3564  "Deleting Data " + name + " is impossible! Data not found",
3565  "OptProblemContainer::DeleteAuxiliaryState");
3566  }
3567  auxiliary_state_.erase(it);
3568  }
3569 
3570  /******************************************************/
3571 
3572  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3573  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3574  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3575  template<int, int> class DH>
3576  void
3577  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3578  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryConstraint(
3579  std::string name)
3580  {
3581  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
3582  auxiliary_constraints_.find(name);
3583  if (it == auxiliary_constraints_.end())
3584  {
3585  throw DOpEException(
3586  "Deleting Data " + name + " is impossible! Data not found",
3587  "OptProblemContainer::DeleteAuxiliaryConstraint");
3588  }
3589  auxiliary_constraints_.erase(it);
3590  }
3591 
3592  /******************************************************/
3593 
3594  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3595  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3596  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3597  template<int, int> class DH>
3598  FUNCTIONAL*
3599  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3600  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctional()
3601  {
3602  if (this->GetType() == "aux_functional"
3603  || this->GetType() == "functional_for_ee")
3604  {
3605  //This may no longer happen!
3606  abort();
3607  // return aux_functionals_[this->GetTypeNum()];
3608  }
3609  return functional_;
3610 
3611  }
3612 
3613  /******************************************************/
3614 
3615  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3616  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3617  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3618  template<int, int> class DH>
3619  const FUNCTIONAL*
3620  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3621  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctional() const
3622  {
3623  if (this->GetType() == "aux_functional"
3624  || this->GetType() == "functional_for_ee")
3625  {
3626  //This may no longer happen!
3627  abort();
3628  // return aux_functionals_[this->GetTypeNum()];
3629  }
3630  return functional_;
3631  }
3632 
3633  /******************************************************/
3634 
3635  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3636  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3637  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3638  template<int, int> class DH>
3639  bool
3640  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3641  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasFaces() const
3642  {
3643  if (this->GetType().find("aux_functional") != std::string::npos)
3644  {
3645  return aux_functionals_[this->GetTypeNum()]->HasFaces();
3646  }
3647  else if (this->GetType().find("functional") != std::string::npos)
3648  {
3649  return this->GetFunctional()->HasFaces();
3650  }
3651  else if (this->GetType().find("constraint") != std::string::npos)
3652  {
3653  return this->GetConstraints()->HasFaces();
3654  }
3655  else
3656  {
3657  if ((this->GetType() == "state") || (this->GetType() == "tangent")
3658  || (this->GetType() == "gradient"))
3659  {
3660  return this->GetPDE().HasFaces();
3661  }
3662  else if ((this->GetType() == "adjoint")
3663  || (this->GetType() == "adjoint_hessian")
3664  || (this->GetType() == "hessian"))
3665  {
3666  return this->GetPDE().HasFaces() || this->GetFunctional()->HasFaces();
3667  }
3668  else if (this->GetType() == "adjoint_for_ee")
3669  {
3670  return this->GetPDE().HasFaces()
3671  || aux_functionals_[functional_for_ee_num_]->HasFaces();
3672  }
3673  else
3674  {
3675  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3676  "OptProblemContainer::HasFaces");
3677  }
3678  }
3679  }
3680 
3681  /******************************************************/
3682 
3683  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3684  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3685  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3686  template<int, int> class DH>
3687  bool
3688  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3689  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasPoints() const
3690  {
3691  if (this->GetType().find("constraint") != std::string::npos
3692  || (this->GetType() == "functional")
3693  || this->GetType() == "aux_functional" || (this->GetType() == "state")
3694  || (this->GetType() == "tangent"))
3695  {
3696  // We dont need PointRhs in this cases.
3697  return false;
3698  }
3699  else if ((this->GetType() == "adjoint")
3700  || (this->GetType() == "adjoint_hessian")
3701  || (this->GetType() == "hessian"))
3702  {
3703  return this->GetFunctional()->HasPoints();
3704  }
3705  else if (this->GetType() == "adjoint_for_ee")
3706  {
3707  return aux_functionals_[functional_for_ee_num_]->HasPoints();
3708  }
3709  else if (this->GetType() == "gradient")
3710  {
3711  return this->GetFunctional()->HasPoints();
3712  }
3713  else
3714  {
3715  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3716  "OptProblem::HasPoints");
3717  }
3718  }
3719 // }
3720 
3721  /******************************************************/
3722 
3723  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3724  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3725  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3726  template<int, int> class DH>
3727  bool
3728  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3729  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasInterfaces() const
3730  {
3731  if (this->GetType().find("aux_functional") != std::string::npos)
3732  {
3733  return false;
3734  }
3735  else if (this->GetType().find("functional") != std::string::npos)
3736  {
3737  return false;
3738  }
3739  else if (this->GetType().find("constraint") != std::string::npos)
3740  {
3741  return false;
3742  }
3743  else
3744  {
3745  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee"
3746  || (this->GetType() == "tangent") || (this->GetType() == "gradient")
3747  || (this->GetType() == "adjoint")
3748  || (this->GetType() == "adjoint_hessian")
3749  || (this->GetType() == "hessian"))
3750  {
3751  return this->GetPDE().HasInterfaces();
3752  }
3753  else
3754  {
3755  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3756  "OptProblemContainer::HasInterfaces");
3757  }
3758  }
3759  }
3760 
3761  /******************************************************/
3762 
3763  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3764  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3765  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3766  template<int, int> class DH>
3767  template<typename ELEMENTITERATOR>
3768  bool
3769  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3770  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AtInterface(ELEMENTITERATOR& element, unsigned int face) const
3771  {
3772  if (this->GetType().find("aux_functional") != std::string::npos)
3773  {
3774  return false;
3775  }
3776  else if (this->GetType().find("functional") != std::string::npos)
3777  {
3778  return false;
3779  }
3780  else if (this->GetType().find("constraint") != std::string::npos)
3781  {
3782  return false;
3783  }
3784  else
3785  {
3786  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee"
3787  || (this->GetType() == "tangent") || (this->GetType() == "gradient")
3788  || (this->GetType() == "adjoint")
3789  || (this->GetType() == "adjoint_hessian")
3790  || (this->GetType() == "hessian"))
3791  {
3792  return this->GetPDE().AtInterface(element,face);
3793  }
3794  else
3795  {
3796  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3797  "OptProblemContainer::HasFaces");
3798  }
3799  }
3800  }
3801  /******************************************************/
3802 
3803  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3804  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3805  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3806  template<int, int> class DH>
3807  void
3808  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3809  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetControlDirichletBoundaryColors(
3810  unsigned int color, const std::vector<bool>& comp_mask,
3811  const DOpEWrapper::Function<dealdim>* values)
3812  {
3813  assert(values);
3814 
3815  unsigned int comp = control_dirichlet_colors_.size();
3816  for (unsigned int i = 0; i < control_dirichlet_colors_.size(); ++i)
3817  {
3818  if (control_dirichlet_colors_[i] == color)
3819  {
3820  comp = i;
3821  break;
3822  }
3823  }
3824  if (comp != control_dirichlet_colors_.size())
3825  {
3826  std::stringstream s;
3827  s << "ControlDirichletColor" << color << " has multiple occurences !";
3828  throw DOpEException(s.str(),
3829  "OptProblemContainer::SetControlDirichletBoundary");
3830  }
3831  control_dirichlet_colors_.push_back(color);
3832  control_dirichlet_comps_.push_back(comp_mask);
3833  control_dirichlet_values_.push_back(values);
3834  }
3835 
3836  /******************************************************/
3837 
3838  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3839  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3840  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3841  template<int, int> class DH>
3842  void
3843  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3844  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetDirichletBoundaryColors(
3845  unsigned int color, const std::vector<bool>& comp_mask,
3846  const DD* values)
3847  {
3848  assert(values);
3849  assert(values->n_components() == comp_mask.size());
3850 
3851  unsigned int comp = dirichlet_colors_.size();
3852  for (unsigned int i = 0; i < dirichlet_colors_.size(); ++i)
3853  {
3854  if (dirichlet_colors_[i] == color)
3855  {
3856  comp = i;
3857  break;
3858  }
3859  }
3860  if (comp != dirichlet_colors_.size())
3861  {
3862  std::stringstream s;
3863  s << "DirichletColor" << color << " has multiple occurrences !";
3864  throw DOpEException(s.str(),
3865  "OptProblemContainer::SetDirichletBoundary");
3866  }
3867  dirichlet_colors_.push_back(color);
3868  dirichlet_comps_.push_back(comp_mask);
3869  PrimalDirichletData<DD, VECTOR, dealdim>* data =
3870  new PrimalDirichletData<DD, VECTOR, dealdim>(*values);
3871  primal_dirichlet_values_.push_back(data);
3872  TangentDirichletData<DD, VECTOR, dealdim>* tdata =
3873  new TangentDirichletData<DD, VECTOR, dealdim>(*values);
3874  tangent_dirichlet_values_.push_back(tdata);
3875 
3876  if (values->NeedsControl())
3877  {
3878  control_transposed_dirichlet_colors_.push_back(color);
3879  TransposedGradientDirichletData<DD, VECTOR, dealdim> * gdata =
3880  new TransposedGradientDirichletData<DD, VECTOR, dealdim>(
3881  *values);
3882  transposed_control_gradient_dirichlet_values_.push_back(gdata);
3883  TransposedHessianDirichletData<DD, VECTOR, dealdim> * hdata =
3884  new TransposedHessianDirichletData<DD, VECTOR, dealdim>(
3885  *values);
3886  transposed_control_hessian_dirichlet_values_.push_back(hdata);
3887  }
3888  }
3889 
3890  /******************************************************/
3891 
3892  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3893  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3894  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3895  template<int, int> class DH>
3896  const std::vector<unsigned int>&
3897  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3898  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletColors() const
3899  {
3900  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3901  || this->GetType() == "adjoint_for_ee"
3902  || (this->GetType() == "tangent")
3903  || (this->GetType() == "adjoint_hessian"))
3904  {
3905  return dirichlet_colors_;
3906  }
3907  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
3908  || (this->GetType() == "global_constraint_gradient"))
3909  {
3910  return control_dirichlet_colors_;
3911  }
3912  else
3913  {
3914  throw DOpEException("Unknown Type:" + this->GetType(),
3915  "OptProblemContainer::GetDirichletColors");
3916  }
3917  }
3918  /******************************************************/
3919 
3920  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3921  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3922  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3923  template<int, int> class DH>
3924  const std::vector<unsigned int>&
3925  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3926  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletColors() const
3927  {
3928  if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
3929  {
3930  return control_transposed_dirichlet_colors_;
3931  }
3932  else
3933  {
3934  throw DOpEException("Unknown Type:" + this->GetType(),
3935  "OptProblemContainer::GetTransposedDirichletColors");
3936  }
3937  }
3938 
3939  /******************************************************/
3940 
3941  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3942  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3943  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3944  template<int, int> class DH>
3945  const std::vector<bool>&
3946  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3947  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletCompMask(
3948  unsigned int color) const
3949  {
3950  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3951  || this->GetType() == "adjoint_for_ee"
3952  || (this->GetType() == "tangent")
3953  || (this->GetType() == "adjoint_hessian"))
3954  {
3955  unsigned int comp = dirichlet_colors_.size();
3956  for (unsigned int i = 0; i < dirichlet_colors_.size(); ++i)
3957  {
3958  if (dirichlet_colors_[i] == color)
3959  {
3960  comp = i;
3961  break;
3962  }
3963  }
3964  if (comp == dirichlet_colors_.size())
3965  {
3966  std::stringstream s;
3967  s << "DirichletColor" << color << " has not been found !";
3968  throw DOpEException(s.str(),
3969  "OptProblemContainer::GetDirichletCompMask");
3970  }
3971  return dirichlet_comps_[comp];
3972  }
3973  else if ((this->GetType() == "gradient")
3974  || (this->GetType() == "hessian"))
3975  {
3976  unsigned int comp = control_dirichlet_colors_.size();
3977  for (unsigned int i = 0; i < control_dirichlet_colors_.size(); ++i)
3978  {
3979  if (control_dirichlet_colors_[i] == color)
3980  {
3981  comp = i;
3982  break;
3983  }
3984  }
3985  if (comp == control_dirichlet_colors_.size())
3986  {
3987  std::stringstream s;
3988  s << "ControlDirichletColor" << color << " has not been found !";
3989  throw DOpEException(s.str(),
3990  "OptProblemContainer::GetDirichletCompMask");
3991  }
3992  return control_dirichlet_comps_[comp];
3993  }
3994  else
3995  {
3996  throw DOpEException("Unknown Type:" + this->GetType(),
3997  "OptProblemContainer::GetDirichletCompMask");
3998  }
3999  }
4000  /******************************************************/
4001 
4002  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4003  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4004  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4005  template<int, int> class DH>
4006  const std::vector<bool>&
4007  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4008  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletCompMask(
4009  unsigned int color) const
4010  {
4011  if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
4012  {
4013  unsigned int comp = dirichlet_colors_.size();
4014  for (unsigned int i = 0; i < dirichlet_colors_.size(); ++i)
4015  {
4016  if (dirichlet_colors_[i] == color)
4017  {
4018  comp = i;
4019  break;
4020  }
4021  }
4022  if (comp == dirichlet_colors_.size())
4023  {
4024  std::stringstream s;
4025  s << "DirichletColor" << color << " has not been found !";
4026  throw DOpEException(s.str(),
4027  "OptProblemContainer::GetDirichletCompMask");
4028  }
4029  return dirichlet_comps_[comp];
4030  }
4031  else
4032  {
4033  throw DOpEException("Unknown Type:" + this->GetType(),
4034  "OptProblemContainer::GetTransposedDirichletCompMask");
4035  }
4036  }
4037 
4038  /******************************************************/
4039 
4040  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4041  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4042  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4043  template<int, int> class DH>
4044  const Function<dealdim>&
4045  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4046  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletValues(
4047  unsigned int color,
4048  const std::map<std::string, const dealii::Vector<double>*> &param_values,
4049  const std::map<std::string, const VECTOR*> &domain_values) const
4050  {
4051 
4052  unsigned int col = dirichlet_colors_.size();
4053  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
4054  || this->GetType() == "adjoint_for_ee"
4055  || (this->GetType() == "tangent")
4056  || (this->GetType() == "adjoint_hessian"))
4057  {
4058  for (unsigned int i = 0; i < dirichlet_colors_.size(); ++i)
4059  {
4060  if (dirichlet_colors_[i] == color)
4061  {
4062  col = i;
4063  break;
4064  }
4065  }
4066  if (col == dirichlet_colors_.size())
4067  {
4068  std::stringstream s;
4069  s << "DirichletColor" << color << " has not been found !";
4070  throw DOpEException(s.str(),
4071  "OptProblemContainer::GetDirichletValues");
4072  }
4073  }
4074  else if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
4075  {
4076  col = control_dirichlet_colors_.size();
4077  for (unsigned int i = 0; i < control_dirichlet_colors_.size(); ++i)
4078  {
4079  if (control_dirichlet_colors_[i] == color)
4080  {
4081  col = i;
4082  break;
4083  }
4084  }
4085  if (col == control_dirichlet_colors_.size())
4086  {
4087  std::stringstream s;
4088  s << "ControlDirichletColor" << color << " has not been found !";
4089  throw DOpEException(s.str(),
4090  "OptProblemContainer::GetDirichletValues");
4091  }
4092  }
4093  else
4094  {
4095  throw DOpEException("Unknown Type:" + this->GetType(),
4096  "OptProblemContainer::GetDirichletValues");
4097  }
4098 
4099  if (this->GetType() == "state")
4100  {
4101  primal_dirichlet_values_[col]->ReInit(param_values, domain_values,
4102  color);
4103  return *(primal_dirichlet_values_[col]);
4104  }
4105  else if (this->GetType() == "tangent")
4106  {
4107  tangent_dirichlet_values_[col]->ReInit(param_values, domain_values,
4108  color);
4109  return *(tangent_dirichlet_values_[col]);
4110  }
4111  else if (this->GetType() == "adjoint"
4112  || this->GetType() == "adjoint_for_ee"
4113  || (this->GetType() == "adjoint_hessian"))
4114  {
4115  return *(zero_dirichlet_values_);
4116  }
4117  else if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
4118  {
4119  return *(control_dirichlet_values_[col]);
4120  }
4121  else
4122  {
4123  throw DOpEException("Unknown Type:" + this->GetType(),
4124  "OptProblemContainer::GetDirichletValues");
4125  }
4126  }
4127  /******************************************************/
4128 
4129  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4130  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4131  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4132  template<int, int> class DH>
4133  const TransposedDirichletDataInterface<dealdim>&
4134  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4135  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletValues(
4136  unsigned int color,
4137  const std::map<std::string, const dealii::Vector<double>*> &param_values,
4138  const std::map<std::string, const VECTOR*> &domain_values) const
4139  {
4140  unsigned int col = control_transposed_dirichlet_colors_.size();
4141  if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
4142  {
4143  for (unsigned int i = 0;
4144  i < control_transposed_dirichlet_colors_.size(); ++i)
4145  {
4146  if (control_transposed_dirichlet_colors_[i] == color)
4147  {
4148  col = i;
4149  break;
4150  }
4151  }
4152  if (col == control_transposed_dirichlet_colors_.size())
4153  {
4154  std::stringstream s;
4155  s << "TransposedControlDirichletColor" << color
4156  << " has not been found !";
4157  throw DOpEException(s.str(),
4158  "OptProblemContainer::GetTransposedDirichletValues");
4159  }
4160  }
4161  else
4162  {
4163  throw DOpEException("Unknown Type:" + this->GetType(),
4164  "OptProblemContainer::GetTransposedDirichletValues");
4165  }
4166 
4167  if (this->GetType() == "gradient")
4168  {
4169  transposed_control_gradient_dirichlet_values_[col]->ReInit(param_values,
4170  domain_values, color);
4171  return *(transposed_control_gradient_dirichlet_values_[col]);
4172  }
4173  else if (this->GetType() == "hessian")
4174  {
4175  transposed_control_hessian_dirichlet_values_[col]->ReInit(param_values,
4176  domain_values, color);
4177  return *(transposed_control_hessian_dirichlet_values_[col]);
4178  }
4179  else
4180  {
4181  throw DOpEException("Unknown Type:" + this->GetType(),
4182  "OptProblemContainer::GetTransposedDirichletValues");
4183  }
4184  }
4185 
4186  /******************************************************/
4187 
4188  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4189  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4190  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4191  template<int, int> class DH>
4192  const std::vector<unsigned int>&
4193  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4194  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetBoundaryEquationColors() const
4195  {
4196  if (this->GetType() == "state" || this->GetType() == "tangent")
4197  {
4198  return state_boundary_equation_colors_;
4199  }
4200  else if (this->GetType() == "adjoint"
4201  || this->GetType() == "adjoint_for_ee"
4202  || this->GetType() == "adjoint_hessian")
4203  {
4204  return adjoint_boundary_equation_colors_;
4205  }
4206  else if (this->GetType() == "gradient" || (this->GetType() == "hessian")
4207  || (this->GetType() == "global_constraint_gradient"))
4208  {
4209  return control_boundary_equation_colors_;
4210  }
4211  else
4212  {
4213  throw DOpEException("Unknown Type:" + this->GetType(),
4214  "OptProblemContainer::GetBoundaryEquationColors");
4215  }
4216  }
4217 
4218  /******************************************************/
4219 
4220  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4221  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4222  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4223  template<int, int> class DH>
4224  void
4225  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4226  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetControlBoundaryEquationColors(
4227  unsigned int color)
4228  {
4229  { //Control Boundary Equation colors are simply inserted
4230  unsigned int comp = control_boundary_equation_colors_.size();
4231  for (unsigned int i = 0; i < control_boundary_equation_colors_.size();
4232  ++i)
4233  {
4234  if (control_boundary_equation_colors_[i] == color)
4235  {
4236  comp = i;
4237  break;
4238  }
4239  }
4240  if (comp != control_boundary_equation_colors_.size())
4241  {
4242  std::stringstream s;
4243  s << "Boundary Equation Color" << color
4244  << " has multiple occurences !";
4245  throw DOpEException(s.str(),
4246  "OptProblemContainer::SetControlBoundaryEquationColors");
4247  }
4248  control_boundary_equation_colors_.push_back(color);
4249  }
4250  }
4251  /******************************************************/
4252 
4253  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4254  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4255  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4256  template<int, int> class DH>
4257  void
4258  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4259  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetBoundaryEquationColors(
4260  unsigned int color)
4261  {
4262  { //State Boundary Equation colors are simply inserted
4263  unsigned int comp = state_boundary_equation_colors_.size();
4264  for (unsigned int i = 0; i < state_boundary_equation_colors_.size();
4265  ++i)
4266  {
4267  if (state_boundary_equation_colors_[i] == color)
4268  {
4269  comp = i;
4270  break;
4271  }
4272  }
4273  if (comp != state_boundary_equation_colors_.size())
4274  {
4275  std::stringstream s;
4276  s << "Boundary Equation Color" << color
4277  << " has multiple occurences !";
4278  throw DOpEException(s.str(),
4279  "OptProblemContainer::SetBoundaryEquationColors");
4280  }
4281  state_boundary_equation_colors_.push_back(color);
4282  }
4283  { //For the adjoint they are added with the boundary functional colors
4284  unsigned int comp = adjoint_boundary_equation_colors_.size();
4285  for (unsigned int i = 0; i < adjoint_boundary_equation_colors_.size();
4286  ++i)
4287  {
4288  if (adjoint_boundary_equation_colors_[i] == color)
4289  {
4290  comp = i;
4291  break;
4292  }
4293  }
4294  if (comp != adjoint_boundary_equation_colors_.size())
4295  {
4296  //Seems this color is already added, however it might have been a functional color
4297  //so we don't do anything.
4298  }
4299  else
4300  {
4301  adjoint_boundary_equation_colors_.push_back(color);
4302  }
4303  }
4304  }
4305 
4306  /******************************************************/
4307 
4308  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4309  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4310  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4311  template<int, int> class DH>
4312  const std::vector<unsigned int>&
4313  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4314  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetBoundaryFunctionalColors() const
4315  {
4316  if (this->GetType() == "cost_functional"
4317  || this->GetType() == "aux_functional"
4318  || this->GetType() == "functional_for_ee") //fixme: what about error_evaluation?
4319  {
4320  return boundary_functional_colors_;
4321  }
4322  else
4323  {
4324  throw DOpEException("Unknown Type:" + this->GetType(),
4325  "OptProblemContainer::GetBoundaryFunctionalColors");
4326  }
4327  }
4328 
4329  /******************************************************/
4330 
4331  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4332  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4333  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4334  template<int, int> class DH>
4335  void
4336  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4337  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetBoundaryFunctionalColors(
4338  unsigned int color)
4339  {
4340  { //Boundary Functional colors are simply inserted
4341  unsigned int comp = boundary_functional_colors_.size();
4342  for (unsigned int i = 0; i < boundary_functional_colors_.size(); ++i)
4343  {
4344  if (boundary_functional_colors_[i] == color)
4345  {
4346  comp = i;
4347  break;
4348  }
4349  }
4350  if (comp != boundary_functional_colors_.size())
4351  {
4352  std::stringstream s;
4353  s << "Boundary Functional Color" << color
4354  << " has multiple occurences !";
4355  throw DOpEException(s.str(),
4356  "OptProblemContainer::SetBoundaryFunctionalColors");
4357  }
4358  boundary_functional_colors_.push_back(color);
4359  }
4360  { //For the adjoint they are addeed to the boundary equation colors
4361  unsigned int comp = adjoint_boundary_equation_colors_.size();
4362  for (unsigned int i = 0; i < adjoint_boundary_equation_colors_.size();
4363  ++i)
4364  {
4365  if (adjoint_boundary_equation_colors_[i] == color)
4366  {
4367  comp = i;
4368  break;
4369  }
4370  }
4371  if (comp != adjoint_boundary_equation_colors_.size())
4372  {
4373  //Seems this color is already added, however it might have been a equation color
4374  //so we don't do anything.
4375  }
4376  else
4377  {
4378  adjoint_boundary_equation_colors_.push_back(color);
4379  }
4380  }
4381  }
4382 
4383  /******************************************************/
4384 
4385  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4386  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4387  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4388  template<int, int> class DH>
4389  unsigned int
4390  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4391  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetControlNBlocks() const
4392  {
4393  return this->GetPDE().GetControlNBlocks();
4394  }
4395 
4396  /******************************************************/
4397 
4398  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4399  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4400  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4401  template<int, int> class DH>
4402  unsigned int
4403  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4404  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetStateNBlocks() const
4405  {
4406  return this->GetPDE().GetStateNBlocks();
4407  }
4408 
4409  /******************************************************/
4410 
4411  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4412  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4413  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4414  template<int, int> class DH>
4415  unsigned int
4416  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4417  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetNBlocks() const
4418  {
4419  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee")
4420  || (this->GetType() == "adjoint") || (this->GetType() == "tangent")
4421  || (this->GetType() == "adjoint_hessian"))
4422  {
4423  return this->GetStateNBlocks();
4424  }
4425  else if ((this->GetType() == "gradient")
4426  || (this->GetType() == "hessian"))
4427  {
4428  return this->GetControlNBlocks();
4429  }
4430  else
4431  {
4432  throw DOpEException("Unknown Type:" + this->GetType(),
4433  "OptProblemContainer::GetNBlocks");
4434  }
4435  }
4436 
4437  /******************************************************/
4438 
4439  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4440  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4441  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4442  template<int, int> class DH>
4443  unsigned int
4444  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4445  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFsPerBlock(
4446  unsigned int b) const
4447  {
4448  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
4449  || (this->GetType() == "adjoint_for_ee")
4450  || (this->GetType() == "tangent")
4451  || (this->GetType() == "adjoint_hessian"))
4452  {
4453  return GetSpaceTimeHandler()->GetStateDoFsPerBlock(b);
4454  }
4455  else if ((this->GetType() == "gradient")
4456  || (this->GetType() == "hessian"))
4457  {
4458  return GetSpaceTimeHandler()->GetControlDoFsPerBlock(b);
4459  }
4460  else
4461  {
4462  throw DOpEException("Unknown Type:" + this->GetType(),
4463  "OptProblemContainer::GetDoFsPerBlock");
4464  }
4465  }
4466 
4467  /******************************************************/
4468 
4469  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4470  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4471  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4472  template<int, int> class DH>
4473  const std::vector<unsigned int>&
4474  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4475  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFsPerBlock() const
4476  {
4477  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
4478  || (this->GetType() == "adjoint_for_ee")
4479  || (this->GetType() == "tangent")
4480  || (this->GetType() == "adjoint_hessian"))
4481  {
4482  return GetSpaceTimeHandler()->GetStateDoFsPerBlock();
4483  }
4484  else if ((this->GetType() == "gradient")
4485  || (this->GetType() == "hessian"))
4486  {
4487  return GetSpaceTimeHandler()->GetControlDoFsPerBlock();
4488  }
4489  else
4490  {
4491  throw DOpEException("Unknown Type:" + this->GetType(),
4492  "OptProblemContainer::GetDoFsPerBlock");
4493  }
4494  }
4495 
4496  /******************************************************/
4497 
4498  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4499  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4500  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4501  template<int, int> class DH>
4502  const dealii::ConstraintMatrix&
4503  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4504  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFConstraints() const
4505  {
4506  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
4507  || (this->GetType() == "adjoint_for_ee")
4508  || (this->GetType() == "tangent")
4509  || (this->GetType() == "adjoint_hessian"))
4510  {
4511  return GetSpaceTimeHandler()->GetStateDoFConstraints();
4512  }
4513  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
4514  || (this->GetType() == "global_constraint_gradient"))
4515  {
4516  return GetSpaceTimeHandler()->GetControlDoFConstraints();
4517  }
4518  else
4519  {
4520  throw DOpEException("Unknown Type:" + this->GetType(),
4521  "OptProblemContainer::GetDoFConstraints");
4522  }
4523  }
4524 
4525  /******************************************************/
4526 
4527  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4528  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4529  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4530  template<int, int> class DH>
4531  bool
4532  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4533  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::NeedTimeFunctional() const
4534  {
4535  if (this->GetType() == "cost_functional")
4536  return GetFunctional()->NeedTime();
4537  else if (this->GetType() == "aux_functional")
4538  return aux_functionals_[this->GetTypeNum()]->NeedTime();
4539  else if (this->GetType() == "functional_for_ee")
4540  return aux_functionals_[functional_for_ee_num_]->NeedTime();
4541  else
4542  throw DOpEException("Not implemented",
4543  "OptProblemContainer::NeedTimeFunctional");
4544  }
4545 
4546  /******************************************************/
4547 
4548  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4549  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4550  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4551  template<int, int> class DH>
4552  bool
4553  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4554  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasControlInDirichletData() const
4555  {
4556  return (!control_transposed_dirichlet_colors_.empty());
4557  }
4558 
4559 }
4560 #endif
DOpE::OptProblemContainer::ElementTimeMatrix
void ElementTimeMatrix(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix)
DOpE::OptProblemContainer::ElementMatrix
void ElementMatrix(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)
Definition: optproblemcontainer.h:2849
DOpE::OptProblemContainer::DeleteAuxiliaryState
void DeleteAuxiliaryState(std::string name)
Definition: optproblemcontainer.h:3556
DOpE::OptProblemContainer::GetOutputHandler
DOpEOutputHandler< VECTOR > * GetOutputHandler()
Definition: optproblemcontainer.h:822
DOpE::OptProblemContainer::Init_PointRhs
void Init_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: optproblemcontainer.h:1254
DOpE::OptProblemContainer::OptProblemContainer
OptProblemContainer(FUNCTIONAL &functional, PDE &pde, CONSTRAINTS &constraints, SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > &STH)
Definition: optproblemcontainer.h:1440
DOpE::OptProblemContainer::BoundaryEquation
void BoundaryEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:2084
DOpE::OptProblemContainer::GetUpdateFlags
dealii::UpdateFlags GetUpdateFlags() const
Definition: optproblemcontainer.h:3166
DOpE::StateProblem
Definition: stateproblem.h:48
DOpE::OptProblemContainer::GetAuxiliaryConstraint
const ConstraintVector< VECTOR > * GetAuxiliaryConstraint(std::string name)
Definition: optproblemcontainer.h:914
DOpE::OptProblemContainer::GetNFunctionals
unsigned int GetNFunctionals() const
Definition: optproblemcontainer.h:754
DOpE::OptProblemContainer::FaceFunctional
double FaceFunctional(const FACEDATACONTAINER &fdc)
Definition: optproblemcontainer.h:1763
DOpE::OptProblemContainer::Init_ElementEquation
void Init_ElementEquation(const DATACONTAINER &edc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1189
DOpE::OptProblemContainer::RegisterOutputHandler
void RegisterOutputHandler(DOpEOutputHandler< VECTOR > *OH)
Definition: optproblemcontainer.h:164
DOpE::OptProblemContainer::BoundaryRhs
void BoundaryRhs(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:2675
DOpE::ConstraintVector
Definition: constraintvector.h:48
DOpE::OptProblemContainer::AddPreviousAuxiliaryToIntegrator
void AddPreviousAuxiliaryToIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:1030
DOpE::OptProblemContainer::AlgebraicFunctional
double AlgebraicFunctional(const std::map< std::string, const dealii::Vector< double > * > &values, const std::map< std::string, const VECTOR * > &block_values)
Definition: optproblemcontainer.h:1796
DOpE::OptProblemContainer::ComputeLocalControlConstraints
void ComputeLocalControlConstraints(VECTOR &constraints, const std::map< std::string, const dealii::Vector< double > * > &values, const std::map< std::string, const VECTOR * > &block_values)
Definition: optproblemcontainer.h:3074
DOpE::OptProblemContainer::ElementTimeMatrixExplicit
void ElementTimeMatrixExplicit(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix)
Definition: optproblemcontainer.h:2929
DOpE::OptProblemContainer::ElementTimeEquationExplicit
void ElementTimeEquationExplicit(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:1954
DOpE::OptProblemContainer::GetBoundaryFunctionalColors
const std::vector< unsigned int > & GetBoundaryFunctionalColors() const
Definition: optproblemcontainer.h:4314
DOpE::OptProblemContainer::ElementRhs
void ElementRhs(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:2135
DOpE::OptProblemContainer::AddNextAuxiliaryToIntegrator
void AddNextAuxiliaryToIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:1003
DOpE::OptProblemContainer::DeletePreviousAuxiliaryFromIntegrator
void DeletePreviousAuxiliaryFromIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:1084
DOpE::OptProblemContainer::ReInit
void ReInit(std::string algo_type)
Definition: optproblemcontainer.h:1508
DOpE::OptProblemContainer::GetBlockVector
const VECTOR * GetBlockVector(const std::map< std::string, const VECTOR * > &values, std::string name)
Definition: optproblemcontainer.h:1353
DOpE::DOpEtypes::initial
Definition: dopetypes.h:106
DOpE::OptProblemContainer::GetInitialValues
const dealii::Function< dealdim > & GetInitialValues() const
Definition: optproblemcontainer.h:662
DOpE::OptProblemContainer::GetFESystem
const FE< dealdim, dealdim > & GetFESystem() const
Definition: optproblemcontainer.h:3131
DOpE::OptProblemContainer::PostProcessConstraints
void PostProcessConstraints(ConstraintVector< VECTOR > &g) const
Definition: optproblemcontainer.h:3418
DOpE::OptProblemContainer::GetFaceUpdateFlags
dealii::UpdateFlags GetFaceUpdateFlags() const
Definition: optproblemcontainer.h:3211
DOpE::OptProblemContainer::GetFunctionalType
std::string GetFunctionalType() const
Definition: optproblemcontainer.h:3255
DOpE::OptProblemContainer::AddFunctional
void AddFunctional(FUNCTIONAL_INTERFACE *F)
Definition: optproblemcontainer.h:692
DOpE::OptProblemContainer::GetFunctional
FUNCTIONAL * GetFunctional()
Definition: optproblemcontainer.h:3600
DOpE::TangentDirichletData
Definition: tangentdirichletdata.h:39
DOpE::OptProblemContainer::HasInterfaces
bool HasInterfaces() const
Definition: optproblemcontainer.h:3729
DOpE::OptProblemContainer::BoundaryFunctional
double BoundaryFunctional(const FACEDATACONTAINER &fdc)
Definition: optproblemcontainer.h:1725
DOpE::OptProblemContainer::ElementFunctional
double ElementFunctional(const DATACONTAINER &edc)
Definition: optproblemcontainer.h:1635
DOpE::ProblemContainerInternal::GetType
std::string GetType() const
Definition: problemcontainer_internal.h:109
DOpE::OptProblemContainer::GetBoundaryEquationColors
const std::vector< unsigned int > & GetBoundaryEquationColors() const
Definition: optproblemcontainer.h:4194
DOpE::OptProblemContainer::SetDirichletBoundaryColors
void SetDirichletBoundaryColors(unsigned int color, const std::vector< bool > &comp_mask, const DD *values)
Definition: optproblemcontainer.h:3844
DOpE::OptProblemContainer::GetDoFsPerBlock
const std::vector< unsigned int > & GetDoFsPerBlock() const
Definition: optproblemcontainer.h:4475
DOpE::OptProblemContainer::SetControlBoundaryEquationColors
void SetControlBoundaryEquationColors(unsigned int color)
Definition: optproblemcontainer.h:4226
DOpE::OptProblemContainer::PointRhs
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: optproblemcontainer.h:2323
DOpE::ProblemContainerInternal
Definition: problemcontainer_internal.h:38
DOpE::DOpEOutputHandler
Definition: optproblemcontainer.h:70
DOpE::OptProblemContainer::GetControlBoxConstraints
void GetControlBoxConstraints(VECTOR &lb, VECTOR &ub) const
Definition: optproblemcontainer.h:566
DOpE::OptProblemContainer::AddAuxiliaryControl
void AddAuxiliaryControl(const ControlVector< VECTOR > *c, std::string name)
Definition: optproblemcontainer.h:3431
DOpE::OptProblemContainer::GetDirichletColors
const std::vector< unsigned int > & GetDirichletColors() const
Definition: optproblemcontainer.h:3898
DOpE::OptProblemContainer::HasPoints
bool HasPoints() const
Definition: optproblemcontainer.h:3689
DOpE::OptProblemContainer::GetDoFType
std::string GetDoFType() const
Definition: optproblemcontainer.h:3103
DOpE::OptProblemContainer::GetName
virtual std::string GetName() const
Definition: optproblemcontainer.h:115
DOpE::OptProblemContainer::ElementEquation
void ElementEquation(const DATACONTAINER &edc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1833
DOpE::OptProblemContainer::ElementTimeEquation
void ElementTimeEquation(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:1910
DOpE::OptProblemContainer::GetDirichletCompMask
const std::vector< bool > & GetDirichletCompMask(unsigned int color) const
Definition: optproblemcontainer.h:3947
DOpE::OptProblemContainer::SetType
void SetType(std::string type, unsigned int num=0)
Definition: optproblemcontainer.h:1550
DOpE::OptProblemContainer::PointFunctional
double PointFunctional(const std::map< std::string, const dealii::Vector< double > * > &param_values, const std::map< std::string, const VECTOR * > &domain_values)
Definition: optproblemcontainer.h:1671
DOpE::TimeIterator
Definition: timeiterator.h:62
DOpE::DOpEtypes::stationary
Definition: dopetypes.h:105
DOpE::OptProblemContainer::GetNBlocks
unsigned int GetNBlocks() const
Definition: optproblemcontainer.h:4417
DOpE::OptProblemContainer::GetAuxiliaryState
const StateVector< VECTOR > * GetAuxiliaryState(std::string name) const
Definition: optproblemcontainer.h:3513
DOpE::OptProblemContainer::GetDirichletValues
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: optproblemcontainer.h:4046
DOpE::OptProblemContainer::SetTime
void SetTime(double time, unsigned int time_dof_number, const TimeIterator &interval, bool initial=false)
Definition: optproblemcontainer.h:3310
DOpE::ControlVector
Definition: controlvector.h:49
DOpE::OptProblemContainer::Init_ElementMatrix
void Init_ElementMatrix(const DATACONTAINER &edc, dealii::FullMatrix< double > &local_entry_matrix, double scale, double scale_ico)
Definition: optproblemcontainer.h:1304
DOpE::OptProblemContainer::GetStateProblem
StateProblem< OptProblemContainer< FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH >, PDE, DD, SPARSITYPATTERN, VECTOR, dealdim > & GetStateProblem()
Definition: optproblemcontainer.h:128
DOpE::ProblemContainerInternal::GetPDE
const PDE & GetPDE() const
Definition: problemcontainer_internal.h:103
DOpE::OptProblemContainer::EEFunctionalIsCost
bool EEFunctionalIsCost() const
Definition: optproblemcontainer.h:1328
DOpE::OptProblemContainer::GetDoFConstraints
const dealii::ConstraintMatrix & GetDoFConstraints() const
Definition: optproblemcontainer.h:4504
DOpE::OptProblemContainer::GetConstraintType
std::string GetConstraintType() const
Definition: optproblemcontainer.h:3297
DOpE::OptProblemContainer
Definition: optproblemcontainer.h:101
problemcontainer_internal.h
DOpE::OptProblemContainer::GetStateNBlocks
unsigned int GetStateNBlocks()
Definition: optproblemcontainer.h:1164
DOpE::OptProblemContainer::FaceEquation
void FaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:2002
DOpE::OptProblemContainer::AddAuxiliaryConstraint
void AddAuxiliaryConstraint(const ConstraintVector< VECTOR > *c, std::string name)
Definition: optproblemcontainer.h:3472
DOpE::OptProblemContainer::GetAuxiliaryControl
const ControlVector< VECTOR > * GetAuxiliaryControl(std::string name) const
Definition: optproblemcontainer.h:3492
DOpE::OptProblemContainer::GetExceptionHandler
DOpEExceptionHandler< VECTOR > * GetExceptionHandler()
Definition: optproblemcontainer.h:813
DOpE::OptProblemContainer::FaceRhs
void FaceRhs(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:2506
DOpE::OptProblemContainer::GetSpaceTimeHandler
const SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > * GetSpaceTimeHandler() const
Definition: optproblemcontainer.h:847
DOpE::OptProblemContainer::RegisterExceptionHandler
void RegisterExceptionHandler(DOpEExceptionHandler< VECTOR > *OH)
Definition: optproblemcontainer.h:172
DOpE::TransposedHessianDirichletData
Definition: transposedhessiandirichletdata.h:41
DOpE::OptProblemContainer::GetConstraints
CONSTRAINTS * GetConstraints()
Definition: optproblemcontainer.h:1342
DOpE::OptProblemContainer::FaceMatrix
void FaceMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)
DOpE::OptProblemContainer::DeleteAuxiliaryFromIntegrator
void DeleteAuxiliaryFromIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:1108
DOpE::OptProblemContainer::GetVector
const dealii::Vector< double > * GetVector(const std::map< std::string, const Vector< double > * > &values, std::string name)
Definition: optproblemcontainer.h:1366
DOpE::OptProblemContainer::GetSpaceTimeHandler
SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > * GetSpaceTimeHandler()
Definition: optproblemcontainer.h:855
DOpE::OptProblemContainer::GetConstraints
const CONSTRAINTS * GetConstraints() const
Definition: optproblemcontainer.h:1347
DOpE::OptProblemContainer::HasFaces
bool HasFaces() const
Definition: optproblemcontainer.h:3641
DOpE::PrimalDirichletData
Definition: primaldirichletdata.h:41
DOpE::OptProblemContainer::GetTransposedDirichletColors
const std::vector< unsigned int > & GetTransposedDirichletColors() const
Definition: optproblemcontainer.h:3926
DOpE::StateVector
Definition: statevector.h:50
DOpE::OptProblemContainer::GetTransposedDirichletValues
const TransposedDirichletDataInterface< dealdim > & GetTransposedDirichletValues(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: optproblemcontainer.h:4135
DOpE::OptProblemContainer::DeleteNextAuxiliaryFromIntegrator
void DeleteNextAuxiliaryFromIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:1058
DOpE::OptProblemContainer::ComputeSparsityPattern
void ComputeSparsityPattern(SPARSITYPATTERN &sparsity) const
Definition: optproblemcontainer.h:3384
DOpE::OptProblemContainer::AtInterface
bool AtInterface(ELEMENTITERATOR &element, unsigned int face) const
Definition: optproblemcontainer.h:3770
DOpE::OptProblemContainer::AlgebraicResidual
void AlgebraicResidual(VECTOR &residual, const std::map< std::string, const dealii::Vector< double > * > &values, const std::map< std::string, const VECTOR * > &block_values)
Definition: optproblemcontainer.h:1883
DOpE::OptProblemContainer::InterfaceEquation
void InterfaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:2052
DOpE::SpaceTimeHandler
Definition: spacetimehandler.h:71
fevalues_wrapper.h
DOpE::OptProblemContainer::GetBaseProblem
OptProblemContainer< FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH > & GetBaseProblem()
Definition: optproblemcontainer.h:145
DOpE::OptProblemContainer::SetFunctionalForErrorEstimation
void SetFunctionalForErrorEstimation(std::string functional_name)
Definition: optproblemcontainer.h:718
DOpE::OptProblemContainer::GetStateNBlocks
unsigned int GetStateNBlocks() const
Definition: optproblemcontainer.h:4404
DOpE::OptProblemContainer::GetControlBlockComponent
std::vector< unsigned int > & GetControlBlockComponent()
Definition: optproblemcontainer.h:1172
DOpE::OptProblemContainer::AddAuxiliaryState
void AddAuxiliaryState(const StateVector< VECTOR > *c, std::string name)
Definition: optproblemcontainer.h:3452
DOpE::OptProblemContainer::NeedTimeFunctional
bool NeedTimeFunctional() const
Definition: optproblemcontainer.h:4533
dofhandler_wrapper.h
DOpE::OptProblemContainer::GetFunctionalName
std::string GetFunctionalName() const
Definition: optproblemcontainer.h:3276
function_wrapper.h
DOpE::TransposedDirichletDataInterface
Definition: transposeddirichletdatainterface.h:37
DOpE::OptProblemContainer::GetControlNBlocks
unsigned int GetControlNBlocks() const
Definition: optproblemcontainer.h:4391
DOpE::OptProblemContainer::GetTransposedDirichletCompMask
const std::vector< bool > & GetTransposedDirichletCompMask(unsigned int color) const
Definition: optproblemcontainer.h:4008
DOpE::OptProblemContainer::~OptProblemContainer
virtual ~OptProblemContainer()
Definition: optproblemcontainer.h:1465
DOpE::OptProblemContainer::HasControlInDirichletData
bool HasControlInDirichletData() const
Definition: optproblemcontainer.h:4554
DOpE::OptProblemContainer::GetStateBlockComponent
std::vector< unsigned int > & GetStateBlockComponent()
Definition: optproblemcontainer.h:1179
DOpE::OptProblemContainer::SetInitialValues
void SetInitialValues(const dealii::Function< dealdim > *values)
Definition: optproblemcontainer.h:657
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::TransposedGradientDirichletData
Definition: transposedgradientdirichletdata.h:41
DOpE::OptProblemContainer::AddAuxiliaryToIntegrator
void AddAuxiliaryToIntegrator(INTEGRATOR &integrator)
Definition: optproblemcontainer.h:936
DOpE::OptProblemContainer::SetBoundaryEquationColors
void SetBoundaryEquationColors(unsigned int color)
Definition: optproblemcontainer.h:4259
DOpE::OptProblemContainer::SetControlDirichletBoundaryColors
void SetControlDirichletBoundaryColors(unsigned int color, const std::vector< bool > &comp_mask, const DOpEWrapper::Function< dealdim > *values)
Definition: optproblemcontainer.h:3809
DOpE::OptProblemContainer::DeleteAuxiliaryConstraint
void DeleteAuxiliaryConstraint(std::string name)
Definition: optproblemcontainer.h:3578
DOpE::OptProblemContainer::GetFunctionalPosition
const std::map< std::string, unsigned int > & GetFunctionalPosition() const
Definition: optproblemcontainer.h:1158
DOpE::OptProblemContainer::SetBoundaryFunctionalColors
void SetBoundaryFunctionalColors(unsigned int color)
Definition: optproblemcontainer.h:4337
DOpE::OptProblemContainer::DeleteAuxiliaryControl
void DeleteAuxiliaryControl(std::string name)
Definition: optproblemcontainer.h:3534
DOpE::OptProblemContainer::BoundaryMatrix
void BoundaryMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_matrix, double scale=1., double scale_ico=1.)
DOpE::OptProblemContainer::Init_ElementRhs
void Init_ElementRhs(const DATACONTAINER &edc, dealii::Vector< double > &local_vector, double scale)
Definition: optproblemcontainer.h:1208
DOpE::OptProblemContainer::InterfaceMatrix
void InterfaceMatrix(const FACEDATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix, double scale=1., double scale_ico=1.)
DOpE::DOpEExceptionHandler
Definition: optproblemcontainer.h:72
Generated on Mon Jan 25 2016 11:32:18 for DOpE by   doxygen 1.8.6