24 #ifndef STAT_PDE_PROBLEM_H_
25 #define STAT_PDE_PROBLEM_H_
33 #include <lac/vector.h>
35 #include <lac/block_sparsity_pattern.h>
36 #include <lac/block_sparse_matrix.h>
49 #include <base/data_out_base.h>
50 #include <base/utilities.h>
51 #include <numerics/data_out.h>
52 #include <numerics/matrix_tools.h>
53 #include <numerics/vector_tools.h>
54 #include <base/function.h>
55 #include <lac/sparse_matrix.h>
56 #include <lac/compressed_simple_sparsity_pattern.h>
57 #include <lac/block_sparsity_pattern.h>
58 #include <lac/sparse_direct.h>
77 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
78 typename VECTOR,
int dealdim>
94 template<
typename INTEGRATORDATACONT>
97 int base_priority = 0);
111 template<
typename INTEGRATORDATACONT>
114 INTEGRATORDATACONT& idc2,
int base_priority = 0);
178 template<
class DWRC,
class PDE>
192 GetU().PrintInfos(out);
208 WriteToFile(
const VECTOR &v, std::string name, std::string outfile,
209 std::string dof_type, std::string filetype);
228 std::string outfile, std::string dof_type, std::string filetype);
317 return build_state_matrix_;
323 return build_adjoint_matrix_;
329 return state_reinit_;
335 return adjoint_reinit_;
342 return build_state_matrix_;
348 return build_adjoint_matrix_;
354 return state_reinit_;
360 return adjoint_reinit_;
393 StateVector<VECTOR> u_;
394 StateVector<VECTOR> z_for_ee_;
396 INTEGRATOR integrator_;
397 NONLINEARSOLVER nonlinear_state_solver_;
398 NONLINEARSOLVER nonlinear_adjoint_solver_;
400 bool build_state_matrix_;
401 bool build_adjoint_matrix_;
403 bool adjoint_reinit_;
407 friend class SolutionExtractor<
408 StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>,
415 using namespace dealii;
418 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
419 typename VECTOR,
int dealdim>
424 NONLINEARSOLVER::declare_params(param_reader);
427 Patterns::Integer(0));
432 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
433 typename VECTOR,
int dealdim>
434 template<
typename INTEGRATORDATACONT>
440 OP->GetSpaceTimeHandler(), state_behavior, param_reader), z_for_ee_(
441 OP->GetSpaceTimeHandler(), state_behavior, param_reader), integrator_(
442 idc), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
443 integrator_, param_reader)
446 n_patches_ = param_reader.
get_integer(
"number of patches");
449 state_reinit_ =
true;
450 adjoint_reinit_ =
true;
454 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
455 typename VECTOR,
int dealdim>
456 template<
typename INTEGRATORDATACONT>
460 INTEGRATORDATACONT& idc2,
int base_priority)
462 OP->GetSpaceTimeHandler(), state_behavior, param_reader), z_for_ee_(
463 OP->GetSpaceTimeHandler(), state_behavior, param_reader), integrator_(
464 idc, idc2), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
465 integrator_, param_reader)
469 state_reinit_ =
true;
470 adjoint_reinit_ =
true;
476 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
477 typename VECTOR,
int dealdim>
484 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
485 typename VECTOR,
int dealdim>
492 return nonlinear_state_solver_;
494 else if (type ==
"adjoint_for_ee")
496 return nonlinear_adjoint_solver_;
500 throw DOpEException(
"No Solver for Problem type:`" + type +
"' found",
501 "StatReducedProblem::GetNonlinearSolver");
508 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
509 typename VECTOR,
int dealdim>
518 state_reinit_ =
true;
519 adjoint_reinit_ =
true;
522 build_state_matrix_ =
true;
523 build_adjoint_matrix_ =
true;
526 GetZForEE().ReInit();
531 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
532 typename VECTOR,
int dealdim>
536 this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
537 this->GetOutputHandler()->Write(
"Computing State Solution:",
538 4 + this->GetBasePriority());
540 this->SetProblemType(
"state");
541 auto& problem = this->GetProblem()->GetStateProblem();
542 if (state_reinit_ ==
true)
544 GetNonlinearSolver(
"state").ReInit(problem);
545 state_reinit_ =
false;
548 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
550 build_state_matrix_ = this->GetNonlinearSolver(
"state").NonlinearSolve(
551 problem, (GetU().GetSpacialVector()),
true, build_state_matrix_);
554 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
556 this->GetOutputHandler()->Write((GetU().GetSpacialVector()),
557 "State" + this->GetPostIndex(), problem.GetDoFType());
563 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
564 typename VECTOR,
int dealdim>
569 this->GetOutputHandler()->Write(
"Computing Dual for Error Estimation:",
570 4 + this->GetBasePriority());
573 this->SetProblemType(
"adjoint_for_ee");
578 "StatPDEProblem::ComputeDualForErrorEstimation");
582 auto& problem = *(this->GetProblem());
583 if (adjoint_reinit_ ==
true)
585 GetNonlinearSolver(
"adjoint_for_ee").ReInit(problem);
586 adjoint_reinit_ =
false;
589 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
591 this->GetIntegrator().AddDomainData(
"state",
592 &(GetU().GetSpacialVector()));
595 build_adjoint_matrix_ =
596 this->GetNonlinearSolver(
"adjoint_for_ee").NonlinearSolve(problem,
597 (GetZForEE().GetSpacialVector()),
true, build_adjoint_matrix_);
599 this->GetIntegrator().DeleteDomainData(
"state");
602 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
604 this->GetOutputHandler()->Write((GetZForEE().GetSpacialVector()),
605 "Adjoint_for_ee" + this->GetPostIndex(), problem.GetDoFType());
611 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
612 typename VECTOR,
int dealdim>
616 this->ComputeReducedState();
618 this->GetOutputHandler()->Write(
"Computing Functionals:",
619 4 + this->GetBasePriority());
621 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
623 this->GetIntegrator().AddDomainData(
"state",
624 &(GetU().GetSpacialVector()));
627 for (
unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
632 this->SetProblemType(
"aux_functional", i);
633 if (this->GetProblem()->GetFunctionalType().find(
"domain")
634 != std::string::npos)
637 ret += this->GetIntegrator().ComputeDomainScalar(
638 *(this->GetProblem()));
640 if (this->GetProblem()->GetFunctionalType().find(
"point")
641 != std::string::npos)
644 ret += this->GetIntegrator().ComputePointScalar(
645 *(this->GetProblem()));
647 if (this->GetProblem()->GetFunctionalType().find(
"boundary")
648 != std::string::npos)
651 ret += this->GetIntegrator().ComputeBoundaryScalar(
652 *(this->GetProblem()));
654 if (this->GetProblem()->GetFunctionalType().find(
"face")
655 != std::string::npos)
658 ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
664 "Unknown Functional Type: "
665 + this->GetProblem()->GetFunctionalType(),
666 "StatPDEProblem::ComputeReducedFunctionals");
668 this->GetFunctionalValues()[i].push_back(ret);
669 std::stringstream out;
670 this->GetOutputHandler()->InitOut(out);
671 out << this->GetProblem()->GetFunctionalName() <<
": " << ret;
672 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
675 this->GetIntegrator().DeleteDomainData(
"state");
678 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
683 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
684 typename VECTOR,
int dealdim>
689 this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
692 this->GetOutputHandler()->Write(
"Computing Functionals:",
693 4 + this->GetBasePriority());
695 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
697 this->GetIntegrator().AddDomainData(
"state",
698 &(GetU().GetSpacialVector()));
701 for (
unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
706 this->SetProblemType(
"aux_functional", i);
707 if (this->GetProblem()->GetFunctionalType().find(
"domain")
708 != std::string::npos)
711 ret += this->GetIntegrator().ComputeDomainScalar(
712 *(this->GetProblem()));
714 if (this->GetProblem()->GetFunctionalType().find(
"point")
715 != std::string::npos)
718 ret += this->GetIntegrator().ComputePointScalar(
719 *(this->GetProblem()));
721 if (this->GetProblem()->GetFunctionalType().find(
"boundary")
722 != std::string::npos)
725 ret += this->GetIntegrator().ComputeBoundaryScalar(
726 *(this->GetProblem()));
728 if (this->GetProblem()->GetFunctionalType().find(
"face")
729 != std::string::npos)
732 ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
738 "Unknown Functional Type: "
739 + this->GetProblem()->GetFunctionalType(),
740 "StatPDEProblem::ComputeReducedFunctionals");
742 this->GetFunctionalValues()[i].push_back(ret);
743 std::stringstream out;
744 this->GetOutputHandler()->InitOut(out);
745 out << this->GetProblem()->GetFunctionalName() <<
": " << ret;
746 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
749 this->GetIntegrator().DeleteDomainData(
"state");
752 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
757 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
758 typename VECTOR,
int dealdim>
759 template<
class DWRC,
class PDE>
762 DWRC& dwrc, PDE& pde)
765 pde.ResidualModifier = boost::bind<void>(
766 boost::mem_fn(&DWRC::ResidualModifier), boost::ref(dwrc), _1);
767 pde.VectorResidualModifier = boost::bind<void>(
768 boost::mem_fn(&DWRC::VectorResidualModifier), boost::ref(dwrc), _1);
771 const unsigned int n_elements =
772 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler().get_tria().n_active_cells();
773 dwrc.ReInit(n_elements);
776 this->ComputeDualForErrorEstimation(dwrc.GetWeightComputation());
779 this->GetOutputHandler()->Write(
"Computing Error Indicators:",
780 4 + this->GetBasePriority());
782 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
785 this->GetIntegrator().AddDomainData(
"state",
786 &(GetU().GetSpacialVector()));
788 this->GetIntegrator().AddDomainData(
"adjoint_for_ee",
789 &(GetZForEE().GetSpacialVector()));
792 this->SetProblemType(
"error_evaluation");
795 dwrc.PrepareWeights(GetU(), GetZForEE());
798 this->GetIntegrator().ComputeRefinementIndicators(*this->GetProblem(),
804 const float error = dwrc.GetError();
807 dwrc.ClearWeightData();
808 this->GetIntegrator().DeleteDomainData(
"state");
810 this->GetIntegrator().DeleteDomainData(
"adjoint_for_ee");
812 this->GetProblem()->DeleteAuxiliaryFromIntegrator(
813 this->GetIntegrator());
815 std::stringstream out;
816 this->GetOutputHandler()->InitOut(out);
817 out <<
"Error estimate using " << dwrc.GetName();
819 out <<
" for the " << this->GetProblem()->GetFunctionalName();
820 out <<
": " << error;
821 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
822 this->GetOutputHandler()->WriteElementwise(dwrc.GetErrorIndicators(),
823 "Error_Indicators" + this->GetPostIndex(),
824 this->GetProblem()->GetDoFType());
830 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
831 typename VECTOR,
int dealdim>
834 const Vector<double> &v, std::string name, std::string outfile,
835 std::string dof_type, std::string filetype)
837 if (dof_type ==
"state")
840 this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
841 data_out.attach_dof_handler(
842 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
844 data_out.add_data_vector(v, name);
845 data_out.build_patches(n_patches_);
847 std::ofstream output(outfile.c_str());
849 if (filetype ==
".vtk")
851 data_out.write_vtk(output);
856 "Don't know how to write filetype `" + filetype +
"'!",
857 "StatPDEProblem::WriteToFileElementwise");
863 throw DOpEException(
"No such DoFHandler `" + dof_type +
"'!",
864 "StatPDEProblem::WriteToFileElementwise");
870 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
871 typename VECTOR,
int dealdim>
874 const VECTOR &v, std::string name, std::string outfile,
875 std::string dof_type, std::string filetype)
877 if (dof_type ==
"state")
880 this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
881 data_out.attach_dof_handler(
882 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
884 data_out.add_data_vector(v, name);
887 data_out.build_patches(
888 this->GetProblem()->GetSpaceTimeHandler()->GetMapping()[0]);
890 std::ofstream output(outfile.c_str());
892 if (filetype ==
".vtk")
894 data_out.write_vtk(output);
896 else if (filetype ==
".gpl")
898 data_out.write_gnuplot(output);
903 "Don't know how to write filetype `" + filetype +
"'!",
904 "StatPDEProblem::WriteToFile");
910 throw DOpEException(
"No such DoFHandler `" + dof_type +
"'!",
911 "StatPDEProblem::WriteToFile");
917 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
918 typename VECTOR,
int dealdim>
923 throw DOpEException(
"This Problem does not support ControlVectors",
"StatPDEProblem::WriteToFile");
StateVector< VECTOR > & GetU()
Definition: statpdeproblem.h:288
void StateSizeInfo(std::stringstream &out)
Definition: statpdeproblem.h:190
void WriteToFileElementwise(const Vector< double > &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statpdeproblem.h:833
Definition: statpdeproblem.h:79
void ComputeReducedFunctionals()
Definition: statpdeproblem.h:614
void WriteToFile(const std::vector< double > &, std::string)
Definition: statpdeproblem.h:251
Definition: parameterreader.h:36
const bool & GetBuildAdjointMatrix() const
Definition: statpdeproblem.h:321
INTEGRATOR & GetIntegrator()
Definition: statpdeproblem.h:309
const StateVector< VECTOR > & GetZForEE() const
Definition: statpdeproblem.h:296
const bool & GetStateReinit() const
Definition: statpdeproblem.h:327
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
Definition: parameterreader.h:98
const bool & GetBuildStateMatrix() const
Definition: statpdeproblem.h:315
NONLINEARSOLVER & GetNonlinearSolver(std::string type)
Definition: statpdeproblem.h:487
Definition: controlvector.h:49
Definition: dopetypes.h:82
Definition: pdeprobleminterface.h:53
StateVector< VECTOR > & GetZForEE()
Definition: statpdeproblem.h:301
bool & GetBuildAdjointMatrix()
Definition: statpdeproblem.h:346
static void declare_params(ParameterReader ¶m_reader)
Definition: statpdeproblem.h:421
Definition: statevector.h:50
bool & GetStateReinit()
Definition: statpdeproblem.h:352
const StateVector< VECTOR > & GetU() const
Definition: statpdeproblem.h:283
bool & GetAdjointReinit()
Definition: statpdeproblem.h:358
void WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statpdeproblem.h:873
virtual ~StatPDEProblem()
Definition: statpdeproblem.h:478
void ComputeRefinementIndicators(DWRC &dwrc, PDE &pde)
Definition: statpdeproblem.h:761
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
VectorStorageType
Definition: dopetypes.h:120
void ComputeDualForErrorEstimation(DOpEtypes::WeightComputation)
Definition: statpdeproblem.h:566
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
void ReInit()
Definition: statpdeproblem.h:511
WeightComputation
Definition: dopetypes.h:80
const std::map< std::string, const VECTOR * > & GetUserDomainData() const
Definition: reducedprobleminterface.h:301
bool & GetBuildStateMatrix()
Definition: statpdeproblem.h:340
StatPDEProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior, ParameterReader ¶m_reader, INTEGRATORDATACONT &idc, int base_priority=0)
Definition: statpdeproblem.h:435
const bool & GetAdjointReinit() const
Definition: statpdeproblem.h:333
void ComputeReducedState()
Definition: statpdeproblem.h:534
Definition: dopeexception.h:35