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