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 <lac/vector.h>
50 #include <lac/full_matrix.h>
51 #include <grid/tria_iterator.h>
52 #include <dofs/dof_handler.h>
53 #include <dofs/dof_accessor.h>
54 #include <dofs/dof_tools.h>
55 #include <fe/fe_system.h>
56 #include <fe/fe_values.h>
57 #include <fe/mapping.h>
58 #include <base/quadrature_lib.h>
59 #include <lac/block_sparsity_pattern.h>
60 #include <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  ~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 
837  void
838  SetTime(double time, const TimeIterator& interval);
839 
840  /******************************************************/
841 
842  const SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>*
843  GetSpaceTimeHandler() const
844  {
845  return _STH;
846  }
847 
848  /******************************************************/
849 
850  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>*
851  GetSpaceTimeHandler()
852  {
853  return _STH;
854  }
855 
856  /******************************************************/
857 
858  void
859  ComputeSparsityPattern(SPARSITYPATTERN & sparsity) const;
860 
861  /******************************************************/
862 
863  void
864  PostProcessConstraints(ConstraintVector<VECTOR>& g) const;
865 
866  /******************************************************/
867 
868  void
869  AddAuxiliaryControl(const ControlVector<VECTOR>* c, std::string name);
870 
871  /******************************************************/
872 
873  void
874  AddAuxiliaryState(const StateVector<VECTOR>* c, std::string name);
875 
876  /******************************************************/
877 
878  void
879  AddAuxiliaryConstraint(const ConstraintVector<VECTOR>* c,
880  std::string name);
881 
882  /******************************************************/
883 
884  const ControlVector<VECTOR>*
885  GetAuxiliaryControl(std::string name) const;
886 
887  /******************************************************/
888 
889  const StateVector<VECTOR>*
890  GetAuxiliaryState(std::string name) const;
891 
892  /******************************************************/
893 
894  void
895  DeleteAuxiliaryControl(std::string name);
896 
897  /******************************************************/
898 
899  void
900  DeleteAuxiliaryState(std::string name);
901 
902  /******************************************************/
903 
904  void
905  DeleteAuxiliaryConstraint(std::string name);
906 
907  /*****************************************************************/
908 
909  const ConstraintVector<VECTOR>*
910  GetAuxiliaryConstraint(std::string name)
911  {
912  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
913  _auxiliary_constraints.find(name);
914  if (it == _auxiliary_constraints.end())
915  {
916  throw DOpEException("Could not find data" + name,
917  "OptProblemContainer::GetAuxiliaryConstraint");
918  }
919  return it->second;
920  }
921  /*****************************************************************/
922 
923  template<typename INTEGRATOR>
924  void
925  AddAuxiliaryToIntegrator(INTEGRATOR& integrator)
926  {
927  {
928  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
929  _auxiliary_controls.begin();
930  for (; it != _auxiliary_controls.end(); it++)
931  {
932  if (dopedim == dealdim)
933  {
934  integrator.AddDomainData(it->first,
935  &(it->second->GetSpacialVector()));
936  }
937  else if (dopedim == 0)
938  {
939  integrator.AddParamData(it->first,
940  &(it->second->GetSpacialVectorCopy()));
941  }
942  else
943  {
944  throw DOpEException("dopedim not implemented",
945  "OptProblemContainer::AddAuxiliaryToIntegrator");
946  }
947  }
948  }
949  {
950  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
951  _auxiliary_state.begin();
952  for (; it != _auxiliary_state.end(); it++)
953  {
954  integrator.AddDomainData(it->first,
955  &(it->second->GetSpacialVector()));
956  }
957  }
958  {
959  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
960  _auxiliary_constraints.begin();
961  for (; it != _auxiliary_constraints.end(); it++)
962  {
963  integrator.AddDomainData(it->first + "_local",
964  &(it->second->GetSpacialVector("local")));
965  integrator.AddParamData(it->first + "_global",
966  &(it->second->GetGlobalConstraints()));
967  }
968  }
969  }
970 
971  /******************************************************/
972 
973  template<typename INTEGRATOR>
974  void
975  DeleteAuxiliaryFromIntegrator(INTEGRATOR& integrator)
976  {
977  {
978  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
979  _auxiliary_controls.begin();
980  for (; it != _auxiliary_controls.end(); it++)
981  {
982  if (dopedim == dealdim)
983  {
984  integrator.DeleteDomainData(it->first);
985  }
986  else if (dopedim == 0)
987  {
988  integrator.DeleteParamData(it->first);
989  it->second->UnLockCopy();
990  }
991  else
992  {
993  throw DOpEException("dopedim not implemented",
994  "OptProblemContainer::AddAuxiliaryToIntegrator");
995  }
996  }
997  }
998  {
999  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
1000  _auxiliary_state.begin();
1001  for (; it != _auxiliary_state.end(); it++)
1002  {
1003  integrator.DeleteDomainData(it->first);
1004  }
1005  }
1006  {
1007  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
1008  _auxiliary_constraints.begin();
1009  for (; it != _auxiliary_constraints.end(); it++)
1010  {
1011  integrator.DeleteDomainData(it->first + "_local");
1012  integrator.DeleteParamData(it->first + "_global");
1013  }
1014  }
1015  }
1016 
1017  /******************************************************/
1018  const std::map<std::string, unsigned int>&
1019  GetFunctionalPosition() const
1020  {
1021  return _functional_position;
1022  }
1023 
1024  unsigned int
1025  GetStateNBlocks()
1026  {
1027  return this->GetPDE().GetStateNBlocks();
1028  }
1029 
1030  /******************************************************/
1031 
1032  std::vector<unsigned int>&
1033  GetControlBlockComponent()
1034  {
1035  return this->GetPDE().GetControlBlockComponent();
1036  }
1037  /******************************************************/
1038 
1039  std::vector<unsigned int>&
1040  GetStateBlockComponent()
1041  {
1042  return this->GetPDE().GetStateBlockComponent();
1043  }
1044 
1045  /******************************************************/
1046  /******************************************************/
1047  /****For the initial values ***************/
1048  template<typename DATACONTAINER>
1049  void
1050  Init_ElementEquation(const DATACONTAINER& edc,
1051  dealii::Vector<double> &local_vector, double scale,
1052  double scale_ico)
1053  {
1054  if (this->GetType() == "adjoint" || this->GetType() == "tangent"
1055  || this->GetType() == "adjoint_hessian")
1056  {
1057  this->GetPDE().Init_ElementEquation(edc, local_vector, scale,
1058  scale_ico);
1059  }
1060  else
1061  {
1062  throw DOpEException("Not implemented",
1063  "OptProblemContainer::Init_ElementEquation");
1064  }
1065  }
1066 
1067  template<typename DATACONTAINER>
1068  void
1069  Init_ElementRhs(const DATACONTAINER& edc,
1070  dealii::Vector<double> &local_vector, double scale)
1071  {
1072  //FIXME: We should take care of the cases of boundary and face functionals...
1073  if (this->GetType() == "adjoint")
1074  {
1075  if (GetFunctional()->NeedTime())
1076  {
1077  if (GetFunctional()->GetType().find("timelocal")
1078  != std::string::npos)
1079  {
1080  if (GetFunctional()->GetType().find("domain")
1081  != std::string::npos)
1082  {
1083  GetFunctional()->ElementValue_U(edc, local_vector, scale);
1084  }
1085  }
1086  }
1087  }
1088  else if (this->GetType() == "tangent")
1089  {
1090  this->GetPDE().Init_ElementRhs_QT(edc, local_vector, scale);
1091  }
1092  else if (this->GetType() == "adjoint_hessian")
1093  {
1094  if (GetFunctional()->NeedTime())
1095  {
1096  if (GetFunctional()->GetType().find("timelocal")
1097  != std::string::npos)
1098  {
1099  if (GetFunctional()->GetType().find("domain")
1100  != std::string::npos)
1101  {
1102  GetFunctional()->ElementValue_UU(edc, local_vector, scale);
1103  }
1104  }
1105  }
1106  }
1107  else
1108  {
1109  throw DOpEException("Not implemented",
1110  "OptProblemContainer::Init_ElementRhs");
1111  }
1112  }
1113 
1114  void
1115  Init_PointRhs(
1116  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1117  const std::map<std::string, const VECTOR*> &domain_values,
1118  VECTOR& rhs_vector, double scale = 1.)
1119  {
1120  if (this->GetType() == "adjoint")
1121  {
1122  if (GetFunctional()->NeedTime())
1123  {
1124  if (GetFunctional()->GetType().find("timelocal")
1125  != std::string::npos)
1126  {
1127  if (GetFunctional()->GetType().find("point")
1128  != std::string::npos)
1129  {
1130  GetFunctional()->PointValue_U(param_values, domain_values,
1131  rhs_vector, scale);
1132  }
1133  }
1134  }
1135  }
1136  else if (this->GetType() == "tangent")
1137  {
1138  //Nothing to do for tangent, since no point sources are allowed in the initial condition.
1139  }
1140  else if (this->GetType() == "adjoint_hessian")
1141  {
1142  if (GetFunctional()->NeedTime())
1143  {
1144  if (GetFunctional()->GetType().find("timelocal")
1145  != std::string::npos)
1146  {
1147  if (GetFunctional()->GetType().find("point")
1148  != std::string::npos)
1149  {
1150  GetFunctional()->PointValue_UU(param_values, domain_values,
1151  rhs_vector, scale);
1152  }
1153  }
1154  }
1155  }
1156  else
1157  {
1158  throw DOpEException("Not implemented",
1159  "OptProblemContainer::Init_PointRhs");
1160  }
1161  }
1162 
1163  template<typename DATACONTAINER>
1164  void
1165  Init_ElementMatrix(const DATACONTAINER& edc,
1166  dealii::FullMatrix<double> &local_entry_matrix, double scale,
1167  double scale_ico)
1168  {
1169  if (this->GetType() == "adjoint" || this->GetType() == "tangent"
1170  || this->GetType() == "adjoint_hessian")
1171  {
1172  this->GetPDE().Init_ElementMatrix(edc, local_entry_matrix, scale,
1173  scale_ico);
1174  }
1175  else
1176  {
1177  throw DOpEException("Not implemented",
1178  "OptProblemContainer::Init_ElementMatrix");
1179  }
1180  }
1181 
1182  /******************************************************/
1183 
1188  bool
1189  EEFunctionalIsCost() const
1190  {
1191  return _functional_for_ee_is_cost;
1192  }
1193 
1194  protected:
1195  FUNCTIONAL*
1196  GetFunctional();
1197  const FUNCTIONAL*
1198  GetFunctional() const;
1199  CONSTRAINTS *
1200  GetConstraints()
1201  {
1202  return _constraints;
1203  }
1204  const CONSTRAINTS *
1205  GetConstraints() const
1206  {
1207  return _constraints;
1208  }
1209 
1210  const VECTOR*
1211  GetBlockVector(const std::map<std::string, const VECTOR*>& values,
1212  std::string name)
1213  {
1214  typename std::map<std::string, const VECTOR*>::const_iterator it =
1215  values.find(name);
1216  if (it == values.end())
1217  {
1218  throw DOpEException("Did not find " + name,
1219  "OptProblemContainer::GetBlockVector");
1220  }
1221  return it->second;
1222  }
1223  const dealii::Vector<double>*
1224  GetVector(const std::map<std::string, const Vector<double>*>& values,
1225  std::string name)
1226  {
1227  typename std::map<std::string, const Vector<double>*>::const_iterator it =
1228  values.find(name);
1229  if (it == values.end())
1230  {
1231  throw DOpEException("Did not find " + name,
1232  "OptProblemContainer::GetVector");
1233  }
1234  return it->second;
1235  }
1236  /******************************************************/
1237  private:
1238  DOpEExceptionHandler<VECTOR>* _ExceptionHandler;
1239  DOpEOutputHandler<VECTOR>* _OutputHandler;
1240  std::string _algo_type;
1241 
1242  bool _functional_for_ee_is_cost;
1243  unsigned int _functional_for_ee_num;
1244  std::vector<FUNCTIONAL_INTERFACE*> _aux_functionals;
1245  std::map<std::string, unsigned int> _functional_position;
1246  FUNCTIONAL* _functional;
1247  CONSTRAINTS* _constraints;
1248  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>* _STH;
1249 
1250  std::vector<unsigned int> _control_dirichlet_colors;
1251  std::vector<unsigned int> _control_transposed_dirichlet_colors;
1252  std::vector<std::vector<bool> > _control_dirichlet_comps;
1253  std::vector<const DOpEWrapper::Function<dealdim>*> _control_dirichlet_values;
1254  std::vector<
1255  TransposedGradientDirichletData<DD, VECTOR, dealdim>*> _transposed_control_gradient_dirichlet_values;
1256  std::vector<
1257  TransposedHessianDirichletData<DD, VECTOR, dealdim>*> _transposed_control_hessian_dirichlet_values;
1258 
1259  std::vector<unsigned int> _dirichlet_colors;
1260  std::vector<std::vector<bool> > _dirichlet_comps;
1261  std::vector<PrimalDirichletData<DD, VECTOR, dealdim>*> _primal_dirichlet_values;
1262  std::vector<TangentDirichletData<DD, VECTOR, dealdim>*> _tangent_dirichlet_values;
1263  const dealii::Function<dealdim>* _zero_dirichlet_values;
1264 
1265  const dealii::Function<dealdim>* _initial_values;
1266 
1267  std::vector<unsigned int> _control_boundary_equation_colors;
1268  std::vector<unsigned int> _state_boundary_equation_colors;
1269  std::vector<unsigned int> _adjoint_boundary_equation_colors;
1270 
1271  std::vector<unsigned int> _boundary_functional_colors;
1272 
1273  std::map<std::string, const ControlVector<VECTOR>*> _auxiliary_controls;
1274  std::map<std::string, const StateVector<VECTOR>*> _auxiliary_state;
1275  std::map<std::string, const ConstraintVector<VECTOR>*> _auxiliary_constraints;
1276 
1277  StateProblem<
1278  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1279  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>,
1280  PDE, DD, SPARSITYPATTERN, VECTOR, dealdim> * _state_problem;
1281 
1282  friend class StateProblem<
1283  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1284  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>,
1285  PDE, DD, SPARSITYPATTERN, VECTOR, dealdim> ;
1286  };
1287  /******************************************************/
1288 
1289  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1290  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1291  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1292  template<int, int> class DH>
1293  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1294  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::OptProblemContainer(
1295  FUNCTIONAL& functional, PDE& pde, CONSTRAINTS& constraints,
1296  SpaceTimeHandler<FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim>& STH) :
1297  ProblemContainerInternal<PDE>(pde), _functional(&functional), _constraints(
1298  &constraints), _STH(&STH), _state_problem(NULL)
1299  {
1300  _ExceptionHandler = NULL;
1301  _OutputHandler = NULL;
1302  _zero_dirichlet_values = new ZeroFunction<dealdim>(
1303  this->GetPDE().GetStateNComponents());
1304  _algo_type = "";
1305  _functional_position[_functional->GetName()] = 0;
1306  //remember! At _functional_values[0] we store always the cost functional!
1307  _functional_for_ee_num = dealii::numbers::invalid_unsigned_int;
1308  }
1309 
1310  /******************************************************/
1311 
1312  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1313  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1314  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1315  template<int, int> class DH>
1316  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1317  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::~OptProblemContainer()
1318  {
1319  if (_zero_dirichlet_values != NULL)
1320  {
1321  delete _zero_dirichlet_values;
1322  }
1323 
1324  for (unsigned int i = 0;
1325  i < _transposed_control_gradient_dirichlet_values.size(); i++)
1326  {
1327  if (_transposed_control_gradient_dirichlet_values[i] != NULL)
1328  delete _transposed_control_gradient_dirichlet_values[i];
1329  }
1330  for (unsigned int i = 0;
1331  i < _transposed_control_hessian_dirichlet_values.size(); i++)
1332  {
1333  if (_transposed_control_hessian_dirichlet_values[i] != NULL)
1334  delete _transposed_control_hessian_dirichlet_values[i];
1335  }
1336  for (unsigned int i = 0; i < _primal_dirichlet_values.size(); i++)
1337  {
1338  if (_primal_dirichlet_values[i] != NULL)
1339  delete _primal_dirichlet_values[i];
1340  }
1341  for (unsigned int i = 0; i < _tangent_dirichlet_values.size(); i++)
1342  {
1343  if (_tangent_dirichlet_values[i] != NULL)
1344  delete _tangent_dirichlet_values[i];
1345  }
1346  if (_state_problem != NULL)
1347  {
1348  delete _state_problem;
1349  }
1350  }
1351 
1352  /******************************************************/
1353 
1354  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1355  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1356  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1357  template<int, int> class DH>
1358  void
1359  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1360  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ReInit(
1361  std::string algo_type)
1362  {
1363  if (_state_problem != NULL)
1364  {
1365  delete _state_problem;
1366  _state_problem = NULL;
1367  }
1368 
1369  if (_algo_type != algo_type && _algo_type != "")
1370  {
1371  throw DOpEException("Conflicting Algorithms!",
1372  "OptProblemContainer::ReInit");
1373  }
1374  else
1375  {
1376  _algo_type = algo_type;
1377  this->SetTypeInternal("");
1378 
1379  if (_algo_type == "reduced")
1380  {
1381  GetSpaceTimeHandler()->ReInit(this->GetPDE().GetControlNBlocks(),
1382  this->GetPDE().GetControlBlockComponent(),
1383  this->GetPDE().GetStateNBlocks(),
1384  this->GetPDE().GetStateBlockComponent());
1385  }
1386  else
1387  {
1388  throw DOpEException("Unknown Algorithm " + _algo_type,
1389  "OptProblemContainer::ReInit");
1390  }
1391  }
1392  }
1393 
1394  /******************************************************/
1395 
1396  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1397  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1398  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1399  template<int, int> class DH>
1400  void
1401  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1402  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetType(
1403  std::string type, unsigned int num)
1404  {
1405  if (this->GetType() != type || this->GetTypeNum() != num)
1406  {
1407  this->SetTypeNumInternal(num);
1408  this->SetTypeInternal(type);
1409  this->GetPDE().SetProblemType(type);
1410  if (_functional_for_ee_num != dealii::numbers::invalid_unsigned_int)
1411  _aux_functionals[_functional_for_ee_num]->SetProblemType(type);
1412  this->GetConstraints()->SetProblemType(type, num);
1413 
1414 #if dope_dimension > 0
1415  if(dealdim == dopedim)
1416  {
1417  //Prepare DoFHandlerPointer
1418 
1419  {
1420  if(this->GetType() == "state" ||this->GetType() == "adjoint"
1421  || this->GetType() == "adjoint_for_ee" || this->GetType() == "cost_functional"
1422  || this->GetType() == "aux_functional" || this->GetType() == "functional_for_ee"
1423  || this->GetType() == "tangent" || this->GetType() == "adjoint_hessian"
1424  || this->GetType() == "error_evaluation"
1425  || this->GetType().find("constraints") != std::string::npos)
1426  {
1427  GetSpaceTimeHandler()->SetDoFHandlerOrdering(1,0);
1428  }
1429  else if (this->GetType() == "gradient"||this->GetType() == "hessian"||this->GetType() == "hessian_inverse" || this->GetType() == "global_constraint_gradient"|| this->GetType() == "global_constraint_hessian")
1430  {
1431  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0,1);
1432  }
1433  else
1434  {
1435  throw DOpEException("_problem_type : "+this->GetType()+" not implemented!", "OptProblemContainer::SetType");
1436  }
1437  }
1438  }
1439  else
1440  {
1441  throw DOpEException("dopedim not implemented", "OptProblemContainer::SetType");
1442  }
1443 #else
1444  //dopedim ==0
1445  {
1446  //Prepare DoFHandlerPointer
1447  {
1448 
1449  if (this->GetType() == "state" || this->GetType() == "adjoint"
1450  || this->GetType() == "adjoint_for_ee"
1451  || this->GetType() == "functional_for_ee"
1452  || this->GetType() == "cost_functional"
1453  || this->GetType() == "aux_functional"
1454  || this->GetType() == "tangent"
1455  || this->GetType() == "error_evaluation"
1456  || this->GetType() == "adjoint_hessian")
1457  {
1458  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0, 0);
1459  }
1460  else if (this->GetType() == "gradient"
1461  || this->GetType() == "hessian_inverse"
1462  || this->GetType() == "hessian")
1463  {
1464  GetSpaceTimeHandler()->SetDoFHandlerOrdering(0, 0);
1465  }
1466  else
1467  {
1468  throw DOpEException(
1469  "_problem_type : " + this->GetType() + " not implemented!",
1470  "OptProblemContainer::SetType");
1471  }
1472  }
1473  }
1474 #endif
1475  }
1476  }
1477 
1478  /******************************************************/
1479 
1480  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1481  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1482  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1483  template<int, int> class DH>
1484  template<typename DATACONTAINER>
1485  double
1486  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1487  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementFunctional(
1488  const DATACONTAINER& edc)
1489  {
1490 
1491  if (this->GetType() == "cost_functional")
1492  {
1493  // state values in quadrature points
1494  return GetFunctional()->ElementValue(edc);
1495  }
1496  else if (this->GetType() == "aux_functional")
1497  {
1498  // state values in quadrature points
1499  return _aux_functionals[this->GetTypeNum()]->ElementValue(edc);
1500  }
1501  else if (this->GetType() == "functional_for_ee")
1502  { //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1503  return _aux_functionals[_functional_for_ee_num]->ElementValue(edc);
1504  }
1505  else if (this->GetType().find("constraints") != std::string::npos)
1506  {
1507  return GetConstraints()->ElementValue(edc);
1508  }
1509  else
1510  {
1511  throw DOpEException("Not implemented",
1512  "OptProblemContainer::ElementFunctional");
1513  }
1514  }
1515 
1516  /******************************************************/
1517  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1518  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1519  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1520  template<int, int> class DH>
1521  double
1522  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1523  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PointFunctional(
1524  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1525  const std::map<std::string, const VECTOR*> &domain_values)
1526  {
1527  if (this->GetType() == "cost_functional")
1528  {
1529  // state values in quadrature points
1530  return GetFunctional()->PointValue(
1531  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1532  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1533  domain_values);
1534 
1535  } //endif cost_functional
1536  else if (this->GetType() == "aux_functional")
1537  {
1538  // state values in quadrature points
1539  return _aux_functionals[this->GetTypeNum()]->PointValue(
1540  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1541  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1542  domain_values);
1543 
1544  } //endif aux_functional
1545  else if (this->GetType() == "functional_for_ee")
1546  {
1547  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1548  return _aux_functionals[_functional_for_ee_num]->PointValue(
1549  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1550  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1551  domain_values);
1552  } //endif functional_for_ee
1553  else if (this->GetType().find("constraints") != std::string::npos)
1554  {
1555  return GetConstraints()->PointValue(
1556  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
1557  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
1558  domain_values);
1559 
1560  } //endif constraints
1561  else
1562  {
1563  throw DOpEException("Not implemented",
1564  "OptProblemContainer::PointFunctional");
1565  }
1566  }
1567 
1568  /******************************************************/
1569 
1570  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1571  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1572  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1573  template<int, int> class DH>
1574  template<typename FACEDATACONTAINER>
1575  double
1576  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1577  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryFunctional(
1578  const FACEDATACONTAINER& fdc)
1579  {
1580  if (this->GetType() == "cost_functional")
1581  {
1582  // state values in quadrature points
1583  return GetFunctional()->BoundaryValue(fdc);
1584  }
1585  else if (this->GetType() == "aux_functional")
1586  {
1587  // state values in quadrature points
1588  return _aux_functionals[this->GetTypeNum()]->BoundaryValue(fdc);
1589  }
1590  else if (this->GetType() == "functional_for_ee")
1591  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1592  {
1593  return _aux_functionals[_functional_for_ee_num]->BoundaryValue(fdc);
1594  }
1595  else if (this->GetType().find("constraints") != std::string::npos)
1596  {
1597  return GetConstraints()->BoundaryValue(fdc);
1598  }
1599  else
1600  {
1601  throw DOpEException("Not implemented",
1602  "OptProblemContainer::BoundaryFunctional");
1603  }
1604  }
1605 
1606  /******************************************************/
1607 
1608  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1609  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1610  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1611  template<int, int> class DH>
1612  template<typename FACEDATACONTAINER>
1613  double
1614  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1615  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceFunctional(
1616  const FACEDATACONTAINER& fdc)
1617  {
1618  if (this->GetType() == "cost_functional")
1619  {
1620  // state values in quadrature points
1621  return GetFunctional()->FaceValue(fdc);
1622  }
1623  else if (this->GetType() == "aux_functional")
1624  {
1625  // state values in quadrature points
1626  return _aux_functionals[this->GetTypeNum()]->FaceValue(fdc);
1627  }
1628  else if (this->GetType() == "functional_for_ee")
1629  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1630  {
1631  return _aux_functionals[_functional_for_ee_num]->FaceValue(fdc);
1632  }
1633  else
1634  {
1635  throw DOpEException("Not implemented",
1636  "OptProblemContainer::FaceFunctional");
1637  }
1638  }
1639 
1640  /******************************************************/
1641 
1642  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1643  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1644  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1645  template<int, int> class DH>
1646  double
1647  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1648  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AlgebraicFunctional(
1649  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1650  const std::map<std::string, const VECTOR*> &domain_values)
1651  {
1652  if (this->GetType() == "cost_functional")
1653  {
1654  // state values in quadrature points
1655  return GetFunctional()->AlgebraicValue(param_values, domain_values);
1656  }
1657  else if (this->GetType() == "aux_functional")
1658  {
1659  // state values in quadrature points
1660  return _aux_functionals[this->GetTypeNum()]->AlgebraicValue(
1661  param_values, domain_values);
1662  }
1663  else if (this->GetType() == "functional_for_ee")
1664  //TODO ist das hier korrekt? Sollten wir eigentlich nicht benoetigen.
1665  {
1666  return _aux_functionals[_functional_for_ee_num]->AlgebraicValue(
1667  param_values, domain_values);
1668  }
1669  else
1670  {
1671  throw DOpEException("Not implemented",
1672  "OptProblemContainer::AlgebraicFunctional");
1673  }
1674  }
1675 
1676  /******************************************************/
1677 
1678  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1679  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1680  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1681  template<int, int> class DH>
1682  template<typename DATACONTAINER>
1683  void
1684  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1685  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementEquation(
1686  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1687  double scale, double scale_ico)
1688  {
1689 
1690  if (this->GetType() == "state")
1691  {
1692  // state values in quadrature points
1693  this->GetPDE().ElementEquation(edc, local_vector, scale, scale_ico);
1694  }
1695  else if ((this->GetType() == "adjoint")
1696  || (this->GetType() == "adjoint_for_ee"))
1697  {
1698  // state values in quadrature points
1699  this->GetPDE().ElementEquation_U(edc, local_vector, scale,
1700  scale_ico);
1701  }
1702  else if (this->GetType() == "adjoint_hessian")
1703  {
1704  // state values in quadrature points
1705  this->GetPDE().ElementEquation_UTT(edc, local_vector, scale,
1706  scale_ico);
1707  }
1708  else if (this->GetType() == "tangent")
1709  {
1710  // state values in quadrature points
1711  this->GetPDE().ElementEquation_UT(edc, local_vector, scale,
1712  scale_ico);
1713  }
1714  else if ((this->GetType() == "gradient")
1715  || (this->GetType() == "hessian"))
1716  {
1717  // control values in quadrature points
1718  this->GetPDE().ControlElementEquation(edc, local_vector, scale);
1719  }
1720  else
1721  {
1722  throw DOpEException("Not implemented",
1723  "OptProblemContainer::ElementEquation");
1724  }
1725  }
1726 
1727  /******************************************************/
1728 
1729  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1730  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1731  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1732  template<int, int> class DH>
1733  void
1734  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
1735  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AlgebraicResidual(
1736  VECTOR& residual,
1737  const std::map<std::string, const dealii::Vector<double>*> &param_values,
1738  const std::map<std::string, const VECTOR*> &domain_values)
1739  {
1740  if (this->GetType() == "gradient")
1741  {
1742  // state values in quadrature points
1743  return GetFunctional()->AlgebraicGradient_Q(residual, param_values,
1744  domain_values);
1745  }
1746  else
1747  {
1748  throw DOpEException("Not implemented",
1749  "OptProblemContainer::AlgebraicFunctional");
1750  }
1751  }
1752 
1753  /******************************************************/
1754 
1755  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1756  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1757  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1758  template<int, int> class DH>
1759  template<typename DATACONTAINER>
1760  void
1761  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1762  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeEquation(
1763  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1764  double scale)
1765  {
1766 
1767  if (this->GetType() == "state")
1768  {
1769  this->GetPDE().ElementTimeEquation(edc, local_vector, scale);
1770  }
1771  else if (this->GetType() == "adjoint"
1772  || this->GetType() == "adjoint_for_ee")
1773  {
1774  this->GetPDE().ElementTimeEquation_U(edc, local_vector, scale);
1775  }
1776  else if (this->GetType() == "adjoint_hessian")
1777  {
1778  this->GetPDE().ElementTimeEquation_UTT(edc, local_vector, scale);
1779  }
1780  else if (this->GetType() == "tangent")
1781  {
1782  this->GetPDE().ElementTimeEquation_UT(edc, local_vector, scale);
1783  }
1784  else if ((this->GetType() == "gradient")
1785  || (this->GetType() == "hessian"))
1786  {
1787  throw DOpEException("Not implemented",
1788  "OptProblemContainer::ElementTimeEquation");
1789  }
1790  else
1791  {
1792  throw DOpEException("Not implemented",
1793  "OptProblemContainer::ElementTimeEquation");
1794  }
1795  }
1796 
1797  /******************************************************/
1798 
1799  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1800  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1801  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1802  template<int, int> class DH>
1803  template<typename DATACONTAINER>
1804  void
1805  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1806  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeEquationExplicit(
1807  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1808  double scale)
1809  {
1810 
1811  if (this->GetType() == "state")
1812  {
1813  this->GetPDE().ElementTimeEquationExplicit(edc, local_vector,
1814  scale);
1815  }
1816  else if (this->GetType() == "adjoint"
1817  || this->GetType() == "adjoint_for_ee")
1818  {
1819  this->GetPDE().ElementTimeEquationExplicit_U(edc, local_vector,
1820  scale);
1821  }
1822  else if (this->GetType() == "adjoint_hessian")
1823  {
1824  this->GetPDE().ElementTimeEquationExplicit_UTT(edc, local_vector,
1825  scale);
1826  }
1827  else if (this->GetType() == "tangent")
1828  {
1829  this->GetPDE().ElementTimeEquationExplicit_UT(edc, local_vector,
1830  scale);
1831  }
1832  else if ((this->GetType() == "gradient")
1833  || (this->GetType() == "hessian"))
1834  {
1835  throw DOpEException("Not implemented",
1836  "OptProblemContainer::ElementTimeEquationExplicit");
1837  }
1838  else
1839  {
1840  throw DOpEException("Not implemented",
1841  "OptProblemContainer::ElementTimeEquationExplicit");
1842  }
1843  }
1844 
1845  /******************************************************/
1846 
1847  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1848  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1849  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1850  template<int, int> class DH>
1851  template<typename FACEDATACONTAINER>
1852  void
1853  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1854  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceEquation(
1855  const FACEDATACONTAINER& fdc,
1856  dealii::Vector<double> &local_vector, double scale,
1857  double scale_ico)
1858  {
1859  if (this->GetType() == "state")
1860  {
1861  // state values in quadrature points
1862  this->GetPDE().FaceEquation(fdc, local_vector, scale, scale_ico);
1863  }
1864  else if (this->GetType() == "adjoint"
1865  || this->GetType() == "adjoint_for_ee")
1866  {
1867  // state values in quadrature points
1868  this->GetPDE().FaceEquation_U(fdc, local_vector, scale,
1869  scale_ico);
1870  }
1871  else if (this->GetType() == "adjoint_hessian")
1872  {
1873  // state values in quadrature points
1874  this->GetPDE().FaceEquation_UTT(fdc, local_vector, scale,
1875  scale_ico);
1876  }
1877  else if (this->GetType() == "tangent")
1878  {
1879  // state values in quadrature points
1880  this->GetPDE().FaceEquation_UT(fdc, local_vector, scale,
1881  scale_ico);
1882  }
1883 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
1884 // {
1885 // // control values in quadrature points
1886 // this->GetPDE().ControlFaceEquation(fdc, local_vector, scale);
1887 // }
1888  else
1889  {
1890  throw DOpEException("Not implemented",
1891  "OptProblemContainer::FaceEquation");
1892  }
1893  }
1894 
1895  /******************************************************/
1896 
1897  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1898  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1899  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1900  template<int, int> class DH>
1901  template<typename FACEDATACONTAINER>
1902  void
1903  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1904  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::InterfaceEquation(
1905  const FACEDATACONTAINER& fdc,
1906  dealii::Vector<double> &local_vector, double scale,
1907  double scale_ico)
1908  {
1909  if (this->GetType() == "state")
1910  {
1911  this->GetPDE().InterfaceEquation(fdc, local_vector, scale,
1912  scale_ico);
1913  }
1914  else if (this->GetType() == "adjoint"
1915  || this->GetType() == "adjoint_for_ee")
1916  {
1917  // state values in quadrature points
1918  this->GetPDE().InterfaceEquation_U(fdc, local_vector, scale,
1919  scale_ico);
1920  }
1921  else
1922  {
1923  throw DOpEException("Not implemented",
1924  "OptProblemContainer::InterfaceEquation");
1925  }
1926  }
1927  /******************************************************/
1928 
1929  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1930  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1931  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1932  template<int, int> class DH>
1933  template<typename FACEDATACONTAINER>
1934  void
1935  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1936  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryEquation(
1937  const FACEDATACONTAINER& fdc,
1938  dealii::Vector<double> &local_vector, double scale,
1939  double scale_ico)
1940  {
1941  if (this->GetType() == "state")
1942  {
1943  // state values in quadrature points
1944  this->GetPDE().BoundaryEquation(fdc, local_vector, scale,
1945  scale_ico);
1946  }
1947  else if (this->GetType() == "adjoint"
1948  || this->GetType() == "adjoint_for_ee")
1949  {
1950  // state values in quadrature points
1951  this->GetPDE().BoundaryEquation_U(fdc, local_vector, scale,
1952  scale_ico);
1953  }
1954  else if (this->GetType() == "adjoint_hessian")
1955  {
1956  // state values in quadrature points
1957  this->GetPDE().BoundaryEquation_UTT(fdc, local_vector, scale,
1958  scale_ico);
1959  }
1960  else if (this->GetType() == "tangent")
1961  {
1962  // state values in quadrature points
1963  this->GetPDE().BoundaryEquation_UT(fdc, local_vector, scale,
1964  scale_ico);
1965  }
1966 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
1967 // {
1968 // // control values in quadrature points
1969 // this->GetPDE().ControlBoundaryEquation(fdc, local_vector, scale);
1970 // }
1971  else
1972  {
1973  throw DOpEException("Not implemented",
1974  "OptProblemContainer::ElementBoundaryEquation");
1975  }
1976  }
1977 
1978  /******************************************************/
1979 
1980  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
1981  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
1982  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
1983  template<int, int> class DH>
1984  template<typename DATACONTAINER>
1985  void
1986  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
1987  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementRhs(
1988  const DATACONTAINER& edc, dealii::Vector<double> &local_vector,
1989  double scale)
1990  {
1991  //FIXME: In timedependent Problems one should evaluate the
1992  // Functional terms here only if they are time distributed,
1993  //otherwise the scaling should be different (e.g. +-1 for local in time
1994  //Functionals)
1995  //The same applies to all ...Rhs functions, except init_...
1996 
1997  if (this->GetType() == "state")
1998  {
1999  // state values in quadrature points
2000  this->GetPDE().ElementRightHandSide(edc, local_vector, scale);
2001  }
2002  else if (this->GetType() == "adjoint")
2003  {
2004  // state values in quadrature points
2005  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2006  GetFunctional()->ElementValue_U(edc, local_vector, scale);
2007  }
2008  else if (this->GetType() == "adjoint_for_ee")
2009  {
2010  //values of the derivative of the functional for error estimation
2011  _aux_functionals[_functional_for_ee_num]->ElementValue_U(edc,
2012  local_vector, scale);
2013  }
2014  else if (this->GetType() == "tangent")
2015  {
2016  // state values in quadrature points
2017  scale *= -1;
2018  this->GetPDE().ElementEquation_QT(edc, local_vector, scale, scale);
2019  }
2020  else if (this->GetType() == "adjoint_hessian")
2021  {
2022  // state values in quadrature points
2023  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2024  {
2025  GetFunctional()->ElementValue_UU(edc, local_vector, scale);
2026  GetFunctional()->ElementValue_QU(edc, local_vector, scale);
2027  }
2028  scale *= -1;
2029  this->GetPDE().ElementEquation_UU(edc, local_vector, scale, scale);
2030  this->GetPDE().ElementEquation_QU(edc, local_vector, scale, scale);
2031  //TODO: make some example where this realy matters to check if this is right
2032  this->GetPDE().ElementTimeEquationExplicit_UU(edc, local_vector,
2033  scale);
2034  }
2035  else if (this->GetType() == "gradient")
2036  {
2037  if (GetSpaceTimeHandler()->GetControlType()
2038  == DOpEtypes::ControlType::initial)
2039  {
2040  this->GetPDE().Init_ElementRhs_Q(edc, local_vector, scale);
2041  }
2042  // state values in quadrature points
2043  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2044  GetFunctional()->ElementValue_Q(edc, local_vector, scale);
2045 
2046  scale *= -1;
2047  this->GetPDE().ElementEquation_Q(edc, local_vector, scale, scale);
2048  }
2049  else if (this->GetType() == "hessian")
2050  {
2051  if (GetSpaceTimeHandler()->GetControlType()
2052  == DOpEtypes::ControlType::initial)
2053  {
2054  this->GetPDE().Init_ElementRhs_QTT(edc, local_vector, scale);
2055  this->GetPDE().Init_ElementRhs_QQ(edc, local_vector, scale);
2056  }
2057 
2058  if (GetFunctional()->GetType().find("domain") != std::string::npos)
2059  {
2060  GetFunctional()->ElementValue_QQ(edc, local_vector, scale);
2061  GetFunctional()->ElementValue_UQ(edc, local_vector, scale);
2062  }
2063 
2064  scale *= -1;
2065  this->GetPDE().ElementEquation_QTT(edc, local_vector, scale, scale);
2066  this->GetPDE().ElementEquation_UQ(edc, local_vector, scale, scale);
2067  this->GetPDE().ElementEquation_QQ(edc, local_vector, scale, scale);
2068  }
2069  else if (this->GetType() == "global_constraint_gradient")
2070  {
2071  GetConstraints()->ElementValue_Q(edc, local_vector, scale);
2072  }
2073  else if (this->GetType() == "global_constraint_hessian")
2074  {
2075  GetConstraints()->ElementValue_QQ(edc, local_vector, scale);
2076  }
2077  else
2078  {
2079  throw DOpEException("Not implemented",
2080  "OptProblemContainer::ElementRhs");
2081  }
2082 
2083  }
2084 
2085  /******************************************************/
2086 
2087  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2088  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2089  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2090  template<int, int> class DH>
2091  void
2092  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2093  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PointRhs(
2094  const std::map<std::string, const dealii::Vector<double>*> &param_values,
2095  const std::map<std::string, const VECTOR*> &domain_values,
2096  VECTOR& rhs_vector, double scale)
2097  {
2098 
2099  if (this->GetType() == "adjoint")
2100  {
2101  // state values in quadrature points
2102  if (GetFunctional()->GetType().find("point") != std::string::npos)
2103  {
2104  GetFunctional()->PointValue_U(
2105  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2106  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2107  domain_values, rhs_vector, scale);
2108  }
2109  }
2110  else if (this->GetType() == "adjoint_for_ee")
2111  {
2112  //values of the derivative of the functional for error estimation
2113  _aux_functionals[_functional_for_ee_num]->PointValue_U(
2114  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2115  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2116  domain_values, rhs_vector, scale);
2117  }
2118  else if (this->GetType() == "adjoint_hessian")
2119  {
2120  // state values in quadrature points
2121  if (GetFunctional()->GetType().find("point") != std::string::npos)
2122  {
2123  GetFunctional()->PointValue_UU(
2124  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2125  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2126  domain_values, rhs_vector, scale);
2127  GetFunctional()->PointValue_QU(
2128  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2129  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2130  domain_values, rhs_vector, scale);
2131  }
2132  }
2133  else if (this->GetType() == "gradient")
2134  {
2135  // state values in quadrature points
2136  if (GetFunctional()->GetType().find("point") != std::string::npos)
2137  {
2138  GetFunctional()->PointValue_Q(
2139  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2140  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2141  domain_values, rhs_vector, scale);
2142  }
2143  }
2144  else if (this->GetType() == "hessian")
2145  {
2146  // state values in quadrature points
2147  if (GetFunctional()->GetType().find("point") != std::string::npos)
2148  {
2149  GetFunctional()->PointValue_QQ(
2150  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2151  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2152  domain_values, rhs_vector, scale);
2153  GetFunctional()->PointValue_UQ(
2154  this->GetSpaceTimeHandler()->GetControlDoFHandler(),
2155  this->GetSpaceTimeHandler()->GetStateDoFHandler(), param_values,
2156  domain_values, rhs_vector, scale);
2157  }
2158  }
2159  else
2160  {
2161  throw DOpEException("Not implemented", "OptProblem::PointRhs");
2162  }
2163 
2164  }
2165 
2166  /******************************************************/
2167 
2168  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2169  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2170  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2171  template<int, int> class DH>
2172  template<typename FACEDATACONTAINER>
2173  void
2174  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2175  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceRhs(
2176  const FACEDATACONTAINER& fdc,
2177  dealii::Vector<double> &local_vector, double scale)
2178  {
2179 
2180  if (this->GetType() == "state")
2181  {
2182  // state values in face quadrature points
2183  this->GetPDE().FaceRightHandSide(fdc, local_vector, scale);
2184  }
2185  else if (this->GetType() == "adjoint")
2186  {
2187  // state values in quadrature points
2188  if (GetFunctional()->GetType().find("face") != std::string::npos)
2189  GetFunctional()->FaceValue_U(fdc, local_vector, scale);
2190  }
2191  else if (this->GetType() == "adjoint_for_ee")
2192  {
2193  //values of the derivative of the functional for error estimation
2194  if (_aux_functionals[_functional_for_ee_num]->GetType().find("face")
2195  != std::string::npos)
2196  _aux_functionals[_functional_for_ee_num]->FaceValue_U(fdc,
2197  local_vector, scale);
2198  }
2199  else if (this->GetType() == "tangent")
2200  {
2201  // state values in quadrature points
2202  scale *= -1;
2203  this->GetPDE().FaceEquation_QT(fdc, local_vector, scale, scale);
2204  }
2205  else if (this->GetType() == "adjoint_hessian")
2206  {
2207  // state values in quadrature points
2208  if (GetFunctional()->GetType().find("face") != std::string::npos)
2209  {
2210  GetFunctional()->FaceValue_UU(fdc, local_vector, scale);
2211  GetFunctional()->FaceValue_QU(fdc, local_vector, scale);
2212  }
2213 
2214  scale *= -1;
2215  this->GetPDE().FaceEquation_UU(fdc, local_vector, scale, scale);
2216 
2217  this->GetPDE().FaceEquation_QU(fdc, local_vector, scale, scale);
2218  }
2219  else if (this->GetType() == "gradient")
2220  {
2221  // state values in quadrature points
2222  if (GetFunctional()->GetType().find("face") != std::string::npos)
2223  {
2224  GetFunctional()->FaceValue_Q(fdc, local_vector, scale);
2225  }
2226 
2227  scale *= -1;
2228  this->GetPDE().FaceEquation_Q(fdc, local_vector, scale, scale);
2229  }
2230  else if (this->GetType() == "hessian")
2231  {
2232  // state values in quadrature points
2233  if (GetFunctional()->GetType().find("face") != std::string::npos)
2234  {
2235  GetFunctional()->FaceValue_QQ(fdc, local_vector, scale);
2236  GetFunctional()->FaceValue_UQ(fdc, local_vector, scale);
2237  }
2238 
2239  scale *= -1;
2240  this->GetPDE().FaceEquation_QTT(fdc, local_vector, scale, scale);
2241 
2242  this->GetPDE().FaceEquation_UQ(fdc, local_vector, scale, scale);
2243 
2244  this->GetPDE().FaceEquation_QQ(fdc, local_vector, scale, scale);
2245  }
2246  else
2247  {
2248  throw DOpEException("Not implemented",
2249  "OptProblemContainer::FaceRhs");
2250  }
2251 
2252  }
2253 
2254  /******************************************************/
2255 
2256  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2257  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2258  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2259  template<int, int> class DH>
2260  template<typename FACEDATACONTAINER>
2261  void
2262  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2263  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryRhs(
2264  const FACEDATACONTAINER& fdc,
2265  dealii::Vector<double> &local_vector, double scale)
2266  {
2267 
2268  if (this->GetType() == "state")
2269  {
2270  // state values in face quadrature points
2271  this->GetPDE().BoundaryRightHandSide(fdc, local_vector, scale);
2272  }
2273  else if (this->GetType() == "adjoint")
2274  {
2275  // state values in quadrature points
2276  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2277  GetFunctional()->BoundaryValue_U(fdc, local_vector, scale);
2278  }
2279  else if (this->GetType() == "adjoint_for_ee")
2280  {
2281  //values of the derivative of the functional for error estimation
2282  if (_aux_functionals[_functional_for_ee_num]->GetType().find(
2283  "boundary") != std::string::npos)
2284  _aux_functionals[_functional_for_ee_num]->BoundaryValue_U(fdc,
2285  local_vector, scale);
2286  }
2287  else if (this->GetType() == "tangent")
2288  {
2289  // state values in quadrature points
2290  scale *= -1;
2291  this->GetPDE().BoundaryEquation_QT(fdc, local_vector, scale,
2292  scale);
2293  }
2294  else if (this->GetType() == "adjoint_hessian")
2295  {
2296  // state values in quadrature points
2297  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2298  {
2299  GetFunctional()->BoundaryValue_UU(fdc, local_vector, scale);
2300  GetFunctional()->BoundaryValue_QU(fdc, local_vector, scale);
2301  }
2302 
2303  scale *= -1;
2304  this->GetPDE().BoundaryEquation_UU(fdc, local_vector, scale,
2305  scale);
2306  this->GetPDE().BoundaryEquation_QU(fdc, local_vector, scale,
2307  scale);
2308  }
2309  else if (this->GetType() == "gradient")
2310  {
2311  // state values in quadrature points
2312  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2313  GetFunctional()->BoundaryValue_Q(fdc, local_vector, scale);
2314 
2315  scale *= -1;
2316  this->GetPDE().BoundaryEquation_Q(fdc, local_vector, scale,
2317  scale);
2318  }
2319  else if (this->GetType() == "hessian")
2320  {
2321  // state values in quadrature points
2322  if (GetFunctional()->GetType().find("boundary") != std::string::npos)
2323  {
2324  GetFunctional()->BoundaryValue_QQ(fdc, local_vector, scale);
2325  GetFunctional()->BoundaryValue_UQ(fdc, local_vector, scale);
2326  }
2327 
2328  scale *= -1;
2329  this->GetPDE().BoundaryEquation_QTT(fdc, local_vector, scale,
2330  scale);
2331  this->GetPDE().BoundaryEquation_UQ(fdc, local_vector, scale,
2332  scale);
2333  this->GetPDE().BoundaryEquation_QQ(fdc, local_vector, scale,
2334  scale);
2335  }
2336  else
2337  {
2338  throw DOpEException("Not implemented",
2339  "OptProblemContainer::BoundaryRhs");
2340  }
2341 
2342  }
2343 
2344  /******************************************************/
2345 
2346  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2347  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2348  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2349  template<int, int> class DH>
2350  template<typename DATACONTAINER>
2351  void
2352  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2353  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementMatrix(
2354  const DATACONTAINER& edc,
2355  dealii::FullMatrix<double> &local_entry_matrix, double scale,
2356  double scale_ico)
2357  {
2358 
2359  if (this->GetType() == "state" || this->GetType() == "tangent")
2360  {
2361  // state values in quadrature points
2362  this->GetPDE().ElementMatrix(edc, local_entry_matrix, scale, scale_ico);
2363  }
2364  else if (this->GetType() == "adjoint"
2365  || this->GetType() == "adjoint_for_ee"
2366  || this->GetType() == "adjoint_hessian")
2367  {
2368  // state values in quadrature points
2369  this->GetPDE().ElementMatrix_T(edc, local_entry_matrix, scale,
2370  scale_ico);
2371  }
2372  else if ((this->GetType() == "gradient")
2373  || (this->GetType() == "hessian"))
2374  {
2375  // control values in quadrature points
2376  this->GetPDE().ControlElementMatrix(edc, local_entry_matrix);
2377  }
2378  else
2379  {
2380  throw DOpEException("Not implemented",
2381  "OptProblemContainer::ElementMatrix");
2382  }
2383 
2384  }
2385 
2386  /******************************************************/
2387 
2388  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2389  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2390  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2391  template<int, int> class DH>
2392  template<typename DATACONTAINER>
2393  void
2394  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2395  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeMatrix(
2396  const DATACONTAINER& edc, FullMatrix<double> &local_entry_matrix)
2397  {
2398 
2399  if (this->GetType() == "state" || this->GetType() == "tangent")
2400  {
2401  // state values in quadrature points
2402  this->GetPDE().ElementTimeMatrix(edc, local_entry_matrix);
2403  }
2404  else if (this->GetType() == "adjoint"
2405  || this->GetType() == "adjoint_for_ee"
2406  || this->GetType() == "adjoint_hessian")
2407  {
2408  this->GetPDE().ElementTimeMatrix_T(edc, local_entry_matrix);
2409  }
2410  else if ((this->GetType() == "gradient")
2411  || (this->GetType() == "hessian"))
2412  {
2413  throw DOpEException("Not implemented",
2414  "OptProblemContainer::ElementTimeMatrix");
2415  }
2416  else
2417  {
2418  throw DOpEException("Not implemented",
2419  "OptProblemContainer::ElementTimeMatrix");
2420  }
2421 
2422  }
2423 
2424  /******************************************************/
2425 
2426  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2427  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2428  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2429  template<int, int> class DH>
2430  template<typename DATACONTAINER>
2431  void
2432  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2433  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ElementTimeMatrixExplicit(
2434  const DATACONTAINER& edc,
2435  dealii::FullMatrix<double> &local_entry_matrix)
2436  {
2437 
2438  if (this->GetType() == "state" || this->GetType() == "tangent")
2439  {
2440  this->GetPDE().ElementTimeMatrixExplicit(edc, local_entry_matrix);
2441  }
2442  else if (this->GetType() == "adjoint"
2443  || this->GetType() == "adjoint_hessian")
2444  {
2445  this->GetPDE().ElementTimeMatrixExplicit_T(edc, local_entry_matrix);
2446  }
2447  else if ((this->GetType() == "gradient")
2448  || (this->GetType() == "hessian"))
2449  {
2450  throw DOpEException("Not implemented",
2451  "OptProblemContainer::ElementTimeMatrixExplicit");
2452  }
2453  else
2454  {
2455  throw DOpEException("Not implemented",
2456  "OptProblemContainer::ElementTimeMatrixExplicit");
2457  }
2458 
2459  }
2460 
2461  /******************************************************/
2462 
2463  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2464  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2465  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2466  template<int, int> class DH>
2467  template<typename FACEDATACONTAINER>
2468  void
2469  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2470  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::FaceMatrix(
2471  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
2472  double scale, double scale_ico)
2473  {
2474  if (this->GetType() == "state" || this->GetType() == "tangent")
2475  {
2476  // state values in face quadrature points
2477  this->GetPDE().FaceMatrix(fdc, local_entry_matrix, scale, scale_ico);
2478  }
2479  else if (this->GetType() == "adjoint"
2480  || this->GetType() == "adjoint_for_ee"
2481  || this->GetType() == "adjoint_hessian")
2482  {
2483  // state values in quadrature points
2484  this->GetPDE().FaceMatrix_T(fdc, local_entry_matrix, scale,
2485  scale_ico);
2486  }
2487 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
2488 // {
2489 // // control values in quadrature points
2490 // this->GetPDE().ControlFaceMatrix(fdc, local_entry_matrix);
2491 // }
2492  else
2493  {
2494  throw DOpEException("Not implemented",
2495  "OptProblemContainer::NewtonFaceMatrix");
2496  }
2497 
2498  }
2499 
2500  /******************************************************/
2501 
2502  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2503  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2504  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2505  template<int, int> class DH>
2506  template<typename FACEDATACONTAINER>
2507  void
2508  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2509  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::InterfaceMatrix(
2510  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_entry_matrix,
2511  double scale, double scale_ico)
2512  {
2513  if (this->GetType() == "state")
2514  {
2515  this->GetPDE().InterfaceMatrix(fdc, local_entry_matrix, scale,
2516  scale_ico);
2517  }
2518  else if (this->GetType() == "adjoint"
2519  || this->GetType() == "adjoint_for_ee"
2520  || this->GetType() == "adjoint_hessian")
2521  {
2522  this->GetPDE().InterfaceMatrix_T(fdc, local_entry_matrix, scale,
2523  scale_ico);
2524  }
2525  else
2526  {
2527  throw DOpEException("Not implemented",
2528  "OptProblemContainer::NewtonInterfaceMatrix");
2529  }
2530  }
2531 
2532  /******************************************************/
2533 
2534  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2535  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2536  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2537  template<int, int> class DH>
2538  template<typename FACEDATACONTAINER>
2539  void
2540  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD,
2541  CONSTRAINTS, SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::BoundaryMatrix(
2542  const FACEDATACONTAINER& fdc, FullMatrix<double> &local_matrix,
2543  double scale, double scale_ico)
2544  {
2545  if (this->GetType() == "state" || this->GetType() == "tangent")
2546  {
2547  // state values in face quadrature points
2548  this->GetPDE().BoundaryMatrix(fdc, local_matrix, scale,
2549  scale_ico);
2550  }
2551  else if (this->GetType() == "adjoint"
2552  || this->GetType() == "adjoint_for_ee"
2553  || this->GetType() == "adjoint_hessian")
2554  {
2555  // state values in quadrature points
2556  this->GetPDE().BoundaryMatrix_T(fdc, local_matrix, scale,
2557  scale_ico);
2558  }
2559 // else if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
2560 // {
2561 // // control values in quadrature points
2562 // this->GetPDE().ControlBoundaryMatrix(fdc, local_matrix);
2563 // }
2564  else
2565  {
2566  throw DOpEException("Not implemented",
2567  "OptProblemContainer::ElementBoundaryMatrix");
2568  }
2569 
2570  }
2571 
2572  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2573  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2574  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2575  template<int, int> class DH>
2576  void
2577  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2578  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ComputeLocalControlConstraints(
2579  VECTOR& constraints,
2580  const std::map<std::string, const dealii::Vector<double>*> &/*values*/,
2581  const std::map<std::string, const VECTOR*> &block_values)
2582  {
2583  if (this->GetType() == "constraints")
2584  {
2585  if (this->GetSpaceTimeHandler()->GetNLocalConstraints() != 0)
2586  {
2587  const VECTOR& control = *GetBlockVector(block_values, "control");
2588  this->GetConstraints()->EvaluateLocalControlConstraints(control,
2589  constraints);
2590  }
2591  }
2592  else
2593  {
2594  throw DOpEException("Wrong problem type" + this->GetType(),
2595  "OptProblemContainer::ComputeLocalConstraints");
2596  }
2597  }
2598 
2599  /******************************************************/
2600 
2601  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2602  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2603  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2604  template<int, int> class DH>
2605  std::string
2606  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2607  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFType() const
2608  {
2609  if (this->GetType() == "state" || this->GetType() == "adjoint"
2610  || this->GetType() == "adjoint_for_ee" || this->GetType() == "tangent"
2611  || this->GetType() == "adjoint_hessian")
2612  {
2613  return "state";
2614  }
2615  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
2616  || (this->GetType() == "hessian_inverse"))
2617  {
2618  return "control";
2619  }
2620  else
2621  {
2622  throw DOpEException("Unknown Type:" + this->GetType(),
2623  "OptProblemContainer::GetDoFType");
2624  }
2625  }
2626 
2627  /******************************************************/
2628 
2629  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2630  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2631  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2632  template<int, int> class DH>
2633  const FE<dealdim, dealdim>&
2634  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2635  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFESystem() const
2636  {
2637  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
2638  || this->GetType() == "adjoint_for_ee" || this->GetType() == "tangent"
2639  || this->GetType() == "adjoint_hessian")
2640  {
2641  return this->GetSpaceTimeHandler()->GetFESystem("state");
2642  }
2643  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
2644  || (this->GetType() == "global_constraint_gradient"))
2645  {
2646 #if dope_dimension > 0
2647  if(dopedim == dealdim)
2648  return this->GetSpaceTimeHandler()->GetFESystem("control");
2649  else
2650  throw DOpEException("Non matching dimensions!","OptProblemContainer::GetFESystem");
2651 #else
2652  return this->GetSpaceTimeHandler()->GetFESystem("state");
2653 #endif
2654  }
2655  else
2656  {
2657  throw DOpEException("Unknown Type:" + this->GetType(),
2658  "OptProblemContainer::GetFESystem");
2659  }
2660  }
2661 
2662  /******************************************************/
2663 
2664  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2665  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2666  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2667  template<int, int> class DH>
2668  UpdateFlags
2669  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2670  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetUpdateFlags() const
2671  {
2672 
2673  UpdateFlags r;
2674  if (this->GetType().find("aux_functional") != std::string::npos)
2675  {
2676  r = _aux_functionals[this->GetTypeNum()]->GetUpdateFlags();
2677  }
2678  else if (this->GetType().find("functional") != std::string::npos)
2679  {
2680  r = this->GetFunctional()->GetUpdateFlags();
2681  }
2682  else if (this->GetType().find("constraints") != std::string::npos)
2683  {
2684  r = this->GetConstraints()->GetUpdateFlags();
2685  }
2686  else if (this->GetType() == "functional_for_ee")
2687  {
2688  r = _aux_functionals[_functional_for_ee_num]->GetUpdateFlags();
2689  }
2690  else
2691  {
2692  r = this->GetPDE().GetUpdateFlags();
2693  if (this->GetType() == "adjoint_hessian" || this->GetType() == "adjoint"
2694  || (this->GetType() == "hessian")
2695  || (this->GetType() == "gradient"))
2696  {
2697  r = r | this->GetFunctional()->GetUpdateFlags();
2698  }
2699  else if (this->GetType() == "adjoint_for_ee")
2700  {
2701  r = r | _aux_functionals[_functional_for_ee_num]->GetUpdateFlags();
2702  }
2703  }
2704  return r | update_JxW_values;
2705  }
2706 
2707  /******************************************************/
2708 
2709  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2710  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2711  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2712  template<int, int> class DH>
2713  UpdateFlags
2714  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2715  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFaceUpdateFlags() const
2716  {
2717  UpdateFlags r;
2718  if (this->GetType().find("aux_functional") != std::string::npos)
2719  {
2720  r = _aux_functionals[this->GetTypeNum()]->GetFaceUpdateFlags();
2721  }
2722  else if (this->GetType().find("functional") != std::string::npos)
2723  {
2724  r = this->GetFunctional()->GetFaceUpdateFlags();
2725  }
2726  else if (this->GetType().find("constraints") != std::string::npos)
2727  {
2728  r = this->GetConstraints()->GetFaceUpdateFlags();
2729  }
2730  else if (this->GetType() == "functional_for_ee")
2731  {
2732  r = _aux_functionals[_functional_for_ee_num]->GetUpdateFlags();
2733  }
2734  else
2735  {
2736  r = this->GetPDE().GetFaceUpdateFlags();
2737  if (this->GetType() == "adjoint_hessian" || this->GetType() == "adjoint"
2738  || (this->GetType() == "hessian"))
2739  {
2740  r = r | this->GetFunctional()->GetFaceUpdateFlags();
2741  }
2742  else if (this->GetType() == "adjoint_for_ee")
2743  {
2744  r = r
2745  | _aux_functionals[_functional_for_ee_num]->GetFaceUpdateFlags();
2746  }
2747  }
2748  return r | update_JxW_values;
2749  }
2750 
2751  /******************************************************/
2752 
2753  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2754  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2755  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2756  template<int, int> class DH>
2757  std::string
2758  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2759  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctionalType() const
2760  {
2761  if (this->GetType() == "aux_functional")
2762  {
2763  return _aux_functionals[this->GetTypeNum()]->GetType();
2764  }
2765  else if (this->GetType() == "functional_for_ee")
2766  {
2767  return _aux_functionals[_functional_for_ee_num]->GetType();
2768  }
2769  return GetFunctional()->GetType();
2770  }
2771 
2772  /******************************************************/
2773 
2774  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2775  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2776  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2777  template<int, int> class DH>
2778  std::string
2779  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2780  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctionalName() const
2781  {
2782  if (this->GetType() == "aux_functional")
2783  {
2784  return _aux_functionals[this->GetTypeNum()]->GetName();
2785  }
2786  else if (this->GetType() == "functional_for_ee")
2787  {
2788  return _aux_functionals[_functional_for_ee_num]->GetName();
2789  }
2790  return GetFunctional()->GetName();
2791  }
2792 
2793  /******************************************************/
2794 
2795  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2796  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2797  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2798  template<int, int> class DH>
2799  std::string
2800  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2801  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetConstraintType() const
2802  {
2803  return GetConstraints()->GetType();
2804  }
2805 
2806  /******************************************************/
2807 
2808  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2809  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2810  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2811  template<int, int> class DH>
2812  void
2813  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2814  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetTime(double time,
2815  const TimeIterator& interval)
2816  {
2817  GetSpaceTimeHandler()->SetInterval(interval);
2818  // GetSpaceTimeHandler()->SetTimeDoFNumber(time_point);
2819 
2820  { //Zeit an Dirichlet Werte uebermitteln
2821  for (unsigned int i = 0;
2822  i < _transposed_control_gradient_dirichlet_values.size(); i++)
2823  _transposed_control_gradient_dirichlet_values[i]->SetTime(time);
2824  for (unsigned int i = 0;
2825  i < _transposed_control_hessian_dirichlet_values.size(); i++)
2826  _transposed_control_hessian_dirichlet_values[i]->SetTime(time);
2827  for (unsigned int i = 0; i < _primal_dirichlet_values.size(); i++)
2828  _primal_dirichlet_values[i]->SetTime(time);
2829  for (unsigned int i = 0; i < _tangent_dirichlet_values.size(); i++)
2830  _tangent_dirichlet_values[i]->SetTime(time);
2831  for (unsigned int i = 0; i < _control_dirichlet_values.size(); i++)
2832  _control_dirichlet_values[i]->SetTime(time);
2833  //Functionals
2834  GetFunctional()->SetTime(time);
2835  for (unsigned int i = 0; i < _aux_functionals.size(); i++)
2836  _aux_functionals[i]->SetTime(time);
2837  //PDE
2838  this->GetPDE().SetTime(time);
2839  }
2840  //Update Auxiliary Control, State and Constraint Vectors
2841  {
2842  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
2843  _auxiliary_controls.begin();
2844  for (; it != _auxiliary_controls.end(); it++)
2845  {
2846  it->second->SetTime(time, interval);
2847  }
2848  }
2849  {
2850  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
2851  _auxiliary_state.begin();
2852  for (; it != _auxiliary_state.end(); it++)
2853  {
2854  it->second->SetTime(time, interval);
2855  }
2856  }
2857  {
2858  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
2859  _auxiliary_constraints.begin();
2860  for (; it != _auxiliary_constraints.end(); it++)
2861  {
2862  it->second->SetTime(time, interval);
2863  }
2864  }
2865 
2866  }
2867 
2868  /******************************************************/
2869 
2870  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2871  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2872  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2873  template<int, int> class DH>
2874  void
2875  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2876  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::ComputeSparsityPattern(
2877  SPARSITYPATTERN & sparsity) const
2878  {
2879  if (this->GetType() == "state" || this->GetType() == "tangent"
2880  || this->GetType() == "adjoint_for_ee" || this->GetType() == "adjoint"
2881  || this->GetType() == "adjoint_hessian")
2882  {
2883  this->GetSpaceTimeHandler()->ComputeStateSparsityPattern(sparsity);
2884  }
2885  else if ((this->GetType() == "gradient")
2886  || (this->GetType() == "hessian"))
2887  {
2888 #if dope_dimension > 0
2889  this->GetSpaceTimeHandler()->ComputeControlSparsityPattern(sparsity);
2890 #else
2891  throw DOpEException("Wrong dimension",
2892  "OptProblemContainer::ComputeSparsityPattern");
2893 #endif
2894  }
2895  else
2896  {
2897  throw DOpEException("Unknown type " + this->GetType(),
2898  "OptProblemContainer::ComputeSparsityPattern");
2899  }
2900  }
2901 
2902  /******************************************************/
2903 
2904  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2905  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2906  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2907  template<int, int> class DH>
2908  void
2909  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2910  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::PostProcessConstraints(
2911  ConstraintVector<VECTOR>& g) const
2912  {
2913  return this->GetConstraints()->PostProcessConstraints(g);
2914  }
2915  /******************************************************/
2916 
2917  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2918  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2919  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2920  template<int, int> class DH>
2921  void
2922  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2923  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryControl(
2924  const ControlVector<VECTOR>* c, std::string name)
2925  {
2926  if (_auxiliary_controls.find(name) != _auxiliary_controls.end())
2927  {
2928  throw DOpEException(
2929  "Adding multiple Data with name " + name + " is prohibited!",
2930  "OptProblemContainer::AddAuxiliaryControl");
2931  }
2932  _auxiliary_controls.insert(
2933  std::pair<std::string, const ControlVector<VECTOR>*>(name, c));
2934  }
2935 
2936  /******************************************************/
2937 
2938  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2939  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2940  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2941  template<int, int> class DH>
2942  void
2943  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2944  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryState(
2945  const StateVector<VECTOR>* c, std::string name)
2946  {
2947  if (_auxiliary_state.find(name) != _auxiliary_state.end())
2948  {
2949  throw DOpEException(
2950  "Adding multiple Data with name " + name + " is prohibited!",
2951  "OptProblemContainer::AddAuxiliaryState");
2952  }
2953  _auxiliary_state.insert(
2954  std::pair<std::string, const StateVector<VECTOR>*>(name, c));
2955  }
2956  /******************************************************/
2957 
2958  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2959  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2960  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2961  template<int, int> class DH>
2962  void
2963  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2964  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::AddAuxiliaryConstraint(
2965  const ConstraintVector<VECTOR>* c, std::string name)
2966  {
2967  if (_auxiliary_constraints.find(name) != _auxiliary_constraints.end())
2968  {
2969  throw DOpEException(
2970  "Adding multiple Data with name " + name + " is prohibited!",
2971  "OptProblemContainer::AddAuxiliaryConstraint");
2972  }
2973  _auxiliary_constraints.insert(
2974  std::pair<std::string, const ConstraintVector<VECTOR>*>(name, c));
2975  }
2976  /******************************************************/
2977 
2978  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
2979  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
2980  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
2981  template<int, int> class DH>
2982  const ControlVector<VECTOR>*
2983  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
2984  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetAuxiliaryControl(
2985  std::string name) const
2986  {
2987  typename std::map<std::string, const ControlVector<VECTOR> *>::const_iterator it =
2988  _auxiliary_controls.find(name);
2989  if (it == _auxiliary_controls.end())
2990  {
2991  throw DOpEException("Could not find Data with name " + name,
2992  "OptProblemContainer::GetAuxiliaryControl");
2993  }
2994  return it->second;
2995  }
2996 
2997  /******************************************************/
2998 
2999  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3000  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3001  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3002  template<int, int> class DH>
3003  const StateVector<VECTOR>*
3004  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3005  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetAuxiliaryState(
3006  std::string name) const
3007  {
3008  typename std::map<std::string, const StateVector<VECTOR> *>::const_iterator it =
3009  _auxiliary_state.find(name);
3010  if (it == _auxiliary_state.end())
3011  {
3012  throw DOpEException("Could not find Data with name " + name,
3013  "OptProblemContainer::GetAuxiliaryState");
3014  }
3015  return it->second;
3016  }
3017 
3018  /******************************************************/
3019 
3020  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3021  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3022  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3023  template<int, int> class DH>
3024  void
3025  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3026  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryControl(
3027  std::string name)
3028  {
3029  typename std::map<std::string, const ControlVector<VECTOR> *>::iterator it =
3030  _auxiliary_controls.find(name);
3031  if (it == _auxiliary_controls.end())
3032  {
3033  throw DOpEException(
3034  "Deleting Data " + name + " is impossible! Data not found",
3035  "OptProblemContainer::DeleteAuxiliaryControl");
3036  }
3037  _auxiliary_controls.erase(it);
3038  }
3039 
3040  /******************************************************/
3041 
3042  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3043  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3044  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3045  template<int, int> class DH>
3046  void
3047  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3048  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryState(
3049  std::string name)
3050  {
3051  typename std::map<std::string, const StateVector<VECTOR> *>::iterator it =
3052  _auxiliary_state.find(name);
3053  if (it == _auxiliary_state.end())
3054  {
3055  throw DOpEException(
3056  "Deleting Data " + name + " is impossible! Data not found",
3057  "OptProblemContainer::DeleteAuxiliaryState");
3058  }
3059  _auxiliary_state.erase(it);
3060  }
3061 
3062  /******************************************************/
3063 
3064  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3065  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3066  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3067  template<int, int> class DH>
3068  void
3069  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3070  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::DeleteAuxiliaryConstraint(
3071  std::string name)
3072  {
3073  typename std::map<std::string, const ConstraintVector<VECTOR> *>::iterator it =
3074  _auxiliary_constraints.find(name);
3075  if (it == _auxiliary_constraints.end())
3076  {
3077  throw DOpEException(
3078  "Deleting Data " + name + " is impossible! Data not found",
3079  "OptProblemContainer::DeleteAuxiliaryConstraint");
3080  }
3081  _auxiliary_constraints.erase(it);
3082  }
3083 
3084  /******************************************************/
3085 
3086  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3087  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3088  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3089  template<int, int> class DH>
3090  FUNCTIONAL*
3091  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3092  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctional()
3093  {
3094  if (this->GetType() == "aux_functional"
3095  || this->GetType() == "functional_for_ee")
3096  {
3097  //This may no longer happen!
3098  abort();
3099  // return _aux_functionals[this->GetTypeNum()];
3100  }
3101  return _functional;
3102 
3103  }
3104 
3105  /******************************************************/
3106 
3107  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3108  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3109  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3110  template<int, int> class DH>
3111  const FUNCTIONAL*
3112  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3113  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetFunctional() const
3114  {
3115  if (this->GetType() == "aux_functional"
3116  || this->GetType() == "functional_for_ee")
3117  {
3118  //This may no longer happen!
3119  abort();
3120  // return _aux_functionals[this->GetTypeNum()];
3121  }
3122  return _functional;
3123  }
3124 
3125  /******************************************************/
3126 
3127  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3128  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3129  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3130  template<int, int> class DH>
3131  bool
3132  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3133  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasFaces() const
3134  {
3135  if (this->GetType().find("aux_functional") != std::string::npos)
3136  {
3137  return _aux_functionals[this->GetTypeNum()]->HasFaces();
3138  }
3139  else if (this->GetType().find("functional") != std::string::npos)
3140  {
3141  return this->GetFunctional()->HasFaces();
3142  }
3143  else if (this->GetType().find("constraint") != std::string::npos)
3144  {
3145  return this->GetConstraints()->HasFaces();
3146  }
3147  else
3148  {
3149  if ((this->GetType() == "state") || (this->GetType() == "tangent")
3150  || (this->GetType() == "gradient"))
3151  {
3152  return this->GetPDE().HasFaces();
3153  }
3154  else if ((this->GetType() == "adjoint")
3155  || (this->GetType() == "adjoint_hessian")
3156  || (this->GetType() == "hessian"))
3157  {
3158  return this->GetPDE().HasFaces() || this->GetFunctional()->HasFaces();
3159  }
3160  else if (this->GetType() == "adjoint_for_ee")
3161  {
3162  return this->GetPDE().HasFaces()
3163  || _aux_functionals[_functional_for_ee_num]->HasFaces();
3164  }
3165  else
3166  {
3167  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3168  "OptProblemContainer::HasFaces");
3169  }
3170  }
3171  }
3172 
3173  /******************************************************/
3174 
3175  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3176  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3177  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3178  template<int, int> class DH>
3179  bool
3180  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3181  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasPoints() const
3182  {
3183  if (this->GetType().find("constraint") != std::string::npos
3184  || (this->GetType() == "functional")
3185  || this->GetType() == "aux_functional" || (this->GetType() == "state")
3186  || (this->GetType() == "tangent"))
3187  {
3188  // We dont need PointRhs in this cases.
3189  return false;
3190  }
3191  else if ((this->GetType() == "adjoint")
3192  || (this->GetType() == "adjoint_hessian")
3193  || (this->GetType() == "hessian"))
3194  {
3195  return this->GetFunctional()->HasPoints();
3196  }
3197  else if (this->GetType() == "adjoint_for_ee")
3198  {
3199  return _aux_functionals[_functional_for_ee_num]->HasPoints();
3200  }
3201  else if (this->GetType() == "gradient")
3202  {
3203  return this->GetFunctional()->HasPoints();
3204  }
3205  else
3206  {
3207  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3208  "OptProblem::HasPoints");
3209  }
3210  }
3211 // }
3212 
3213  /******************************************************/
3214 
3215  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3216  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3217  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3218  template<int, int> class DH>
3219  bool
3220  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3221  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasInterfaces() const
3222  {
3223  if (this->GetType().find("aux_functional") != std::string::npos)
3224  {
3225  return false;
3226  }
3227  else if (this->GetType().find("functional") != std::string::npos)
3228  {
3229  return false;
3230  }
3231  else if (this->GetType().find("constraint") != std::string::npos)
3232  {
3233  return false;
3234  }
3235  else
3236  {
3237  if ((this->GetType() == "state") || this->GetType() == "adjoint_for_ee"
3238  || (this->GetType() == "tangent") || (this->GetType() == "gradient")
3239  || (this->GetType() == "adjoint")
3240  || (this->GetType() == "adjoint_hessian")
3241  || (this->GetType() == "hessian"))
3242  {
3243  return this->GetPDE().HasInterfaces();
3244  }
3245  else
3246  {
3247  throw DOpEException("Unknown Type: '" + this->GetType() + "'!",
3248  "OptProblemContainer::HasFaces");
3249  }
3250  }
3251  }
3252 
3253  /******************************************************/
3254 
3255  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3256  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3257  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3258  template<int, int> class DH>
3259  void
3260  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3261  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetControlDirichletBoundaryColors(
3262  unsigned int color, const std::vector<bool>& comp_mask,
3263  const DOpEWrapper::Function<dealdim>* values)
3264  {
3265  assert(values);
3266 
3267  unsigned int comp = _control_dirichlet_colors.size();
3268  for (unsigned int i = 0; i < _control_dirichlet_colors.size(); ++i)
3269  {
3270  if (_control_dirichlet_colors[i] == color)
3271  {
3272  comp = i;
3273  break;
3274  }
3275  }
3276  if (comp != _control_dirichlet_colors.size())
3277  {
3278  std::stringstream s;
3279  s << "ControlDirichletColor" << color << " has multiple occurences !";
3280  throw DOpEException(s.str(),
3281  "OptProblemContainer::SetControlDirichletBoundary");
3282  }
3283  _control_dirichlet_colors.push_back(color);
3284  _control_dirichlet_comps.push_back(comp_mask);
3285  _control_dirichlet_values.push_back(values);
3286  }
3287 
3288  /******************************************************/
3289 
3290  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3291  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3292  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3293  template<int, int> class DH>
3294  void
3295  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3296  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetDirichletBoundaryColors(
3297  unsigned int color, const std::vector<bool>& comp_mask,
3298  const DD* values)
3299  {
3300  assert(values);
3301  assert(values->n_components() == comp_mask.size());
3302 
3303  unsigned int comp = _dirichlet_colors.size();
3304  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
3305  {
3306  if (_dirichlet_colors[i] == color)
3307  {
3308  comp = i;
3309  break;
3310  }
3311  }
3312  if (comp != _dirichlet_colors.size())
3313  {
3314  std::stringstream s;
3315  s << "DirichletColor" << color << " has multiple occurrences !";
3316  throw DOpEException(s.str(),
3317  "OptProblemContainer::SetDirichletBoundary");
3318  }
3319  _dirichlet_colors.push_back(color);
3320  _dirichlet_comps.push_back(comp_mask);
3321  PrimalDirichletData<DD, VECTOR, dealdim>* data =
3322  new PrimalDirichletData<DD, VECTOR, dealdim>(*values);
3323  _primal_dirichlet_values.push_back(data);
3324  TangentDirichletData<DD, VECTOR, dealdim>* tdata =
3325  new TangentDirichletData<DD, VECTOR, dealdim>(*values);
3326  _tangent_dirichlet_values.push_back(tdata);
3327 
3328  if (values->NeedsControl())
3329  {
3330  _control_transposed_dirichlet_colors.push_back(color);
3331  TransposedGradientDirichletData<DD, VECTOR, dealdim> * gdata =
3332  new TransposedGradientDirichletData<DD, VECTOR, dealdim>(
3333  *values);
3334  _transposed_control_gradient_dirichlet_values.push_back(gdata);
3335  TransposedHessianDirichletData<DD, VECTOR, dealdim> * hdata =
3336  new TransposedHessianDirichletData<DD, VECTOR, dealdim>(
3337  *values);
3338  _transposed_control_hessian_dirichlet_values.push_back(hdata);
3339  }
3340  }
3341 
3342  /******************************************************/
3343 
3344  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3345  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3346  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3347  template<int, int> class DH>
3348  const std::vector<unsigned int>&
3349  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3350  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletColors() const
3351  {
3352  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3353  || this->GetType() == "adjoint_for_ee"
3354  || (this->GetType() == "tangent")
3355  || (this->GetType() == "adjoint_hessian"))
3356  {
3357  return _dirichlet_colors;
3358  }
3359  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
3360  || (this->GetType() == "global_constraint_gradient"))
3361  {
3362  return _control_dirichlet_colors;
3363  }
3364  else
3365  {
3366  throw DOpEException("Unknown Type:" + this->GetType(),
3367  "OptProblemContainer::GetDirichletColors");
3368  }
3369  }
3370  /******************************************************/
3371 
3372  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3373  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3374  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3375  template<int, int> class DH>
3376  const std::vector<unsigned int>&
3377  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3378  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletColors() const
3379  {
3380  if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
3381  {
3382  return _control_transposed_dirichlet_colors;
3383  }
3384  else
3385  {
3386  throw DOpEException("Unknown Type:" + this->GetType(),
3387  "OptProblemContainer::GetTransposedDirichletColors");
3388  }
3389  }
3390 
3391  /******************************************************/
3392 
3393  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3394  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3395  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3396  template<int, int> class DH>
3397  const std::vector<bool>&
3398  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3399  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletCompMask(
3400  unsigned int color) const
3401  {
3402  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3403  || this->GetType() == "adjoint_for_ee"
3404  || (this->GetType() == "tangent")
3405  || (this->GetType() == "adjoint_hessian"))
3406  {
3407  unsigned int comp = _dirichlet_colors.size();
3408  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
3409  {
3410  if (_dirichlet_colors[i] == color)
3411  {
3412  comp = i;
3413  break;
3414  }
3415  }
3416  if (comp == _dirichlet_colors.size())
3417  {
3418  std::stringstream s;
3419  s << "DirichletColor" << color << " has not been found !";
3420  throw DOpEException(s.str(),
3421  "OptProblemContainer::GetDirichletCompMask");
3422  }
3423  return _dirichlet_comps[comp];
3424  }
3425  else if ((this->GetType() == "gradient")
3426  || (this->GetType() == "hessian"))
3427  {
3428  unsigned int comp = _control_dirichlet_colors.size();
3429  for (unsigned int i = 0; i < _control_dirichlet_colors.size(); ++i)
3430  {
3431  if (_control_dirichlet_colors[i] == color)
3432  {
3433  comp = i;
3434  break;
3435  }
3436  }
3437  if (comp == _control_dirichlet_colors.size())
3438  {
3439  std::stringstream s;
3440  s << "ControlDirichletColor" << color << " has not been found !";
3441  throw DOpEException(s.str(),
3442  "OptProblemContainer::GetDirichletCompMask");
3443  }
3444  return _control_dirichlet_comps[comp];
3445  }
3446  else
3447  {
3448  throw DOpEException("Unknown Type:" + this->GetType(),
3449  "OptProblemContainer::GetDirichletCompMask");
3450  }
3451  }
3452  /******************************************************/
3453 
3454  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3455  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3456  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3457  template<int, int> class DH>
3458  const std::vector<bool>&
3459  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3460  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletCompMask(
3461  unsigned int color) const
3462  {
3463  if ((this->GetType() == "gradient") || (this->GetType() == "hessian"))
3464  {
3465  unsigned int comp = _dirichlet_colors.size();
3466  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
3467  {
3468  if (_dirichlet_colors[i] == color)
3469  {
3470  comp = i;
3471  break;
3472  }
3473  }
3474  if (comp == _dirichlet_colors.size())
3475  {
3476  std::stringstream s;
3477  s << "DirichletColor" << color << " has not been found !";
3478  throw DOpEException(s.str(),
3479  "OptProblemContainer::GetDirichletCompMask");
3480  }
3481  return _dirichlet_comps[comp];
3482  }
3483  else
3484  {
3485  throw DOpEException("Unknown Type:" + this->GetType(),
3486  "OptProblemContainer::GetTransposedDirichletCompMask");
3487  }
3488  }
3489 
3490  /******************************************************/
3491 
3492  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3493  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3494  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3495  template<int, int> class DH>
3496  const Function<dealdim>&
3497  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3498  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDirichletValues(
3499  unsigned int color,
3500  const std::map<std::string, const dealii::Vector<double>*> &param_values,
3501  const std::map<std::string, const VECTOR*> &domain_values) const
3502  {
3503 
3504  unsigned int col = _dirichlet_colors.size();
3505  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3506  || this->GetType() == "adjoint_for_ee"
3507  || (this->GetType() == "tangent")
3508  || (this->GetType() == "adjoint_hessian"))
3509  {
3510  for (unsigned int i = 0; i < _dirichlet_colors.size(); ++i)
3511  {
3512  if (_dirichlet_colors[i] == color)
3513  {
3514  col = i;
3515  break;
3516  }
3517  }
3518  if (col == _dirichlet_colors.size())
3519  {
3520  std::stringstream s;
3521  s << "DirichletColor" << color << " has not been found !";
3522  throw DOpEException(s.str(),
3523  "OptProblemContainer::GetDirichletValues");
3524  }
3525  }
3526  else if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
3527  {
3528  col = _control_dirichlet_colors.size();
3529  for (unsigned int i = 0; i < _control_dirichlet_colors.size(); ++i)
3530  {
3531  if (_control_dirichlet_colors[i] == color)
3532  {
3533  col = i;
3534  break;
3535  }
3536  }
3537  if (col == _control_dirichlet_colors.size())
3538  {
3539  std::stringstream s;
3540  s << "ControlDirichletColor" << color << " has not been found !";
3541  throw DOpEException(s.str(),
3542  "OptProblemContainer::GetDirichletValues");
3543  }
3544  }
3545  else
3546  {
3547  throw DOpEException("Unknown Type:" + this->GetType(),
3548  "OptProblemContainer::GetDirichletValues");
3549  }
3550 
3551  if (this->GetType() == "state")
3552  {
3553  _primal_dirichlet_values[col]->ReInit(param_values, domain_values,
3554  color);
3555  return *(_primal_dirichlet_values[col]);
3556  }
3557  else if (this->GetType() == "tangent")
3558  {
3559  _tangent_dirichlet_values[col]->ReInit(param_values, domain_values,
3560  color);
3561  return *(_tangent_dirichlet_values[col]);
3562  }
3563  else if (this->GetType() == "adjoint"
3564  || this->GetType() == "adjoint_for_ee"
3565  || (this->GetType() == "adjoint_hessian"))
3566  {
3567  return *(_zero_dirichlet_values);
3568  }
3569  else if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
3570  {
3571  return *(_control_dirichlet_values[col]);
3572  }
3573  else
3574  {
3575  throw DOpEException("Unknown Type:" + this->GetType(),
3576  "OptProblemContainer::GetDirichletValues");
3577  }
3578  }
3579  /******************************************************/
3580 
3581  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3582  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3583  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3584  template<int, int> class DH>
3585  const TransposedDirichletDataInterface<dealdim>&
3586  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3587  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetTransposedDirichletValues(
3588  unsigned int color,
3589  const std::map<std::string, const dealii::Vector<double>*> &param_values,
3590  const std::map<std::string, const VECTOR*> &domain_values) const
3591  {
3592  unsigned int col = _control_transposed_dirichlet_colors.size();
3593  if (this->GetType() == "gradient" || (this->GetType() == "hessian"))
3594  {
3595  for (unsigned int i = 0;
3596  i < _control_transposed_dirichlet_colors.size(); ++i)
3597  {
3598  if (_control_transposed_dirichlet_colors[i] == color)
3599  {
3600  col = i;
3601  break;
3602  }
3603  }
3604  if (col == _control_transposed_dirichlet_colors.size())
3605  {
3606  std::stringstream s;
3607  s << "TransposedControlDirichletColor" << color
3608  << " has not been found !";
3609  throw DOpEException(s.str(),
3610  "OptProblemContainer::GetTransposedDirichletValues");
3611  }
3612  }
3613  else
3614  {
3615  throw DOpEException("Unknown Type:" + this->GetType(),
3616  "OptProblemContainer::GetTransposedDirichletValues");
3617  }
3618 
3619  if (this->GetType() == "gradient")
3620  {
3621  _transposed_control_gradient_dirichlet_values[col]->ReInit(param_values,
3622  domain_values, color);
3623  return *(_transposed_control_gradient_dirichlet_values[col]);
3624  }
3625  else if (this->GetType() == "hessian")
3626  {
3627  _transposed_control_hessian_dirichlet_values[col]->ReInit(param_values,
3628  domain_values, color);
3629  return *(_transposed_control_hessian_dirichlet_values[col]);
3630  }
3631  else
3632  {
3633  throw DOpEException("Unknown Type:" + this->GetType(),
3634  "OptProblemContainer::GetTransposedDirichletValues");
3635  }
3636  }
3637 
3638  /******************************************************/
3639 
3640  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3641  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3642  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3643  template<int, int> class DH>
3644  const std::vector<unsigned int>&
3645  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3646  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetBoundaryEquationColors() const
3647  {
3648  if (this->GetType() == "state" || this->GetType() == "tangent")
3649  {
3650  return _state_boundary_equation_colors;
3651  }
3652  else if (this->GetType() == "adjoint"
3653  || this->GetType() == "adjoint_for_ee"
3654  || this->GetType() == "adjoint_hessian")
3655  {
3656  return _adjoint_boundary_equation_colors;
3657  }
3658  else if (this->GetType() == "gradient" || (this->GetType() == "hessian")
3659  || (this->GetType() == "global_constraint_gradient"))
3660  {
3661  return _control_boundary_equation_colors;
3662  }
3663  else
3664  {
3665  throw DOpEException("Unknown Type:" + this->GetType(),
3666  "OptProblemContainer::GetBoundaryEquationColors");
3667  }
3668  }
3669 
3670  /******************************************************/
3671 
3672  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3673  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3674  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3675  template<int, int> class DH>
3676  void
3677  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3678  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetControlBoundaryEquationColors(
3679  unsigned int color)
3680  {
3681  { //Control Boundary Equation colors are simply inserted
3682  unsigned int comp = _control_boundary_equation_colors.size();
3683  for (unsigned int i = 0; i < _control_boundary_equation_colors.size();
3684  ++i)
3685  {
3686  if (_control_boundary_equation_colors[i] == color)
3687  {
3688  comp = i;
3689  break;
3690  }
3691  }
3692  if (comp != _control_boundary_equation_colors.size())
3693  {
3694  std::stringstream s;
3695  s << "Boundary Equation Color" << color
3696  << " has multiple occurences !";
3697  throw DOpEException(s.str(),
3698  "OptProblemContainer::SetControlBoundaryEquationColors");
3699  }
3700  _control_boundary_equation_colors.push_back(color);
3701  }
3702  }
3703  /******************************************************/
3704 
3705  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3706  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3707  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3708  template<int, int> class DH>
3709  void
3710  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3711  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetBoundaryEquationColors(
3712  unsigned int color)
3713  {
3714  { //State Boundary Equation colors are simply inserted
3715  unsigned int comp = _state_boundary_equation_colors.size();
3716  for (unsigned int i = 0; i < _state_boundary_equation_colors.size();
3717  ++i)
3718  {
3719  if (_state_boundary_equation_colors[i] == color)
3720  {
3721  comp = i;
3722  break;
3723  }
3724  }
3725  if (comp != _state_boundary_equation_colors.size())
3726  {
3727  std::stringstream s;
3728  s << "Boundary Equation Color" << color
3729  << " has multiple occurences !";
3730  throw DOpEException(s.str(),
3731  "OptProblemContainer::SetBoundaryEquationColors");
3732  }
3733  _state_boundary_equation_colors.push_back(color);
3734  }
3735  { //For the adjoint they are added with the boundary functional colors
3736  unsigned int comp = _adjoint_boundary_equation_colors.size();
3737  for (unsigned int i = 0; i < _adjoint_boundary_equation_colors.size();
3738  ++i)
3739  {
3740  if (_adjoint_boundary_equation_colors[i] == color)
3741  {
3742  comp = i;
3743  break;
3744  }
3745  }
3746  if (comp != _adjoint_boundary_equation_colors.size())
3747  {
3748  //Seems this color is already added, however it might have been a functional color
3749  //so we don't do anything.
3750  }
3751  else
3752  {
3753  _adjoint_boundary_equation_colors.push_back(color);
3754  }
3755  }
3756  }
3757 
3758  /******************************************************/
3759 
3760  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3761  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3762  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3763  template<int, int> class DH>
3764  const std::vector<unsigned int>&
3765  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3766  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetBoundaryFunctionalColors() const
3767  {
3768  if (this->GetType() == "cost_functional"
3769  || this->GetType() == "aux_functional"
3770  || this->GetType() == "functional_for_ee") //fixme: what about error_evaluation?
3771  {
3772  return _boundary_functional_colors;
3773  }
3774  else
3775  {
3776  throw DOpEException("Unknown Type:" + this->GetType(),
3777  "OptProblemContainer::GetBoundaryFunctionalColors");
3778  }
3779  }
3780 
3781  /******************************************************/
3782 
3783  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3784  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3785  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3786  template<int, int> class DH>
3787  void
3788  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3789  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::SetBoundaryFunctionalColors(
3790  unsigned int color)
3791  {
3792  { //Boundary Functional colors are simply inserted
3793  unsigned int comp = _boundary_functional_colors.size();
3794  for (unsigned int i = 0; i < _boundary_functional_colors.size(); ++i)
3795  {
3796  if (_boundary_functional_colors[i] == color)
3797  {
3798  comp = i;
3799  break;
3800  }
3801  }
3802  if (comp != _boundary_functional_colors.size())
3803  {
3804  std::stringstream s;
3805  s << "Boundary Functional Color" << color
3806  << " has multiple occurences !";
3807  throw DOpEException(s.str(),
3808  "OptProblemContainer::SetBoundaryFunctionalColors");
3809  }
3810  _boundary_functional_colors.push_back(color);
3811  }
3812  { //For the adjoint they are addeed to the boundary equation colors
3813  unsigned int comp = _adjoint_boundary_equation_colors.size();
3814  for (unsigned int i = 0; i < _adjoint_boundary_equation_colors.size();
3815  ++i)
3816  {
3817  if (_adjoint_boundary_equation_colors[i] == color)
3818  {
3819  comp = i;
3820  break;
3821  }
3822  }
3823  if (comp != _adjoint_boundary_equation_colors.size())
3824  {
3825  //Seems this color is already added, however it might have been a equation color
3826  //so we don't do anything.
3827  }
3828  else
3829  {
3830  _adjoint_boundary_equation_colors.push_back(color);
3831  }
3832  }
3833  }
3834 
3835  /******************************************************/
3836 
3837  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3838  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3839  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3840  template<int, int> class DH>
3841  unsigned int
3842  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3843  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetControlNBlocks() const
3844  {
3845  return this->GetPDE().GetControlNBlocks();
3846  }
3847 
3848  /******************************************************/
3849 
3850  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3851  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3852  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3853  template<int, int> class DH>
3854  unsigned int
3855  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3856  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetStateNBlocks() const
3857  {
3858  return this->GetPDE().GetStateNBlocks();
3859  }
3860 
3861  /******************************************************/
3862 
3863  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3864  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3865  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3866  template<int, int> class DH>
3867  unsigned int
3868  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3869  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetNBlocks() const
3870  {
3871  if ((this->GetType() == "state") || (this->GetType() == "adjoint_for_ee")
3872  || (this->GetType() == "adjoint") || (this->GetType() == "tangent")
3873  || (this->GetType() == "adjoint_hessian"))
3874  {
3875  return this->GetStateNBlocks();
3876  }
3877  else if ((this->GetType() == "gradient")
3878  || (this->GetType() == "hessian"))
3879  {
3880  return this->GetControlNBlocks();
3881  }
3882  else
3883  {
3884  throw DOpEException("Unknown Type:" + this->GetType(),
3885  "OptProblemContainer::GetNBlocks");
3886  }
3887  }
3888 
3889  /******************************************************/
3890 
3891  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3892  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3893  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3894  template<int, int> class DH>
3895  unsigned int
3896  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3897  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFsPerBlock(
3898  unsigned int b) 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 GetSpaceTimeHandler()->GetStateDoFsPerBlock(b);
3906  }
3907  else if ((this->GetType() == "gradient")
3908  || (this->GetType() == "hessian"))
3909  {
3910  return GetSpaceTimeHandler()->GetControlDoFsPerBlock(b);
3911  }
3912  else
3913  {
3914  throw DOpEException("Unknown Type:" + this->GetType(),
3915  "OptProblemContainer::GetDoFsPerBlock");
3916  }
3917  }
3918 
3919  /******************************************************/
3920 
3921  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3922  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3923  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3924  template<int, int> class DH>
3925  const std::vector<unsigned int>&
3926  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3927  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFsPerBlock() const
3928  {
3929  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3930  || (this->GetType() == "adjoint_for_ee")
3931  || (this->GetType() == "tangent")
3932  || (this->GetType() == "adjoint_hessian"))
3933  {
3934  return GetSpaceTimeHandler()->GetStateDoFsPerBlock();
3935  }
3936  else if ((this->GetType() == "gradient")
3937  || (this->GetType() == "hessian"))
3938  {
3939  return GetSpaceTimeHandler()->GetControlDoFsPerBlock();
3940  }
3941  else
3942  {
3943  throw DOpEException("Unknown Type:" + this->GetType(),
3944  "OptProblemContainer::GetDoFsPerBlock");
3945  }
3946  }
3947 
3948  /******************************************************/
3949 
3950  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3951  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3952  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3953  template<int, int> class DH>
3954  const dealii::ConstraintMatrix&
3955  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3956  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::GetDoFConstraints() const
3957  {
3958  if ((this->GetType() == "state") || (this->GetType() == "adjoint")
3959  || (this->GetType() == "adjoint_for_ee")
3960  || (this->GetType() == "tangent")
3961  || (this->GetType() == "adjoint_hessian"))
3962  {
3963  return GetSpaceTimeHandler()->GetStateDoFConstraints();
3964  }
3965  else if ((this->GetType() == "gradient") || (this->GetType() == "hessian")
3966  || (this->GetType() == "global_constraint_gradient"))
3967  {
3968  return GetSpaceTimeHandler()->GetControlDoFConstraints();
3969  }
3970  else
3971  {
3972  throw DOpEException("Unknown Type:" + this->GetType(),
3973  "OptProblemContainer::GetDoFConstraints");
3974  }
3975  }
3976 
3977  /******************************************************/
3978 
3979  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
3980  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
3981  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
3982  template<int, int> class DH>
3983  bool
3984  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
3985  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::NeedTimeFunctional() const
3986  {
3987  if (this->GetType() == "cost_functional")
3988  return GetFunctional()->NeedTime();
3989  else if (this->GetType() == "aux_functional")
3990  return _aux_functionals[this->GetTypeNum()]->NeedTime();
3991  else if (this->GetType() == "functional_for_ee")
3992  return _aux_functionals[_functional_for_ee_num]->NeedTime();
3993  else
3994  throw DOpEException("Not implemented",
3995  "OptProblemContainer::NeedTimeFunctional");
3996  }
3997 
3998  /******************************************************/
3999 
4000  template<typename FUNCTIONAL_INTERFACE, typename FUNCTIONAL, typename PDE,
4001  typename DD, typename CONSTRAINTS, typename SPARSITYPATTERN,
4002  typename VECTOR, int dopedim, int dealdim, template<int, int> class FE,
4003  template<int, int> class DH>
4004  bool
4005  OptProblemContainer<FUNCTIONAL_INTERFACE, FUNCTIONAL, PDE, DD, CONSTRAINTS,
4006  SPARSITYPATTERN, VECTOR, dopedim, dealdim, FE, DH>::HasControlInDirichletData() const
4007  {
4008  return (!_control_transposed_dirichlet_colors.empty());
4009  }
4010 
4011 }
4012 #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:2353
DOpE::OptProblemContainer::DeleteAuxiliaryState
void DeleteAuxiliaryState(std::string name)
Definition: optproblemcontainer.h:3048
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:1115
DOpE::OptProblemContainer::OptProblemContainer
OptProblemContainer(FUNCTIONAL &functional, PDE &pde, CONSTRAINTS &constraints, SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > &STH)
Definition: optproblemcontainer.h:1294
DOpE::OptProblemContainer::BoundaryEquation
void BoundaryEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1936
DOpE::OptProblemContainer::GetUpdateFlags
dealii::UpdateFlags GetUpdateFlags() const
Definition: optproblemcontainer.h:2670
DOpE::StateProblem
Definition: stateproblem.h:48
DOpE::OptProblemContainer::GetAuxiliaryConstraint
const ConstraintVector< VECTOR > * GetAuxiliaryConstraint(std::string name)
Definition: optproblemcontainer.h:910
DOpE::OptProblemContainer::SetTime
void SetTime(double time, const TimeIterator &interval)
Definition: optproblemcontainer.h:2814
DOpE::OptProblemContainer::GetNFunctionals
unsigned int GetNFunctionals() const
Definition: optproblemcontainer.h:754
DOpE::OptProblemContainer::FaceFunctional
double FaceFunctional(const FACEDATACONTAINER &fdc)
Definition: optproblemcontainer.h:1615
DOpE::OptProblemContainer::Init_ElementEquation
void Init_ElementEquation(const DATACONTAINER &edc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1050
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:2263
DOpE::ConstraintVector
Definition: constraintvector.h:47
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:1648
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:2578
DOpE::OptProblemContainer::ElementTimeMatrixExplicit
void ElementTimeMatrixExplicit(const DATACONTAINER &dc, dealii::FullMatrix< double > &local_entry_matrix)
Definition: optproblemcontainer.h:2433
DOpE::OptProblemContainer::ElementTimeEquationExplicit
void ElementTimeEquationExplicit(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:1806
DOpE::OptProblemContainer::GetBoundaryFunctionalColors
const std::vector< unsigned int > & GetBoundaryFunctionalColors() const
Definition: optproblemcontainer.h:3766
DOpE::OptProblemContainer::ElementRhs
void ElementRhs(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:1987
DOpE::OptProblemContainer::ReInit
void ReInit(std::string algo_type)
Definition: optproblemcontainer.h:1360
DOpE::OptProblemContainer::GetBlockVector
const VECTOR * GetBlockVector(const std::map< std::string, const VECTOR * > &values, std::string name)
Definition: optproblemcontainer.h:1211
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:2635
DOpE::OptProblemContainer::PostProcessConstraints
void PostProcessConstraints(ConstraintVector< VECTOR > &g) const
Definition: optproblemcontainer.h:2910
DOpE::OptProblemContainer::GetFaceUpdateFlags
dealii::UpdateFlags GetFaceUpdateFlags() const
Definition: optproblemcontainer.h:2715
DOpE::OptProblemContainer::GetFunctionalType
std::string GetFunctionalType() const
Definition: optproblemcontainer.h:2759
DOpE::OptProblemContainer::AddFunctional
void AddFunctional(FUNCTIONAL_INTERFACE *F)
Definition: optproblemcontainer.h:692
DOpE::OptProblemContainer::GetFunctional
FUNCTIONAL * GetFunctional()
Definition: optproblemcontainer.h:3092
DOpE::TangentDirichletData
Definition: tangentdirichletdata.h:39
DOpE::OptProblemContainer::HasInterfaces
bool HasInterfaces() const
Definition: optproblemcontainer.h:3221
DOpE::OptProblemContainer::BoundaryFunctional
double BoundaryFunctional(const FACEDATACONTAINER &fdc)
Definition: optproblemcontainer.h:1577
DOpE::OptProblemContainer::ElementFunctional
double ElementFunctional(const DATACONTAINER &edc)
Definition: optproblemcontainer.h:1487
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:3646
DOpE::OptProblemContainer::SetDirichletBoundaryColors
void SetDirichletBoundaryColors(unsigned int color, const std::vector< bool > &comp_mask, const DD *values)
Definition: optproblemcontainer.h:3296
DOpE::OptProblemContainer::GetDoFsPerBlock
const std::vector< unsigned int > & GetDoFsPerBlock() const
Definition: optproblemcontainer.h:3927
DOpE::OptProblemContainer::SetControlBoundaryEquationColors
void SetControlBoundaryEquationColors(unsigned int color)
Definition: optproblemcontainer.h:3678
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:2093
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:2923
DOpE::OptProblemContainer::GetDirichletColors
const std::vector< unsigned int > & GetDirichletColors() const
Definition: optproblemcontainer.h:3350
DOpE::OptProblemContainer::HasPoints
bool HasPoints() const
Definition: optproblemcontainer.h:3181
DOpE::OptProblemContainer::GetDoFType
std::string GetDoFType() const
Definition: optproblemcontainer.h:2607
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:1685
DOpE::OptProblemContainer::ElementTimeEquation
void ElementTimeEquation(const DATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale=1.)
Definition: optproblemcontainer.h:1762
DOpE::OptProblemContainer::GetDirichletCompMask
const std::vector< bool > & GetDirichletCompMask(unsigned int color) const
Definition: optproblemcontainer.h:3399
DOpE::OptProblemContainer::SetType
void SetType(std::string type, unsigned int num=0)
Definition: optproblemcontainer.h:1402
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:1523
DOpE::TimeIterator
Definition: timeiterator.h:63
DOpE::OptProblemContainer::GetNBlocks
unsigned int GetNBlocks() const
Definition: optproblemcontainer.h:3869
DOpE::OptProblemContainer::GetAuxiliaryState
const StateVector< VECTOR > * GetAuxiliaryState(std::string name) const
Definition: optproblemcontainer.h:3005
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:3498
DOpE::ControlVector
Definition: controlvector.h:48
DOpE::OptProblemContainer::Init_ElementMatrix
void Init_ElementMatrix(const DATACONTAINER &edc, dealii::FullMatrix< double > &local_entry_matrix, double scale, double scale_ico)
Definition: optproblemcontainer.h:1165
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:1189
DOpE::OptProblemContainer::GetDoFConstraints
const dealii::ConstraintMatrix & GetDoFConstraints() const
Definition: optproblemcontainer.h:3956
DOpE::OptProblemContainer::GetConstraintType
std::string GetConstraintType() const
Definition: optproblemcontainer.h:2801
DOpE::OptProblemContainer
Definition: optproblemcontainer.h:101
problemcontainer_internal.h
DOpE::OptProblemContainer::GetStateNBlocks
unsigned int GetStateNBlocks()
Definition: optproblemcontainer.h:1025
DOpE::OptProblemContainer::FaceEquation
void FaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1854
DOpE::OptProblemContainer::AddAuxiliaryConstraint
void AddAuxiliaryConstraint(const ConstraintVector< VECTOR > *c, std::string name)
Definition: optproblemcontainer.h:2964
DOpE::OptProblemContainer::GetAuxiliaryControl
const ControlVector< VECTOR > * GetAuxiliaryControl(std::string name) const
Definition: optproblemcontainer.h:2984
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:2175
DOpE::OptProblemContainer::GetSpaceTimeHandler
const SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > * GetSpaceTimeHandler() const
Definition: optproblemcontainer.h:843
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:1200
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:975
DOpE::OptProblemContainer::GetVector
const dealii::Vector< double > * GetVector(const std::map< std::string, const Vector< double > * > &values, std::string name)
Definition: optproblemcontainer.h:1224
DOpE::OptProblemContainer::GetSpaceTimeHandler
SpaceTimeHandler< FE, DH, SPARSITYPATTERN, VECTOR, dopedim, dealdim > * GetSpaceTimeHandler()
Definition: optproblemcontainer.h:851
DOpE::OptProblemContainer::GetConstraints
const CONSTRAINTS * GetConstraints() const
Definition: optproblemcontainer.h:1205
DOpE::OptProblemContainer::HasFaces
bool HasFaces() const
Definition: optproblemcontainer.h:3133
DOpE::PrimalDirichletData
Definition: primaldirichletdata.h:41
DOpE::OptProblemContainer::GetTransposedDirichletColors
const std::vector< unsigned int > & GetTransposedDirichletColors() const
Definition: optproblemcontainer.h:3378
DOpE::StateVector
Definition: statevector.h:49
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:3587
DOpE::OptProblemContainer::ComputeSparsityPattern
void ComputeSparsityPattern(SPARSITYPATTERN &sparsity) const
Definition: optproblemcontainer.h:2876
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:1735
DOpE::OptProblemContainer::InterfaceEquation
void InterfaceEquation(const FACEDATACONTAINER &dc, dealii::Vector< double > &local_vector, double scale, double scale_ico)
Definition: optproblemcontainer.h:1904
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:3856
DOpE::OptProblemContainer::GetControlBlockComponent
std::vector< unsigned int > & GetControlBlockComponent()
Definition: optproblemcontainer.h:1033
DOpE::OptProblemContainer::AddAuxiliaryState
void AddAuxiliaryState(const StateVector< VECTOR > *c, std::string name)
Definition: optproblemcontainer.h:2944
DOpE::OptProblemContainer::NeedTimeFunctional
bool NeedTimeFunctional() const
Definition: optproblemcontainer.h:3985
dofhandler_wrapper.h
DOpE::OptProblemContainer::GetFunctionalName
std::string GetFunctionalName() const
Definition: optproblemcontainer.h:2780
function_wrapper.h
DOpE::TransposedDirichletDataInterface
Definition: transposeddirichletdatainterface.h:37
DOpE::OptProblemContainer::GetControlNBlocks
unsigned int GetControlNBlocks() const
Definition: optproblemcontainer.h:3843
DOpE::OptProblemContainer::GetTransposedDirichletCompMask
const std::vector< bool > & GetTransposedDirichletCompMask(unsigned int color) const
Definition: optproblemcontainer.h:3460
DOpE::OptProblemContainer::~OptProblemContainer
~OptProblemContainer()
Definition: optproblemcontainer.h:1317
DOpE::OptProblemContainer::HasControlInDirichletData
bool HasControlInDirichletData() const
Definition: optproblemcontainer.h:4006
DOpE::OptProblemContainer::GetStateBlockComponent
std::vector< unsigned int > & GetStateBlockComponent()
Definition: optproblemcontainer.h:1040
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:925
DOpE::OptProblemContainer::SetBoundaryEquationColors
void SetBoundaryEquationColors(unsigned int color)
Definition: optproblemcontainer.h:3711
DOpE::OptProblemContainer::SetControlDirichletBoundaryColors
void SetControlDirichletBoundaryColors(unsigned int color, const std::vector< bool > &comp_mask, const DOpEWrapper::Function< dealdim > *values)
Definition: optproblemcontainer.h:3261
DOpE::OptProblemContainer::DeleteAuxiliaryConstraint
void DeleteAuxiliaryConstraint(std::string name)
Definition: optproblemcontainer.h:3070
DOpE::OptProblemContainer::GetFunctionalPosition
const std::map< std::string, unsigned int > & GetFunctionalPosition() const
Definition: optproblemcontainer.h:1019
DOpE::OptProblemContainer::SetBoundaryFunctionalColors
void SetBoundaryFunctionalColors(unsigned int color)
Definition: optproblemcontainer.h:3789
DOpE::OptProblemContainer::DeleteAuxiliaryControl
void DeleteAuxiliaryControl(std::string name)
Definition: optproblemcontainer.h:3026
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:1069
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 Fri Oct 24 2014 17:56:39 for DOpE by   doxygen 1.8.6