DOpE
statreducedproblem.h
Go to the documentation of this file.
1 
24 #ifndef _STAT_REDUCED_PROBLEM_H_
25 #define _STAT_REDUCED_PROBLEM_H_
26 
27 #include "reducedprobleminterface.h"
28 #include "integrator.h"
29 #include "parameterreader.h"
30 #include "statevector.h"
31 #include "stateproblem.h"
32 
33 #include <lac/vector.h>
34 
35 #include <lac/block_sparsity_pattern.h>
36 #include <lac/block_sparse_matrix.h>
37 
38 #include "optproblemcontainer.h"
39 #include "pdeinterface.h"
40 #include "functionalinterface.h"
41 #include "dirichletdatainterface.h"
42 #include "dopeexception.h"
43 #include "newtonsolver.h"
44 #include "newtonsolvermixeddims.h"
45 #include "cglinearsolver.h"
46 #include "gmreslinearsolver.h"
47 #include "directlinearsolver.h"
48 #include "voidlinearsolver.h"
49 #include "constraintinterface.h"
50 #include "solutionextractor.h"
51 
52 #include <base/data_out_base.h>
53 #include <numerics/data_out.h>
54 #include <numerics/matrix_tools.h>
55 #include <numerics/vector_tools.h>
56 #include <base/function.h>
57 #include <lac/sparse_matrix.h>
58 #include <lac/compressed_simple_sparsity_pattern.h>
59 #include <lac/block_sparsity_pattern.h>
60 #include <lac/sparse_direct.h>
61 
62 #include <fstream>
63 namespace DOpE
64 {
80  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
81  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
82  typename VECTOR, int dopedim, int dealdim>
83  class StatReducedProblem : public ReducedProblemInterface<PROBLEM, VECTOR>
84  {
85  public:
98  template<typename INTEGRATORDATACONT>
99  StatReducedProblem(PROBLEM *OP, std::string state_behavior,
100  ParameterReader &param_reader, INTEGRATORDATACONT& idc,
101  int base_priority = 0);
102 
116  template<typename STATEINTEGRATORDATACONT,
117  typename CONTROLINTEGRATORCONT>
118  StatReducedProblem(PROBLEM *OP, std::string state_behavior,
119  ParameterReader &param_reader, CONTROLINTEGRATORCONT& c_idc,
120  STATEINTEGRATORDATACONT & s_idc, int base_priority = 0);
121 
122  ~StatReducedProblem();
123 
124  /******************************************************/
125 
131  static void
132  declare_params(ParameterReader &param_reader);
133 
134  /******************************************************/
135 
141  void
142  ReInit();
143 
144  /******************************************************/
145 
151  bool
152  ComputeReducedConstraints(const ControlVector<VECTOR>& q,
153  ConstraintVector<VECTOR>& g);
154 
155  /******************************************************/
156 
162  void
163  GetControlBoxConstraints(ControlVector<VECTOR>& lb,
164  ControlVector<VECTOR>& ub);
165 
166  /******************************************************/
167 
173  void
174  ComputeReducedGradient(const ControlVector<VECTOR>& q,
175  ControlVector<VECTOR>& gradient,
176  ControlVector<VECTOR>& gradient_transposed);
177 
178  /******************************************************/
179 
185  double
186  ComputeReducedCostFunctional(const ControlVector<VECTOR>& q);
187 
188  /******************************************************/
189 
195  void
196  ComputeReducedFunctionals(const ControlVector<VECTOR>& q);
197 
198  /******************************************************/
199 
219  template<class DWRC,class PDE>
220  void
221  ComputeRefinementIndicators(const ControlVector<VECTOR>& q,
222  DWRC& dwrc, PDE& pde);
223 
224  /******************************************************/
225 
231  void
232  ComputeReducedHessianVector(const ControlVector<VECTOR>& q,
233  const ControlVector<VECTOR>& direction,
234  ControlVector<VECTOR>& hessian_direction,
235  ControlVector<VECTOR>& hessian_direction_transposed);
236 
237  /******************************************************/
238 
244  void
245  ComputeReducedGradientOfGlobalConstraints(unsigned int num,
246  const ControlVector<VECTOR>& q, const ConstraintVector<VECTOR>& g,
247  ControlVector<VECTOR>& gradient,
248  ControlVector<VECTOR>& gradient_transposed);
249 
250  /******************************************************/
251 
257  void
258  StateSizeInfo(std::stringstream& out)
259  {
260  GetU().PrintInfos(out);
261  }
262 
263  /******************************************************/
264 
275  void
276  WriteToFile(const VECTOR &v, std::string name, std::string outfile,
277  std::string dof_type, std::string filetype);
278 
279  /******************************************************/
280 
291  void
292  WriteToFile(const ControlVector<VECTOR> &v, std::string name,
293  std::string outfile, std::string dof_type, std::string filetype);
294 
302  void
303  WriteToFile(const std::vector<double> &/*v*/,
304  std::string /*outfile*/)
305  {
306  abort();
307  }
308 
309  protected:
317  void
318  ComputeReducedState(const ControlVector<VECTOR>& q);
319 
320  /******************************************************/
329  void
330  ComputeReducedAdjoint(const ControlVector<VECTOR>& q);
331 
332  /******************************************************/
333 
344  void
345  ComputeDualForErrorEstimation(const ControlVector<VECTOR>& q,
346  DOpEtypes::WeightComputation weight_comp);
347 
348  const StateVector<VECTOR> &
349  GetU() const
350  {
351  return _u;
352  }
353  StateVector<VECTOR> &
354  GetU()
355  {
356  return _u;
357  }
358  StateVector<VECTOR> &
359  GetZ()
360  {
361  return _z;
362  }
363  StateVector<VECTOR> &
364  GetDU()
365  {
366  return _du;
367  }
368  StateVector<VECTOR> &
369  GetDZ()
370  {
371  return _dz;
372  }
376  const StateVector<VECTOR> &
377  GetZForEE() const
378  {
379  return _z_for_ee;
380  }
381  StateVector<VECTOR> &
382  GetZForEE()
383  {
384  return _z_for_ee;
385  }
386 
387  NONLINEARSOLVER&
388  GetNonlinearSolver(std::string type);
389  CONTROLNONLINEARSOLVER&
390  GetControlNonlinearSolver();
391  INTEGRATOR&
392  GetIntegrator()
393  {
394  return _integrator;
395  }
396  CONTROLINTEGRATOR&
397  GetControlIntegrator()
398  {
399  return _control_integrator;
400  }
401 
402  private:
403  StateVector<VECTOR> _u;
404  StateVector<VECTOR> _z;
405  StateVector<VECTOR> _du;
406  StateVector<VECTOR> _dz;
407  StateVector<VECTOR> _z_for_ee;
408 
409  INTEGRATOR _integrator;
410  CONTROLINTEGRATOR _control_integrator;
411  NONLINEARSOLVER _nonlinear_state_solver;
412  NONLINEARSOLVER _nonlinear_adjoint_solver;
413  CONTROLNONLINEARSOLVER _nonlinear_gradient_solver;
414 
415  bool _build_state_matrix, _build_adjoint_matrix, _build_control_matrix;
416  bool _state_reinit, _adjoint_reinit, _gradient_reinit;
417 
418  friend class SolutionExtractor<
419  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
420  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>,
421  VECTOR> ;
422  };
423 
424  /*************************************************************************/
425  /*****************************IMPLEMENTATION******************************/
426  /*************************************************************************/
427  using namespace dealii;
428 
429  /******************************************************/
430  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
431  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
432  typename VECTOR, int dopedim, int dealdim>
433  void
434  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
435  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::declare_params(
436  ParameterReader &param_reader)
437  {
438  NONLINEARSOLVER::declare_params(param_reader);
439  }
440  /******************************************************/
441 
442  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
443  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
444  typename VECTOR, int dopedim, int dealdim>
445  template<typename INTEGRATORDATACONT>
446  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
447  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::StatReducedProblem(
448  PROBLEM *OP, std::string state_behavior,
449  ParameterReader &param_reader, INTEGRATORDATACONT& idc,
450  int base_priority)
451  : ReducedProblemInterface<PROBLEM, VECTOR>(OP,
452  base_priority), _u(OP->GetSpaceTimeHandler(), state_behavior,
453  param_reader), _z(OP->GetSpaceTimeHandler(), state_behavior,
454  param_reader), _du(OP->GetSpaceTimeHandler(), state_behavior,
455  param_reader), _dz(OP->GetSpaceTimeHandler(), state_behavior,
456  param_reader), _z_for_ee(OP->GetSpaceTimeHandler(),
457  state_behavior, param_reader), _integrator(idc), _control_integrator(
458  idc), _nonlinear_state_solver(_integrator, param_reader), _nonlinear_adjoint_solver(
459  _integrator, param_reader), _nonlinear_gradient_solver(
460  _control_integrator, param_reader)
461 
462  {
463  //ReducedProblems should be ReInited
464  {
465  _state_reinit = true;
466  _adjoint_reinit = true;
467  _gradient_reinit = true;
468  }
469  }
470 
471  /******************************************************/
472 
473  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
474  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
475  typename VECTOR, int dopedim, int dealdim>
476  template<typename STATEINTEGRATORDATACONT, typename CONTROLINTEGRATORCONT>
477  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
478  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::StatReducedProblem(
479  PROBLEM *OP, std::string state_behavior,
480  ParameterReader &param_reader, CONTROLINTEGRATORCONT& c_idc,
481  STATEINTEGRATORDATACONT & s_idc, int base_priority)
482  : ReducedProblemInterface<PROBLEM, VECTOR>(OP,
483  base_priority), _u(OP->GetSpaceTimeHandler(), state_behavior,
484  param_reader), _z(OP->GetSpaceTimeHandler(), state_behavior,
485  param_reader), _du(OP->GetSpaceTimeHandler(), state_behavior,
486  param_reader), _dz(OP->GetSpaceTimeHandler(), state_behavior,
487  param_reader), _z_for_ee(OP->GetSpaceTimeHandler(),
488  state_behavior, param_reader), _integrator(s_idc), _control_integrator(
489  c_idc), _nonlinear_state_solver(_integrator, param_reader), _nonlinear_adjoint_solver(
490  _integrator, param_reader), _nonlinear_gradient_solver(
491  _control_integrator, param_reader)
492 
493  {
494  //ReducedProblems should be ReInited
495  {
496  _state_reinit = true;
497  _adjoint_reinit = true;
498  _gradient_reinit = true;
499  }
500  }
501 
502  /******************************************************/
503 
504  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
505  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
506  typename VECTOR, int dopedim, int dealdim>
507  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
508  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::~StatReducedProblem()
509  {
510  }
511 
512  /******************************************************/
513 
514  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
515  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
516  typename VECTOR, int dopedim, int dealdim>
517  NONLINEARSOLVER&
518  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
519  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetNonlinearSolver(
520  std::string type)
521  {
522  if ((type == "state") || (type == "tangent"))
523  {
524  return _nonlinear_state_solver;
525  }
526  else if ((type == "adjoint") || (type == "adjoint_hessian")
527  || (type == "adjoint_for_ee"))
528  {
529  return _nonlinear_adjoint_solver;
530  }
531  else
532  {
533  throw DOpEException("No Solver for Problem type:`" + type + "' found",
534  "StatReducedProblem::GetNonlinearSolver");
535 
536  }
537  }
538  /******************************************************/
539 
540  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
541  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
542  typename VECTOR, int dopedim, int dealdim>
543  CONTROLNONLINEARSOLVER&
544  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
545  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetControlNonlinearSolver()
546  {
547  if ((this->GetProblem()->GetType() == "gradient")
548  || (this->GetProblem()->GetType() == "hessian"))
549  {
550  return _nonlinear_gradient_solver;
551  }
552  else
553  {
554  throw DOpEException(
555  "No Solver for Problem type:`" + this->GetProblem()->GetType()
556  + "' found", "StatReducedProblem::GetControlNonlinearSolver");
557 
558  }
559  }
560  /******************************************************/
561 
562  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
563  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
564  typename VECTOR, int dopedim, int dealdim>
565  void
566  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
567  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ReInit()
568  {
569  ReducedProblemInterface<PROBLEM, VECTOR>::ReInit();
570 
571  //Some Solvers must be reinited when called
572  // Better have subproblems, so that solver can be reinited here
573  {
574  _state_reinit = true;
575  _adjoint_reinit = true;
576  _gradient_reinit = true;
577  }
578 
579  _build_state_matrix = true;
580  _build_adjoint_matrix = true;
581 
582  GetU().ReInit();
583  GetZ().ReInit();
584  GetDU().ReInit();
585  GetDZ().ReInit();
586  GetZForEE().ReInit();
587 
588  _build_control_matrix = true;
589  }
590 
591  /******************************************************/
592 
593  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
594  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
595  typename VECTOR, int dopedim, int dealdim>
596  void
597  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
598  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedState(
599  const ControlVector<VECTOR>& q)
600  {
601  this->InitializeFunctionalValues(
602  this->GetProblem()->GetNFunctionals() + 1);
603 
604  this->GetOutputHandler()->Write("Computing State Solution:",
605  4 + this->GetBasePriority());
606 
607  this->SetProblemType("state");
608  auto& problem = this->GetProblem()->GetStateProblem();
609  if (_state_reinit == true)
610  {
611  GetNonlinearSolver("state").ReInit(problem);
612  _state_reinit = false;
613  }
614 
615  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
616  if (dopedim == dealdim)
617  {
618  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
619  }
620  else if (dopedim == 0)
621  {
622  this->GetIntegrator().AddParamData("control",
623  &(q.GetSpacialVectorCopy()));
624  }
625  else
626  {
627  throw DOpEException("dopedim not implemented",
628  "StatReducedProblem::ComputeReducedState");
629  }
630 
631  _build_state_matrix = this->GetNonlinearSolver("state").NonlinearSolve(
632  problem, (GetU().GetSpacialVector()), true, _build_state_matrix);
633 
634  if (dopedim == dealdim)
635  {
636  this->GetIntegrator().DeleteDomainData("control");
637  }
638  else if (dopedim == 0)
639  {
640  this->GetIntegrator().DeleteParamData("control");
641  q.UnLockCopy();
642  }
643  else
644  {
645  throw DOpEException("dopedim not implemented",
646  "StatReducedProblem::ComputeReducedState");
647  }
648  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
649 
650  this->GetOutputHandler()->Write((GetU().GetSpacialVector()),
651  "State" + this->GetPostIndex(), problem.GetDoFType());
652 
653  }
654  /******************************************************/
655 
656  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
657  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
658  typename VECTOR, int dopedim, int dealdim>
659  bool
660  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
661  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedConstraints(
662  const ControlVector<VECTOR>& q, ConstraintVector<VECTOR>& g)
663  {
664  this->GetOutputHandler()->Write("Evaluating Constraints:",
665  4 + this->GetBasePriority());
666 
667  this->SetProblemType("constraints");
668 
669  g = 0;
670  //Local constraints
671  // this->GetProblem()->ComputeLocalConstraints(q.GetSpacialVector(), GetU().GetSpacialVector(),
672  // g.GetSpacialVector("local"));
673  if (dopedim == dealdim)
674  {
675  this->GetControlIntegrator().AddDomainData("control",
676  &(q.GetSpacialVector()));
677  }
678  else if (dopedim == 0)
679  {
680  this->GetControlIntegrator().AddParamData("control",
681  &(q.GetSpacialVectorCopy()));
682  }
683  else
684  {
685  throw DOpEException("dopedim not implemented",
686  "StatReducedProblem::ComputeReducedConstraints");
687  }
688  this->GetControlIntegrator().ComputeLocalControlConstraints(
689  *(this->GetProblem()), g.GetSpacialVector("local"));
690  if (dopedim == dealdim)
691  {
692  this->GetControlIntegrator().DeleteDomainData("control");
693  }
694  else if (dopedim == 0)
695  {
696  this->GetControlIntegrator().DeleteParamData("control");
697  q.UnLockCopy();
698  }
699  else
700  {
701  throw DOpEException("dopedim not implemented",
702  "StatReducedProblem::ComputeReducedConstraints");
703  }
704  //Global in Space-Time Constraints
705  dealii::Vector<double>& gc = g.GetGlobalConstraints();
706  //dealii::Vector<double> global_values(gc.size());
707 
708  unsigned int nglobal = gc.size(); //global_values.size();
709 
710  if (nglobal > 0)
711  {
712  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
713 
714  for (unsigned int i = 0; i < nglobal; i++)
715  {
716  //this->SetProblemType("local_global_constraints", i);
717  this->SetProblemType("global_constraints", i);
718  this->GetIntegrator().AddDomainData("state",
719  &(GetU().GetSpacialVector()));
720  if (dopedim == dealdim)
721  {
722  this->GetIntegrator().AddDomainData("control",
723  &(q.GetSpacialVector()));
724  }
725  else if (dopedim == 0)
726  {
727  this->GetIntegrator().AddParamData("control",
728  &(q.GetSpacialVectorCopy()));
729  }
730  else
731  {
732  throw DOpEException("dopedim not implemented",
733  "StatReducedProblem::ComputeReducedConstraints");
734  }
735 
736  double ret = 0;
737  bool found = false;
738 
739  if (this->GetProblem()->GetConstraintType().find("domain")
740  != std::string::npos)
741  {
742  found = true;
743  ret += this->GetIntegrator().ComputeDomainScalar(
744  *(this->GetProblem()));
745  }
746  if (this->GetProblem()->GetConstraintType().find("point")
747  != std::string::npos)
748  {
749  found = true;
750  ret += this->GetIntegrator().ComputePointScalar(
751  *(this->GetProblem()));
752  }
753  if (this->GetProblem()->GetConstraintType().find("boundary")
754  != std::string::npos)
755  {
756  found = true;
757  ret += this->GetIntegrator().ComputeBoundaryScalar(
758  *(this->GetProblem()));
759  }
760  if (this->GetProblem()->GetConstraintType().find("face")
761  != std::string::npos)
762  {
763  found = true;
764  ret += this->GetIntegrator().ComputeFaceScalar(
765  *(this->GetProblem()));
766  }
767 
768  if (!found)
769  {
770  throw DOpEException(
771  "Unknown Constraint Type: "
772  + this->GetProblem()->GetConstraintType(),
773  "StatReducedProblem::ComputeReducedConstraints");
774  }
775  // global_values(i) = ret;
776  gc(i) = ret;
777 
778  if (dopedim == dealdim)
779  {
780  this->GetIntegrator().DeleteDomainData("control");
781  }
782  else if (dopedim == 0)
783  {
784  this->GetIntegrator().DeleteParamData("control");
785  q.UnLockCopy();
786  }
787  else
788  {
789  throw DOpEException("dopedim not implemented",
790  "StatReducedProblem::ComputeReducedConstraints");
791  }
792  this->GetIntegrator().DeleteDomainData("state");
793  }
794 
795  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
796  this->GetIntegrator());
797  //gc = global_values;
798  }
799 
800  //Check that no global in space, local in time constraints are given!
801  if (g.HasType("local_global_control") || g.HasType("local_global_state"))
802  {
803  throw DOpEException(
804  "There are global in space, local in time constraints given. In the stationary case they should be moved to global in space and time!",
805  "StatReducedProblem::ComputeReducedConstraints");
806  }
807 
808  //this->GetProblem()->PostProcessConstraints(g, true);
809  this->GetProblem()->PostProcessConstraints(g);
810 
811 
812  return g.IsFeasible();
813  }
814 
815  /******************************************************/
816 
817  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
818  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
819  typename VECTOR, int dopedim, int dealdim>
820  void
821  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
822  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetControlBoxConstraints(
823  ControlVector<VECTOR>& lb, ControlVector<VECTOR>& ub)
824  {
825  this->GetProblem()->GetControlBoxConstraints(lb.GetSpacialVector(),
826  ub.GetSpacialVector());
827  }
828 
829  /******************************************************/
830 
831  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
832  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
833  typename VECTOR, int dopedim, int dealdim>
834  void
835  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
836  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedAdjoint(
837  const ControlVector<VECTOR>& q)
838  {
839  this->GetOutputHandler()->Write("Computing Reduced Adjoint:",
840  4 + this->GetBasePriority());
841 
842  this->SetProblemType("adjoint");
843  if (_adjoint_reinit == true)
844  {
845  GetNonlinearSolver("adjoint").ReInit(*(this->GetProblem()));
846  _adjoint_reinit = false;
847  }
848 
849  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
850 
851  this->GetIntegrator().AddDomainData("state",
852  &(GetU().GetSpacialVector())); //&(GetU().GetSpacialVector())
853 
854  if (dopedim == dealdim)
855  {
856  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
857  }
858  else if (dopedim == 0)
859  {
860  this->GetIntegrator().AddParamData("control",
861  &(q.GetSpacialVectorCopy()));
862  }
863  else
864  {
865  throw DOpEException("dopedim not implemented",
866  "StatReducedProblem::ComputeReducedAdjoint");
867  }
868 
869  _build_adjoint_matrix =
870  this->GetNonlinearSolver("adjoint").NonlinearSolve(
871  *(this->GetProblem()), (GetZ().GetSpacialVector()), true,
872  _build_adjoint_matrix);
873 
874  if (dopedim == dealdim)
875  {
876  this->GetIntegrator().DeleteDomainData("control");
877  }
878  else if (dopedim == 0)
879  {
880  this->GetIntegrator().DeleteParamData("control");
881  q.UnLockCopy();
882  }
883  else
884  {
885  throw DOpEException("dopedim not implemented",
886  "StatReducedProblem::ComputeReducedAdjoint");
887  }
888  this->GetIntegrator().DeleteDomainData("state");
889 
890  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
891 
892  this->GetOutputHandler()->Write((GetZ().GetSpacialVector()),
893  "Adjoint" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
894  }
895 
896  /******************************************************/
897 
898  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
899  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
900  typename VECTOR, int dopedim, int dealdim>
901  void
902  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
903  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeDualForErrorEstimation(
904  const ControlVector<VECTOR>& q,
905  DOpEtypes::WeightComputation weight_comp)
906  {
907  this->GetOutputHandler()->Write("Computing Dual for Error Estimation:",
908  4 + this->GetBasePriority());
909 
910  if (weight_comp == DOpEtypes::higher_order_interpolation)
911  {
912  this->SetProblemType("adjoint_for_ee");
913  }
914  else
915  {
916  throw DOpEException("Unknown WeightComputation",
917  "StatPDEProblem::ComputeDualForErrorEstimation");
918  }
919 
920  // auto& problem = this->GetProblem()->GetStateProblem();//Hier ist adjoint problem einzufuegen
921  auto& problem = *(this->GetProblem());
922  if (_adjoint_reinit == true)
923  {
924  GetNonlinearSolver("adjoint_for_ee").ReInit(problem);
925  _adjoint_reinit = false;
926  }
927 
928  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
929 
930  this->GetIntegrator().AddDomainData("state",
931  &(GetU().GetSpacialVector())); //&(GetU().GetSpacialVector())
932 
933  if (dopedim == dealdim)
934  {
935  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
936  }
937  else if (dopedim == 0)
938  {
939  this->GetIntegrator().AddParamData("control",
940  &(q.GetSpacialVectorCopy()));
941  }
942  else
943  {
944  throw DOpEException("dopedim not implemented",
945  "StatReducedProblem::ComputeReducedAdjoint");
946  }
947 
948  _build_adjoint_matrix =
949  this->GetNonlinearSolver("adjoint_for_ee").NonlinearSolve(problem,
950  (GetZForEE().GetSpacialVector()), true, _build_adjoint_matrix);
951 
952  if (dopedim == dealdim)
953  {
954  this->GetIntegrator().DeleteDomainData("control");
955  }
956  else if (dopedim == 0)
957  {
958  this->GetIntegrator().DeleteParamData("control");
959  q.UnLockCopy();
960  }
961  else
962  {
963  throw DOpEException("dopedim not implemented",
964  "StatReducedProblem::ComputeReducedAdjoint");
965  }
966 
967  this->GetIntegrator().DeleteDomainData("state");
968 
969  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
970 
971  this->GetOutputHandler()->Write((GetZForEE().GetSpacialVector()),
972  "Adjoint_for_ee" + this->GetPostIndex(), problem.GetDoFType());
973 
974  }
975 
976  /******************************************************/
977 
978  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
979  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
980  typename VECTOR, int dopedim, int dealdim>
981  void
982  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
983  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedGradient(
984  const ControlVector<VECTOR>& q, ControlVector<VECTOR>& gradient,
985  ControlVector<VECTOR>& gradient_transposed)
986  {
987  this->ComputeReducedAdjoint(q);
988 
989  this->GetOutputHandler()->Write("Computing Reduced Gradient:",
990  4 + this->GetBasePriority());
991 
992  //Preparations for ControlInTheDirichletData
993  VECTOR tmp;
994  if (this->GetProblem()->HasControlInDirichletData())
995  {
996  tmp.reinit(GetU().GetSpacialVector());
997  this->SetProblemType("adjoint");
998  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
999 
1000  if (dopedim == dealdim)
1001  {
1002  this->GetIntegrator().AddDomainData("control",
1003  &(q.GetSpacialVector()));
1004  }
1005  else if (dopedim == 0)
1006  {
1007  this->GetIntegrator().AddParamData("control",
1008  &(q.GetSpacialVectorCopy()));
1009  }
1010  else
1011  {
1012  throw DOpEException("dopedim not implemented",
1013  "StatReducedProblem::ComputeReducedGradient");
1014  }
1015 
1016  this->GetIntegrator().AddDomainData("state",
1017  &(GetU().GetSpacialVector()));
1018  this->GetIntegrator().AddDomainData("last_newton_solution",
1019  &(GetZ().GetSpacialVector()));
1020 
1021  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1022  tmp, false);
1023  tmp *= -1.;
1024 
1025  if (dopedim == dealdim)
1026  {
1027  this->GetIntegrator().DeleteDomainData("control");
1028  }
1029  else if (dopedim == 0)
1030  {
1031  this->GetIntegrator().DeleteParamData("control");
1032  q.UnLockCopy();
1033  }
1034  else
1035  {
1036  throw DOpEException("dopedim not implemented",
1037  "StatReducedProblem::ComputeReducedGradient");
1038  }
1039 
1040  this->GetIntegrator().DeleteDomainData("state");
1041  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1042  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1043  this->GetIntegrator());
1044  }
1045  //Endof Dirichletdata Preparations
1046  this->SetProblemType("gradient");
1047  if (_gradient_reinit == true)
1048  {
1049  GetControlNonlinearSolver().ReInit(*(this->GetProblem()));
1050  _gradient_reinit = false;
1051  }
1052 
1053  this->GetProblem()->AddAuxiliaryToIntegrator(
1054  this->GetControlIntegrator());
1055 
1056  if (dopedim == dealdim)
1057  {
1058  this->GetControlIntegrator().AddDomainData("control",
1059  &(q.GetSpacialVector()));
1060  }
1061  else if (dopedim == 0)
1062  {
1063  this->GetControlIntegrator().AddParamData("control",
1064  &(q.GetSpacialVectorCopy()));
1065  }
1066  else
1067  {
1068  throw DOpEException("dopedim not implemented",
1069  "StatReducedProblem::ComputeReducedGradient");
1070  }
1071 
1072  this->GetControlIntegrator().AddDomainData("state",
1073  &(GetU().GetSpacialVector()));
1074  this->GetControlIntegrator().AddDomainData("adjoint",
1075  &(GetZ().GetSpacialVector()));
1076  if (this->GetProblem()->HasControlInDirichletData())
1077  this->GetControlIntegrator().AddDomainData("adjoint_residual", &tmp);
1078 
1079  gradient_transposed = 0.;
1080  if (dopedim == dealdim)
1081  {
1082  this->GetControlIntegrator().AddDomainData("last_newton_solution",
1083  &(gradient_transposed.GetSpacialVector()));
1084  this->GetControlIntegrator().ComputeNonlinearResidual(
1085  *(this->GetProblem()), gradient.GetSpacialVector(), true);
1086  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
1087  }
1088  else if (dopedim == 0)
1089  {
1090  this->GetControlIntegrator().AddParamData("last_newton_solution",
1091  &(gradient_transposed.GetSpacialVectorCopy()));
1092  this->GetControlIntegrator().ComputeNonlinearResidual(
1093  *(this->GetProblem()), gradient.GetSpacialVector(), true);
1094 
1095  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
1096  gradient_transposed.UnLockCopy();
1097 
1098  }
1099 
1100  gradient *= -1.;
1101  gradient_transposed = gradient;
1102 
1103  //Compute l^2 representation of the Gradient
1104 
1105  _build_control_matrix = this->GetControlNonlinearSolver().NonlinearSolve(
1106  *(this->GetProblem()), gradient_transposed.GetSpacialVector(), true,
1107  _build_control_matrix);
1108  if (dopedim == dealdim)
1109  {
1110  this->GetControlIntegrator().DeleteDomainData("control");
1111  }
1112  else if (dopedim == 0)
1113  {
1114  this->GetControlIntegrator().DeleteParamData("control");
1115  q.UnLockCopy();
1116  }
1117  else
1118  {
1119  throw DOpEException("dopedim not implemented",
1120  "StatReducedProblem::ComputeReducedGradient");
1121  }
1122  this->GetControlIntegrator().DeleteDomainData("state");
1123  this->GetControlIntegrator().DeleteDomainData("adjoint");
1124  if (this->GetProblem()->HasControlInDirichletData())
1125  this->GetControlIntegrator().DeleteDomainData("adjoint_residual");
1126 
1127  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1128  this->GetControlIntegrator());
1129 
1130  this->GetOutputHandler()->Write(gradient,
1131  "Gradient" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1132  this->GetOutputHandler()->Write(gradient_transposed,
1133  "Gradient_Transposed" + this->GetPostIndex(),
1134  this->GetProblem()->GetDoFType());
1135 
1136  }
1137 
1138  /******************************************************/
1139 
1140  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1141  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1142  typename VECTOR, int dopedim, int dealdim>
1143  double
1144  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1145  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedCostFunctional(
1146  const ControlVector<VECTOR>& q)
1147  {
1148  this->ComputeReducedState(q);
1149 
1150  double ret = 0;
1151  bool found = false;
1152 
1153  this->GetOutputHandler()->Write("Computing Cost Functional:",
1154  4 + this->GetBasePriority());
1155 
1156  this->SetProblemType("cost_functional");
1157 
1158  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1159 
1160  if (dopedim == dealdim)
1161  {
1162  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1163  }
1164  else if (dopedim == 0)
1165  {
1166  this->GetIntegrator().AddParamData("control",
1167  &(q.GetSpacialVectorCopy()));
1168  }
1169  else
1170  {
1171  throw DOpEException("dopedim not implemented",
1172  "StatReducedProblem::ComputeReducedCostFunctional");
1173  }
1174  this->GetIntegrator().AddDomainData("state",
1175  &(GetU().GetSpacialVector()));
1176 
1177  if (this->GetProblem()->GetFunctionalType().find("domain")
1178  != std::string::npos)
1179  {
1180  found = true;
1181  ret += this->GetIntegrator().ComputeDomainScalar(*(this->GetProblem()));
1182  }
1183  if (this->GetProblem()->GetFunctionalType().find("point")
1184  != std::string::npos)
1185  {
1186  found = true;
1187  ret += this->GetIntegrator().ComputePointScalar(*(this->GetProblem()));
1188  }
1189  if (this->GetProblem()->GetFunctionalType().find("boundary")
1190  != std::string::npos)
1191  {
1192  found = true;
1193  ret += this->GetIntegrator().ComputeBoundaryScalar(
1194  *(this->GetProblem()));
1195  }
1196  if (this->GetProblem()->GetFunctionalType().find("face")
1197  != std::string::npos)
1198  {
1199  found = true;
1200  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1201  }
1202 
1203  if (!found)
1204  {
1205  throw DOpEException(
1206  "Unknown Functional Type: "
1207  + this->GetProblem()->GetFunctionalType(),
1208  "StatReducedProblem::ComputeReducedCostFunctional");
1209  }
1210 
1211  if (dopedim == dealdim)
1212  {
1213  this->GetIntegrator().DeleteDomainData("control");
1214  }
1215  else if (dopedim == 0)
1216  {
1217  this->GetIntegrator().DeleteParamData("control");
1218  q.UnLockCopy();
1219  }
1220  else
1221  {
1222  throw DOpEException("dopedim not implemented",
1223  "StatReducedProblem::ComputeReducedCostFunctional");
1224  }
1225  this->GetIntegrator().DeleteDomainData("state");
1226  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1227 
1228  this->GetFunctionalValues()[0].push_back(ret);
1229  return ret;
1230  }
1231 
1232  /******************************************************/
1233 
1234  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1235  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1236  typename VECTOR, int dopedim, int dealdim>
1237  void
1238  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1239  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedFunctionals(
1240  const ControlVector<VECTOR>& q)
1241  {
1242  this->GetOutputHandler()->Write("Computing Functionals:",
1243  4 + this->GetBasePriority());
1244 
1245  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1246 
1247  if (dopedim == dealdim)
1248  {
1249  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1250  }
1251  else if (dopedim == 0)
1252  {
1253  this->GetIntegrator().AddParamData("control",
1254  &(q.GetSpacialVectorCopy()));
1255  }
1256  else
1257  {
1258  throw DOpEException("dopedim not implemented",
1259  "StatReducedProblem::ComputeReducedFunctionals");
1260  }
1261  this->GetIntegrator().AddDomainData("state",
1262  &(GetU().GetSpacialVector()));
1263 
1264  for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
1265  {
1266  double ret = 0;
1267  bool found = false;
1268 
1269  this->SetProblemType("aux_functional", i);
1270  if (this->GetProblem()->GetFunctionalType().find("domain")
1271  != std::string::npos)
1272  {
1273  found = true;
1274  ret += this->GetIntegrator().ComputeDomainScalar(
1275  *(this->GetProblem()));
1276  }
1277  if (this->GetProblem()->GetFunctionalType().find("point")
1278  != std::string::npos)
1279  {
1280  found = true;
1281  ret += this->GetIntegrator().ComputePointScalar(
1282  *(this->GetProblem()));
1283  }
1284  if (this->GetProblem()->GetFunctionalType().find("boundary")
1285  != std::string::npos)
1286  {
1287  found = true;
1288  ret += this->GetIntegrator().ComputeBoundaryScalar(
1289  *(this->GetProblem()));
1290  }
1291  if (this->GetProblem()->GetFunctionalType().find("face")
1292  != std::string::npos)
1293  {
1294  found = true;
1295  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1296  }
1297 
1298  if (!found)
1299  {
1300  throw DOpEException(
1301  "Unknown Functional Type: "
1302  + this->GetProblem()->GetFunctionalType(),
1303  "StatReducedProblem::ComputeReducedFunctionals");
1304  }
1305  this->GetFunctionalValues()[i + 1].push_back(ret);
1306  std::stringstream out;
1307  this->GetOutputHandler()->InitOut(out);
1308  out << this->GetProblem()->GetFunctionalName() << ": " << ret;
1309  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
1310  }
1311 
1312  if (dopedim == dealdim)
1313  {
1314  this->GetIntegrator().DeleteDomainData("control");
1315  }
1316  else if (dopedim == 0)
1317  {
1318  this->GetIntegrator().DeleteParamData("control");
1319  q.UnLockCopy();
1320  }
1321  else
1322  {
1323  throw DOpEException("dopedim not implemented",
1324  "StatReducedProblem::ComputeReducedFunctionals");
1325  }
1326  this->GetIntegrator().DeleteDomainData("state");
1327  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1328 
1329  }
1330 
1331  /******************************************************/
1332 
1333  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1334  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1335  typename VECTOR, int dopedim, int dealdim>
1336  template<class DWRC,class PDE>
1337  void
1338  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1339  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeRefinementIndicators(
1340  const ControlVector<VECTOR>& q, DWRC& dwrc, PDE& pde)
1341 // StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1342 // CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeRefinementIndicators(
1343 // const ControlVector<VECTOR>& q, DWRDataContainerBase<VECTOR>& dwrc)
1344  {
1345  //Attach the ResidualModifier to the PDE.
1346  pde.ResidualModifier = boost::bind<void>(boost::mem_fn(&DWRC::ResidualModifier),boost::ref(dwrc),_1);
1347  pde.VectorResidualModifier = boost::bind<void>(boost::mem_fn(&DWRC::VectorResidualModifier),boost::ref(dwrc),_1);
1348 
1349  //first we reinit the dwrdatacontainer (this
1350  //sets the weight-vectors to their correct length)
1351  const unsigned int n_elements =
1352  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler().get_tria().n_active_cells();
1353  dwrc.ReInit(n_elements);
1354 
1355  //Estimation for Costfunctional or if no dual is needed
1356  if(this->GetProblem()->EEFunctionalIsCost() || !dwrc.NeedDual())
1357  {
1358  this->GetOutputHandler()->Write("Computing Error Indicators:",
1359  4 + this->GetBasePriority());
1360  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1361 
1362  //add the primal and dual solution to the integrator
1363  this->GetIntegrator().AddDomainData("state",
1364  &(GetU().GetSpacialVector()));
1365  this->GetIntegrator().AddDomainData("adjoint_for_ee",
1366  &(GetZ().GetSpacialVector()));
1367 
1368  if (dopedim == dealdim)
1369  {
1370  this->GetIntegrator().AddDomainData("control",
1371  &(q.GetSpacialVector()));
1372  }
1373  else if (dopedim == 0)
1374  {
1375  this->GetIntegrator().AddParamData("control",
1376  &(q.GetSpacialVectorCopy()));
1377  }
1378  else
1379  {
1380  throw DOpEException("dopedim not implemented",
1381  "StatReducedProblem::ComputeRefinementIndicators");
1382  }
1383 
1384  this->SetProblemType("error_evaluation");
1385 
1386  //prepare the weights...
1387  dwrc.PrepareWeights(GetU(), GetZ());
1388  dwrc.PrepareWeights(q);
1389 
1390  //now we finally compute the refinement indicators
1391  this->GetIntegrator().ComputeRefinementIndicators(*this->GetProblem(),
1392  dwrc);
1393  // release the lock on the refinement indicators (see dwrcontainer.h)
1394  dwrc.ReleaseLock();
1395  dwrc.ClearWeightData();
1396 
1397  // clear the data
1398  if (dopedim == dealdim)
1399  {
1400  this->GetIntegrator().DeleteDomainData("control");
1401  }
1402  else if (dopedim == 0)
1403  {
1404  this->GetIntegrator().DeleteParamData("control");
1405  q.UnLockCopy();
1406  }
1407  else
1408  {
1409  throw DOpEException("dopedim not implemented",
1410  "StatReducedProblem::ComputeRefinementIndicators");
1411  }
1412  this->GetIntegrator().DeleteDomainData("state");
1413  this->GetIntegrator().DeleteDomainData("adjoint_for_ee");
1414  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1415  this->GetIntegrator());
1416  }
1417  else //Estimation for other (not the cost) functional
1418  {
1419  throw DOpEException("Estimating the error in other functionals than cost is not implemented",
1420  "StatReducedProblem::ComputeRefinementIndicators");
1421 
1422  }
1423 
1424  std::stringstream out;
1425  this->GetOutputHandler()->InitOut(out);
1426  out << "Error estimate using "<<dwrc.GetName();
1427  if(dwrc.NeedDual())
1428  out<<" For the computation of "<<this->GetProblem()->GetFunctionalName();
1429  out<< ": "<< dwrc.GetError();
1430  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
1431  }
1432 
1433  /******************************************************/
1434 
1435  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1436  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1437  typename VECTOR, int dopedim, int dealdim>
1438  void
1439  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1440  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedHessianVector(
1441  const ControlVector<VECTOR>& q, const ControlVector<VECTOR>& direction,
1442  ControlVector<VECTOR>& hessian_direction,
1443  ControlVector<VECTOR>& hessian_direction_transposed)
1444  {
1445  this->GetOutputHandler()->Write("Computing ReducedHessianVector:",
1446  4 + this->GetBasePriority());
1447  this->GetOutputHandler()->Write("\tSolving Tangent:",
1448  5 + this->GetBasePriority());
1449 
1450  this->SetProblemType("tangent");
1451 
1452  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1453 
1454  this->GetIntegrator().AddDomainData("state",
1455  &(GetU().GetSpacialVector()));
1456  this->GetControlIntegrator().AddDomainData("state",
1457  &(GetU().GetSpacialVector()));
1458 
1459  if (dopedim == dealdim)
1460  {
1461  this->GetIntegrator().AddDomainData("dq",
1462  &(direction.GetSpacialVector()));
1463  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1464  }
1465  else if (dopedim == 0)
1466  {
1467  this->GetIntegrator().AddParamData("dq",
1468  &(direction.GetSpacialVectorCopy()));
1469  this->GetIntegrator().AddParamData("control",
1470  &(q.GetSpacialVectorCopy()));
1471  }
1472  else
1473  {
1474  throw DOpEException("dopedim not implemented",
1475  "StatReducedProblem::ComputeReducedHessianVector");
1476  }
1477 
1478  //tangent Matrix is the same as state matrix
1479  _build_state_matrix = this->GetNonlinearSolver("tangent").NonlinearSolve(
1480  *(this->GetProblem()), (GetDU().GetSpacialVector()), true,
1481  _build_state_matrix);
1482 
1483  this->GetOutputHandler()->Write((GetDU().GetSpacialVector()),
1484  "Tangent" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1485 
1486  this->GetOutputHandler()->Write("\tSolving Adjoint Hessian:",
1487  5 + this->GetBasePriority());
1488  this->SetProblemType("adjoint_hessian");
1489  this->GetIntegrator().AddDomainData("adjoint",
1490  &(GetZ().GetSpacialVector()));
1491  this->GetIntegrator().AddDomainData("tangent",
1492  &(GetDU().GetSpacialVector()));
1493  this->GetControlIntegrator().AddDomainData("adjoint",
1494  &(GetZ().GetSpacialVector()));
1495  this->GetControlIntegrator().AddDomainData("tangent",
1496  &(GetDU().GetSpacialVector()));
1497 
1498  //adjoint_hessian Matrix is the same as adjoint matrix
1499  _build_adjoint_matrix =
1500  this->GetNonlinearSolver("adjoint_hessian").NonlinearSolve(
1501  *(this->GetProblem()), (GetDZ().GetSpacialVector()), true,
1502  _build_adjoint_matrix);
1503 
1504  this->GetOutputHandler()->Write((GetDZ().GetSpacialVector()),
1505  "Hessian" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1506 
1507  this->GetIntegrator().AddDomainData("adjoint_hessian",
1508  &(GetDZ().GetSpacialVector()));
1509  this->GetControlIntegrator().AddDomainData("adjoint_hessian",
1510  &(GetDZ().GetSpacialVector()));
1511 
1512  this->GetOutputHandler()->Write(
1513  "\tComputing Representation of the Hessian:",
1514  5 + this->GetBasePriority());
1515  //Preparations for Control In The Dirichlet Data
1516  VECTOR tmp;
1517  VECTOR tmp_second;
1518  if (this->GetProblem()->HasControlInDirichletData())
1519  {
1520  tmp.reinit(GetU().GetSpacialVector());
1521  tmp_second.reinit(GetU().GetSpacialVector());
1522  this->SetProblemType("adjoint");
1523  this->GetIntegrator().AddDomainData("last_newton_solution",
1524  &(GetZ().GetSpacialVector()));
1525 
1526  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1527  tmp_second, false);
1528  tmp_second *= -1.;
1529 
1530  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1531  this->SetProblemType("adjoint_hessian");
1532  this->GetIntegrator().AddDomainData("last_newton_solution",
1533  &(GetDZ().GetSpacialVector()));
1534 
1535  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1536  tmp, false);
1537  tmp *= -1.;
1538 
1539  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1540  }
1541  //Endof Dirichletdata Preparations
1542  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1543 
1544  this->SetProblemType("hessian");
1545  this->GetProblem()->AddAuxiliaryToIntegrator(
1546  this->GetControlIntegrator());
1547  if (dopedim == dealdim)
1548  {
1549  this->GetIntegrator().DeleteDomainData("dq");
1550  this->GetIntegrator().DeleteDomainData("control");
1551  this->GetControlIntegrator().AddDomainData("dq",
1552  &(direction.GetSpacialVector()));
1553  this->GetControlIntegrator().AddDomainData("control",
1554  &(q.GetSpacialVector()));
1555  }
1556  else if (dopedim == 0)
1557  {
1558  this->GetIntegrator().DeleteParamData("dq");
1559  this->GetIntegrator().DeleteParamData("control");
1560  direction.UnLockCopy();
1561  q.UnLockCopy();
1562  this->GetControlIntegrator().AddParamData("dq",
1563  &(direction.GetSpacialVectorCopy()));
1564  this->GetControlIntegrator().AddParamData("control",
1565  &(q.GetSpacialVectorCopy()));
1566  }
1567  else
1568  {
1569  throw DOpEException("dopedim not implemented",
1570  "StatReducedProblem::ComputeReducedHessianVector");
1571  }
1572  if (this->GetProblem()->HasControlInDirichletData())
1573  {
1574  this->GetControlIntegrator().AddDomainData("adjoint_residual", &tmp);
1575  this->GetControlIntegrator().AddDomainData("hessian_residual",
1576  &tmp_second);
1577  }
1578 
1579  {
1580  hessian_direction_transposed = 0.;
1581  if (dopedim == dealdim)
1582  {
1583  this->GetControlIntegrator().AddDomainData("last_newton_solution",
1584  &(hessian_direction_transposed.GetSpacialVector()));
1585  this->GetControlIntegrator().ComputeNonlinearResidual(
1586  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
1587  true);
1588  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
1589  }
1590  else if (dopedim == 0)
1591  {
1592  this->GetControlIntegrator().AddParamData("last_newton_solution",
1593  &(hessian_direction_transposed.GetSpacialVectorCopy()));
1594  this->GetControlIntegrator().ComputeNonlinearResidual(
1595  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
1596  true);
1597  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
1598  hessian_direction_transposed.UnLockCopy();
1599  }
1600  hessian_direction *= -1.;
1601  hessian_direction_transposed = hessian_direction;
1602  //Compute l^2 representation of the HessianVector
1603  //hessian Matrix is the same as control matrix
1604  _build_control_matrix =
1605  this->GetControlNonlinearSolver().NonlinearSolve(
1606  *(this->GetProblem()),
1607  hessian_direction_transposed.GetSpacialVector(), true,
1608  _build_control_matrix);
1609 
1610  this->GetOutputHandler()->Write(hessian_direction,
1611  "HessianDirection" + this->GetPostIndex(),
1612  this->GetProblem()->GetDoFType());
1613  this->GetOutputHandler()->Write(hessian_direction_transposed,
1614  "HessianDirection_Transposed" + this->GetPostIndex(),
1615  this->GetProblem()->GetDoFType());
1616  }
1617 
1618  if (dopedim == dealdim)
1619  {
1620  this->GetControlIntegrator().DeleteDomainData("dq");
1621  this->GetControlIntegrator().DeleteDomainData("control");
1622  }
1623  else if (dopedim == 0)
1624  {
1625  this->GetControlIntegrator().DeleteParamData("dq");
1626  this->GetControlIntegrator().DeleteParamData("control");
1627  direction.UnLockCopy();
1628  q.UnLockCopy();
1629  }
1630  else
1631  {
1632  throw DOpEException("dopedim not implemented",
1633  "StatReducedProblem::ComputeReducedHessianVector");
1634  }
1635  this->GetIntegrator().DeleteDomainData("state");
1636  this->GetIntegrator().DeleteDomainData("adjoint");
1637  this->GetIntegrator().DeleteDomainData("tangent");
1638  this->GetIntegrator().DeleteDomainData("adjoint_hessian");
1639  this->GetControlIntegrator().DeleteDomainData("state");
1640  this->GetControlIntegrator().DeleteDomainData("adjoint");
1641  this->GetControlIntegrator().DeleteDomainData("tangent");
1642  this->GetControlIntegrator().DeleteDomainData("adjoint_hessian");
1643  if (this->GetProblem()->HasControlInDirichletData())
1644  {
1645  this->GetControlIntegrator().DeleteDomainData("adjoint_residual");
1646  this->GetControlIntegrator().DeleteDomainData("hessian_residual");
1647  }
1648  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1649  this->GetControlIntegrator());
1650 
1651  }
1652 
1653  /******************************************************/
1654 
1655  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1656  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1657  typename VECTOR, int dopedim, int dealdim>
1658  void
1659  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1660  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedGradientOfGlobalConstraints(
1661  unsigned int num, const ControlVector<VECTOR>& q,
1662  const ConstraintVector<VECTOR>& g, ControlVector<VECTOR>& gradient,
1663  ControlVector<VECTOR>& gradient_transposed)
1664  {
1665  //FIXME: If the global constraints depend on u we need to calculate a corresponding
1666  // dual solution before we can calculate the gradient.
1667  std::stringstream out;
1668  out << "Computing Reduced Gradient of global constraint " << num << " :";
1669  this->GetOutputHandler()->Write(out, 4 + this->GetBasePriority());
1670  //Compute derivatives of global constraints
1671  this->SetProblemType("global_constraint_gradient", num);
1672 
1673  if (dopedim == dealdim)
1674  {
1675  this->GetControlIntegrator().AddDomainData("control",
1676  &(q.GetSpacialVector()));
1677  }
1678  else if (dopedim == 0)
1679  {
1680  this->GetControlIntegrator().AddParamData("control",
1681  &(q.GetSpacialVectorCopy()));
1682  }
1683  else
1684  {
1685  throw DOpEException("dopedim not implemented",
1686  "StatReducedProblem::ComputeReducedGradient");
1687  }
1688  this->GetProblem()->AddAuxiliaryToIntegrator(
1689  this->GetControlIntegrator());
1690  this->GetControlIntegrator().AddDomainData("constraints_local",
1691  &g.GetSpacialVector("local"));
1692  this->GetControlIntegrator().AddParamData("constraints_global",
1693  &g.GetGlobalConstraints());
1694 
1695  //Compute
1696  this->GetControlIntegrator().ComputeNonlinearRhs(*(this->GetProblem()),
1697  gradient.GetSpacialVector(), true);
1698  gradient_transposed = gradient;
1699 
1700  this->GetControlIntegrator().DeleteDomainData("constraints_local");
1701  this->GetControlIntegrator().DeleteParamData("constraints_global");
1702  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1703  this->GetControlIntegrator());
1704  if (dopedim == dealdim)
1705  {
1706  this->GetControlIntegrator().DeleteDomainData("control");
1707  }
1708  else if (dopedim == 0)
1709  {
1710  this->GetControlIntegrator().DeleteParamData("control");
1711  q.UnLockCopy();
1712  }
1713  else
1714  {
1715  throw DOpEException("dopedim not implemented",
1716  "StatReducedProblem::ComputeReducedGradient");
1717  }
1718  }
1719 
1720  /******************************************************/
1721 
1722  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1723  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1724  typename VECTOR, int dopedim, int dealdim>
1725  void
1726  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1727  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(
1728  const VECTOR &v, std::string name, std::string outfile,
1729  std::string dof_type, std::string filetype)
1730  {
1731  if (dof_type == "state")
1732  {
1733  auto& data_out =
1734  this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1735  data_out.attach_dof_handler(
1736  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
1737 
1738  data_out.add_data_vector(v, name);
1739  data_out.build_patches();
1740 
1741  std::ofstream output(outfile.c_str());
1742 
1743  if (filetype == ".vtk")
1744  {
1745  data_out.write_vtk(output);
1746  }
1747  else if (filetype == ".gpl")
1748  {
1749  data_out.write_gnuplot(output);
1750  }
1751  else
1752  {
1753  throw DOpEException(
1754  "Don't know how to write filetype `" + filetype + "'!",
1755  "StatReducedProblem::WriteToFile");
1756  }
1757  data_out.clear();
1758  }
1759  else if (dof_type == "control")
1760  {
1761 #if dope_dimension >0
1762  auto& data_out = this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1763  data_out.attach_dof_handler (this->GetProblem()->GetSpaceTimeHandler()->GetControlDoFHandler());
1764 
1765  data_out.add_data_vector (v,name);
1766  data_out.build_patches ();
1767 
1768  std::ofstream output(outfile.c_str());
1769 
1770  if(filetype == ".vtk")
1771  {
1772  data_out.write_vtk (output);
1773  }
1774  else if (filetype == ".gpl")
1775  {
1776  data_out.write_gnuplot(output);
1777  }
1778  else
1779  {
1780  throw DOpEException("Don't know how to write filetype `" + filetype + "'!",
1781  "StatReducedProblem::WriteToFile");
1782  }
1783  data_out.clear();
1784 #else
1785  if (filetype == ".txt")
1786  {
1787  std::ofstream output(outfile.c_str());
1788  Vector<double> off;
1789  off = v;
1790  for (unsigned int i = 0; i < off.size(); i++)
1791  {
1792  output << off(i) << std::endl;
1793  }
1794  }
1795  else
1796  {
1797  throw DOpEException(
1798  "Don't know how to write filetype `" + filetype + "'!",
1799  "StatReducedProblem::WriteToFile");
1800  }
1801 #endif
1802  }
1803  else
1804  {
1805  throw DOpEException("No such DoFHandler `" + dof_type + "'!",
1806  "StatReducedProblem::WriteToFile");
1807  }
1808  }
1809 
1810  /******************************************************/
1811 
1812  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1813  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1814  typename VECTOR, int dopedim, int dealdim>
1815  void
1816  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1817  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(
1818  const ControlVector<VECTOR> &v, std::string name, std::string outfile,
1819  std::string dof_type, std::string filetype)
1820  {
1821  WriteToFile(v.GetSpacialVector(), name, outfile, dof_type, filetype);
1822  }
1823 
1824 /******************************************************/
1825 }
1826 #endif
DOpE::StatReducedProblem::ComputeRefinementIndicators
void ComputeRefinementIndicators(const ControlVector< VECTOR > &q, DWRC &dwrc, PDE &pde)
Definition: statreducedproblem.h:1339
DOpE::StatReducedProblem::GetControlIntegrator
CONTROLINTEGRATOR & GetControlIntegrator()
Definition: statreducedproblem.h:397
DOpE::ControlVector::UnLockCopy
void UnLockCopy() const
Definition: controlvector.h:196
DOpE::ConstraintVector
Definition: constraintvector.h:47
DOpE::ControlVector::GetSpacialVectorCopy
const dealii::Vector< double > & GetSpacialVectorCopy() const
Definition: controlvector.cc:132
DOpE::StatReducedProblem::GetNonlinearSolver
NONLINEARSOLVER & GetNonlinearSolver(std::string type)
Definition: statreducedproblem.h:519
DOpE::StatReducedProblem::GetControlNonlinearSolver
CONTROLNONLINEARSOLVER & GetControlNonlinearSolver()
Definition: statreducedproblem.h:545
DOpE::StatReducedProblem::ComputeReducedFunctionals
void ComputeReducedFunctionals(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:1239
DOpE::StatReducedProblem::ReInit
void ReInit()
Definition: statreducedproblem.h:567
DOpE::StatReducedProblem::GetIntegrator
INTEGRATOR & GetIntegrator()
Definition: statreducedproblem.h:392
DOpE::StatReducedProblem::StatReducedProblem
StatReducedProblem(PROBLEM *OP, std::string state_behavior, ParameterReader &param_reader, INTEGRATORDATACONT &idc, int base_priority=0)
Definition: statreducedproblem.h:447
DOpE::StatReducedProblem::ComputeDualForErrorEstimation
void ComputeDualForErrorEstimation(const ControlVector< VECTOR > &q, DOpEtypes::WeightComputation weight_comp)
Definition: statreducedproblem.h:903
DOpE::StatReducedProblem::GetU
StateVector< VECTOR > & GetU()
Definition: statreducedproblem.h:354
DOpE::ParameterReader
Definition: parameterreader.h:36
DOpE::StatReducedProblem::WriteToFile
void WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statreducedproblem.h:1727
DOpE::StatReducedProblem::ComputeReducedAdjoint
void ComputeReducedAdjoint(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:836
DOpE::ConstraintVector::GetGlobalConstraints
const dealii::Vector< double > & GetGlobalConstraints() const
Definition: constraintvector.cc:174
DOpE::StatReducedProblem
Definition: statreducedproblem.h:83
DOpE::StatReducedProblem::GetZ
StateVector< VECTOR > & GetZ()
Definition: statreducedproblem.h:359
DOpE::StatReducedProblem::GetZForEE
StateVector< VECTOR > & GetZForEE()
Definition: statreducedproblem.h:382
DOpE::ControlVector
Definition: controlvector.h:48
DOpE::StatReducedProblem::declare_params
static void declare_params(ParameterReader &param_reader)
Definition: statreducedproblem.h:435
DOpE::StatReducedProblem::GetControlBoxConstraints
void GetControlBoxConstraints(ControlVector< VECTOR > &lb, ControlVector< VECTOR > &ub)
Definition: statreducedproblem.h:822
DOpE::StatReducedProblem::ComputeReducedCostFunctional
double ComputeReducedCostFunctional(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:1145
DOpE::SolutionExtractor
Definition: solutionextractor.h:40
DOpE::StatReducedProblem::GetZForEE
const StateVector< VECTOR > & GetZForEE() const
Definition: statreducedproblem.h:377
DOpE::StatReducedProblem::GetDZ
StateVector< VECTOR > & GetDZ()
Definition: statreducedproblem.h:369
DOpE::StatReducedProblem::ComputeReducedGradient
void ComputeReducedGradient(const ControlVector< VECTOR > &q, ControlVector< VECTOR > &gradient, ControlVector< VECTOR > &gradient_transposed)
Definition: statreducedproblem.h:983
DOpE::StateVector
Definition: statevector.h:49
DOpE::ControlVector::GetSpacialVector
VECTOR & GetSpacialVector()
Definition: controlvector.cc:118
DOpE::ConstraintVector::HasType
bool HasType(std::string name) const
Definition: constraintvector.cc:119
DOpE::DOpEtypes::WeightComputation
WeightComputation
Definition: dopetypes.h:78
DOpE::StatReducedProblem::ComputeReducedState
void ComputeReducedState(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:598
DOpE::ReducedProblemInterface
Definition: reducedprobleminterface.h:338
DOpE::StatReducedProblem::ComputeReducedHessianVector
void ComputeReducedHessianVector(const ControlVector< VECTOR > &q, const ControlVector< VECTOR > &direction, ControlVector< VECTOR > &hessian_direction, ControlVector< VECTOR > &hessian_direction_transposed)
Definition: statreducedproblem.h:1440
DOpE::StatReducedProblem::ComputeReducedConstraints
bool ComputeReducedConstraints(const ControlVector< VECTOR > &q, ConstraintVector< VECTOR > &g)
Definition: statreducedproblem.h:661
DOpE::StatReducedProblem::StateSizeInfo
void StateSizeInfo(std::stringstream &out)
Definition: statreducedproblem.h:258
DOpE::ConstraintVector::GetSpacialVector
VECTOR & GetSpacialVector(std::string name)
Definition: constraintvector.cc:130
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::StatReducedProblem::WriteToFile
void WriteToFile(const std::vector< double > &, std::string)
Definition: statreducedproblem.h:303
DOpE::StatReducedProblem::GetU
const StateVector< VECTOR > & GetU() const
Definition: statreducedproblem.h:349
DOpE::ConstraintVector::IsFeasible
virtual bool IsFeasible() const
Definition: constraintvector.cc:574
DOpE::StatReducedProblem::GetDU
StateVector< VECTOR > & GetDU()
Definition: statreducedproblem.h:364
DOpE::ReducedProblemInterface::ReInit
virtual void ReInit()
Definition: reducedprobleminterface.h:360
DOpE::StatReducedProblem::~StatReducedProblem
~StatReducedProblem()
Definition: statreducedproblem.h:508
DOpE::StatReducedProblem::ComputeReducedGradientOfGlobalConstraints
void ComputeReducedGradientOfGlobalConstraints(unsigned int num, const ControlVector< VECTOR > &q, const ConstraintVector< VECTOR > &g, ControlVector< VECTOR > &gradient, ControlVector< VECTOR > &gradient_transposed)
Definition: statreducedproblem.h:1660
Generated on Fri Oct 24 2014 17:56:40 for DOpE by   doxygen 1.8.6