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);
244 std::string outfile, std::string dof_type, std::string filetype);
320 return _build_state_matrix;
326 return _build_adjoint_matrix;
332 return _state_reinit;
338 return _adjoint_reinit;
345 return _build_state_matrix;
351 return _build_adjoint_matrix;
357 return _state_reinit;
363 return _adjoint_reinit;
396 StateVector<VECTOR> _u;
397 StateVector<VECTOR> _z_for_ee;
399 INTEGRATOR _integrator;
400 NONLINEARSOLVER _nonlinear_state_solver;
401 NONLINEARSOLVER _nonlinear_adjoint_solver;
403 bool _build_state_matrix;
404 bool _build_adjoint_matrix;
406 bool _adjoint_reinit;
410 friend class SolutionExtractor<
411 StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>,
418 using namespace dealii;
421 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
422 typename VECTOR,
int dealdim>
427 NONLINEARSOLVER::declare_params(param_reader);
430 Patterns::Integer(0));
435 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
436 typename VECTOR,
int dealdim>
437 template<
typename INTEGRATORDATACONT>
439 PROBLEM *OP, std::string state_behavior,
443 OP->GetSpaceTimeHandler(), state_behavior, param_reader), _z_for_ee(
444 OP->GetSpaceTimeHandler(), state_behavior, param_reader), _integrator(
445 idc), _nonlinear_state_solver(_integrator, param_reader), _nonlinear_adjoint_solver(
446 _integrator, param_reader)
449 _n_patches = param_reader.
get_integer(
"number of patches");
452 _state_reinit =
true;
453 _adjoint_reinit =
true;
457 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
458 typename VECTOR,
int dealdim>
459 template<
typename INTEGRATORDATACONT>
461 PROBLEM *OP, std::string state_behavior,
463 INTEGRATORDATACONT& idc2,
int base_priority)
465 OP->GetSpaceTimeHandler(), state_behavior, param_reader), _z_for_ee(
466 OP->GetSpaceTimeHandler(), state_behavior, param_reader), _integrator(
467 idc, idc2), _nonlinear_state_solver(_integrator, param_reader), _nonlinear_adjoint_solver(
468 _integrator, param_reader)
472 _state_reinit =
true;
473 _adjoint_reinit =
true;
479 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
480 typename VECTOR,
int dealdim>
487 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
488 typename VECTOR,
int dealdim>
495 return _nonlinear_state_solver;
497 else if (type ==
"adjoint_for_ee")
499 return _nonlinear_adjoint_solver;
503 throw DOpEException(
"No Solver for Problem type:`" + type +
"' found",
504 "StatReducedProblem::GetNonlinearSolver");
511 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
512 typename VECTOR,
int dealdim>
521 _state_reinit =
true;
522 _adjoint_reinit =
true;
525 _build_state_matrix =
true;
526 _build_adjoint_matrix =
true;
529 GetZForEE().ReInit();
534 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
535 typename VECTOR,
int dealdim>
539 this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
540 this->GetOutputHandler()->Write(
"Computing State Solution:",
541 4 + this->GetBasePriority());
543 this->SetProblemType(
"state");
544 auto& problem = this->GetProblem()->GetStateProblem();
545 if (_state_reinit ==
true)
547 GetNonlinearSolver(
"state").ReInit(problem);
548 _state_reinit =
false;
551 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
553 _build_state_matrix = this->GetNonlinearSolver(
"state").NonlinearSolve(
554 problem, (GetU().GetSpacialVector()),
true, _build_state_matrix);
557 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
559 this->GetOutputHandler()->Write((GetU().GetSpacialVector()),
560 "State" + this->GetPostIndex(), problem.GetDoFType());
566 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
567 typename VECTOR,
int dealdim>
572 this->GetOutputHandler()->Write(
"Computing Dual for Error Estimation:",
573 4 + this->GetBasePriority());
576 this->SetProblemType(
"adjoint_for_ee");
581 "StatPDEProblem::ComputeDualForErrorEstimation");
585 auto& problem = *(this->GetProblem());
586 if (_adjoint_reinit ==
true)
588 GetNonlinearSolver(
"adjoint_for_ee").ReInit(problem);
589 _adjoint_reinit =
false;
592 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
594 this->GetIntegrator().AddDomainData(
"state",
595 &(GetU().GetSpacialVector()));
598 _build_adjoint_matrix =
599 this->GetNonlinearSolver(
"adjoint_for_ee").NonlinearSolve(problem,
600 (GetZForEE().GetSpacialVector()),
true, _build_adjoint_matrix);
602 this->GetIntegrator().DeleteDomainData(
"state");
605 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
607 this->GetOutputHandler()->Write((GetZForEE().GetSpacialVector()),
608 "Adjoint_for_ee" + this->GetPostIndex(), problem.GetDoFType());
614 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
615 typename VECTOR,
int dealdim>
619 this->ComputeReducedState();
621 this->GetOutputHandler()->Write(
"Computing Functionals:",
622 4 + this->GetBasePriority());
624 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
626 this->GetIntegrator().AddDomainData(
"state",
627 &(GetU().GetSpacialVector()));
630 for (
unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
635 this->SetProblemType(
"aux_functional", i);
636 if (this->GetProblem()->GetFunctionalType().find(
"domain")
637 != std::string::npos)
640 ret += this->GetIntegrator().ComputeDomainScalar(
641 *(this->GetProblem()));
643 if (this->GetProblem()->GetFunctionalType().find(
"point")
644 != std::string::npos)
647 ret += this->GetIntegrator().ComputePointScalar(
648 *(this->GetProblem()));
650 if (this->GetProblem()->GetFunctionalType().find(
"boundary")
651 != std::string::npos)
654 ret += this->GetIntegrator().ComputeBoundaryScalar(
655 *(this->GetProblem()));
657 if (this->GetProblem()->GetFunctionalType().find(
"face")
658 != std::string::npos)
661 ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
667 "Unknown Functional Type: "
668 + this->GetProblem()->GetFunctionalType(),
669 "StatPDEProblem::ComputeReducedFunctionals");
671 this->GetFunctionalValues()[i].push_back(ret);
672 std::stringstream out;
673 this->GetOutputHandler()->InitOut(out);
674 out << this->GetProblem()->GetFunctionalName() <<
": " << ret;
675 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
678 this->GetIntegrator().DeleteDomainData(
"state");
681 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
686 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
687 typename VECTOR,
int dealdim>
692 this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
695 this->GetOutputHandler()->Write(
"Computing Functionals:",
696 4 + this->GetBasePriority());
698 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
700 this->GetIntegrator().AddDomainData(
"state",
701 &(GetU().GetSpacialVector()));
704 for (
unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
709 this->SetProblemType(
"aux_functional", i);
710 if (this->GetProblem()->GetFunctionalType().find(
"domain")
711 != std::string::npos)
714 ret += this->GetIntegrator().ComputeDomainScalar(
715 *(this->GetProblem()));
717 if (this->GetProblem()->GetFunctionalType().find(
"point")
718 != std::string::npos)
721 ret += this->GetIntegrator().ComputePointScalar(
722 *(this->GetProblem()));
724 if (this->GetProblem()->GetFunctionalType().find(
"boundary")
725 != std::string::npos)
728 ret += this->GetIntegrator().ComputeBoundaryScalar(
729 *(this->GetProblem()));
731 if (this->GetProblem()->GetFunctionalType().find(
"face")
732 != std::string::npos)
735 ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
741 "Unknown Functional Type: "
742 + this->GetProblem()->GetFunctionalType(),
743 "StatPDEProblem::ComputeReducedFunctionals");
745 this->GetFunctionalValues()[i].push_back(ret);
746 std::stringstream out;
747 this->GetOutputHandler()->InitOut(out);
748 out << this->GetProblem()->GetFunctionalName() <<
": " << ret;
749 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
752 this->GetIntegrator().DeleteDomainData(
"state");
755 this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
760 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
761 typename VECTOR,
int dealdim>
762 template<
class DWRC,
class PDE>
765 DWRC& dwrc, PDE& pde)
768 pde.ResidualModifier = boost::bind<void>(
769 boost::mem_fn(&DWRC::ResidualModifier), boost::ref(dwrc), _1);
770 pde.VectorResidualModifier = boost::bind<void>(
771 boost::mem_fn(&DWRC::VectorResidualModifier), boost::ref(dwrc), _1);
774 const unsigned int n_elements =
775 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler().get_tria().n_active_cells();
776 dwrc.ReInit(n_elements);
779 this->ComputeDualForErrorEstimation(dwrc.GetWeightComputation());
782 this->GetOutputHandler()->Write(
"Computing Error Indicators:",
783 4 + this->GetBasePriority());
785 this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
788 this->GetIntegrator().AddDomainData(
"state",
789 &(GetU().GetSpacialVector()));
791 this->GetIntegrator().AddDomainData(
"adjoint_for_ee",
792 &(GetZForEE().GetSpacialVector()));
795 this->SetProblemType(
"error_evaluation");
798 dwrc.PrepareWeights(GetU(), GetZForEE());
801 this->GetIntegrator().ComputeRefinementIndicators(*this->GetProblem(),
807 const float error = dwrc.GetError();
810 dwrc.ClearWeightData();
811 this->GetIntegrator().DeleteDomainData(
"state");
813 this->GetIntegrator().DeleteDomainData(
"adjoint_for_ee");
815 this->GetProblem()->DeleteAuxiliaryFromIntegrator(
816 this->GetIntegrator());
818 std::stringstream out;
819 this->GetOutputHandler()->InitOut(out);
820 out <<
"Error estimate using " << dwrc.GetName();
822 out <<
" for the " << this->GetProblem()->GetFunctionalName();
823 out <<
": " << error;
824 this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
825 this->GetOutputHandler()->WriteElementwise(dwrc.GetErrorIndicators(),
826 "Error_Indicators" + this->GetPostIndex(),
827 this->GetProblem()->GetDoFType());
833 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
834 typename VECTOR,
int dealdim>
837 const Vector<double> &v, std::string name, std::string outfile,
838 std::string dof_type, std::string filetype)
840 if (dof_type ==
"state")
843 this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
844 data_out.attach_dof_handler(
845 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
847 data_out.add_data_vector(v, name);
848 data_out.build_patches(_n_patches);
850 std::ofstream output(outfile.c_str());
852 if (filetype ==
".vtk")
854 data_out.write_vtk(output);
859 "Don't know how to write filetype `" + filetype +
"'!",
860 "StatPDEProblem::WriteToFileElementwise");
866 throw DOpEException(
"No such DoFHandler `" + dof_type +
"'!",
867 "StatPDEProblem::WriteToFileElementwise");
873 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
874 typename VECTOR,
int dealdim>
877 const VECTOR &v, std::string name, std::string outfile,
878 std::string dof_type, std::string filetype)
880 if (dof_type ==
"state")
883 this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
884 data_out.attach_dof_handler(
885 this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
887 data_out.add_data_vector(v, name);
890 data_out.build_patches(
891 this->GetProblem()->GetSpaceTimeHandler()->GetMapping()[0]);
893 std::ofstream output(outfile.c_str());
895 if (filetype ==
".vtk")
897 data_out.write_vtk(output);
899 else if (filetype ==
".gpl")
901 data_out.write_gnuplot(output);
906 "Don't know how to write filetype `" + filetype +
"'!",
907 "StatPDEProblem::WriteToFile");
913 throw DOpEException(
"No such DoFHandler `" + dof_type +
"'!",
914 "StatPDEProblem::WriteToFile");
920 template<
typename NONLINEARSOLVER,
typename INTEGRATOR,
typename PROBLEM,
921 typename VECTOR,
int dealdim>
925 std::string dof_type, std::string filetype)
StateVector< VECTOR > & GetU()
Definition: statpdeproblem.h:291
void StateSizeInfo(std::stringstream &out)
Definition: statpdeproblem.h:190
StatPDEProblem(PROBLEM *OP, std::string state_behavior, ParameterReader ¶m_reader, INTEGRATORDATACONT &idc, int base_priority=0)
Definition: statpdeproblem.h:438
void WriteToFileElementwise(const Vector< double > &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statpdeproblem.h:836
Definition: statpdeproblem.h:79
void ComputeReducedFunctionals()
Definition: statpdeproblem.h:617
void WriteToFile(const std::vector< double > &, std::string)
Definition: statpdeproblem.h:254
Definition: parameterreader.h:36
const bool & GetBuildAdjointMatrix() const
Definition: statpdeproblem.h:324
INTEGRATOR & GetIntegrator()
Definition: statpdeproblem.h:312
const StateVector< VECTOR > & GetZForEE() const
Definition: statpdeproblem.h:299
const bool & GetStateReinit() const
Definition: statpdeproblem.h:330
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:318
NONLINEARSOLVER & GetNonlinearSolver(std::string type)
Definition: statpdeproblem.h:490
Definition: controlvector.h:48
Definition: dopetypes.h:80
Definition: pdeprobleminterface.h:53
StateVector< VECTOR > & GetZForEE()
Definition: statpdeproblem.h:304
bool & GetBuildAdjointMatrix()
Definition: statpdeproblem.h:349
static void declare_params(ParameterReader ¶m_reader)
Definition: statpdeproblem.h:424
Definition: statevector.h:49
VECTOR & GetSpacialVector()
Definition: controlvector.cc:118
bool & GetStateReinit()
Definition: statpdeproblem.h:355
const StateVector< VECTOR > & GetU() const
Definition: statpdeproblem.h:286
bool & GetAdjointReinit()
Definition: statpdeproblem.h:361
void WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statpdeproblem.h:876
virtual ~StatPDEProblem()
Definition: statpdeproblem.h:481
void ComputeRefinementIndicators(DWRC &dwrc, PDE &pde)
Definition: statpdeproblem.h:764
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
void ComputeDualForErrorEstimation(DOpEtypes::WeightComputation)
Definition: statpdeproblem.h:569
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
void ReInit()
Definition: statpdeproblem.h:514
WeightComputation
Definition: dopetypes.h:78
const std::map< std::string, const VECTOR * > & GetUserDomainData() const
Definition: reducedprobleminterface.h:303
bool & GetBuildStateMatrix()
Definition: statpdeproblem.h:343
const bool & GetAdjointReinit() const
Definition: statpdeproblem.h:336
void ComputeReducedState()
Definition: statpdeproblem.h:537
Definition: dopeexception.h:35