diff --git a/applications/precipitateEvolution_pfunction/PLibrary/PLibrary.hh b/applications/precipitateEvolution_pfunction/PLibrary/PLibrary.hh index 4dbfcfee3..0d4a6c4cb 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/PLibrary.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/PLibrary.hh @@ -6,37 +6,36 @@ #ifndef PLIBRARY_HH #define PLIBRARY_HH -#include #include "IntegrationTools/PFunction.hh" #include "IntegrationTools/PPieceWise.hh" +#include + namespace PRISMS { - /// Library where you can find functions and basis sets - /// - namespace PLibrary - { - // Use these functions to checkout objects which manage their own memory - - void checkout( std::string name, PSimpleFunction< double*, double > &simplefunc); - - void checkout( std::string name, PFunction< double*, double > &func); - - - - - // Use these functions to checkout new 'Base' objects which the user must delete + /// Library where you can find functions and basis sets + /// + namespace PLibrary + { + // Use these functions to checkout objects which manage their own memory - void checkout( std::string name, PSimpleBase< double*, double > *&simplefunc); + void + checkout(std::string name, PSimpleFunction &simplefunc); - void checkout( std::string name, PFuncBase< double*, double > *&func); + void + checkout(std::string name, PFunction &func); + // Use these functions to checkout new 'Base' objects which the user must delete + void + checkout(std::string name, PSimpleBase *&simplefunc); - } + void + checkout(std::string name, PFuncBase *&func); -} + } // namespace PLibrary +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh index d267db05b..1e95e3061 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh @@ -6,128 +6,136 @@ #ifndef pfunct_McV_HH #define pfunct_McV_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_McV_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+00; - } - - public: - - pfunct_McV_f() - { - this->_name = "pfunct_McV_f"; - } - - std::string csrc() const - { - return "1.0000000000000000e+00"; - } - - std::string sym() const - { - return "1.0"; - } - - std::string latex() const - { - return "1.0"; - } - - pfunct_McV_f* clone() const - { - return new pfunct_McV_f(*this); - } - }; - - template - class pfunct_McV : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_McV() - { - construct(); - } - - pfunct_McV(const pfunct_McV &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - - } - - pfunct_McV& operator=( pfunct_McV RHS ) - { - using std::swap; - - swap(_val, RHS._val); - - return *this; - } - - ~pfunct_McV() - { - delete _val; - - } - - pfunct_McV* clone() const - { - return new pfunct_McV(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - double operator()() const - { - return (*_val)(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_McV"; - this->_var_name.clear(); - this->_var_name.push_back("c"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - if(!allocate) return; - - _val = new pfunct_McV_f(); - } - - }; - - -} + template + class pfunct_McV_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+00; + } + + public: + pfunct_McV_f() + { + this->_name = "pfunct_McV_f"; + } + + std::string + csrc() const + { + return "1.0000000000000000e+00"; + } + + std::string + sym() const + { + return "1.0"; + } + + std::string + latex() const + { + return "1.0"; + } + + pfunct_McV_f * + clone() const + { + return new pfunct_McV_f(*this); + } + }; + + template + class pfunct_McV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_McV() + { + construct(); + } + + pfunct_McV(const pfunct_McV &RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_McV & + operator=(pfunct_McV RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_McV() + { + delete _val; + } + + pfunct_McV * + clone() const + { + return new pfunct_McV(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + double + operator()() const + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_McV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_McV_f(); + } + }; + +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh index eea85cc73..c39a38a2f 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh @@ -6,128 +6,136 @@ #ifndef pfunct_Mn1V_HH #define pfunct_Mn1V_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_Mn1V_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+02; - } - - public: - - pfunct_Mn1V_f() - { - this->_name = "pfunct_Mn1V_f"; - } - - std::string csrc() const - { - return "1.0000000000000000e+02"; - } - - std::string sym() const - { - return "100.0"; - } - - std::string latex() const - { - return "100.0"; - } - - pfunct_Mn1V_f* clone() const - { - return new pfunct_Mn1V_f(*this); - } - }; - - template - class pfunct_Mn1V : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_Mn1V() - { - construct(); - } - - pfunct_Mn1V(const pfunct_Mn1V &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - - } - - pfunct_Mn1V& operator=( pfunct_Mn1V RHS ) - { - using std::swap; - - swap(_val, RHS._val); - - return *this; - } - - ~pfunct_Mn1V() - { - delete _val; - - } - - pfunct_Mn1V* clone() const - { - return new pfunct_Mn1V(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - double operator()() const - { - return (*_val)(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_Mn1V"; - this->_var_name.clear(); - this->_var_name.push_back("n1"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - if(!allocate) return; - - _val = new pfunct_Mn1V_f(); - } - - }; - - -} + template + class pfunct_Mn1V_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn1V_f() + { + this->_name = "pfunct_Mn1V_f"; + } + + std::string + csrc() const + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const + { + return "100.0"; + } + + std::string + latex() const + { + return "100.0"; + } + + pfunct_Mn1V_f * + clone() const + { + return new pfunct_Mn1V_f(*this); + } + }; + + template + class pfunct_Mn1V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn1V() + { + construct(); + } + + pfunct_Mn1V(const pfunct_Mn1V &RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn1V & + operator=(pfunct_Mn1V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn1V() + { + delete _val; + } + + pfunct_Mn1V * + clone() const + { + return new pfunct_Mn1V(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + double + operator()() const + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn1V"; + this->_var_name.clear(); + this->_var_name.push_back("n1"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn1V_f(); + } + }; + +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh index 7d0d2d531..a4a596b84 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh @@ -6,128 +6,136 @@ #ifndef pfunct_Mn2V_HH #define pfunct_Mn2V_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_Mn2V_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+02; - } - - public: - - pfunct_Mn2V_f() - { - this->_name = "pfunct_Mn2V_f"; - } - - std::string csrc() const - { - return "1.0000000000000000e+02"; - } - - std::string sym() const - { - return "100.0"; - } - - std::string latex() const - { - return "100.0"; - } - - pfunct_Mn2V_f* clone() const - { - return new pfunct_Mn2V_f(*this); - } - }; - - template - class pfunct_Mn2V : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_Mn2V() - { - construct(); - } - - pfunct_Mn2V(const pfunct_Mn2V &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - - } - - pfunct_Mn2V& operator=( pfunct_Mn2V RHS ) - { - using std::swap; - - swap(_val, RHS._val); - - return *this; - } - - ~pfunct_Mn2V() - { - delete _val; - - } - - pfunct_Mn2V* clone() const - { - return new pfunct_Mn2V(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - double operator()() const - { - return (*_val)(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_Mn2V"; - this->_var_name.clear(); - this->_var_name.push_back("n2"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - if(!allocate) return; - - _val = new pfunct_Mn2V_f(); - } - - }; - - -} + template + class pfunct_Mn2V_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn2V_f() + { + this->_name = "pfunct_Mn2V_f"; + } + + std::string + csrc() const + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const + { + return "100.0"; + } + + std::string + latex() const + { + return "100.0"; + } + + pfunct_Mn2V_f * + clone() const + { + return new pfunct_Mn2V_f(*this); + } + }; + + template + class pfunct_Mn2V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn2V() + { + construct(); + } + + pfunct_Mn2V(const pfunct_Mn2V &RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn2V & + operator=(pfunct_Mn2V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn2V() + { + delete _val; + } + + pfunct_Mn2V * + clone() const + { + return new pfunct_Mn2V(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + double + operator()() const + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn2V"; + this->_var_name.clear(); + this->_var_name.push_back("n2"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn2V_f(); + } + }; + +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh index 6b2b87f23..a785980b2 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh @@ -6,128 +6,136 @@ #ifndef pfunct_Mn3V_HH #define pfunct_Mn3V_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_Mn3V_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+02; - } - - public: - - pfunct_Mn3V_f() - { - this->_name = "pfunct_Mn3V_f"; - } - - std::string csrc() const - { - return "1.0000000000000000e+02"; - } - - std::string sym() const - { - return "100.0"; - } - - std::string latex() const - { - return "100.0"; - } - - pfunct_Mn3V_f* clone() const - { - return new pfunct_Mn3V_f(*this); - } - }; - - template - class pfunct_Mn3V : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_Mn3V() - { - construct(); - } - - pfunct_Mn3V(const pfunct_Mn3V &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - - } - - pfunct_Mn3V& operator=( pfunct_Mn3V RHS ) - { - using std::swap; - - swap(_val, RHS._val); - - return *this; - } - - ~pfunct_Mn3V() - { - delete _val; - - } - - pfunct_Mn3V* clone() const - { - return new pfunct_Mn3V(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - double operator()() const - { - return (*_val)(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_Mn3V"; - this->_var_name.clear(); - this->_var_name.push_back("n3"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - if(!allocate) return; - - _val = new pfunct_Mn3V_f(); - } - - }; - - -} + template + class pfunct_Mn3V_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn3V_f() + { + this->_name = "pfunct_Mn3V_f"; + } + + std::string + csrc() const + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const + { + return "100.0"; + } + + std::string + latex() const + { + return "100.0"; + } + + pfunct_Mn3V_f * + clone() const + { + return new pfunct_Mn3V_f(*this); + } + }; + + template + class pfunct_Mn3V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn3V() + { + construct(); + } + + pfunct_Mn3V(const pfunct_Mn3V &RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn3V & + operator=(pfunct_Mn3V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn3V() + { + delete _val; + } + + pfunct_Mn3V * + clone() const + { + return new pfunct_Mn3V(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + double + operator()() const + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn3V"; + this->_var_name.clear(); + this->_var_name.push_back("n3"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn3V_f(); + } + }; + +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh index 57557931b..9946e3d2a 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh @@ -6,259 +6,295 @@ #ifndef pfunct_faV_HH #define pfunct_faV_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_faV_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 5.1622000000000003e+00*(var[0]*var[0])+-2.7374999999999998e+00*(var[0]*var[0]*var[0])+-4.7759999999999998e+00*var[0]+1.3687000000000000e+00*((var[0]*var[0])*(var[0]*var[0]))-1.6704000000000001e+00; - } - - public: - - pfunct_faV_f() - { - this->_name = "pfunct_faV_f"; - } - - std::string csrc() const - { - return " 5.1622000000000003e+00*(var[0]*var[0])+-2.7374999999999998e+00*(var[0]*var[0]*var[0])+-4.7759999999999998e+00*var[0]+1.3687000000000000e+00*((var[0]*var[0])*(var[0]*var[0]))-1.6704000000000001e+00"; - } - - std::string sym() const - { - return "-1.6704+(1.3687)*c^4+(5.1622)*c^2-(2.7375)*c^3-(4.776)*c"; - } - - std::string latex() const - { - return "-1.6704+{(5.1622)} c^{2}-{(4.776)} c-{(2.7375)} c^{3}+{(1.3687)} c^{4}"; - } - - pfunct_faV_f* clone() const - { - return new pfunct_faV_f(*this); - } - }; - - template< class VarContainer> - class pfunct_faV_grad_0 : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 5.4748000000000001e+00*(var[0]*var[0]*var[0])+-8.2125000000000004e+00*(var[0]*var[0])+1.0324400000000001e+01*var[0]-4.7759999999999998e+00; - } - - public: - - pfunct_faV_grad_0() - { - this->_name = "pfunct_faV_grad_0"; - } - - std::string csrc() const - { - return " 5.4748000000000001e+00*(var[0]*var[0]*var[0])+-8.2125000000000004e+00*(var[0]*var[0])+1.0324400000000001e+01*var[0]-4.7759999999999998e+00"; - } - - std::string sym() const - { - return "-4.776+(10.3244)*c+(5.4748)*c^3-(8.2125)*c^2"; - } - - std::string latex() const - { - return "-4.776+{(5.4748)} c^{3}-{(8.2125)} c^{2}+{(10.3244)} c"; - } - - pfunct_faV_grad_0* clone() const - { - return new pfunct_faV_grad_0(*this); - } - }; - - template< class VarContainer> - class pfunct_faV_hess_0_0 : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.6424399999999999e+01*(var[0]*var[0])+-1.6425000000000001e+01*var[0]+1.0324400000000001e+01; - } - - public: - - pfunct_faV_hess_0_0() - { - this->_name = "pfunct_faV_hess_0_0"; - } - - std::string csrc() const - { - return " 1.6424399999999999e+01*(var[0]*var[0])+-1.6425000000000001e+01*var[0]+1.0324400000000001e+01"; - } - - std::string sym() const - { - return "10.3244-(16.425)*c+(16.4244)*c^2"; - } - - std::string latex() const - { - return "10.3244+{(16.4244)} c^{2}-{(16.425)} c"; - } - - pfunct_faV_hess_0_0* clone() const - { - return new pfunct_faV_hess_0_0(*this); - } - }; - - template - class pfunct_faV : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_faV() - { - construct(); - } - - pfunct_faV(const pfunct_faV &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - _grad_val[0] = RHS._grad_val[0]->clone(); - _hess_val[0][0] = RHS._hess_val[0][0]->clone(); - - } - - pfunct_faV& operator=( pfunct_faV RHS ) - { - using std::swap; - - swap(_val, RHS._val); - swap(_grad_val[0], RHS._grad_val[0]); - swap(_hess_val[0][0], RHS._hess_val[0][0]); - - return *this; - } - - ~pfunct_faV() - { - delete _val; - - delete _grad_val[0]; - delete [] _grad_val; - - delete _hess_val[0][0]; - delete [] _hess_val[0]; - delete [] _hess_val; - } - - pfunct_faV* clone() const - { - return new pfunct_faV(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - PSimpleFunction< VarContainer, double> grad_simplefunction(size_type di) const - { - return PSimpleFunction< VarContainer, double>( *_grad_val[di] ); - } - - PSimpleFunction< VarContainer, double> hess_simplefunction(size_type di, size_type dj) const - { - return PSimpleFunction< VarContainer, double>( *_hess_val[di][dj] ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - double grad(const VarContainer &var, size_type di) - { - return (*_grad_val[di])(var); - } - - double hess(const VarContainer &var, size_type di, size_type dj) - { - return (*_hess_val[di][dj])(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - void eval_grad(const VarContainer &var) - { - (*_grad_val[0])(var); - } - - void eval_hess(const VarContainer &var) - { - (*_hess_val[0][0])(var); - } - - double operator()() const - { - return (*_val)(); - } - - double grad(size_type di) const - { - return (*_grad_val[di])(); - } - - double hess(size_type di, size_type dj) const - { - return (*_hess_val[di][dj])(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_faV"; - this->_var_name.clear(); - this->_var_name.push_back("c"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - _grad_val = new PSimpleBase< VarContainer, double>*[1]; - - _hess_val = new PSimpleBase< VarContainer, double>**[1]; - _hess_val[0] = new PSimpleBase< VarContainer, double>*[1]; - - if(!allocate) return; - - _val = new pfunct_faV_f(); - - _grad_val[0] = new pfunct_faV_grad_0(); - - _hess_val[0][0] = new pfunct_faV_hess_0_0(); - } - - }; - - -} + template + class pfunct_faV_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 5.1622000000000003e+00 * (var[0] * var[0]) + + -2.7374999999999998e+00 * (var[0] * var[0] * var[0]) + + -4.7759999999999998e+00 * var[0] + + 1.3687000000000000e+00 * ((var[0] * var[0]) * (var[0] * var[0])) - + 1.6704000000000001e+00; + } + + public: + pfunct_faV_f() + { + this->_name = "pfunct_faV_f"; + } + + std::string + csrc() const + { + return " 5.1622000000000003e+00*(var[0]*var[0])+-2.7374999999999998e+00*(var[0]*" + "var[0]*var[0])+-4.7759999999999998e+00*var[0]+1.3687000000000000e+00*((var[" + "0]*var[0])*(var[0]*var[0]))-1.6704000000000001e+00"; + } + + std::string + sym() const + { + return "-1.6704+(1.3687)*c^4+(5.1622)*c^2-(2.7375)*c^3-(4.776)*c"; + } + + std::string + latex() const + { + return "-1.6704+{(5.1622)} c^{2}-{(4.776)} c-{(2.7375)} c^{3}+{(1.3687)} c^{4}"; + } + + pfunct_faV_f * + clone() const + { + return new pfunct_faV_f(*this); + } + }; + + template + class pfunct_faV_grad_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 5.4748000000000001e+00 * (var[0] * var[0] * var[0]) + + -8.2125000000000004e+00 * (var[0] * var[0]) + + 1.0324400000000001e+01 * var[0] - 4.7759999999999998e+00; + } + + public: + pfunct_faV_grad_0() + { + this->_name = "pfunct_faV_grad_0"; + } + + std::string + csrc() const + { + return " 5.4748000000000001e+00*(var[0]*var[0]*var[0])+-8.2125000000000004e+00*(" + "var[0]*var[0])+1.0324400000000001e+01*var[0]-4.7759999999999998e+00"; + } + + std::string + sym() const + { + return "-4.776+(10.3244)*c+(5.4748)*c^3-(8.2125)*c^2"; + } + + std::string + latex() const + { + return "-4.776+{(5.4748)} c^{3}-{(8.2125)} c^{2}+{(10.3244)} c"; + } + + pfunct_faV_grad_0 * + clone() const + { + return new pfunct_faV_grad_0(*this); + } + }; + + template + class pfunct_faV_hess_0_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.6424399999999999e+01 * (var[0] * var[0]) + + -1.6425000000000001e+01 * var[0] + 1.0324400000000001e+01; + } + + public: + pfunct_faV_hess_0_0() + { + this->_name = "pfunct_faV_hess_0_0"; + } + + std::string + csrc() const + { + return " 1.6424399999999999e+01*(var[0]*var[0])+-1.6425000000000001e+01*var[0]+1." + "0324400000000001e+01"; + } + + std::string + sym() const + { + return "10.3244-(16.425)*c+(16.4244)*c^2"; + } + + std::string + latex() const + { + return "10.3244+{(16.4244)} c^{2}-{(16.425)} c"; + } + + pfunct_faV_hess_0_0 * + clone() const + { + return new pfunct_faV_hess_0_0(*this); + } + }; + + template + class pfunct_faV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_faV() + { + construct(); + } + + pfunct_faV(const pfunct_faV &RHS) + { + construct(false); + + _val = RHS._val->clone(); + _grad_val[0] = RHS._grad_val[0]->clone(); + _hess_val[0][0] = RHS._hess_val[0][0]->clone(); + } + + pfunct_faV & + operator=(pfunct_faV RHS) + { + using std::swap; + + swap(_val, RHS._val); + swap(_grad_val[0], RHS._grad_val[0]); + swap(_hess_val[0][0], RHS._hess_val[0][0]); + + return *this; + } + + ~pfunct_faV() + { + delete _val; + + delete _grad_val[0]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete[] _hess_val[0]; + delete[] _hess_val; + } + + pfunct_faV * + clone() const + { + return new pfunct_faV(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + PSimpleFunction + grad_simplefunction(size_type di) const + { + return PSimpleFunction(*_grad_val[di]); + } + + PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + return PSimpleFunction(*_hess_val[di][dj]); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + double + grad(const VarContainer &var, size_type di) + { + return (*_grad_val[di])(var); + } + + double + hess(const VarContainer &var, size_type di, size_type dj) + { + return (*_hess_val[di][dj])(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + void + eval_grad(const VarContainer &var) + { + (*_grad_val[0])(var); + } + + void + eval_hess(const VarContainer &var) + { + (*_hess_val[0][0])(var); + } + + double + operator()() const + { + return (*_val)(); + } + + double + grad(size_type di) const + { + return (*_grad_val[di])(); + } + + double + hess(size_type di, size_type dj) const + { + return (*_hess_val[di][dj])(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_faV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + _grad_val = new PSimpleBase *[1]; + + _hess_val = new PSimpleBase **[1]; + _hess_val[0] = new PSimpleBase *[1]; + + if (!allocate) + return; + + _val = new pfunct_faV_f(); + + _grad_val[0] = new pfunct_faV_grad_0(); + + _hess_val[0][0] = new pfunct_faV_hess_0_0(); + } + }; + +} // namespace PRISMS #endif diff --git a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh index a9bafc7f3..ce40da99c 100644 --- a/applications/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh +++ b/applications/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh @@ -6,259 +6,286 @@ #ifndef pfunct_fbV_HH #define pfunct_fbV_HH +#include "IntegrationTools/PFunction.hh" + #include #include -#include "IntegrationTools/PFunction.hh" namespace PRISMS { - template< class VarContainer> - class pfunct_fbV_f : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return -5.9745999999999997e+00*var[0]+5.0000000000000000e+00*(var[0]*var[0])-1.5924000000000000e+00; - } - - public: - - pfunct_fbV_f() - { - this->_name = "pfunct_fbV_f"; - } - - std::string csrc() const - { - return " -5.9745999999999997e+00*var[0]+5.0000000000000000e+00*(var[0]*var[0])-1.5924000000000000e+00"; - } - - std::string sym() const - { - return "-1.5924-(5.9746)*c+(5.0)*c^2"; - } - - std::string latex() const - { - return "-1.5924+{(5.0)} c^{2}-{(5.9746)} c"; - } - - pfunct_fbV_f* clone() const - { - return new pfunct_fbV_f(*this); - } - }; - - template< class VarContainer> - class pfunct_fbV_grad_0 : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+01*var[0]-5.9745999999999997e+00; - } - - public: - - pfunct_fbV_grad_0() - { - this->_name = "pfunct_fbV_grad_0"; - } - - std::string csrc() const - { - return " 1.0000000000000000e+01*var[0]-5.9745999999999997e+00"; - } - - std::string sym() const - { - return "-5.9746+(10.0)*c"; - } - - std::string latex() const - { - return "-5.9746+{(10.0)} c"; - } - - pfunct_fbV_grad_0* clone() const - { - return new pfunct_fbV_grad_0(*this); - } - }; - - template< class VarContainer> - class pfunct_fbV_hess_0_0 : public PSimpleBase< VarContainer, double> - { - double eval( const VarContainer &var) const - { - return 1.0000000000000000e+01; - } - - public: - - pfunct_fbV_hess_0_0() - { - this->_name = "pfunct_fbV_hess_0_0"; - } - - std::string csrc() const - { - return "1.0000000000000000e+01"; - } - - std::string sym() const - { - return "10.0"; - } - - std::string latex() const - { - return "10.0"; - } - - pfunct_fbV_hess_0_0* clone() const - { - return new pfunct_fbV_hess_0_0(*this); - } - }; - - template - class pfunct_fbV : public PFuncBase< VarContainer, double> - { - public: - - typedef typename PFuncBase< VarContainer, double>::size_type size_type; - - PSimpleBase< VarContainer, double> *_val; - PSimpleBase< VarContainer, double> **_grad_val; - PSimpleBase< VarContainer, double> ***_hess_val; - - pfunct_fbV() - { - construct(); - } - - pfunct_fbV(const pfunct_fbV &RHS ) - { - construct(false); - - _val = RHS._val->clone(); - _grad_val[0] = RHS._grad_val[0]->clone(); - _hess_val[0][0] = RHS._hess_val[0][0]->clone(); - - } - - pfunct_fbV& operator=( pfunct_fbV RHS ) - { - using std::swap; - - swap(_val, RHS._val); - swap(_grad_val[0], RHS._grad_val[0]); - swap(_hess_val[0][0], RHS._hess_val[0][0]); - - return *this; - } - - ~pfunct_fbV() - { - delete _val; - - delete _grad_val[0]; - delete [] _grad_val; - - delete _hess_val[0][0]; - delete [] _hess_val[0]; - delete [] _hess_val; - } - - pfunct_fbV* clone() const - { - return new pfunct_fbV(*this); - } - - PSimpleFunction< VarContainer, double> simplefunction() const - { - return PSimpleFunction< VarContainer, double>( *_val ); - } - - PSimpleFunction< VarContainer, double> grad_simplefunction(size_type di) const - { - return PSimpleFunction< VarContainer, double>( *_grad_val[di] ); - } - - PSimpleFunction< VarContainer, double> hess_simplefunction(size_type di, size_type dj) const - { - return PSimpleFunction< VarContainer, double>( *_hess_val[di][dj] ); - } - - double operator()(const VarContainer &var) - { - return (*_val)(var); - } - - double grad(const VarContainer &var, size_type di) - { - return (*_grad_val[di])(var); - } - - double hess(const VarContainer &var, size_type di, size_type dj) - { - return (*_hess_val[di][dj])(var); - } - - void eval(const VarContainer &var) - { - (*_val)(var); - } - - void eval_grad(const VarContainer &var) - { - (*_grad_val[0])(var); - } - - void eval_hess(const VarContainer &var) - { - (*_hess_val[0][0])(var); - } - - double operator()() const - { - return (*_val)(); - } - - double grad(size_type di) const - { - return (*_grad_val[di])(); - } - - double hess(size_type di, size_type dj) const - { - return (*_hess_val[di][dj])(); - } - - private: - void construct(bool allocate = true) - { - this->_name = "pfunct_fbV"; - this->_var_name.clear(); - this->_var_name.push_back("c"); - this->_var_description.clear(); - this->_var_description.push_back("concentration"); - - _grad_val = new PSimpleBase< VarContainer, double>*[1]; - - _hess_val = new PSimpleBase< VarContainer, double>**[1]; - _hess_val[0] = new PSimpleBase< VarContainer, double>*[1]; - - if(!allocate) return; - - _val = new pfunct_fbV_f(); - - _grad_val[0] = new pfunct_fbV_grad_0(); - - _hess_val[0][0] = new pfunct_fbV_hess_0_0(); - } - - }; - - -} + template + class pfunct_fbV_f : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return -5.9745999999999997e+00 * var[0] + + 5.0000000000000000e+00 * (var[0] * var[0]) - 1.5924000000000000e+00; + } + + public: + pfunct_fbV_f() + { + this->_name = "pfunct_fbV_f"; + } + + std::string + csrc() const + { + return " -5.9745999999999997e+00*var[0]+5.0000000000000000e+00*(var[0]*var[0])-1." + "5924000000000000e+00"; + } + + std::string + sym() const + { + return "-1.5924-(5.9746)*c+(5.0)*c^2"; + } + + std::string + latex() const + { + return "-1.5924+{(5.0)} c^{2}-{(5.9746)} c"; + } + + pfunct_fbV_f * + clone() const + { + return new pfunct_fbV_f(*this); + } + }; + + template + class pfunct_fbV_grad_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+01 * var[0] - 5.9745999999999997e+00; + } + + public: + pfunct_fbV_grad_0() + { + this->_name = "pfunct_fbV_grad_0"; + } + + std::string + csrc() const + { + return " 1.0000000000000000e+01*var[0]-5.9745999999999997e+00"; + } + + std::string + sym() const + { + return "-5.9746+(10.0)*c"; + } + + std::string + latex() const + { + return "-5.9746+{(10.0)} c"; + } + + pfunct_fbV_grad_0 * + clone() const + { + return new pfunct_fbV_grad_0(*this); + } + }; + + template + class pfunct_fbV_hess_0_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const + { + return 1.0000000000000000e+01; + } + + public: + pfunct_fbV_hess_0_0() + { + this->_name = "pfunct_fbV_hess_0_0"; + } + + std::string + csrc() const + { + return "1.0000000000000000e+01"; + } + + std::string + sym() const + { + return "10.0"; + } + + std::string + latex() const + { + return "10.0"; + } + + pfunct_fbV_hess_0_0 * + clone() const + { + return new pfunct_fbV_hess_0_0(*this); + } + }; + + template + class pfunct_fbV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_fbV() + { + construct(); + } + + pfunct_fbV(const pfunct_fbV &RHS) + { + construct(false); + + _val = RHS._val->clone(); + _grad_val[0] = RHS._grad_val[0]->clone(); + _hess_val[0][0] = RHS._hess_val[0][0]->clone(); + } + + pfunct_fbV & + operator=(pfunct_fbV RHS) + { + using std::swap; + + swap(_val, RHS._val); + swap(_grad_val[0], RHS._grad_val[0]); + swap(_hess_val[0][0], RHS._hess_val[0][0]); + + return *this; + } + + ~pfunct_fbV() + { + delete _val; + + delete _grad_val[0]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete[] _hess_val[0]; + delete[] _hess_val; + } + + pfunct_fbV * + clone() const + { + return new pfunct_fbV(*this); + } + + PSimpleFunction + simplefunction() const + { + return PSimpleFunction(*_val); + } + + PSimpleFunction + grad_simplefunction(size_type di) const + { + return PSimpleFunction(*_grad_val[di]); + } + + PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + return PSimpleFunction(*_hess_val[di][dj]); + } + + double + operator()(const VarContainer &var) + { + return (*_val)(var); + } + + double + grad(const VarContainer &var, size_type di) + { + return (*_grad_val[di])(var); + } + + double + hess(const VarContainer &var, size_type di, size_type dj) + { + return (*_hess_val[di][dj])(var); + } + + void + eval(const VarContainer &var) + { + (*_val)(var); + } + + void + eval_grad(const VarContainer &var) + { + (*_grad_val[0])(var); + } + + void + eval_hess(const VarContainer &var) + { + (*_hess_val[0][0])(var); + } + + double + operator()() const + { + return (*_val)(); + } + + double + grad(size_type di) const + { + return (*_grad_val[di])(); + } + + double + hess(size_type di, size_type dj) const + { + return (*_hess_val[di][dj])(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_fbV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + _grad_val = new PSimpleBase *[1]; + + _hess_val = new PSimpleBase **[1]; + _hess_val[0] = new PSimpleBase *[1]; + + if (!allocate) + return; + + _val = new pfunct_fbV_f(); + + _grad_val[0] = new pfunct_fbV_grad_0(); + + _hess_val[0][0] = new pfunct_fbV_hess_0_0(); + } + }; + +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/PPieceWise.hh b/include/IntegrationTools/PPieceWise.hh index 3ba618892..23c7288be 100644 --- a/include/IntegrationTools/PPieceWise.hh +++ b/include/IntegrationTools/PPieceWise.hh @@ -1,4 +1,4 @@ -#include "piecewise/SimplePiece.hh" +#include "piecewise/PPieceWiseFuncBase.hh" #include "piecewise/PPieceWiseSimpleBase.hh" #include "piecewise/Piece.hh" -#include "piecewise/PPieceWiseFuncBase.hh" +#include "piecewise/SimplePiece.hh" diff --git a/include/IntegrationTools/datastruc/Bin.hh b/include/IntegrationTools/datastruc/Bin.hh index a595474a5..06a5e4d7b 100644 --- a/include/IntegrationTools/datastruc/Bin.hh +++ b/include/IntegrationTools/datastruc/Bin.hh @@ -2,226 +2,239 @@ #ifndef Bin_HH #define Bin_HH +#include "PNDArray.hh" + #include #include -#include "PNDArray.hh" namespace PRISMS { - - /// A class for binning things - /// template class T: the type of thing being binned - /// - should have T::operator==() defined to use 'Bin::add_once' - /// template class Coordinate: - /// - expected to be 'vector like' - /// - should have Coordinate::operator[]() defined + + /// A class for binning things + /// template class T: the type of thing being binned + /// - should have T::operator==() defined to use 'Bin::add_once' + /// template class Coordinate: + /// - expected to be 'vector like' + /// - should have Coordinate::operator[]() defined + /// + template + class Bin + { + PNDArray> _item; + + std::vector _min; // min coord + std::vector _incr; // histogram spacing along each direction (uniform along + // each dimension) + std::vector _N; // number of bins along each direction + std::vector _max; // max coord + std::vector _indices; + + public: + Bin() {}; + + /// Construct a Bin + /// 'min': minimum value of each coordinate component + /// 'incr': bin spacing along each direction (this is uniform along each dimension) + /// 'N': number of bins along each direction /// - template< class T, class Coordinate> - class Bin + /// For example, to bin the range (0->10, 0->20, 10->100), with size 1 bin spacing: + /// Bin( {0,0,10}, {1,1,1}, {10, 20, 90}) + /// + Bin(const std::vector &min, std::vector &incr, std::vector &N) { - PNDArray > _item; - - std::vector _min; // min coord - std::vector _incr; // histogram spacing along each direction (uniform along each dimension) - std::vector _N; // number of bins along each direction - std::vector _max; // max coord - std::vector _indices; - - public: - - Bin(){}; - - /// Construct a Bin - /// 'min': minimum value of each coordinate component - /// 'incr': bin spacing along each direction (this is uniform along each dimension) - /// 'N': number of bins along each direction - /// - /// For example, to bin the range (0->10, 0->20, 10->100), with size 1 bin spacing: - /// Bin( {0,0,10}, {1,1,1}, {10, 20, 90}) - /// - Bin(const std::vector &min, std::vector &incr, std::vector &N) - { - _min = min; - _incr = incr; - _N = N; - _max = std::vector(_min.size()); - for( unsigned int i=0; i<_min.size(); i++) - _max[i] = _min[i] + _incr[i]*_N[i]; - _indices = _N; - _item.resize(_N); - } - - void clear() - { - (*this) = Bin(); - } - - std::vector& min() - { - return _min; - } - - std::vector& max() - { - return _max; - } - - /// Add a new item to the bin containing a given coordinate - void add( const T &newitem, const Coordinate &coord) - { - indices( coord, _indices); - _item(_indices).push_back(newitem); - } - - /// Add a new item to the bin containing a given coordinate, - /// if an equivalent item is not already present - /// return 'true' if added succesfully, 'false' if not - void add_once( const T &newitem, const Coordinate &coord) + _min = min; + _incr = incr; + _N = N; + _max = std::vector(_min.size()); + for (unsigned int i = 0; i < _min.size(); i++) + _max[i] = _min[i] + _incr[i] * _N[i]; + _indices = _N; + _item.resize(_N); + } + + void + clear() + { + (*this) = Bin(); + } + + std::vector & + min() + { + return _min; + } + + std::vector & + max() + { + return _max; + } + + /// Add a new item to the bin containing a given coordinate + void + add(const T &newitem, const Coordinate &coord) + { + indices(coord, _indices); + _item(_indices).push_back(newitem); + } + + /// Add a new item to the bin containing a given coordinate, + /// if an equivalent item is not already present + /// return 'true' if added succesfully, 'false' if not + void + add_once(const T &newitem, const Coordinate &coord) + { + indices(coord, _indices); + std::vector &singlebin = _item(_indices); + for (int i = 0; i < singlebin.size(); i++) { - indices( coord, _indices); - std::vector &singlebin = _item(_indices); - for( int i=0; i - void add_range( const T &newitem, const PCoord &min, const PCoord &max) + + singlebin.push_back(newitem); + } + + /// Add a new item to all bins that fall in cuboid defined by 'min' and 'max' + /// Coordinates + template + void + add_range(const T &newitem, const PCoord &min, const PCoord &max) + { + // std::cout << "begin add_range()" << std::endl; + + // std::cout << "min: " << min << " max: " << max << std::endl; + + int i; + std::vector index_min(_item.order()); + std::vector index_max(_item.order()); + indices(min, index_min); + indices(max, index_max); + + // std::cout << "MIN: "; + // for( j=0; j index_min(_item.order()); - std::vector index_max(_item.order()); - indices( min, index_min); - indices( max, index_max); - - //std::cout << "MIN: "; - //for( j=0; j=0; i--) + i++; + if (i == _item.order()) { - _indices[i]=index_min[i]; + return; } - /////////////// } - - //std::cout << "finish add_range()" << std::endl; - - } - - std::vector& contents( const Coordinate &coord) - { - indices( coord, _indices); - return _item(_indices); + + _indices[i]++; + + for (i = i - 1; i >= 0; i--) + { + _indices[i] = index_min[i]; + } + /////////////// } - - // maximum size of any bin - int max_size() + + // std::cout << "finish add_range()" << std::endl; + } + + std::vector & + contents(const Coordinate &coord) + { + indices(coord, _indices); + return _item(_indices); + } + + // maximum size of any bin + int + max_size() + { + unsigned int max = 0; + for (int i = 0; i < _item.volume(); i++) + if (_item(i).size() > max) + max = _item(i).size(); + return max; + } + + private: + /// Set 'term' to be the indices into '_item' PNDArray of the bin that contains + /// 'coord' + /// Return 'false' if unsuccesful, 'true' if succesful + template + void + indices(const PCoord &coord, std::vector &term) + { + for (int i = 0; i < _item.order(); i++) { - unsigned int max = 0; - for( int i=0; i<_item.volume(); i++) - if( _item(i).size() > max) - max = _item(i).size(); - return max; + // std::cout << "i: " << i << std::endl; + // std::cout << "coord: " << coord[i] << std::endl; + // std::cout << "_min: " << _min[i] << std::endl; + // std::cout << "_max: " << _max[i] << std::endl; + // std::cout << "_incr: " << _incr[i] << std::endl; + + if ((coord[i] < _min[i]) || (coord[i] > _max[i])) + throw std::domain_error("Invalid coord, out of bin range"); + term[i] = std::floor((coord[i] - _min[i]) / _incr[i]); } - - private: - - /// Set 'term' to be the indices into '_item' PNDArray of the bin that contains 'coord' - /// Return 'false' if unsuccesful, 'true' if succesful - template< class PCoord> - void indices(const PCoord &coord, std::vector &term) + } + + /// Set 'term' to be the indices into '_item' PNDArray of the bin that contains + /// 'coord' + /// Return 'false' if unsuccesful, 'true' if succesful + template + void + max_indices(const PCoord &coord, std::vector &term) + { + for (int i = 0; i < _item.order(); i++) { - for( int i=0; i<_item.order(); i++) + // std::cout << "i: " << i << std::endl; + // std::cout << "coord: " << coord[i] << std::endl; + // std::cout << "_min: " << _min[i] << std::endl; + // std::cout << "_max: " << _max[i] << std::endl; + // std::cout << "_incr: " << _incr[i] << std::endl; + + if ((coord[i] < _min[i]) || (coord[i] > _max[i])) + throw std::domain_error("Invalid coord, out of bin range"); + + // std::cout << "coord: " << coord[i] << " _min: " << _min[i] << " _incr: " << + // _incr[i] << std::endl; std::cout << std::floor( (coord[i] - + // _min[i])/_incr[i]) << " " << (coord[i] - _min[i])/_incr[i] << std::endl; + + if (std::floor((coord[i] - _min[i]) / _incr[i]) == + (coord[i] - _min[i]) / _incr[i]) { - //std::cout << "i: " << i << std::endl; - //std::cout << "coord: " << coord[i] << std::endl; - //std::cout << "_min: " << _min[i] << std::endl; - //std::cout << "_max: " << _max[i] << std::endl; - //std::cout << "_incr: " << _incr[i] << std::endl; - - if( (coord[i] < _min[i]) || (coord[i] > _max[i]) ) - throw std::domain_error("Invalid coord, out of bin range"); - term[i] = std::floor( (coord[i] - _min[i])/_incr[i]); + // std::cout << "reduce max" << std::endl; + term[i] = std::floor((coord[i] - _min[i]) / _incr[i]); + term[i]--; } - } - - /// Set 'term' to be the indices into '_item' PNDArray of the bin that contains 'coord' - /// Return 'false' if unsuccesful, 'true' if succesful - template< class PCoord> - void max_indices(const PCoord &coord, std::vector &term) - { - for( int i=0; i<_item.order(); i++) + else { - //std::cout << "i: " << i << std::endl; - //std::cout << "coord: " << coord[i] << std::endl; - //std::cout << "_min: " << _min[i] << std::endl; - //std::cout << "_max: " << _max[i] << std::endl; - //std::cout << "_incr: " << _incr[i] << std::endl; - - if( (coord[i] < _min[i]) || (coord[i] > _max[i]) ) - throw std::domain_error("Invalid coord, out of bin range"); - - //std::cout << "coord: " << coord[i] << " _min: " << _min[i] << " _incr: " << _incr[i] << std::endl; - //std::cout << std::floor( (coord[i] - _min[i])/_incr[i]) << " " << (coord[i] - _min[i])/_incr[i] << std::endl; - - if( std::floor( (coord[i] - _min[i])/_incr[i]) == (coord[i] - _min[i])/_incr[i]) - { - //std::cout << "reduce max" << std::endl; - term[i] = std::floor( (coord[i] - _min[i])/_incr[i]); - term[i]--; - } - else - { - term[i] = std::floor( (coord[i] - _min[i])/_incr[i]); - } + term[i] = std::floor((coord[i] - _min[i]) / _incr[i]); } } - }; - -} + } + }; +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/datastruc/PNDArray.hh b/include/IntegrationTools/datastruc/PNDArray.hh index 9efc6190c..1d16163cd 100644 --- a/include/IntegrationTools/datastruc/PNDArray.hh +++ b/include/IntegrationTools/datastruc/PNDArray.hh @@ -2,159 +2,175 @@ #ifndef PNDArray_HH #define PNDArray_HH -#include -#include #include +#include +#include namespace PRISMS { - - /// A Tensor class, - /// takes int valued IndexContainer and returns OutType value - /// - template< class OutType> - class PNDArray + + /// A Tensor class, + /// takes int valued IndexContainer and returns OutType value + /// + template + class PNDArray + { + std::vector _dim; // dimension along each compenent + std::vector _unroll; // used to translate from tensor indices to linear index + std::vector _val; // unrolled list of coefficients (first index is outer + // loop) + int _order; // _dim.size() + int _volume; // _coeff_tensor.size() = product_i(_dim[i]) + + public: + PNDArray() + : _order(0) + , _volume(0) + {} + + PNDArray(const std::vector &dim) { - std::vector< int > _dim; // dimension along each compenent - std::vector< int > _unroll; // used to translate from tensor indices to linear index - std::vector< OutType > _val; // unrolled list of coefficients (first index is outer loop) - int _order; // _dim.size() - int _volume; // _coeff_tensor.size() = product_i(_dim[i]) - - public: - - PNDArray() : _order(0), _volume(0) {} - PNDArray(const std::vector &dim) - { - resize(dim); - } - - PNDArray(const std::vector &dim, const std::vector &value) - { - resize(dim); - if( _volume != value.size()) - { - std::cerr << "Error in PNDArray(const std::vector &dim, const std::vector &value)." << std::endl; - std::cerr << " value.size() does not match volume based on dim." << std::endl; - exit(1); - } - _val = value; - } - - int order() const - { - return _order; - } - - int volume() const - { - return _volume; - } - - void resize(const std::vector &dim) - { - _dim = dim; - _order = _dim.size(); - _volume = calc_volume(_dim); - _val.resize(_volume); - generate_unroll(); - } - - void reshape(const std::vector &dim) - { - if( calc_volume(dim) == _volume) - { - _dim = dim; - _order = _dim.size(); - generate_unroll(); - } - else - { - std::cerr << "Error in PNDArray::reshape. Volume is not equivalent." << std::endl; - exit(1); - } - } - - void clear() - { - _val.clear(); - _dim.clear(); - _unroll.clear(); - _order = 0; - _volume = 0; - - } - - const std::vector& dim() const - { - return _dim; - } - - int dim(int i) const + resize(dim); + } + + PNDArray(const std::vector &dim, const std::vector &value) + { + resize(dim); + if (_volume != value.size()) { - return _dim[i]; + std::cerr << "Error in PNDArray(const std::vector &dim, const " + "std::vector &value)." + << std::endl; + std::cerr << " value.size() does not match volume based on dim." << std::endl; + exit(1); } - - OutType& operator()( int i) + _val = value; + } + + int + order() const + { + return _order; + } + + int + volume() const + { + return _volume; + } + + void + resize(const std::vector &dim) + { + _dim = dim; + _order = _dim.size(); + _volume = calc_volume(_dim); + _val.resize(_volume); + generate_unroll(); + } + + void + reshape(const std::vector &dim) + { + if (calc_volume(dim) == _volume) { - return _val[i]; + _dim = dim; + _order = _dim.size(); + generate_unroll(); } - - template - OutType& operator()( const IndexContainer &term) + else { - return _val[linear_index(term)]; + std::cerr << "Error in PNDArray::reshape. Volume is not equivalent." + << std::endl; + exit(1); } - - template - int linear_index( const IndexContainer &term) const + } + + void + clear() + { + _val.clear(); + _dim.clear(); + _unroll.clear(); + _order = 0; + _volume = 0; + } + + const std::vector & + dim() const + { + return _dim; + } + + int + dim(int i) const + { + return _dim[i]; + } + + OutType & + operator()(int i) + { + return _val[i]; + } + + template + OutType & + operator()(const IndexContainer &term) + { + return _val[linear_index(term)]; + } + + template + int + linear_index(const IndexContainer &term) const + { + int lindex = 0; + for (unsigned int i = 0; i < _unroll.size(); i++) + lindex += term[i] * _unroll[i]; + return lindex; + } + + template + void + tensor_indices(int lindex, IndexContainer &term) const + { // assumes term.size() == order() (the tensor order) + // not sure if this is how we want to do it, but it avoids assuming push_back() or + // resize() + + for (int i = 0; i < _unroll.size(); i++) { - int lindex = 0; - for( unsigned int i=0; i<_unroll.size(); i++) - lindex += term[i]*_unroll[i]; - return lindex; - } - - template - void tensor_indices( int lindex, IndexContainer &term) const - { // assumes term.size() == order() (the tensor order) - // not sure if this is how we want to do it, but it avoids assuming push_back() or resize() - - for( int i=0; i<_unroll.size(); i++) - { - term[i] = lindex/_unroll[i]; - lindex -= term[i]*_unroll[i]; - } + term[i] = lindex / _unroll[i]; + lindex -= term[i] * _unroll[i]; } - - - private: - - int calc_volume( const std::vector &dim) + } + + private: + int + calc_volume(const std::vector &dim) + { + if (dim.size() == 0) { - if( dim.size() == 0) - { - return 0; - } - - int vol = 1; - for( unsigned int i=0; i=0; i--) - _unroll[i] = _unroll[i+1]*_dim[i+1]; + vol *= dim[i]; } - - }; -} + return vol; + } + void + generate_unroll() + { + _unroll.resize(_dim.size()); + _unroll[_dim.size() - 1] = 1; + for (int i = _dim.size() - 2; i >= 0; i--) + _unroll[i] = _unroll[i + 1] * _dim[i + 1]; + } + }; +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/extern/PExtern.hh b/include/IntegrationTools/extern/PExtern.hh index 8fe5e57d9..e193a6002 100644 --- a/include/IntegrationTools/extern/PExtern.hh +++ b/include/IntegrationTools/extern/PExtern.hh @@ -2,284 +2,510 @@ #ifndef PExtern_HH #define PExtern_HH -#include -#include -#include -#include - -#include "../PFunction.hh" #include "../PField.hh" +#include "../PFunction.hh" +#include +#include +#include +#include -// In future, might have more complicated OutType, +// In future, might have more complicated OutType, // so make all have 'void' return and pass everything by reference - extern "C" { - // Functions for using a PSimpleBase externally (say Python or Fortran) - // written for VarContainer=double*, OutType=double, hence 'dsd' in function names - - void PSimpleFunction_dsd_new(char* name, PRISMS::PSimpleBase* &f); - - void PSimpleFunction_dsd_delete(PRISMS::PSimpleBase* &f); - - void PSimpleFunction_dsd_name(PRISMS::PSimpleBase* f, char* name); - - void PSimpleFunction_dsd_calc( PRISMS::PSimpleBase* f, double* var, double &val); - - void PSimpleFunction_dsd_get( PRISMS::PSimpleBase* f, double &val); - - - - // Functions for using a PSimpleBase externally (say Python or Fortran) - // written for VarContainer=double, OutType=double, hence 'dd' in function names - - void PSimpleFunction_dd_new(char* name, PRISMS::PSimpleBase* &f); - - void PSimpleFunction_dd_delete(PRISMS::PSimpleBase* &f); - - void PSimpleFunction_dd_name(PRISMS::PSimpleBase* f, char* name); - - void PSimpleFunction_dd_calc( PRISMS::PSimpleBase* f, double var, double &val); - - void PSimpleFunction_dd_get( PRISMS::PSimpleBase* f, double &val); - - - - - - // Functions for using a PFuncBase externally (say Python or Fortran) - // written for VarContainer=double*, OutType=double, hence 'dsd' in function names - - void PFunction_dsd_new(char* name, PRISMS::PFuncBase* &f); - - void PFunction_dsd_delete(PRISMS::PFuncBase* &f); - - void PFunction_dsd_name(PRISMS::PFuncBase* f, char* name); - - void PFunction_dsd_size(PRISMS::PFuncBase* f, int &size); - - void PFunction_dsd_var_name(PRISMS::PFuncBase* f, int i, char* var_name); - - void PFunction_dsd_var_description(PRISMS::PFuncBase* f, int i, char* var_description); - - //void PFunction_dsd_simplefunc(PRISMS::PFuncBase* f, PSimpleBase* &simplefunc); - //void PFunction_dsd_grad_simplefunc(PRISMS::PFuncBase* f, int di, PSimpleBase* &simplefunc); - //void PFunction_dsd_hess_simplefunc(PRISMS::PFuncBase* f, int di, int dj, PSimpleBase* &simplefunc); - - void PFunction_dsd_calc(PRISMS::PFuncBase* f, double* var, double &val); - - void PFunction_dsd_calc_grad(PRISMS::PFuncBase* f, double* var, int di, double &val); - - void PFunction_dsd_calc_hess(PRISMS::PFuncBase* f, double* var, int di, int dj, double &val); - - void PFunction_dsd_eval(PRISMS::PFuncBase* f, double* var); - - void PFunction_dsd_eval_grad(PRISMS::PFuncBase* f, double* var, int di); - - void PFunction_dsd_eval_hess(PRISMS::PFuncBase* f, double* var, int di, int dj); - - void PFunction_dsd_get(PRISMS::PFuncBase* f, double &val); - - void PFunction_dsd_get_grad(PRISMS::PFuncBase* f, int di, double &val); - - void PFunction_dsd_get_hess(PRISMS::PFuncBase* f, int di, int dj, double &val); - - - - - - // Functions for using a PBasisSetBase externally (say Python or Fortran) - // written for InType=double, OutType=double, hence 'dd' in function names - - void PBasisSet_dd_new(char* name, PRISMS::PBasisSetBase* &b, int N); - - void PBasisSet_dd_delete(PRISMS::PBasisSetBase* &b); - - void PBasisSet_dd_name(PRISMS::PBasisSetBase* b, char* name); - - void PBasisSet_dd_description(PRISMS::PBasisSetBase* b, char* description); - - void PBasisSet_dd_size(PRISMS::PBasisSetBase* b, int &size); - - void PBasisSet_dd_resize(PRISMS::PBasisSetBase* b, int N); - - void PBasisSet_dd_max_size(PRISMS::PBasisSetBase* b, int &max_size); - - //void PBasisSet_dd_basis_function(PRISMS::PFuncBase* b, int term, PFuncBase* &f); - - void PBasisSet_dd_calc(PRISMS::PBasisSetBase* b, int term, double var, double &val); - - void PBasisSet_dd_calc_grad(PRISMS::PBasisSetBase* b, int term, double var, double &val); - - void PBasisSet_dd_calc_hess(PRISMS::PBasisSetBase* b, int term, double var, double &val); - - void PBasisSet_dd_eval(PRISMS::PBasisSetBase* b, double var); - - void PBasisSet_dd_eval_grad(PRISMS::PBasisSetBase* b, double var); - - void PBasisSet_dd_eval_hess(PRISMS::PBasisSetBase* b, double var); - - void PBasisSet_dd_get(PRISMS::PBasisSetBase* b, int term, double &val); - - void PBasisSet_dd_get_grad(PRISMS::PBasisSetBase* b, int term, double &val); - - void PBasisSet_dd_get_hess(PRISMS::PBasisSetBase* b, int term, double &val); - - //void PBasisSet_dd_getall(PRISMS::PBasisSetBase* b, const double* &val); - //void PBasisSet_dd_getall_grad(PRISMS::PBasisSetBase* b, const double* &val); - //void PBasisSet_dd_getall_hess(PRISMS::PBasisSetBase* b, const double* &val); - - - - - - // Functions for using a PSeriesFunction externally (say Python or Fortran) - // written for InType=double, OutType=double, VarContainer=double*, IndexContainer=int*, hence 'dsis' - - // Functions for using a PSeriesFunction externally (say Python or Fortran) - // written for InType=double, OutType=double, VarContainer=double*, IndexContainer=int*, hence 'dsis' - - void PSeriesFunction_dsis_new(PRISMS::PSeriesFunction* &f); - - void PSeriesFunction_dsis_setnew(PRISMS::PSeriesFunction* &f, PRISMS::PBasisSetBase** basis_set, int order); - - void PSeriesFunction_dsis_delete(PRISMS::PSeriesFunction* &f); - - void PSeriesFunction_dsis_clear(PRISMS::PSeriesFunction* f); - - void PSeriesFunction_dsis_set(PRISMS::PSeriesFunction* f, PRISMS::PBasisSetBase** basis_set, int order); - - void PSeriesFunction_dsis_order(PRISMS::PSeriesFunction* f, int &order); - - void PSeriesFunction_dsis_volume(PRISMS::PSeriesFunction* f, int &volume); - - void PSeriesFunction_dsis_dim(PRISMS::PSeriesFunction* f, int i, int &dim); - - void PSeriesFunction_dsis_get_linear_coeff(PRISMS::PSeriesFunction* f, int i, double &coeff); - - void PSeriesFunction_dsis_get_tensor_coeff(PRISMS::PSeriesFunction* f, int* term, double &coeff); - - void PSeriesFunction_dsis_set_linear_coeff(PRISMS::PSeriesFunction* f, int i, double coeff); - - void PSeriesFunction_dsis_set_tensor_coeff(PRISMS::PSeriesFunction* f, int* term, double coeff); - - void PSeriesFunction_dsis_linear_index(PRISMS::PSeriesFunction* f, int* term, int &linear_index); - - void PSeriesFunction_dsis_tensor_indices(PRISMS::PSeriesFunction* f, int linear_index, int* term); - - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - - void PSeriesFunction_dsis_calc(PRISMS::PSeriesFunction* f, double* var, double &val); - - void PSeriesFunction_dsis_calc_grad(PRISMS::PSeriesFunction* f, double* var, int di, double &val); - - void PSeriesFunction_dsis_calc_hess(PRISMS::PSeriesFunction* f, double* var, int di, int dj, double &val); - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - - void PSeriesFunction_dsis_eval(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_eval_grad(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_eval_hess(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_get(PRISMS::PSeriesFunction* f, double &val); - - void PSeriesFunction_dsis_get_grad(PRISMS::PSeriesFunction* f, int di, double &val); - - void PSeriesFunction_dsis_get_hess(PRISMS::PSeriesFunction* f, int di, int dj, double &val); - - - // ---------------------------------------------------------- - // Functions for evaluating basis functions & their derivatives: - - // Use these functions if you want to evaluate a single value - - // use basis index and term index for individual basis function - - void PSeriesFunction_dsis_calc_basis(PRISMS::PSeriesFunction* f, int bindex, int term, double* var, double &val); - - void PSeriesFunction_dsis_calc_basis_grad(PRISMS::PSeriesFunction* f, int bindex, int term, double* var, double &val); - - void PSeriesFunction_dsis_calc_basis_hess(PRISMS::PSeriesFunction* f, int bindex, int term, double* var, double &val); - - // or use tensor indices to evaluate basis function product - void PSeriesFunction_dsis_calc_tensor_basis(PRISMS::PSeriesFunction* f, int* term, double* var, double &val); - - void PSeriesFunction_dsis_calc_tensor_basis_grad(PRISMS::PSeriesFunction* f, int* term, double* var, int di, double &val); - - void PSeriesFunction_dsis_calc_tensor_basis_hess(PRISMS::PSeriesFunction* f, int* term, double* var, int di, int dj, double &val); - - // ---------------------------------------------------------- - // Use these functions to evaluate all basis functions, - // then use following methods to access results. - - void PSeriesFunction_dsis_eval_basis_all(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_eval_basis(PRISMS::PSeriesFunction* f, double* var, int i); - - void PSeriesFunction_dsis_eval_basis_grad_all(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_eval_basis_grad(PRISMS::PSeriesFunction* f, double* var, int i); - - void PSeriesFunction_dsis_eval_basis_hess_all(PRISMS::PSeriesFunction* f, double* var); - - void PSeriesFunction_dsis_eval_basis_hess(PRISMS::PSeriesFunction* f, double* var, int i); - - - // use basis index and term index for individual basis function - void PSeriesFunction_dsis_get_basis(PRISMS::PSeriesFunction* f, int bindex, int term, double &val); - - void PSeriesFunction_dsis_get_basis_grad(PRISMS::PSeriesFunction* f, int bindex, int term, double &val); - - void PSeriesFunction_dsis_get_basis_hess(PRISMS::PSeriesFunction* f, int bindex, int term, double &val); - - // or use tensor indices to evaluate basis function product - void PSeriesFunction_dsis_get_tensor_basis(PRISMS::PSeriesFunction* f, int* term, double &val); - - void PSeriesFunction_dsis_get_tensor_basis_grad(PRISMS::PSeriesFunction* f, int* term, int di, double &val); - - void PSeriesFunction_dsis_get_tensor_basis_hess(PRISMS::PSeriesFunction* f, int* term, int di, int dj, double &val); - - - - // Functions for using constructing a 2D PRISMS::Body externally (say Python or Fortran), - // allowing access to PFields - // written for Coordinate=double*, DIM=2 - - void Body2D_new(char* vtkfile, PRISMS::Body* &b); - - void Body2D_delete(PRISMS::Body* &b); - - - // Functions for using a 2D scalar PField externally (say Python or Fortran), as a PFunction. - // From a Body pointer, returns a pointer to a PFuncBase - // written for Coordinate=double*, OutType=double, DIM=2 - - void ScalarField2D(char* name, PRISMS::Body* b, PRISMS::PFuncBase* &f); - - - - // Functions for using constructing a 3D PRISMS::Body externally (say Python or Fortran), - // allowing access to PFields - // written for Coordinate=double*, DIM=3 - - void Body3D_new(char* vtkfile, PRISMS::Body* &b); - - void Body3D_delete(PRISMS::Body* &b); - - - // Functions for using a 3D scalar PField externally (say Python or Fortran), as a PFunction. - // From a Body pointer, returns a pointer to a PFuncBase - // written for Coordinate=double*, OutType=double, DIM=3 - - void ScalarField3D(char* name, PRISMS::Body* b, PRISMS::PFuncBase* &f); - - -} + // Functions for using a PSimpleBase externally (say Python or Fortran) + // written for VarContainer=double*, OutType=double, hence 'dsd' in function names + + void + PSimpleFunction_dsd_new(char *name, PRISMS::PSimpleBase *&f); + + void + PSimpleFunction_dsd_delete(PRISMS::PSimpleBase *&f); + + void + PSimpleFunction_dsd_name(PRISMS::PSimpleBase *f, char *name); + + void + PSimpleFunction_dsd_calc(PRISMS::PSimpleBase *f, + double *var, + double &val); + + void + PSimpleFunction_dsd_get(PRISMS::PSimpleBase *f, double &val); + + // Functions for using a PSimpleBase externally (say Python or Fortran) + // written for VarContainer=double, OutType=double, hence 'dd' in function names + + void + PSimpleFunction_dd_new(char *name, PRISMS::PSimpleBase *&f); + + void + PSimpleFunction_dd_delete(PRISMS::PSimpleBase *&f); + + void + PSimpleFunction_dd_name(PRISMS::PSimpleBase *f, char *name); + + void + PSimpleFunction_dd_calc(PRISMS::PSimpleBase *f, + double var, + double &val); + + void + PSimpleFunction_dd_get(PRISMS::PSimpleBase *f, double &val); + + // Functions for using a PFuncBase externally (say Python or Fortran) + // written for VarContainer=double*, OutType=double, hence 'dsd' in function names + + void + PFunction_dsd_new(char *name, PRISMS::PFuncBase *&f); + + void + PFunction_dsd_delete(PRISMS::PFuncBase *&f); + + void + PFunction_dsd_name(PRISMS::PFuncBase *f, char *name); + + void + PFunction_dsd_size(PRISMS::PFuncBase *f, int &size); + + void + PFunction_dsd_var_name(PRISMS::PFuncBase *f, int i, char *var_name); + + void + PFunction_dsd_var_description(PRISMS::PFuncBase *f, + int i, + char *var_description); + + // void PFunction_dsd_simplefunc(PRISMS::PFuncBase* f, + // PSimpleBase* &simplefunc); void + // PFunction_dsd_grad_simplefunc(PRISMS::PFuncBase* f, int di, + // PSimpleBase* &simplefunc); void + // PFunction_dsd_hess_simplefunc(PRISMS::PFuncBase* f, int di, int dj, + // PSimpleBase* &simplefunc); + + void + PFunction_dsd_calc(PRISMS::PFuncBase *f, double *var, double &val); + + void + PFunction_dsd_calc_grad(PRISMS::PFuncBase *f, + double *var, + int di, + double &val); + + void + PFunction_dsd_calc_hess(PRISMS::PFuncBase *f, + double *var, + int di, + int dj, + double &val); + + void + PFunction_dsd_eval(PRISMS::PFuncBase *f, double *var); + + void + PFunction_dsd_eval_grad(PRISMS::PFuncBase *f, double *var, int di); + + void + PFunction_dsd_eval_hess(PRISMS::PFuncBase *f, + double *var, + int di, + int dj); + + void + PFunction_dsd_get(PRISMS::PFuncBase *f, double &val); + + void + PFunction_dsd_get_grad(PRISMS::PFuncBase *f, int di, double &val); + + void + PFunction_dsd_get_hess(PRISMS::PFuncBase *f, + int di, + int dj, + double &val); + + // Functions for using a PBasisSetBase externally (say Python or Fortran) + // written for InType=double, OutType=double, hence 'dd' in function names + + void + PBasisSet_dd_new(char *name, PRISMS::PBasisSetBase *&b, int N); + + void + PBasisSet_dd_delete(PRISMS::PBasisSetBase *&b); + + void + PBasisSet_dd_name(PRISMS::PBasisSetBase *b, char *name); + + void + PBasisSet_dd_description(PRISMS::PBasisSetBase *b, char *description); + + void + PBasisSet_dd_size(PRISMS::PBasisSetBase *b, int &size); + + void + PBasisSet_dd_resize(PRISMS::PBasisSetBase *b, int N); + + void + PBasisSet_dd_max_size(PRISMS::PBasisSetBase *b, int &max_size); + + // void PBasisSet_dd_basis_function(PRISMS::PFuncBase* b, int term, + // PFuncBase* &f); + + void + PBasisSet_dd_calc(PRISMS::PBasisSetBase *b, + int term, + double var, + double &val); + + void + PBasisSet_dd_calc_grad(PRISMS::PBasisSetBase *b, + int term, + double var, + double &val); + + void + PBasisSet_dd_calc_hess(PRISMS::PBasisSetBase *b, + int term, + double var, + double &val); + + void + PBasisSet_dd_eval(PRISMS::PBasisSetBase *b, double var); + + void + PBasisSet_dd_eval_grad(PRISMS::PBasisSetBase *b, double var); + + void + PBasisSet_dd_eval_hess(PRISMS::PBasisSetBase *b, double var); + + void + PBasisSet_dd_get(PRISMS::PBasisSetBase *b, int term, double &val); + + void + PBasisSet_dd_get_grad(PRISMS::PBasisSetBase *b, int term, double &val); + + void + PBasisSet_dd_get_hess(PRISMS::PBasisSetBase *b, int term, double &val); + + // void PBasisSet_dd_getall(PRISMS::PBasisSetBase* b, const double* + // &val); void PBasisSet_dd_getall_grad(PRISMS::PBasisSetBase* b, const + // double* &val); void PBasisSet_dd_getall_hess(PRISMS::PBasisSetBase* b, + // const double* &val); + + // Functions for using a PSeriesFunction externally (say Python or Fortran) + // written for InType=double, OutType=double, VarContainer=double*, + // IndexContainer=int*, hence 'dsis' + + // Functions for using a PSeriesFunction externally (say Python or Fortran) + // written for InType=double, OutType=double, VarContainer=double*, + // IndexContainer=int*, hence 'dsis' + + void + PSeriesFunction_dsis_new(PRISMS::PSeriesFunction *&f); + + void + PSeriesFunction_dsis_setnew( + PRISMS::PSeriesFunction *&f, + PRISMS::PBasisSetBase **basis_set, + int order); + void + PSeriesFunction_dsis_delete( + PRISMS::PSeriesFunction *&f); + + void + PSeriesFunction_dsis_clear(PRISMS::PSeriesFunction *f); + + void + PSeriesFunction_dsis_set(PRISMS::PSeriesFunction *f, + PRISMS::PBasisSetBase **basis_set, + int order); + + void + PSeriesFunction_dsis_order(PRISMS::PSeriesFunction *f, + int &order); + + void + PSeriesFunction_dsis_volume(PRISMS::PSeriesFunction *f, + int &volume); + + void + PSeriesFunction_dsis_dim(PRISMS::PSeriesFunction *f, + int i, + int &dim); + + void + PSeriesFunction_dsis_get_linear_coeff( + PRISMS::PSeriesFunction *f, + int i, + double &coeff); + + void + PSeriesFunction_dsis_get_tensor_coeff( + PRISMS::PSeriesFunction *f, + int *term, + double &coeff); + + void + PSeriesFunction_dsis_set_linear_coeff( + PRISMS::PSeriesFunction *f, + int i, + double coeff); + + void + PSeriesFunction_dsis_set_tensor_coeff( + PRISMS::PSeriesFunction *f, + int *term, + double coeff); + + void + PSeriesFunction_dsis_linear_index( + PRISMS::PSeriesFunction *f, + int *term, + int &linear_index); + + void + PSeriesFunction_dsis_tensor_indices( + PRISMS::PSeriesFunction *f, + int linear_index, + int *term); + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + + void + PSeriesFunction_dsis_calc(PRISMS::PSeriesFunction *f, + double *var, + double &val); + + void + PSeriesFunction_dsis_calc_grad( + PRISMS::PSeriesFunction *f, + double *var, + int di, + double &val); + + void + PSeriesFunction_dsis_calc_hess( + PRISMS::PSeriesFunction *f, + double *var, + int di, + int dj, + double &val); + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + + void + PSeriesFunction_dsis_eval(PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_eval_grad( + PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_eval_hess( + PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_get(PRISMS::PSeriesFunction *f, + double &val); + + void + PSeriesFunction_dsis_get_grad( + PRISMS::PSeriesFunction *f, + int di, + double &val); + + void + PSeriesFunction_dsis_get_hess( + PRISMS::PSeriesFunction *f, + int di, + int dj, + double &val); + + // ---------------------------------------------------------- + // Functions for evaluating basis functions & their derivatives: + + // Use these functions if you want to evaluate a single value + + // use basis index and term index for individual basis function + + void + PSeriesFunction_dsis_calc_basis( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double *var, + double &val); + + void + PSeriesFunction_dsis_calc_basis_grad( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double *var, + double &val); + + void + PSeriesFunction_dsis_calc_basis_hess( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double *var, + double &val); + + // or use tensor indices to evaluate basis function product + void + PSeriesFunction_dsis_calc_tensor_basis( + PRISMS::PSeriesFunction *f, + int *term, + double *var, + double &val); + + void + PSeriesFunction_dsis_calc_tensor_basis_grad( + PRISMS::PSeriesFunction *f, + int *term, + double *var, + int di, + double &val); + + void + PSeriesFunction_dsis_calc_tensor_basis_hess( + PRISMS::PSeriesFunction *f, + int *term, + double *var, + int di, + int dj, + double &val); + + // ---------------------------------------------------------- + // Use these functions to evaluate all basis functions, + // then use following methods to access results. + + void + PSeriesFunction_dsis_eval_basis_all( + PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_eval_basis( + PRISMS::PSeriesFunction *f, + double *var, + int i); + + void + PSeriesFunction_dsis_eval_basis_grad_all( + PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_eval_basis_grad( + PRISMS::PSeriesFunction *f, + double *var, + int i); + + void + PSeriesFunction_dsis_eval_basis_hess_all( + PRISMS::PSeriesFunction *f, + double *var); + + void + PSeriesFunction_dsis_eval_basis_hess( + PRISMS::PSeriesFunction *f, + double *var, + int i); + + // use basis index and term index for individual basis function + void + PSeriesFunction_dsis_get_basis( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double &val); + + void + PSeriesFunction_dsis_get_basis_grad( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double &val); + + void + PSeriesFunction_dsis_get_basis_hess( + PRISMS::PSeriesFunction *f, + int bindex, + int term, + double &val); + + // or use tensor indices to evaluate basis function product + void + PSeriesFunction_dsis_get_tensor_basis( + PRISMS::PSeriesFunction *f, + int *term, + double &val); + + void + PSeriesFunction_dsis_get_tensor_basis_grad( + PRISMS::PSeriesFunction *f, + int *term, + int di, + double &val); + + void + PSeriesFunction_dsis_get_tensor_basis_hess( + PRISMS::PSeriesFunction *f, + int *term, + int di, + int dj, + double &val); + + // Functions for using constructing a 2D PRISMS::Body externally (say Python or + // Fortran), + // allowing access to PFields + // written for Coordinate=double*, DIM=2 + + void + Body2D_new(char *vtkfile, PRISMS::Body *&b); + + void + Body2D_delete(PRISMS::Body *&b); + + // Functions for using a 2D scalar PField externally (say Python or Fortran), as a + // PFunction. + // From a Body pointer, returns a pointer to a PFuncBase + // written for Coordinate=double*, OutType=double, DIM=2 + + void + ScalarField2D(char *name, + PRISMS::Body *b, + PRISMS::PFuncBase *&f); + + // Functions for using constructing a 3D PRISMS::Body externally (say Python or + // Fortran), + // allowing access to PFields + // written for Coordinate=double*, DIM=3 + + void + Body3D_new(char *vtkfile, PRISMS::Body *&b); + + void + Body3D_delete(PRISMS::Body *&b); + + // Functions for using a 3D scalar PField externally (say Python or Fortran), as a + // PFunction. + // From a Body pointer, returns a pointer to a PFuncBase + // written for Coordinate=double*, OutType=double, DIM=3 + + void + ScalarField3D(char *name, + PRISMS::Body *b, + PRISMS::PFuncBase *&f); +} #endif \ No newline at end of file diff --git a/include/IntegrationTools/extern/PLibraryExtern.hh b/include/IntegrationTools/extern/PLibraryExtern.hh index 1d442ab19..d2a27768a 100644 --- a/include/IntegrationTools/extern/PLibraryExtern.hh +++ b/include/IntegrationTools/extern/PLibraryExtern.hh @@ -1,43 +1,52 @@ #ifndef PLIBRARY_HH #define PLIBRARY_HH -#include -#include #include "../PFunction.hh" +#include +#include namespace PRISMS { - /// Library where you can find functions and basis sets - /// - namespace PLibrary - { - // Use these functions to checkout objects which manage their own memory - - void checkout( std::string name, PSimpleFunction< double, double > &simplefunc); - void checkout( std::string name, PSimpleFunction< std::vector, double > &simplefunc); - void checkout( std::string name, PSimpleFunction< double*, double > &simplefunc); - - void checkout( std::string name, PFunction< std::vector, double > &func); - void checkout( std::string name, PFunction< double*, double > &func); - - void checkout( std::string name, PBasisSet< double, double > &basis_set, int N); - - - - // Use these functions to checkout new 'Base' objects which the user must delete - - void checkout( std::string name, PSimpleBase< double, double > *&simplefunc); - void checkout( std::string name, PSimpleBase< std::vector, double > *&simplefunc); - void checkout( std::string name, PSimpleBase< double*, double > *&simplefunc); - - void checkout( std::string name, PFuncBase< std::vector, double > *&func); - void checkout( std::string name, PFuncBase< double*, double > *&func); - - void checkout( std::string name, PBasisSetBase< double, double > *&basis_set, int N); - } - -} - + /// Library where you can find functions and basis sets + /// + namespace PLibrary + { + // Use these functions to checkout objects which manage their own memory + + void + checkout(std::string name, PSimpleFunction &simplefunc); + void + checkout(std::string name, PSimpleFunction, double> &simplefunc); + void + checkout(std::string name, PSimpleFunction &simplefunc); + + void + checkout(std::string name, PFunction, double> &func); + void + checkout(std::string name, PFunction &func); + + void + checkout(std::string name, PBasisSet &basis_set, int N); + + // Use these functions to checkout new 'Base' objects which the user must delete + + void + checkout(std::string name, PSimpleBase *&simplefunc); + void + checkout(std::string name, PSimpleBase, double> *&simplefunc); + void + checkout(std::string name, PSimpleBase *&simplefunc); + + void + checkout(std::string name, PFuncBase, double> *&func); + void + checkout(std::string name, PFuncBase *&func); + + void + checkout(std::string name, PBasisSetBase *&basis_set, int N); + } // namespace PLibrary + +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/pfield/Body.hh b/include/IntegrationTools/pfield/Body.hh index 5cb94a668..2fb4ea0bd 100644 --- a/include/IntegrationTools/pfield/Body.hh +++ b/include/IntegrationTools/pfield/Body.hh @@ -4,192 +4,188 @@ #include "./Mesh.hh" #include "./PField.hh" -#include #include -#include +#include #include +#include namespace PRISMS { - /// A class for a Body: a combination of Mesh and Field(s)) - /// - template< class Coordinate, int DIM> - class Body - { - public: - - Mesh mesh; + /// A class for a Body: a combination of Mesh and Field(s)) + /// + template + class Body + { + public: + Mesh mesh; - std::vector< PField > scalar_field; + std::vector> scalar_field; - //std::vector< PField, DIM > > vector_field; + // std::vector< PField, DIM > > vector_field; - //std::vector< PField, DIM > > tensor_field; + // std::vector< PField, DIM > > tensor_field; + // ---------------------------------------------------------- + // Constructors + Body() {}; - // ---------------------------------------------------------- - // Constructors - Body(){}; - - /// Read from a 2D vtk file - /// For now: - /// only ASCII files - /// only rectilinear grids (though output as UNSTRUCTURED_GRID) - /// only (2d) Quad elements - /// - void read_vtk( const std::string &vtkfile) - { - std::cout << "Begin reading vtk file" << std::endl; + /// Read from a 2D vtk file + /// For now: + /// only ASCII files + /// only rectilinear grids (though output as UNSTRUCTURED_GRID) + /// only (2d) Quad elements + /// + void + read_vtk(const std::string &vtkfile) + { + std::cout << "Begin reading vtk file" << std::endl; + // read in vtk file here + std::ifstream infile_mesh(vtkfile.c_str()); - // read in vtk file here - std::ifstream infile_mesh(vtkfile.c_str()); + // read mesh info + mesh.read_vtk(infile_mesh); - // read mesh info - mesh.read_vtk(infile_mesh); + std::ifstream infile(vtkfile.c_str()); - std::ifstream infile(vtkfile.c_str()); + // read point data + std::istringstream ss; + std::string str, name, type, line; + int numcomp; + unsigned long int Npoints; - // read point data - std::istringstream ss; - std::string str, name, type, line; - int numcomp; - unsigned long int Npoints; + while (!infile.eof()) + { + std::getline(infile, line); - while(!infile.eof()) + if (line[0] == 'P') { - - std::getline( infile, line); - - if( line[0] == 'P') + if (line.size() > 9 && line.substr(0, 10) == "POINT_DATA") { - if( line.size() > 9 && line.substr(0,10) == "POINT_DATA") - { - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); - ss >> str >> Npoints; - - } + // std::cout << line << "\n"; + ss.clear(); + ss.str(line); + ss >> str >> Npoints; } - else if( line[0] == 'S') + } + else if (line[0] == 'S') + { + if (line.size() > 6 && line.substr(0, 7) == "SCALARS") { - if( line.size() > 6 && line.substr(0,7) == "SCALARS") - { - ss.clear(); - ss.str(line); - ss >> str >> name >> type >> numcomp; - - // read LOOKUP_TABLE line - std::getline( infile, line); + ss.clear(); + ss.str(line); + ss >> str >> name >> type >> numcomp; - // read data - std::cout << "begin reading data" << std::endl; - std::vector data(Npoints); - for( unsigned int i=0; i> data[i]; - //std::cout << data[i] << std::endl; - } - std::cout << " done" << std::endl; - - - // construct field - std::vector var_name(DIM); - std::vector var_description(DIM); - - if( DIM >= 2) - { - var_name[0] = "x"; - var_description[0] = "x coordinate"; - var_name[1] = "y"; - var_description[1] = "y coordinate"; - - } - if( DIM >= 3) - { - var_name[2] = "z"; - var_description[2] = "z coordinate"; - - } + // read LOOKUP_TABLE line + std::getline(infile, line); + // read data + std::cout << "begin reading data" << std::endl; + std::vector data(Npoints); + for (unsigned int i = 0; i < Npoints; i++) + { + infile >> data[i]; + // std::cout << data[i] << std::endl; + } + std::cout << " done" << std::endl; - std::cout << "Construct PField '" << name << "'" << std::endl; - scalar_field.push_back( PField( name, var_name, var_description, mesh, data, 0.0) ); - std::cout << " done" << std::endl; + // construct field + std::vector var_name(DIM); + std::vector var_description(DIM); + if (DIM >= 2) + { + var_name[0] = "x"; + var_description[0] = "x coordinate"; + var_name[1] = "y"; + var_description[1] = "y coordinate"; } - } - // Alternative field descriptor used by ParaView (holds the same information as the "SCALAR" line above) - else if( line[0] == 'F') - { - if( line.size() > 14 && line.substr(0,15) == "FIELD FieldData") + if (DIM >= 3) { - ss.clear(); - ss.str(line); - ss >> str >> numcomp; - - // read LOOKUP_TABLE line - std::getline( infile, line); - - ss.clear(); - ss.str(line); - ss >> name >> numcomp >> Npoints >> type; - - - // read data - std::cout << "begin reading data" << std::endl; - std::vector data(Npoints); - for( unsigned int i=0; i> data[i]; - } - std::cout << " done" << std::endl; - - - // construct field - std::vector var_name(DIM); - std::vector var_description(DIM); + var_name[2] = "z"; + var_description[2] = "z coordinate"; + } - if( DIM >= 2) - { - var_name[0] = "x"; - var_description[0] = "x coordinate"; - var_name[1] = "y"; - var_description[1] = "y coordinate"; + std::cout << "Construct PField '" << name << "'" << std::endl; + scalar_field.push_back(PField(name, + var_name, + var_description, + mesh, + data, + 0.0)); + std::cout << " done" << std::endl; + } + } + // Alternative field descriptor used by ParaView (holds the same information as + // the "SCALAR" line above) + else if (line[0] == 'F') + { + if (line.size() > 14 && line.substr(0, 15) == "FIELD FieldData") + { + ss.clear(); + ss.str(line); + ss >> str >> numcomp; - } - if( DIM >= 3) - { - var_name[2] = "z"; - var_description[2] = "z coordinate"; + // read LOOKUP_TABLE line + std::getline(infile, line); - } + ss.clear(); + ss.str(line); + ss >> name >> numcomp >> Npoints >> type; + // read data + std::cout << "begin reading data" << std::endl; + std::vector data(Npoints); + for (unsigned int i = 0; i < Npoints; i++) + { + infile >> data[i]; + } + std::cout << " done" << std::endl; - std::cout << "Construct PField '" << name << "'" << std::endl; - scalar_field.push_back( PField( name, var_name, var_description, mesh, data, 0.0) ); - std::cout << " done" << std::endl; + // construct field + std::vector var_name(DIM); + std::vector var_description(DIM); + if (DIM >= 2) + { + var_name[0] = "x"; + var_description[0] = "x coordinate"; + var_name[1] = "y"; + var_description[1] = "y coordinate"; } + if (DIM >= 3) + { + var_name[2] = "z"; + var_description[2] = "z coordinate"; + } + + std::cout << "Construct PField '" << name << "'" << std::endl; + scalar_field.push_back(PField(name, + var_name, + var_description, + mesh, + data, + 0.0)); + std::cout << " done" << std::endl; } } - - infile.close(); } - PField& find_scalar_field(std::string name) + infile.close(); + } + + PField & + find_scalar_field(std::string name) + { + for (unsigned int i = 0; i < scalar_field.size(); i++) { - for( unsigned int i=0; i - class Coordinate + + /// A class for a coordinate, templated by dimension + /// This is a possible option anyplace 'Coordinate' class template is used + /// but it is not the only option. Any class that implements + /// 'Coordinate::operator[]()' should work + + template + class Coordinate + { + double _coord[DIM]; + + public: + int + size() const { - double _coord[DIM]; - public: - - - int size() const - { - return DIM; - } - - double& operator[](int i) - { - return _coord[i]; - } - - const double& operator[](int i) const - { - return _coord[i]; - } - - template< int D> friend std::ostream &operator<<(std::ostream &outstream, const Coordinate &coord); - - }; - - template std::ostream &operator<<(std::ostream &outstream, const Coordinate &coord) + return DIM; + } + + double & + operator[](int i) { - for( int i=0; i + friend std::ostream & + operator<<(std::ostream &outstream, const Coordinate &coord); + }; + + template + std::ostream & + operator<<(std::ostream &outstream, const Coordinate &coord) + { + for (int i = 0; i < coord.size(); i++) + { + outstream << coord[i]; + if (i < coord.size() - 1) + outstream << " "; + } + return outstream; + } +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfield/Mesh.hh b/include/IntegrationTools/pfield/Mesh.hh index 173da9e4c..016972cbc 100644 --- a/include/IntegrationTools/pfield/Mesh.hh +++ b/include/IntegrationTools/pfield/Mesh.hh @@ -2,799 +2,855 @@ #ifndef Mesh_HH #define Mesh_HH -#include -#include -#include -#include - #include "../datastruc/Bin.hh" #include "../pfunction/PFuncBase.hh" +#include "./interpolation/Hexahedron.hh" #include "./interpolation/Interpolator.hh" #include "./interpolation/Quad.hh" -#include "./interpolation/Hexahedron.hh" +#include +#include +#include +#include namespace PRISMS { - inline void construct_basis_function( PFuncBase >, double>* &bfunc, const std::string &name) - { - if( name == "Quad") - { - bfunc = new Quad(); - } - else - { - std::cout << "Error in construct_basis_function (2D): unknown name: " << name << std::endl; - exit(1); - } - } - - inline void construct_basis_function( PFuncBase >, double>* &bfunc, const std::string &name) - { - if( name == "Hexahedron") - { - bfunc = new Hexahedron(); - } - else - { - std::cout << "Error in construct_basis_function (3D): unknown name: " << name << std::endl; - exit(1); - } - } + inline void + construct_basis_function(PFuncBase>, double> *&bfunc, + const std::string &name) + { + if (name == "Quad") + { + bfunc = new Quad(); + } + else + { + std::cout << "Error in construct_basis_function (2D): unknown name: " << name + << std::endl; + exit(1); + } + } + + inline void + construct_basis_function(PFuncBase>, double> *&bfunc, + const std::string &name) + { + if (name == "Hexahedron") + { + bfunc = new Hexahedron(); + } + else + { + std::cout << "Error in construct_basis_function (3D): unknown name: " << name + << std::endl; + exit(1); + } + } + + template + void + construct_interpolating_functions( + std::vector *> &interp, + const std::string &name, + unsigned long int cell, + PFuncBase>, double> *bfunc_ptr, + const std::vector &cell_node, + const std::vector> &node) + { + if (name == "Quad") + { + Interpolator *interp_ptr; + + // std::cout << "cell nodes: " << cell_node[0] << " " << cell_node[2] << + // std::endl; + + PRISMS::Coordinate<2> dim; + dim[0] = node[cell_node[2]][0] - node[cell_node[0]][0]; + dim[1] = node[cell_node[2]][1] - node[cell_node[0]][1]; + + // QuadValues(const Coordinate &node, const Coordinate &dim, int node_index) + for (int j = 0; j < 4; j++) + { + interp.push_back(interp_ptr); + interp.back() = new PRISMS::QuadValues(cell_node[j], + cell, + bfunc_ptr, + node[cell_node[j]], + dim, + j); + } + } + else + { + std::cout << "Error in construct_interpolating_function (2D): unknown name: " + << name << std::endl; + exit(1); + } + } + + template + void + construct_interpolating_functions( + std::vector *> &interp, + const std::string &name, + unsigned long int cell, + PFuncBase>, double> *bfunc_ptr, + const std::vector &cell_node, + const std::vector> &node) + { + if (name == "Hexahedron") + { + Interpolator *interp_ptr; + + PRISMS::Coordinate<3> dim; + dim[0] = node[cell_node[6]][0] - node[cell_node[0]][0]; + dim[1] = node[cell_node[6]][1] - node[cell_node[0]][1]; + dim[2] = node[cell_node[6]][2] - node[cell_node[0]][2]; + + // QuadValues(const Coordinate &node, const Coordinate &dim, int node_index) + for (int j = 0; j < 8; j++) + { + interp.push_back(interp_ptr); + interp.back() = new PRISMS::HexahedronValues(cell_node[j], + cell, + bfunc_ptr, + node[cell_node[j]], + dim, + j); + } + } + else + { + std::cout << "Error in construct_interpolating_function (3D): unknown name: " + << name << std::endl; + exit(1); + } + } + + /// A template class for a finite element mesh + /// Needs: Coordinate::operator[]() for use in Bin + /// + template + class Mesh + { + // min and max coordinate of cuboid surrounding the body + PRISMS::Coordinate _min; + PRISMS::Coordinate _max; + + /// Vector of nodal coordinates + /// nodal values live in 'Field' class + /// + std::vector> _node; - template< class Coordinate> - void construct_interpolating_functions( std::vector* > &interp, - const std::string &name, - unsigned long int cell, - PFuncBase >, double>* bfunc_ptr, - const std::vector &cell_node, - const std::vector > &node) - { - if( name == "Quad") - { - Interpolator* interp_ptr; + /// array containing interpolating functions: + /// owns the interpolating functions + /// interpolating functions contain basis function / element info, + /// these point to _bfunc pfunctions which are used to evaluate + /// + std::vector *> _interp; - //std::cout << "cell nodes: " << cell_node[0] << " " << cell_node[2] << std::endl; + /// array containing PFunctions evaluated by interpolating functions + /// owns the pfunctions, which are pointed to by the interpolating functions + /// !!! do not modify after initial construction or pointers will be messed up !!! + /// + std::vector>, double> *> _bfunc; + /// bin of interpolating functions (this might be updated to be either Element or + /// Spline Bins) + /// + Bin *, Coordinate> _bin; - PRISMS::Coordinate<2> dim; - dim[0] = node[ cell_node[2]][0] - node[ cell_node[0]][0]; - dim[1] = node[ cell_node[2]][1] - node[ cell_node[0]][1]; + public: + // still need a constructor + Mesh() {}; - // QuadValues(const Coordinate &node, const Coordinate &dim, int node_index) - for( int j=0; j<4; j++) - { - interp.push_back(interp_ptr); - interp.back() = new PRISMS::QuadValues(cell_node[j], cell, bfunc_ptr, node[ cell_node[j]], dim, j); - } - } - else - { - std::cout << "Error in construct_interpolating_function (2D): unknown name: " << name << std::endl; - exit(1); - } - } - - template< class Coordinate> - void construct_interpolating_functions( std::vector* > &interp, - const std::string &name, - unsigned long int cell, - PFuncBase >, double>* bfunc_ptr, - const std::vector &cell_node, - const std::vector > &node) + ~Mesh() { - if( name == "Hexahedron") + for (unsigned int i = 0; i < _interp.size(); i++) { - Interpolator* interp_ptr; - - PRISMS::Coordinate<3> dim; - dim[0] = node[ cell_node[6]][0] - node[ cell_node[0]][0]; - dim[1] = node[ cell_node[6]][1] - node[ cell_node[0]][1]; - dim[2] = node[ cell_node[6]][2] - node[ cell_node[0]][2]; - - // QuadValues(const Coordinate &node, const Coordinate &dim, int node_index) - for( int j=0; j<8; j++) - { - interp.push_back(interp_ptr); - interp.back() = new PRISMS::HexahedronValues(cell_node[j], cell, bfunc_ptr, node[ cell_node[j]], dim, j); - } + delete _interp[i]; } - else + + for (unsigned int i = 0; i < _bfunc.size(); i++) { - std::cout << "Error in construct_interpolating_function (3D): unknown name: " << name << std::endl; - exit(1); + delete _bfunc[i]; } - } - - + }; - /// A template class for a finite element mesh - /// Needs: Coordinate::operator[]() for use in Bin - /// - template - class Mesh + // reads vtk file through 'CELL_TYPES' and then returns + void + read_vtk(std::ifstream &infile) { - // min and max coordinate of cuboid surrounding the body - PRISMS::Coordinate _min; - PRISMS::Coordinate _max; - - /// Vector of nodal coordinates - /// nodal values live in 'Field' class - /// - std::vector > _node; - - /// array containing interpolating functions: - /// owns the interpolating functions - /// interpolating functions contain basis function / element info, - /// these point to _bfunc pfunctions which are used to evaluate - /// - std::vector* > _interp; - - /// array containing PFunctions evaluated by interpolating functions - /// owns the pfunctions, which are pointed to by the interpolating functions - /// !!! do not modify after initial construction or pointers will be messed up !!! - /// - std::vector< PFuncBase >, double>* > _bfunc; - - /// bin of interpolating functions (this might be updated to be either Element or Spline Bins) - /// - Bin*, Coordinate > _bin; - - public: - - // still need a constructor - Mesh(){}; - - ~Mesh() - { - for( unsigned int i=0; i<_interp.size(); i++) - { - delete _interp[i]; - } - - for( unsigned int i=0; i<_bfunc.size(); i++) - { - delete _bfunc[i]; - } - }; + bool mesh_as_points = true; + std::vector x_coord, y_coord, z_coord; - // reads vtk file through 'CELL_TYPES' and then returns - void read_vtk(std::ifstream &infile) - { - bool mesh_as_points = true; - std::vector x_coord, y_coord, z_coord; + std::istringstream ss; + std::string line, str, type; - std::istringstream ss; - std::string line, str, type; + unsigned long int uli_dummy; + double d_dummy; - unsigned long int uli_dummy; - double d_dummy; + unsigned long int Npoints, Ncells, Ncell_numbers; + std::vector cell_node; - unsigned long int Npoints, Ncells, Ncell_numbers; - std::vector cell_node; + PRISMS::Coordinate _coord; - PRISMS::Coordinate _coord; + while (!infile.eof()) + { + std::getline(infile, line); + // std::cout << "line: " << line << std::endl; - while(!infile.eof()) + if (line[0] == 'P') { - std::getline( infile, line); - //std::cout << "line: " << line << std::endl; - - if( line[0] == 'P') + // read POINTS info: + // POINTS # type + // x y z + // x y z + // ... + if (line.size() > 5 && line.substr(0, 6) == "POINTS") { - // read POINTS info: - // POINTS # type - // x y z - // x y z - // ... - if( line.size() > 5 && line.substr(0,6) == "POINTS") + // read header line + // std::cout << line << "\n"; + ss.clear(); + ss.str(line); + ss >> str >> Npoints >> type; + + // read points + std::vector> value(DIM); + std::vector> hist(DIM); + + std::cout << "Read POINTS: " << Npoints << std::endl; + _node.reserve(Npoints); + std::cout << " reserve OK" << std::endl; + for (unsigned int i = 0; i < Npoints; i++) { - // read header line - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); - ss >> str >> Npoints >> type; - - // read points - std::vector< std::vector > value(DIM); - std::vector< std::vector > hist(DIM); - - std::cout << "Read POINTS: " << Npoints << std::endl; - _node.reserve(Npoints); - std::cout << " reserve OK" << std::endl; - for( unsigned int i=0; i> _coord[0] >> _coord[1] >> d_dummy; - //std::cout << _coord[0] << " " << _coord[1] << " " << d_dummy << std::endl; - } - else if( DIM == 3) - { - infile >> _coord[0] >> _coord[1] >> _coord[2]; - //std::cout << _coord[0] << " " << _coord[1] << " " << _coord[3] << std::endl; - } - - for( int j=0; j> _coord[0] >> _coord[1] >> d_dummy; + // std::cout << _coord[0] << " " << _coord[1] << " " << d_dummy + // << std::endl; + } + else if (DIM == 3) + { + infile >> _coord[0] >> _coord[1] >> _coord[2]; + // std::cout << _coord[0] << " " << _coord[1] << " " << + // _coord[3] << std::endl; } - std::cout << " done" << std::endl; - // create bins - std::vector min; - std::vector N; - std::vector incr; + for (int j = 0; j < DIM; j++) + add_once(value[j], hist[j], _coord[j]); + _node.push_back(_coord); + } + std::cout << " done" << std::endl; - std::cout << "Determine Body size" << std::endl; - for( int j=0; j min; + std::vector N; + std::vector incr; + + std::cout << "Determine Body size" << std::endl; + for (int j = 0; j < DIM; j++) + { + for (unsigned int i = 1; i < hist[j].size(); i++) { - for( unsigned int i=1; i*, Coordinate>(min, incr, N); - std::cout << " done" << std::endl; + std::cout << " done" << std::endl; + std::cout << "Initialize Bin" << std::endl; + _bin = Bin *, Coordinate>(min, incr, N); + std::cout << " done" << std::endl; } - - } - else if( line[0] == 'C') + else if (line[0] == 'C') { - - if( line.size() > 4 && line.substr(0,5) == "CELLS") + if (line.size() > 4 && line.substr(0, 5) == "CELLS") { - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); + // std::cout << line << "\n"; + ss.clear(); + ss.str(line); - ss >> str >> Ncells >> Ncell_numbers; + ss >> str >> Ncells >> Ncell_numbers; - PFuncBase >, double>* bfunc_ptr; - _bfunc.push_back( bfunc_ptr); + PFuncBase>, double> *bfunc_ptr; + _bfunc.push_back(bfunc_ptr); - if( DIM == 2) + if (DIM == 2) { - // add Quad basis function - _interp.reserve(Ncells*4); - construct_basis_function(_bfunc.back(), "Quad"); + // add Quad basis function + _interp.reserve(Ncells * 4); + construct_basis_function(_bfunc.back(), "Quad"); } - else if( DIM == 3) + else if (DIM == 3) { - // add Hexahedron basis function - _interp.reserve(Ncells*8); - construct_basis_function(_bfunc.back(), "Hexahedron"); + // add Hexahedron basis function + _interp.reserve(Ncells * 8); + construct_basis_function(_bfunc.back(), "Hexahedron"); } - bfunc_ptr = _bfunc.back(); + bfunc_ptr = _bfunc.back(); - std::cout << "Read CELLS: " << Ncells << std::endl; - for( unsigned int i=0; i> uli_dummy; + infile >> uli_dummy; - cell_node.resize(uli_dummy); - for( unsigned int j=0; j> cell_node[j]; + infile >> cell_node[j]; } - //std::cout << cell_node[0] << " " << cell_node[1] << " " << cell_node[2] << " " << cell_node[3] << std::endl; + // std::cout << cell_node[0] << " " << cell_node[1] << " " << + // cell_node[2] << " " << cell_node[3] << std::endl; - // create interpolator - if( DIM == 2) + // create interpolator + if (DIM == 2) { - construct_interpolating_functions(_interp, "Quad", i, bfunc_ptr, cell_node, _node); + construct_interpolating_functions(_interp, + "Quad", + i, + bfunc_ptr, + cell_node, + _node); } - else if(DIM == 3) + else if (DIM == 3) { - construct_interpolating_functions(_interp, "Hexahedron", i, bfunc_ptr, cell_node, _node); + construct_interpolating_functions(_interp, + "Hexahedron", + i, + bfunc_ptr, + cell_node, + _node); } - } - std::cout << " done" << std::endl; + std::cout << " done" << std::endl; - // bin interpolators - std::cout << "Bin interpolating functions" << std::endl; + // bin interpolators + std::cout << "Bin interpolating functions" << std::endl; - for( unsigned int i=0; i<_interp.size(); i++) + for (unsigned int i = 0; i < _interp.size(); i++) { - _bin.add_range(_interp[i], _interp[i]->min(), _interp[i]->max()); + _bin.add_range(_interp[i], _interp[i]->min(), _interp[i]->max()); } - std::cout << " done max_bin_size: " << _bin.max_size() << std::endl; - - + std::cout << " done max_bin_size: " << _bin.max_size() << std::endl; } - else if( line.size() > 9 && line.substr(0,10) == "CELL_TYPES") + else if (line.size() > 9 && line.substr(0, 10) == "CELL_TYPES") { - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); + // std::cout << line << "\n"; + ss.clear(); + ss.str(line); - //std::cout << "ss.str()" << ss.str() << std::endl; - ss >> str >> Ncells; + // std::cout << "ss.str()" << ss.str() << std::endl; + ss >> str >> Ncells; - for( unsigned int i=0; i> uli_dummy; + infile >> uli_dummy; - if( uli_dummy != 9 && uli_dummy != 12) + if (uli_dummy != 9 && uli_dummy != 12) { - std::cout << "Error reading CELL_TYPES: CELL TYPE != 9 && != 12" << std::endl; - std::cout << " CELL TYPE: " << uli_dummy << std::endl; - exit(1); + std::cout << "Error reading CELL_TYPES: CELL TYPE != 9 && != 12" + << std::endl; + std::cout << " CELL TYPE: " << uli_dummy << std::endl; + exit(1); } } - return; + return; } - } - if( line[0] == 'X') + if (line[0] == 'X') { - mesh_as_points = false; - - // read X_COORDINATES info: - // X_COORDINATES # type - // x - // x - // ... - if( line.size() > 12 && line.substr(0,13) == "X_COORDINATES") + mesh_as_points = false; + + // read X_COORDINATES info: + // X_COORDINATES # type + // x + // x + // ... + if (line.size() > 12 && line.substr(0, 13) == "X_COORDINATES") { - // read header line - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); - ss >> str >> Npoints >> type; - - std::cout << "Read X_COORDINATES: " << Npoints << std::endl; - _node.reserve(Npoints); - std::cout << " reserve OK" << std::endl; - for( unsigned int i=0; i> str >> Npoints >> type; + + std::cout << "Read X_COORDINATES: " << Npoints << std::endl; + _node.reserve(Npoints); + std::cout << " reserve OK" << std::endl; + for (unsigned int i = 0; i < Npoints; i++) { - double temp_coord; - - infile >> temp_coord; + double temp_coord; - x_coord.push_back(temp_coord); + infile >> temp_coord; + x_coord.push_back(temp_coord); } } } - if( line[0] == 'Y') + if (line[0] == 'Y') { - mesh_as_points = false; - - // read Y_COORDINATES info: - // Y_COORDINATES # type - // y - // y - // ... - if( line.size() > 12 && line.substr(0,13) == "Y_COORDINATES") + mesh_as_points = false; + + // read Y_COORDINATES info: + // Y_COORDINATES # type + // y + // y + // ... + if (line.size() > 12 && line.substr(0, 13) == "Y_COORDINATES") { - - - // read header line - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); - ss >> str >> Npoints >> type; - - // read points - std::vector< std::vector > value(DIM); - std::vector< std::vector > hist(DIM); - - std::cout << "Read Y_COORDINATES: " << Npoints << std::endl; - _node.reserve(Npoints); - std::cout << " reserve OK" << std::endl; - for( unsigned int i=0; i> str >> Npoints >> type; + + // read points + std::vector> value(DIM); + std::vector> hist(DIM); + + std::cout << "Read Y_COORDINATES: " << Npoints << std::endl; + _node.reserve(Npoints); + std::cout << " reserve OK" << std::endl; + for (unsigned int i = 0; i < Npoints; i++) { - double temp_coord; - - infile >> temp_coord; + double temp_coord; - y_coord.push_back(temp_coord); + infile >> temp_coord; + y_coord.push_back(temp_coord); } } } - if( line[0] == 'Z') + if (line[0] == 'Z') { - mesh_as_points = false; - - // read Z_COORDINATES info: - // Z_COORDINATES # type - // z - // z - // ... - if( line.size() > 12 && line.substr(0,13) == "Z_COORDINATES") + mesh_as_points = false; + + // read Z_COORDINATES info: + // Z_COORDINATES # type + // z + // z + // ... + if (line.size() > 12 && line.substr(0, 13) == "Z_COORDINATES") { - - - // read header line - //std::cout << line << "\n"; - ss.clear(); - ss.str(line); - ss >> str >> Npoints >> type; - - // read points - - std::cout << "Read Z_COORDINATES: " << Npoints << std::endl; - _node.reserve(Npoints); - std::cout << " reserve OK" << std::endl; - for( unsigned int i=0; i> str >> Npoints >> type; + + // read points + + std::cout << "Read Z_COORDINATES: " << Npoints << std::endl; + _node.reserve(Npoints); + std::cout << " reserve OK" << std::endl; + for (unsigned int i = 0; i < Npoints; i++) { - double temp_coord; + double temp_coord; - infile >> temp_coord; - - z_coord.push_back(temp_coord); + infile >> temp_coord; + z_coord.push_back(temp_coord); } } } } - if (!mesh_as_points){ - std::vector< std::vector > value(DIM); - std::vector< std::vector > hist(DIM); - - for (unsigned int i=0; i 2){ - _coord[2] = z_coord.at(k); + if (!mesh_as_points) + { + std::vector> value(DIM); + std::vector> hist(DIM); + + for (unsigned int i = 0; i < x_coord.size(); i++) + { + for (unsigned int j = 0; j < y_coord.size(); j++) + { + for (unsigned int k = 0; k < z_coord.size(); k++) + { + _coord[0] = x_coord.at(i); + _coord[1] = y_coord.at(j); + if (DIM > 2) + { + _coord[2] = z_coord.at(k); } - for( int m=0; m min; - std::vector N; - std::vector incr; + // create bins + std::vector min; + std::vector N; + std::vector incr; - std::cout << "Determine Body size" << std::endl; - for( int j=0; j*, Coordinate>(min, incr, N); - std::cout << " done" << std::endl; - - // Now add the cell data - unsigned int Ncells = (x_coord.size()-1) * (y_coord.size()-1); - if (DIM > 2){ - Ncells *= (z_coord.size()-1); - } + std::sort(value[j].begin(), value[j].end()); + // std::cout << "j: " << j << " back(): " << value[j].back() << std::endl; + min.push_back(value[j][0]); + N.push_back(value[j].size()); + incr.push_back((value[j].back() - value[j][0]) / (1.0 * N.back())); - PFuncBase >, double>* bfunc_ptr; - _bfunc.push_back( bfunc_ptr); + // get min and max surrounding coordinates + _min[j] = value[j][0]; + _max[j] = value[j].back(); - if( DIM == 2) - { - // add Quad basis function - _interp.reserve(Ncells*4); - construct_basis_function(_bfunc.back(), "Quad"); - } - else if( DIM == 3) - { - // add Hexahedron basis function - _interp.reserve(Ncells*8); - construct_basis_function(_bfunc.back(), "Hexahedron"); - } - bfunc_ptr = _bfunc.back(); + // for short term, expand bin to avoid edge issues + min[j] -= incr[j]; + N[j] += 2; + } + std::cout << " Min Coordinate: "; + for (int j = 0; j < DIM; j++) + std::cout << _min[j] << " "; + std::cout << std::endl; + std::cout << " Max Coordinate: "; + for (int j = 0; j < DIM; j++) + std::cout << _max[j] << " "; + std::cout << std::endl; + + std::cout << " done" << std::endl; + + std::cout << "Initialize Bin" << std::endl; + _bin = Bin *, Coordinate>(min, incr, N); + std::cout << " done" << std::endl; + + // Now add the cell data + unsigned int Ncells = (x_coord.size() - 1) * (y_coord.size() - 1); + if (DIM > 2) + { + Ncells *= (z_coord.size() - 1); + } - std::cout << "Read CELLS: " << Ncells << std::endl; + PFuncBase>, double> *bfunc_ptr; + _bfunc.push_back(bfunc_ptr); - unsigned int uli_dummy; - if (DIM > 2){ - uli_dummy = 8; - } - else { - uli_dummy = 4; - } - for( unsigned int i=0; i 2) { - construct_interpolating_functions(_interp, "Quad", i, bfunc_ptr, cell_node, _node); + uli_dummy = 8; } - else if(DIM == 3) + else { - construct_interpolating_functions(_interp, "Hexahedron", i, bfunc_ptr, cell_node, _node); + uli_dummy = 4; } + for (unsigned int i = 0; i < Ncells; i++) + { + cell_node.resize(uli_dummy); + for (unsigned int j = 0; j < uli_dummy; j++) + { + cell_node[j] = i * uli_dummy + j; + } - } - std::cout << " done" << std::endl; - - // bin interpolators - std::cout << "Bin interpolating functions" << std::endl; - std::cout << "num nodes: " << _node.size() << std::endl; - for( unsigned int i=0; i<_interp.size(); i++) - { - std::cout << "interp: " << _interp[i] << " " << _interp[i]->min() << " " << _interp[i]->max() << std::endl; - _bin.add_range(_interp[i], _interp[i]->min(), _interp[i]->max()); - } - std::cout << " done max_bin_size: " << _bin.max_size() << std::endl; - - } + if (DIM == 2) + { + double temp = cell_node[2]; + cell_node[2] = cell_node[3]; + cell_node[3] = temp; + } + // std::cout << cell_node[0] << " " << cell_node[1] << " " << cell_node[2] + // << " " << cell_node[3] << std::endl; -} + // create interpolator + if (DIM == 2) + { + construct_interpolating_functions(_interp, + "Quad", + i, + bfunc_ptr, + cell_node, + _node); + } + else if (DIM == 3) + { + construct_interpolating_functions(_interp, + "Hexahedron", + i, + bfunc_ptr, + cell_node, + _node); + } + } + std::cout << " done" << std::endl; - void min( Coordinate &coord) - { - for( int i=0; imin() << " " + << _interp[i]->max() << std::endl; + _bin.add_range(_interp[i], _interp[i]->min(), _interp[i]->max()); + } + std::cout << " done max_bin_size: " << _bin.max_size() << std::endl; } + } - void max( Coordinate &coord) - { - for( int i=0; i &bfunc, std::vector &node_index, int &s) - { - std::vector* > &bin = _bin.contents(coord); - s = bin.size(); + // Set 'bfunc' to evaluated basis functions at 'coord', + // 'node_index' to node indices for each basis function, + // and 's' is the length (number of basis functions) + // - 'bfunc' and 'node_index' are not resized, they must be big enough + // + void + basis_functions(const Coordinate &coord, + std::vector &bfunc, + std::vector &node_index, + int &s) + { + std::vector *> &bin = _bin.contents(coord); + s = bin.size(); - int i=0; - unsigned long int element; + int i = 0; + unsigned long int element; - for( i=0; i &bfunc, std::vector &node_index, int &s) - { - //std::cout << "begin Mesh::grad_basis_functions()" << std::endl; - std::vector* > &bin = _bin.contents(coord); - s = bin.size(); + // Set 'bfunc' to evaluated grad basis functions at coord, and 's' is the length + void + grad_basis_functions(const Coordinate &coord, + int di, + std::vector &bfunc, + std::vector &node_index, + int &s) + { + // std::cout << "begin Mesh::grad_basis_functions()" << std::endl; + std::vector *> &bin = _bin.contents(coord); + s = bin.size(); - int i=0; - unsigned long int element; + int i = 0; + unsigned long int element; - for( i=0; i &bfunc, std::vector &node_index, int &s) - { - std::vector* > &bin = _bin.contents(coord); - s = bin.size(); + // Set 'bfunc' to evaluated hess basis functions at coord, and 's' is the length + void + hess_basis_functions(Coordinate coord, + int di, + int dj, + std::vector &bfunc, + std::vector &node_index, + int &s) + { + std::vector *> &bin = _bin.contents(coord); + s = bin.size(); - int i=0; - unsigned long int element; + int i = 0; + unsigned long int element; - for( i=0; i &list, std::vector &hist, double val) + { + // std::cout << "begin add_once()" << std::endl; - void add_once( std::vector &list, std::vector &hist, double val) + for (unsigned int i = 0; i < list.size(); i++) { - //std::cout << "begin add_once()" << std::endl; - - for( unsigned int i=0; i' for x, y, z - /// FieldType is the datatype for the field, - /// for instance 'double' for temperature or 'std::vector' for vector displacement - /// - /// A field consists of a pointer to a mesh, a list field values (at mesh nodes) - template - class PField : public PFuncBase - { - public: - - typedef typename PFuncBase::size_type size_type; - - // pointer to a Mesh that lives in a Body - Mesh *_mesh; - - // array of field values at mesh nodes - std::vector _field; - - FieldType _zero; - - FieldType _val; - std::vector _grad_val; - std::vector< std::vector > _hess_val; - - // ---------------------------------------------------------- - // Constructors - // PField(); - - PField( const std::string &name, - const std::vector &var_name, - const std::vector &var_description, - Mesh &mesh, - const std::vector &field, - const FieldType &zero) : - PFuncBase(name, var_name, var_description), - _mesh(&mesh), - _field(field), - _zero(zero) - - { - int max = mesh.max_bin_size(); - _bfunc.resize(max); - _node_index.resize(max); - _grad_val.resize(DIM); - _hess_val.resize(DIM); - for( int i=0; i* clone() const - { - return new PField(*this); - } - - - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - FieldType operator()(const Coordinate &coord); - FieldType grad(const Coordinate &coord, size_type di); - FieldType hess(const Coordinate &coord, size_type di, size_type dj); - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - void eval(const Coordinate &coord); - void eval_grad(const Coordinate &coord); - void eval_hess(const Coordinate &coord); - - FieldType operator()() const; - FieldType grad(size_type di) const; - FieldType hess(size_type di, size_type dj) const; - - // PField unique members ------------------------------------------ - - private: - - // temporary vector - std::vector _bfunc; - std::vector _node_index; - int _Nbfunc; - - - }; - - template - FieldType PField::operator()( const Coordinate &coord) - { - //std::cout << "begin calc" << std::endl; - - //std::cout << "coord: "; - //for( int i=0; i - FieldType PField::grad( const Coordinate &coord, size_type di) - { - //std::cout << "begin PField::grad()" << std::endl; - // get evaluated basis functions - (*_mesh).grad_basis_functions(coord, di, _bfunc, _node_index, _Nbfunc); - - // sum them - _grad_val[di] = _zero; - for( int i=0; i<_Nbfunc; i++) - _grad_val[di] += _bfunc[i]*_field[_node_index[i]]; - return _grad_val[di]; - //std::cout << "finish PField::grad()" << std::endl; - - } - - template - FieldType PField::hess( const Coordinate &coord, size_type di, size_type dj) - { - // get evaluated basis functions - (*_mesh).hess_basis_functions(coord, di, dj, _bfunc, _node_index, _Nbfunc); - - // sum them - _hess_val[di][dj] = _zero; - for( int i=0; i<_Nbfunc; i++) - _hess_val[di][dj] += _bfunc[i]*_field[_node_index[i]]; - return _hess_val[di][dj]; - } - - template - void PField::eval( const Coordinate &coord) - { - (*this)(coord); - } + /// A class for a field + /// + /// Varcontainer contains variables, for instance 'std::vector' for x, y, z + /// FieldType is the datatype for the field, + /// for instance 'double' for temperature or 'std::vector' for vector + /// displacement + /// + /// A field consists of a pointer to a mesh, a list field values (at mesh nodes) + template + class PField : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + // pointer to a Mesh that lives in a Body + Mesh *_mesh; + + // array of field values at mesh nodes + std::vector _field; + + FieldType _zero; + + FieldType _val; + std::vector _grad_val; + std::vector> _hess_val; + + // ---------------------------------------------------------- + // Constructors + // PField(); + + PField(const std::string &name, + const std::vector &var_name, + const std::vector &var_description, + Mesh &mesh, + const std::vector &field, + const FieldType &zero) + : PFuncBase(name, var_name, var_description) + , _mesh(&mesh) + , _field(field) + , _zero(zero) - template - void PField::eval_grad( const Coordinate &coord) { - for( int di=0; di - void PField::eval_hess( const Coordinate &coord) - { - for( int di=0; di - FieldType PField::operator()() const + PField * + clone() const { - return _val; + return new PField(*this); } - template - FieldType PField::grad( size_type di) const - { - return _grad_val[di]; - } - - template - FieldType PField::hess( size_type di, size_type dj) const - { - return _hess_val[di][dj]; - } - -} - + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + FieldType + operator()(const Coordinate &coord); + FieldType + grad(const Coordinate &coord, size_type di); + FieldType + hess(const Coordinate &coord, size_type di, size_type dj); + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + void + eval(const Coordinate &coord); + void + eval_grad(const Coordinate &coord); + void + eval_hess(const Coordinate &coord); + + FieldType + operator()() const; + FieldType + grad(size_type di) const; + FieldType + hess(size_type di, size_type dj) const; + + // PField unique members ------------------------------------------ + + private: + // temporary vector + std::vector _bfunc; + std::vector _node_index; + int _Nbfunc; + }; + + template + FieldType + PField::operator()(const Coordinate &coord) + { + // std::cout << "begin calc" << std::endl; + + // std::cout << "coord: "; + // for( int i=0; i + FieldType + PField::grad(const Coordinate &coord, size_type di) + { + // std::cout << "begin PField::grad()" << std::endl; + // get evaluated basis functions + (*_mesh).grad_basis_functions(coord, di, _bfunc, _node_index, _Nbfunc); + + // sum them + _grad_val[di] = _zero; + for (int i = 0; i < _Nbfunc; i++) + _grad_val[di] += _bfunc[i] * _field[_node_index[i]]; + return _grad_val[di]; + // std::cout << "finish PField::grad()" << std::endl; + } + + template + FieldType + PField::hess(const Coordinate &coord, + size_type di, + size_type dj) + { + // get evaluated basis functions + (*_mesh).hess_basis_functions(coord, di, dj, _bfunc, _node_index, _Nbfunc); + + // sum them + _hess_val[di][dj] = _zero; + for (int i = 0; i < _Nbfunc; i++) + _hess_val[di][dj] += _bfunc[i] * _field[_node_index[i]]; + return _hess_val[di][dj]; + } + + template + void + PField::eval(const Coordinate &coord) + { + (*this)(coord); + } + + template + void + PField::eval_grad(const Coordinate &coord) + { + for (int di = 0; di < DIM; di++) + (*this).grad(coord, di); + } + + template + void + PField::eval_hess(const Coordinate &coord) + { + for (int di = 0; di < DIM; di++) + for (int dj = 0; dj < DIM; dj++) + (*this).hess(coord, di, dj); + } + + template + FieldType + PField::operator()() const + { + return _val; + } + + template + FieldType + PField::grad(size_type di) const + { + return _grad_val[di]; + } + + template + FieldType + PField::hess(size_type di, size_type dj) const + { + return _hess_val[di][dj]; + } + +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/pfield/interpolation/Hexahedron.hh b/include/IntegrationTools/pfield/interpolation/Hexahedron.hh index 78118ffda..5b54610f3 100644 --- a/include/IntegrationTools/pfield/interpolation/Hexahedron.hh +++ b/include/IntegrationTools/pfield/interpolation/Hexahedron.hh @@ -2,617 +2,653 @@ #ifndef Hexahedron_HH #define Hexahedron_HH +#include "../../pfunction/PFuncBase.hh" #include "../../pfunction/PSimpleBase.hh" #include "../../pfunction/PSimpleFunction.hh" -#include "../../pfunction/PFuncBase.hh" -#include "./Interpolator.hh" #include "../Coordinate.hh" +#include "./Interpolator.hh" namespace PRISMS { - class Hexahedron_f : public PSimpleBase< std::vector >, double> + class Hexahedron_f : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const { - double eval( const std::vector > &var) const - { - // var[0]: coordinate to be evaluated - // var[1]: nodal coordinate - // var[2]: element dimension - // var[3]: +/- 1 depending on which 'corner' of quad - // - // f = (1.0 - e0)*(1.0 - e1)*(1.0 - e2) - // e = var[3]*(var[0] - var[1])/var[2] - - return (1.0 - var[3][0]*(var[0][0] - var[1][0])/var[2][0])* - (1.0 - var[3][1]*(var[0][1] - var[1][1])/var[2][1])* - (1.0 - var[3][2]*(var[0][2] - var[1][2])/var[2][2]); - } + // var[0]: coordinate to be evaluated + // var[1]: nodal coordinate + // var[2]: element dimension + // var[3]: +/- 1 depending on which 'corner' of quad + // + // f = (1.0 - e0)*(1.0 - e1)*(1.0 - e2) + // e = var[3]*(var[0] - var[1])/var[2] + + return (1.0 - var[3][0] * (var[0][0] - var[1][0]) / var[2][0]) * + (1.0 - var[3][1] * (var[0][1] - var[1][1]) / var[2][1]) * + (1.0 - var[3][2] * (var[0][2] - var[1][2]) / var[2][2]); + } + + public: + Hexahedron_f() + { + this->_name = "Hexahedron_f"; + } - public: + Hexahedron_f * + clone() const + { + return new Hexahedron_f(*this); + } + }; + + class Hexahedron_grad_0 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return -var[3][0] * (1.0 - var[3][1] * (var[0][1] - var[1][1]) / var[2][1]) * + (1.0 - var[3][2] * (var[0][2] - var[1][2]) / var[2][2]) / var[2][0]; + } - Hexahedron_f() - { - this->_name = "Hexahedron_f"; - } + public: + Hexahedron_grad_0() + { + this->_name = "Hexahedron_grad_0"; + } - Hexahedron_f* clone() const - { - return new Hexahedron_f(*this); - } - }; - - class Hexahedron_grad_0 : public PSimpleBase< std::vector >, double> + Hexahedron_grad_0 * + clone() const { - double eval( const std::vector > &var) const - { - return -var[3][0]*(1.0 - var[3][1]*(var[0][1] - var[1][1])/var[2][1])* - (1.0 - var[3][2]*(var[0][2] - var[1][2])/var[2][2])/var[2][0]; - } + return new Hexahedron_grad_0(*this); + } + }; + + class Hexahedron_grad_1 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return -var[3][1] * (1.0 - var[3][0] * (var[0][0] - var[1][0]) / var[2][0]) * + (1.0 - var[3][2] * (var[0][2] - var[1][2]) / var[2][2]) / var[2][1]; + } - public: + public: + Hexahedron_grad_1() + { + this->_name = "Hexahedron_grad_1"; + } - Hexahedron_grad_0() - { - this->_name = "Hexahedron_grad_0"; - } + Hexahedron_grad_1 * + clone() const + { + return new Hexahedron_grad_1(*this); + } + }; + + class Hexahedron_grad_2 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return -var[3][2] * (1.0 - var[3][0] * (var[0][0] - var[1][0]) / var[2][0]) * + (1.0 - var[3][1] * (var[0][1] - var[1][1]) / var[2][1]) / var[2][2]; + } - Hexahedron_grad_0* clone() const - { - return new Hexahedron_grad_0(*this); - } - }; - - class Hexahedron_grad_1 : public PSimpleBase< std::vector >, double> + public: + Hexahedron_grad_2() { - double eval( const std::vector > &var) const - { - return -var[3][1]*(1.0 - var[3][0]*(var[0][0] - var[1][0])/var[2][0])* - (1.0 - var[3][2]*(var[0][2] - var[1][2])/var[2][2])/var[2][1]; - - } + this->_name = "Hexahedron_grad_2"; + } - public: + Hexahedron_grad_2 * + clone() const + { + return new Hexahedron_grad_2(*this); + } + }; + + class Hexahedron_hess_0_0 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return 0.0; + } - Hexahedron_grad_1() - { - this->_name = "Hexahedron_grad_1"; - } + public: + Hexahedron_hess_0_0() + { + this->_name = "Hexahedron_hess_0_0"; + } - Hexahedron_grad_1* clone() const - { - return new Hexahedron_grad_1(*this); - } - }; - - class Hexahedron_grad_2 : public PSimpleBase< std::vector >, double> + Hexahedron_hess_0_0 * + clone() const { - double eval( const std::vector > &var) const - { - return -var[3][2]*(1.0 - var[3][0]*(var[0][0] - var[1][0])/var[2][0])* - (1.0 - var[3][1]*(var[0][1] - var[1][1])/var[2][1])/var[2][2]; - - } + return new Hexahedron_hess_0_0(*this); + } + }; + + class Hexahedron_hess_0_1 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][0] * var[3][1] / var[2][0] / var[2][1]; + } - public: + public: + Hexahedron_hess_0_1() + { + this->_name = "Hexahedron_hess_0_1"; + } - Hexahedron_grad_2() - { - this->_name = "Hexahedron_grad_2"; - } + Hexahedron_hess_0_1 * + clone() const + { + return new Hexahedron_hess_0_1(*this); + } + }; + + class Hexahedron_hess_0_2 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][0] * var[3][2] / var[2][0] / var[2][2]; + } - Hexahedron_grad_2* clone() const - { - return new Hexahedron_grad_2(*this); - } - }; - - class Hexahedron_hess_0_0 : public PSimpleBase< std::vector >, double> + public: + Hexahedron_hess_0_2() { - double eval( const std::vector > &var) const - { - return 0.0; - } + this->_name = "Hexahedron_hess_0_2"; + } - public: + Hexahedron_hess_0_2 * + clone() const + { + return new Hexahedron_hess_0_2(*this); + } + }; + + class Hexahedron_hess_1_0 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][1] * var[3][0] / var[2][1] / var[2][0]; + } - Hexahedron_hess_0_0() - { - this->_name = "Hexahedron_hess_0_0"; - } + public: + Hexahedron_hess_1_0() + { + this->_name = "Hexahedron_hess_1_0"; + } - Hexahedron_hess_0_0* clone() const - { - return new Hexahedron_hess_0_0(*this); - } - }; - - class Hexahedron_hess_0_1 : public PSimpleBase< std::vector >, double> + Hexahedron_hess_1_0 * + clone() const { - double eval( const std::vector > &var) const - { - return var[3][0]*var[3][1]/var[2][0]/var[2][1]; - } + return new Hexahedron_hess_1_0(*this); + } + }; + + class Hexahedron_hess_1_1 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return 0.0; + } - public: + public: + Hexahedron_hess_1_1() + { + this->_name = "Hexahedron_hess_1_1"; + } - Hexahedron_hess_0_1() - { - this->_name = "Hexahedron_hess_0_1"; - } + Hexahedron_hess_1_1 * + clone() const + { + return new Hexahedron_hess_1_1(*this); + } + }; + + class Hexahedron_hess_1_2 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][1] * var[3][2] / var[2][1] / var[2][2]; + } - Hexahedron_hess_0_1* clone() const - { - return new Hexahedron_hess_0_1(*this); - } - }; - - class Hexahedron_hess_0_2 : public PSimpleBase< std::vector >, double> + public: + Hexahedron_hess_1_2() { - double eval( const std::vector > &var) const - { - return var[3][0]*var[3][2]/var[2][0]/var[2][2]; - } + this->_name = "Hexahedron_hess_1_2"; + } - public: + Hexahedron_hess_1_2 * + clone() const + { + return new Hexahedron_hess_1_2(*this); + } + }; + + class Hexahedron_hess_2_0 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][2] * var[3][0] / var[2][2] / var[2][0]; + } - Hexahedron_hess_0_2() - { - this->_name = "Hexahedron_hess_0_2"; - } + public: + Hexahedron_hess_2_0() + { + this->_name = "Hexahedron_hess_2_0"; + } - Hexahedron_hess_0_2* clone() const - { - return new Hexahedron_hess_0_2(*this); - } - }; - - class Hexahedron_hess_1_0 : public PSimpleBase< std::vector >, double> + Hexahedron_hess_2_0 * + clone() const { - double eval( const std::vector > &var) const - { - return var[3][1]*var[3][0]/var[2][1]/var[2][0]; - } + return new Hexahedron_hess_2_0(*this); + } + }; + + class Hexahedron_hess_2_1 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][2] * var[3][1] / var[2][2] / var[2][1]; + } - public: + public: + Hexahedron_hess_2_1() + { + this->_name = "Hexahedron_hess_2_1"; + } - Hexahedron_hess_1_0() - { - this->_name = "Hexahedron_hess_1_0"; - } + Hexahedron_hess_2_1 * + clone() const + { + return new Hexahedron_hess_2_1(*this); + } + }; + + class Hexahedron_hess_2_2 + : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return 0.0; + } - Hexahedron_hess_1_0* clone() const - { - return new Hexahedron_hess_1_0(*this); - } - }; - - class Hexahedron_hess_1_1 : public PSimpleBase< std::vector >, double> + public: + Hexahedron_hess_2_2() { - double eval( const std::vector > &var) const - { - return 0.0; - } + this->_name = "Hexahedron_hess_2_2"; + } - public: + Hexahedron_hess_2_2 * + clone() const + { + return new Hexahedron_hess_2_2(*this); + } + }; + + class Hexahedron : public PFuncBase>, double> + { + PSimpleBase>, double> *_val; + PSimpleBase>, double> **_grad_val; + PSimpleBase>, double> ***_hess_val; + + public: + Hexahedron() + { + construct(); + } - Hexahedron_hess_1_1() - { - this->_name = "Hexahedron_hess_1_1"; - } + Hexahedron(const Hexahedron &RHS) + { + construct(); + } - Hexahedron_hess_1_1* clone() const - { - return new Hexahedron_hess_1_1(*this); - } - }; - - class Hexahedron_hess_1_2 : public PSimpleBase< std::vector >, double> + ~Hexahedron() { - double eval( const std::vector > &var) const - { - return var[3][1]*var[3][2]/var[2][1]/var[2][2]; - } + delete _val; + + delete _grad_val[0]; + delete _grad_val[1]; + delete _grad_val[2]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete _hess_val[0][1]; + delete _hess_val[0][2]; + delete _hess_val[1][0]; + delete _hess_val[1][1]; + delete _hess_val[1][2]; + delete _hess_val[2][0]; + delete _hess_val[2][1]; + delete _hess_val[2][2]; + delete[] _hess_val[0]; + delete[] _hess_val[1]; + delete[] _hess_val[2]; + delete[] _hess_val; + } + + Hexahedron * + clone() const + { + return new Hexahedron(*this); + } - public: + PSimpleFunction>, double> + simplefunction() const + { + return PSimpleFunction>, double>(*_val); + } - Hexahedron_hess_1_2() - { - this->_name = "Hexahedron_hess_1_2"; - } + PSimpleFunction>, double> + grad_simplefunction(size_type di) const + { + return PSimpleFunction>, double>(*_grad_val[di]); + } - Hexahedron_hess_1_2* clone() const - { - return new Hexahedron_hess_1_2(*this); - } - }; - - class Hexahedron_hess_2_0 : public PSimpleBase< std::vector >, double> + PSimpleFunction>, double> + hess_simplefunction(size_type di, size_type dj) const { - double eval( const std::vector > &var) const - { - return var[3][2]*var[3][0]/var[2][2]/var[2][0]; - } + return PSimpleFunction>, double>( + *_hess_val[di][dj]); + } - public: + double + operator()(const std::vector> &var) + { + return (*_val)(var); + } - Hexahedron_hess_2_0() - { - this->_name = "Hexahedron_hess_2_0"; - } + double + grad(const std::vector> &var, size_type di) + { + return (*_grad_val[di])(var); + } - Hexahedron_hess_2_0* clone() const - { - return new Hexahedron_hess_2_0(*this); - } - }; - - class Hexahedron_hess_2_1 : public PSimpleBase< std::vector >, double> + double + hess(const std::vector> &var, size_type di, size_type dj) { - double eval( const std::vector > &var) const - { - return var[3][2]*var[3][1]/var[2][2]/var[2][1]; - } + return (*_hess_val[di][dj])(var); + } - public: + void + eval(const std::vector> &var) + { + (*_val)(var); + } - Hexahedron_hess_2_1() - { - this->_name = "Hexahedron_hess_2_1"; - } + void + eval_grad(const std::vector> &var) + { + (*_grad_val[0])(var); + (*_grad_val[1])(var); + } - Hexahedron_hess_2_1* clone() const - { - return new Hexahedron_hess_2_1(*this); - } - }; - - class Hexahedron_hess_2_2 : public PSimpleBase< std::vector >, double> + void + eval_hess(const std::vector> &var) { - double eval( const std::vector > &var) const - { - return 0.0; - } + (*_hess_val[0][0])(var); + (*_hess_val[0][1])(var); + (*_hess_val[1][0])(var); + (*_hess_val[1][1])(var); + } + + double + operator()() const + { + return (*_val)(); + } - public: + double + grad(size_type di) const + { + return (*_grad_val[di])(); + } - Hexahedron_hess_2_2() - { - this->_name = "Hexahedron_hess_2_2"; - } + double + hess(size_type di, size_type dj) const + { + return (*_hess_val[di][dj])(); + } + + private: + void + construct() + { + this->_name = "Hexahedron"; + this->_var_name.clear(); + this->_var_name.push_back("r"); + this->_var_name.push_back("n"); + this->_var_name.push_back("h"); + this->_var_name.push_back("s"); + this->_var_description.clear(); + this->_var_description.push_back("Coordinate to be evaluated (Cartesian)"); + this->_var_description.push_back("Coordinate of node"); + this->_var_description.push_back("Coordinate containing element dimensions"); + this->_var_description.push_back( + "Coordinate containing +/- 1.0, depending on which corner of quad element"); + + _val = new Hexahedron_f(); + + _grad_val = new PSimpleBase>, double> *[3]; + _grad_val[0] = new Hexahedron_grad_0(); + _grad_val[1] = new Hexahedron_grad_1(); + _grad_val[2] = new Hexahedron_grad_2(); + + _hess_val = new PSimpleBase>, double> **[3]; + _hess_val[0] = new PSimpleBase>, double> *[3]; + _hess_val[1] = new PSimpleBase>, double> *[3]; + _hess_val[2] = new PSimpleBase>, double> *[3]; + _hess_val[0][0] = new Hexahedron_hess_0_0(); + _hess_val[0][1] = new Hexahedron_hess_0_1(); + _hess_val[0][2] = new Hexahedron_hess_0_2(); + _hess_val[1][0] = new Hexahedron_hess_1_0(); + _hess_val[1][1] = new Hexahedron_hess_1_1(); + _hess_val[1][2] = new Hexahedron_hess_1_2(); + _hess_val[2][0] = new Hexahedron_hess_2_0(); + _hess_val[2][1] = new Hexahedron_hess_2_1(); + _hess_val[2][2] = new Hexahedron_hess_2_2(); + } + }; + + /// A base class for interpolating functions + /// + template + class HexahedronValues : public Interpolator + { + //_var[0]: Coordinate _r; // coordinate to evaluate field at + //_var[1]: Coordinate _n; // coordinate of node + //_var[2]: Coordinate _h; // quad dimensions + //_var[3]: Coordinate _s; // +/- 1, depending on orientation of basis function + std::vector> _var; + + public: + typedef typename Interpolator::size_type size_type; + + // node_index: index of node in mesh + // node_index: index of element in mesh + // bfunc: PFuncBase, double>* + // node_coord: Coordinate of node + // dim: Coordinate containing x and y dimension of element + // element_node_index: 0 == bottom left, proceed counter clockwise to 3 == top left of + // element + + HexahedronValues(unsigned long int node_index, + unsigned long int element_index, + PFuncBase>, double> *bfunc, + const PRISMS::Coordinate<3> &node_coord, + const PRISMS::Coordinate<3> &dim, + int element_node_index) + : Interpolator(node_index, element_index, bfunc) + { + _var.resize(4); + + _var[1][0] = node_coord[0]; + _var[1][1] = node_coord[1]; + _var[1][2] = node_coord[2]; + + _var[2][0] = dim[0]; + _var[2][1] = dim[1]; + _var[2][2] = dim[2]; - Hexahedron_hess_2_2* clone() const + if (element_node_index == 0) { - return new Hexahedron_hess_2_2(*this); + _var[3][0] = 1.0; + _var[3][1] = 1.0; + _var[3][2] = 1.0; } - }; - - - class Hexahedron : public PFuncBase >, double> - { - PSimpleBase >, double>* _val; - PSimpleBase >, double>** _grad_val; - PSimpleBase >, double>*** _hess_val; - - public: - Hexahedron() + else if (element_node_index == 1) { - construct(); + _var[3][0] = -1.0; + _var[3][1] = 1.0; + _var[3][2] = 1.0; } - - Hexahedron( const Hexahedron &RHS) + else if (element_node_index == 2) { - construct(); + _var[3][0] = -1.0; + _var[3][1] = -1.0; + _var[3][2] = 1.0; } - - ~Hexahedron() + else if (element_node_index == 3) { - delete _val; - - delete _grad_val[0]; - delete _grad_val[1]; - delete _grad_val[2]; - delete [] _grad_val; - - delete _hess_val[0][0]; - delete _hess_val[0][1]; - delete _hess_val[0][2]; - delete _hess_val[1][0]; - delete _hess_val[1][1]; - delete _hess_val[1][2]; - delete _hess_val[2][0]; - delete _hess_val[2][1]; - delete _hess_val[2][2]; - delete [] _hess_val[0]; - delete [] _hess_val[1]; - delete [] _hess_val[2]; - delete [] _hess_val; + _var[3][0] = 1.0; + _var[3][1] = -1.0; + _var[3][2] = 1.0; } - - Hexahedron* clone() const + else if (element_node_index == 4) { - return new Hexahedron(*this); + _var[3][0] = 1.0; + _var[3][1] = 1.0; + _var[3][2] = -1.0; } - - PSimpleFunction< std::vector >, double> simplefunction() const + else if (element_node_index == 5) { - return PSimpleFunction< std::vector >, double>( *_val ); + _var[3][0] = -1.0; + _var[3][1] = 1.0; + _var[3][2] = -1.0; } - - PSimpleFunction< std::vector >, double> grad_simplefunction(size_type di) const + else if (element_node_index == 6) { - return PSimpleFunction< std::vector >, double>( *_grad_val[di] ); + _var[3][0] = -1.0; + _var[3][1] = -1.0; + _var[3][2] = -1.0; } - - PSimpleFunction< std::vector >, double> hess_simplefunction(size_type di, size_type dj) const + else if (element_node_index == 7) { - return PSimpleFunction< std::vector >, double>( *_hess_val[di][dj] ); + _var[3][0] = 1.0; + _var[3][1] = -1.0; + _var[3][2] = -1.0; } + } - double operator()(const std::vector > &var) - { - return (*_val)(var); - } + PRISMS::Coordinate<3> + min() const + { + PRISMS::Coordinate<3> coord = _var[1]; - double grad(const std::vector > &var, size_type di) - { - return (*_grad_val[di])(var); - } + if (_var[3][0] == -1.0) + coord[0] -= _var[2][0]; + if (_var[3][1] == -1.0) + coord[1] -= _var[2][1]; + if (_var[3][2] == -1.0) + coord[2] -= _var[2][2]; - double hess(const std::vector > &var, size_type di, size_type dj) - { - return (*_hess_val[di][dj])(var); - } + return coord; + } - void eval(const std::vector > &var) - { - (*_val)(var); - } + PRISMS::Coordinate<3> + max() const + { + PRISMS::Coordinate<3> coord = _var[1]; - void eval_grad(const std::vector > &var) - { - (*_grad_val[0])(var); - (*_grad_val[1])(var); - } + if (_var[3][0] == 1.0) + coord[0] += _var[2][0]; + if (_var[3][1] == 1.0) + coord[1] += _var[2][1]; + if (_var[3][2] == 1.0) + coord[2] += _var[2][2]; - void eval_hess(const std::vector > &var) - { - (*_hess_val[0][0])(var); - (*_hess_val[0][1])(var); - (*_hess_val[1][0])(var); - (*_hess_val[1][1])(var); - } + return coord; + } - double operator()() const - { - return (*_val)(); - } + bool + is_in_range(const Coordinate &coord) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + _var[0][2] = coord[2]; + double e; - double grad(size_type di) const + for (int i = 0; i < 3; i++) { - return (*_grad_val[di])(); - } + e = _var[3][i] * (_var[0][i] - _var[1][i]) / _var[2][i]; + if (e < 0.0 || e >= 1.0) + return false; - double hess(size_type di, size_type dj) const - { - return (*_hess_val[di][dj])(); - } - - private: - void construct() - { - this->_name = "Hexahedron"; - this->_var_name.clear(); - this->_var_name.push_back("r"); - this->_var_name.push_back("n"); - this->_var_name.push_back("h"); - this->_var_name.push_back("s"); - this->_var_description.clear(); - this->_var_description.push_back("Coordinate to be evaluated (Cartesian)"); - this->_var_description.push_back("Coordinate of node"); - this->_var_description.push_back("Coordinate containing element dimensions"); - this->_var_description.push_back("Coordinate containing +/- 1.0, depending on which corner of quad element"); - - _val = new Hexahedron_f(); - - _grad_val = new PSimpleBase< std::vector >, double>*[3]; - _grad_val[0] = new Hexahedron_grad_0(); - _grad_val[1] = new Hexahedron_grad_1(); - _grad_val[2] = new Hexahedron_grad_2(); - - _hess_val = new PSimpleBase< std::vector >, double>**[3]; - _hess_val[0] = new PSimpleBase< std::vector >, double>*[3]; - _hess_val[1] = new PSimpleBase< std::vector >, double>*[3]; - _hess_val[2] = new PSimpleBase< std::vector >, double>*[3]; - _hess_val[0][0] = new Hexahedron_hess_0_0(); - _hess_val[0][1] = new Hexahedron_hess_0_1(); - _hess_val[0][2] = new Hexahedron_hess_0_2(); - _hess_val[1][0] = new Hexahedron_hess_1_0(); - _hess_val[1][1] = new Hexahedron_hess_1_1(); - _hess_val[1][2] = new Hexahedron_hess_1_2(); - _hess_val[2][0] = new Hexahedron_hess_2_0(); - _hess_val[2][1] = new Hexahedron_hess_2_1(); - _hess_val[2][2] = new Hexahedron_hess_2_2(); - } - }; - - - /// A base class for interpolating functions - /// - template - class HexahedronValues : public Interpolator - { - //_var[0]: Coordinate _r; // coordinate to evaluate field at - //_var[1]: Coordinate _n; // coordinate of node - //_var[2]: Coordinate _h; // quad dimensions - //_var[3]: Coordinate _s; // +/- 1, depending on orientation of basis function - std::vector< PRISMS::Coordinate<3> > _var; - - public: - - typedef typename Interpolator::size_type size_type; - - // node_index: index of node in mesh - // node_index: index of element in mesh - // bfunc: PFuncBase, double>* - // node_coord: Coordinate of node - // dim: Coordinate containing x and y dimension of element - // element_node_index: 0 == bottom left, proceed counter clockwise to 3 == top left of element - - HexahedronValues( unsigned long int node_index, - unsigned long int element_index, - PFuncBase >, double>* bfunc, - const PRISMS::Coordinate<3> &node_coord, - const PRISMS::Coordinate<3> &dim, - int element_node_index) - : Interpolator(node_index, element_index, bfunc) - { - _var.resize(4); - - _var[1][0] = node_coord[0]; - _var[1][1] = node_coord[1]; - _var[1][2] = node_coord[2]; - - _var[2][0] = dim[0]; - _var[2][1] = dim[1]; - _var[2][2] = dim[2]; - - if( element_node_index == 0) - { - _var[3][0] = 1.0; - _var[3][1] = 1.0; - _var[3][2] = 1.0; - } - else if( element_node_index == 1) - { - _var[3][0] = -1.0; - _var[3][1] = 1.0; - _var[3][2] = 1.0; - } - else if( element_node_index == 2) - { - _var[3][0] = -1.0; - _var[3][1] = -1.0; - _var[3][2] = 1.0; - } - else if( element_node_index == 3) - { - _var[3][0] = 1.0; - _var[3][1] = -1.0; - _var[3][2] = 1.0; - } - else if( element_node_index == 4) - { - _var[3][0] = 1.0; - _var[3][1] = 1.0; - _var[3][2] = -1.0; - } - else if( element_node_index == 5) - { - _var[3][0] = -1.0; - _var[3][1] = 1.0; - _var[3][2] = -1.0; - } - else if( element_node_index == 6) - { - _var[3][0] = -1.0; - _var[3][1] = -1.0; - _var[3][2] = -1.0; - } - else if( element_node_index == 7) - { - _var[3][0] = 1.0; - _var[3][1] = -1.0; - _var[3][2] = -1.0; - } - } - - PRISMS::Coordinate<3> min() const - { - PRISMS::Coordinate<3> coord = _var[1]; - - if( _var[3][0] == -1.0) - coord[0] -= _var[2][0]; - if( _var[3][1] == -1.0) - coord[1] -= _var[2][1]; - if( _var[3][2] == -1.0) - coord[2] -= _var[2][2]; - - return coord; - } - - PRISMS::Coordinate<3> max() const - { - PRISMS::Coordinate<3> coord = _var[1]; - - if( _var[3][0] == 1.0) - coord[0] += _var[2][0]; - if( _var[3][1] == 1.0) - coord[1] += _var[2][1]; - if( _var[3][2] == 1.0) - coord[2] += _var[2][2]; - - return coord; - } - - bool is_in_range(const Coordinate &coord) - { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - _var[0][2] = coord[2]; - double e; - - for( int i=0; i<3; i++) - { - e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; - if( e < 0.0 || e >= 1.0) - return false; - - //if( e == 0.0 && std::signbit(e)) - // return false; - } - - - - //std::cout << "e: " ; - //for( int i=0; i<2; i++) - //{ - // e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; - // std::cout << e << " "; - //} - //std::cout << std::endl; - - return true; - } - - // for the following, - // you are expected to KNOW that the coord is_in_range!!! - - double operator()(const Coordinate &coord) - { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - _var[0][2] = coord[2]; - return (*this->_bfunc)(_var); - } - - double grad(const Coordinate &coord, size_type di) - { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - _var[0][2] = coord[2]; - return (*this->_bfunc).grad(_var,di); - } - - double hess(const Coordinate &coord, size_type di, size_type dj) - { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - _var[0][2] = coord[2]; - return (*this->_bfunc).hess(_var,di,dj); + // if( e == 0.0 && std::signbit(e)) + // return false; } - }; + // std::cout << "e: " ; + // for( int i=0; i<2; i++) + //{ + // e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; + // std::cout << e << " "; + // } + // std::cout << std::endl; + + return true; + } -} + // for the following, + // you are expected to KNOW that the coord is_in_range!!! + double + operator()(const Coordinate &coord) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + _var[0][2] = coord[2]; + return (*this->_bfunc)(_var); + } + + double + grad(const Coordinate &coord, size_type di) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + _var[0][2] = coord[2]; + return (*this->_bfunc).grad(_var, di); + } + + double + hess(const Coordinate &coord, size_type di, size_type dj) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + _var[0][2] = coord[2]; + return (*this->_bfunc).hess(_var, di, dj); + } + }; + +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/pfield/interpolation/Interpolator.hh b/include/IntegrationTools/pfield/interpolation/Interpolator.hh index 325eebb87..8094d6d18 100644 --- a/include/IntegrationTools/pfield/interpolation/Interpolator.hh +++ b/include/IntegrationTools/pfield/interpolation/Interpolator.hh @@ -7,86 +7,96 @@ namespace PRISMS { - - /// A base class for interpolating functions - /// - template - class Interpolator + + /// A base class for interpolating functions + /// + template + class Interpolator + { + protected: + unsigned long int _node; // index of nodal value or control point + unsigned long int _element; // index of element + PFuncBase>, double> *_bfunc; // basis function to + // evaluate + + public: + typedef typename PFuncBase>, double>::size_type + size_type; + + Interpolator(unsigned long int node, + unsigned long int element, + PFuncBase>, double> *bfunc) + : _node(node) + , _element(element) + , _bfunc(bfunc) {}; + + virtual ~Interpolator() {}; + + unsigned long int + node() + { + return _node; + } + + unsigned long int + element() + { + return _element; + } + + virtual PRISMS::Coordinate + min() const + { + undefined("void min(Coordinate &coord) const"); + return PRISMS::Coordinate(); + } + + virtual PRISMS::Coordinate + max() const + { + undefined("void max(Coordinate &coord) const"); + return PRISMS::Coordinate(); + } + + virtual bool + is_in_range(const Coordinate &coord) + { + undefined("bool is_in_range(Coordinate coord) const"); + return false; + } + + virtual double + operator()(const Coordinate &coord) + { + undefined("double operator()(Coordinate coord)"); + return double(); + } + + virtual double + grad(const Coordinate &coord, size_type di) + { + undefined("double grad()(Coordinate coord, size_type di)"); + return double(); + } + + virtual double + hess(const Coordinate &coord, size_type di, size_type dj) + { + undefined("double hess()(Coordinate coord, size_type di, size_type dj)"); + return double(); + } + + private: + void + undefined(std::string fname) const { - protected: - - unsigned long int _node; //index of nodal value or control point - unsigned long int _element; //index of element - PFuncBase >, double>* _bfunc; // basis function to evaluate - - public: - - typedef typename PFuncBase >, double>::size_type size_type; - - Interpolator( unsigned long int node, unsigned long int element, PFuncBase >, double>* bfunc): - _node(node), _element(element), _bfunc(bfunc) - {}; - - virtual ~Interpolator() {}; - - unsigned long int node() - { - return _node; - } - - unsigned long int element() - { - return _element; - } - - virtual PRISMS::Coordinate min() const - { - undefined("void min(Coordinate &coord) const"); - return PRISMS::Coordinate(); - } - - virtual PRISMS::Coordinate max() const - { - undefined("void max(Coordinate &coord) const"); - return PRISMS::Coordinate(); - } - - virtual bool is_in_range(const Coordinate &coord) - { - undefined("bool is_in_range(Coordinate coord) const"); - return false; - } - - virtual double operator()(const Coordinate &coord) - { - undefined("double operator()(Coordinate coord)"); - return double(); - } - - virtual double grad(const Coordinate &coord, size_type di) - { - undefined("double grad()(Coordinate coord, size_type di)"); - return double(); - } - - virtual double hess(const Coordinate &coord, size_type di, size_type dj) - { - undefined("double hess()(Coordinate coord, size_type di, size_type dj)"); - return double(); - } - - private: - - void undefined(std::string fname) const - { - std::cout << "Error in Interpolator." << std::endl; - std::cout << " The member function '" << fname << "' has not been defined." << std::endl; - exit(1); - } - - }; - -} + std::cout << "Error in Interpolator." << std::endl; + std::cout << " The member function '" << fname << "' has not been defined." + << std::endl; + exit(1); + } + }; +} // namespace PRISMS #endif diff --git a/include/IntegrationTools/pfield/interpolation/Quad.hh b/include/IntegrationTools/pfield/interpolation/Quad.hh index 6dc808165..31eb786a6 100644 --- a/include/IntegrationTools/pfield/interpolation/Quad.hh +++ b/include/IntegrationTools/pfield/interpolation/Quad.hh @@ -2,454 +2,478 @@ #ifndef Quad_HH #define Quad_HH +#include "../../pfunction/PFuncBase.hh" #include "../../pfunction/PSimpleBase.hh" #include "../../pfunction/PSimpleFunction.hh" -#include "../../pfunction/PFuncBase.hh" -#include "./Interpolator.hh" #include "../Coordinate.hh" +#include "./Interpolator.hh" namespace PRISMS { - class Quad_f : public PSimpleBase< std::vector >, double> + class Quad_f : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const { - double eval( const std::vector > &var) const - { - // var[0]: coordinate to be evaluated - // var[1]: nodal coordinate - // var[2]: element dimension - // var[3]: +/- 1 depending on which 'corner' of quad - // - // f = (1.0 - e0)*(1.0 - e1) - // e = var[3]*(var[0] - var[1])/var[2] - - return (1.0 - var[3][0]*(var[0][0] - var[1][0])/var[2][0])*(1.0 - var[3][1]*(var[0][1] - var[1][1])/var[2][1]); - } + // var[0]: coordinate to be evaluated + // var[1]: nodal coordinate + // var[2]: element dimension + // var[3]: +/- 1 depending on which 'corner' of quad + // + // f = (1.0 - e0)*(1.0 - e1) + // e = var[3]*(var[0] - var[1])/var[2] + + return (1.0 - var[3][0] * (var[0][0] - var[1][0]) / var[2][0]) * + (1.0 - var[3][1] * (var[0][1] - var[1][1]) / var[2][1]); + } + + public: + Quad_f() + { + this->_name = "Quad_f"; + } - public: + Quad_f * + clone() const + { + return new Quad_f(*this); + } + }; + + class Quad_grad_0 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return -var[3][0] * (1.0 - var[3][1] * (var[0][1] - var[1][1]) / var[2][1]) / + var[2][0]; + } - Quad_f() - { - this->_name = "Quad_f"; - } + public: + Quad_grad_0() + { + this->_name = "Quad_grad_0"; + } - Quad_f* clone() const - { - return new Quad_f(*this); - } - }; - - class Quad_grad_0 : public PSimpleBase< std::vector >, double> + Quad_grad_0 * + clone() const { - double eval( const std::vector > &var) const - { - return -var[3][0]*(1.0 - var[3][1]*(var[0][1] - var[1][1])/var[2][1])/var[2][0]; - } + return new Quad_grad_0(*this); + } + }; + + class Quad_grad_1 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return -var[3][1] * (1.0 - var[3][0] * (var[0][0] - var[1][0]) / var[2][0]) / + var[2][1]; + } - public: + public: + Quad_grad_1() + { + this->_name = "Quad_grad_1"; + } - Quad_grad_0() - { - this->_name = "Quad_grad_0"; - } + Quad_grad_1 * + clone() const + { + return new Quad_grad_1(*this); + } + }; + + class Quad_hess_0_0 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return 0.0; + } - Quad_grad_0* clone() const - { - return new Quad_grad_0(*this); - } - }; - - class Quad_grad_1 : public PSimpleBase< std::vector >, double> + public: + Quad_hess_0_0() { - double eval( const std::vector > &var) const - { - return -var[3][1]*(1.0 - var[3][0]*(var[0][0] - var[1][0])/var[2][0])/var[2][1]; - - } + this->_name = "Quad_hess_0_0"; + } - public: + Quad_hess_0_0 * + clone() const + { + return new Quad_hess_0_0(*this); + } + }; + + class Quad_hess_0_1 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][0] * var[3][1] / var[2][0] / var[2][1]; + } - Quad_grad_1() - { - this->_name = "Quad_grad_1"; - } + public: + Quad_hess_0_1() + { + this->_name = "Quad_hess_0_1"; + } - Quad_grad_1* clone() const - { - return new Quad_grad_1(*this); - } - }; - - class Quad_hess_0_0 : public PSimpleBase< std::vector >, double> + Quad_hess_0_1 * + clone() const { - double eval( const std::vector > &var) const - { - return 0.0; - } + return new Quad_hess_0_1(*this); + } + }; + + class Quad_hess_1_0 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return var[3][0] * var[3][1] / var[2][0] / var[2][1]; + } - public: + public: + Quad_hess_1_0() + { + this->_name = "Quad_hess_1_0"; + } - Quad_hess_0_0() - { - this->_name = "Quad_hess_0_0"; - } + Quad_hess_1_0 * + clone() const + { + return new Quad_hess_1_0(*this); + } + }; + + class Quad_hess_1_1 : public PSimpleBase>, double> + { + double + eval(const std::vector> &var) const + { + return 0.0; + } - Quad_hess_0_0* clone() const - { - return new Quad_hess_0_0(*this); - } - }; - - class Quad_hess_0_1 : public PSimpleBase< std::vector >, double> + public: + Quad_hess_1_1() { - double eval( const std::vector > &var) const - { - return var[3][0]*var[3][1]/var[2][0]/var[2][1]; - } + this->_name = "Quad_hess_1_1"; + } - public: + Quad_hess_1_1 * + clone() const + { + return new Quad_hess_1_1(*this); + } + }; - Quad_hess_0_1() - { - this->_name = "Quad_hess_0_1"; - } + class Quad : public PFuncBase>, double> + { + PSimpleBase>, double> *_val; + PSimpleBase>, double> **_grad_val; + PSimpleBase>, double> ***_hess_val; - Quad_hess_0_1* clone() const - { - return new Quad_hess_0_1(*this); - } - }; - - class Quad_hess_1_0 : public PSimpleBase< std::vector >, double> + public: + typedef PFuncBase>, double>::size_type size_type; + + Quad() { - double eval( const std::vector > &var) const - { - return var[3][0]*var[3][1]/var[2][0]/var[2][1]; - } + construct(); + } - public: + Quad(const Quad &RHS) + { + construct(); + } - Quad_hess_1_0() - { - this->_name = "Quad_hess_1_0"; - } + ~Quad() + { + delete _val; + + delete _grad_val[0]; + delete _grad_val[1]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete _hess_val[0][1]; + delete _hess_val[1][0]; + delete _hess_val[1][1]; + delete[] _hess_val[0]; + delete[] _hess_val[1]; + delete[] _hess_val; + } + + Quad * + clone() const + { + return new Quad(*this); + } - Quad_hess_1_0* clone() const - { - return new Quad_hess_1_0(*this); - } - }; - - class Quad_hess_1_1 : public PSimpleBase< std::vector >, double> + PSimpleFunction>, double> + simplefunction() const { - double eval( const std::vector > &var) const - { - return 0.0; - } + return PSimpleFunction>, double>(*_val); + } - public: + PSimpleFunction>, double> + grad_simplefunction(size_type di) const + { + return PSimpleFunction>, double>(*_grad_val[di]); + } - Quad_hess_1_1() - { - this->_name = "Quad_hess_1_1"; - } + PSimpleFunction>, double> + hess_simplefunction(size_type di, size_type dj) const + { + return PSimpleFunction>, double>( + *_hess_val[di][dj]); + } - Quad_hess_1_1* clone() const - { - return new Quad_hess_1_1(*this); - } - }; - - - class Quad : public PFuncBase >, double> - { - PSimpleBase >, double>* _val; - PSimpleBase >, double>** _grad_val; - PSimpleBase >, double>*** _hess_val; - - public: - - typedef PFuncBase >, double>::size_type size_type; - - Quad() - { - construct(); - } - - Quad( const Quad &RHS) - { - construct(); - } - - ~Quad() - { - delete _val; - - delete _grad_val[0]; - delete _grad_val[1]; - delete [] _grad_val; - - delete _hess_val[0][0]; - delete _hess_val[0][1]; - delete _hess_val[1][0]; - delete _hess_val[1][1]; - delete [] _hess_val[0]; - delete [] _hess_val[1]; - delete [] _hess_val; - } - - Quad* clone() const - { - return new Quad(*this); - } + double + operator()(const std::vector> &var) + { + return (*_val)(var); + } - PSimpleFunction< std::vector >, double> simplefunction() const - { - return PSimpleFunction< std::vector >, double>( *_val ); - } + double + grad(const std::vector> &var, size_type di) + { + return (*_grad_val[di])(var); + } - PSimpleFunction< std::vector >, double> grad_simplefunction(size_type di) const - { - return PSimpleFunction< std::vector >, double>( *_grad_val[di] ); - } + double + hess(const std::vector> &var, size_type di, size_type dj) + { + return (*_hess_val[di][dj])(var); + } - PSimpleFunction< std::vector >, double> hess_simplefunction(size_type di, size_type dj) const - { - return PSimpleFunction< std::vector >, double>( *_hess_val[di][dj] ); - } + void + eval(const std::vector> &var) + { + (*_val)(var); + } - double operator()(const std::vector > &var) - { - return (*_val)(var); - } + void + eval_grad(const std::vector> &var) + { + (*_grad_val[0])(var); + (*_grad_val[1])(var); + } - double grad(const std::vector > &var, size_type di) - { - return (*_grad_val[di])(var); - } + void + eval_hess(const std::vector> &var) + { + (*_hess_val[0][0])(var); + (*_hess_val[0][1])(var); + (*_hess_val[1][0])(var); + (*_hess_val[1][1])(var); + } + + double + operator()() const + { + return (*_val)(); + } - double hess(const std::vector > &var, size_type di, size_type dj) - { - return (*_hess_val[di][dj])(var); - } + double + grad(size_type di) const + { + return (*_grad_val[di])(); + } - void eval(const std::vector > &var) - { - (*_val)(var); - } + double + hess(size_type di, size_type dj) const + { + return (*_hess_val[di][dj])(); + } - void eval_grad(const std::vector > &var) - { - (*_grad_val[0])(var); - (*_grad_val[1])(var); - } + private: + void + construct() + { + this->_name = "Quad"; + this->_var_name.clear(); + this->_var_name.push_back("r"); + this->_var_name.push_back("n"); + this->_var_name.push_back("h"); + this->_var_name.push_back("s"); + this->_var_description.clear(); + this->_var_description.push_back("Coordinate to be evaluated (Cartesian)"); + this->_var_description.push_back("Coordinate of node"); + this->_var_description.push_back("Coordinate containing element dimensions"); + this->_var_description.push_back( + "Coordinate containing +/- 1.0, depending on which corner of quad element"); + + _val = new Quad_f(); + + _grad_val = new PSimpleBase>, double> *[2]; + _grad_val[0] = new Quad_grad_0(); + _grad_val[1] = new Quad_grad_1(); + + _hess_val = new PSimpleBase>, double> **[2]; + _hess_val[0] = new PSimpleBase>, double> *[2]; + _hess_val[1] = new PSimpleBase>, double> *[2]; + _hess_val[0][0] = new Quad_hess_0_0(); + _hess_val[0][1] = new Quad_hess_0_1(); + _hess_val[1][0] = new Quad_hess_1_0(); + _hess_val[1][1] = new Quad_hess_1_1(); + } + }; + + /// A base class for interpolating functions + /// + template + class QuadValues : public Interpolator + { + //_var[0]: Coordinate _r; // coordinate to evaluate field at + //_var[1]: Coordinate _n; // coordinate of node + //_var[2]: Coordinate _h; // quad dimensions + //_var[3]: Coordinate _s; // +/- 1, depending on orientation of basis function + std::vector> _var; + + public: + typedef typename Interpolator::size_type size_type; + + // node_index: index of node in mesh + // node_index: index of element in mesh + // bfunc: PFuncBase, double>* + // node_coord: Coordinate of node + // dim: Coordinate containing x and y dimension of element + // element_node_index: 0 == bottom left, proceed counter clockwise to 3 == top left of + // element + + QuadValues(unsigned long int node_index, + unsigned long int element_index, + PFuncBase>, double> *bfunc, + const PRISMS::Coordinate<2> &node_coord, + const PRISMS::Coordinate<2> &dim, + int element_node_index) + : Interpolator(node_index, element_index, bfunc) + { + _var.resize(4); - void eval_hess(const std::vector > &var) - { - (*_hess_val[0][0])(var); - (*_hess_val[0][1])(var); - (*_hess_val[1][0])(var); - (*_hess_val[1][1])(var); - } + _var[1][0] = node_coord[0]; + _var[1][1] = node_coord[1]; - double operator()() const - { - return (*_val)(); - } + _var[2][0] = dim[0]; + _var[2][1] = dim[1]; - double grad(size_type di) const + if (element_node_index == 0) { - return (*_grad_val[di])(); - } + //_var[1][0] = node_coord[0]; + //_var[1][1] = node_coord[1]; - double hess(size_type di, size_type dj) const - { - return (*_hess_val[di][dj])(); + _var[3][0] = 1.0; + _var[3][1] = 1.0; } - - private: - void construct() + else if (element_node_index == 1) { - this->_name = "Quad"; - this->_var_name.clear(); - this->_var_name.push_back("r"); - this->_var_name.push_back("n"); - this->_var_name.push_back("h"); - this->_var_name.push_back("s"); - this->_var_description.clear(); - this->_var_description.push_back("Coordinate to be evaluated (Cartesian)"); - this->_var_description.push_back("Coordinate of node"); - this->_var_description.push_back("Coordinate containing element dimensions"); - this->_var_description.push_back("Coordinate containing +/- 1.0, depending on which corner of quad element"); - - _val = new Quad_f(); - - _grad_val = new PSimpleBase< std::vector >, double>*[2]; - _grad_val[0] = new Quad_grad_0(); - _grad_val[1] = new Quad_grad_1(); - - _hess_val = new PSimpleBase< std::vector >, double>**[2]; - _hess_val[0] = new PSimpleBase< std::vector >, double>*[2]; - _hess_val[1] = new PSimpleBase< std::vector >, double>*[2]; - _hess_val[0][0] = new Quad_hess_0_0(); - _hess_val[0][1] = new Quad_hess_0_1(); - _hess_val[1][0] = new Quad_hess_1_0(); - _hess_val[1][1] = new Quad_hess_1_1(); - } - }; - - - /// A base class for interpolating functions - /// - template - class QuadValues : public Interpolator - { - //_var[0]: Coordinate _r; // coordinate to evaluate field at - //_var[1]: Coordinate _n; // coordinate of node - //_var[2]: Coordinate _h; // quad dimensions - //_var[3]: Coordinate _s; // +/- 1, depending on orientation of basis function - std::vector< PRISMS::Coordinate<2> > _var; - - public: - - typedef typename Interpolator::size_type size_type; - - // node_index: index of node in mesh - // node_index: index of element in mesh - // bfunc: PFuncBase, double>* - // node_coord: Coordinate of node - // dim: Coordinate containing x and y dimension of element - // element_node_index: 0 == bottom left, proceed counter clockwise to 3 == top left of element - - QuadValues( unsigned long int node_index, - unsigned long int element_index, - PFuncBase >, double>* bfunc, - const PRISMS::Coordinate<2> &node_coord, - const PRISMS::Coordinate<2> &dim, - int element_node_index) - : Interpolator(node_index, element_index, bfunc) - { - _var.resize(4); - - _var[1][0] = node_coord[0]; - _var[1][1] = node_coord[1]; - - _var[2][0] = dim[0]; - _var[2][1] = dim[1]; - - if( element_node_index == 0) - { - //_var[1][0] = node_coord[0]; - //_var[1][1] = node_coord[1]; - - _var[3][0] = 1.0; - _var[3][1] = 1.0; - } - else if( element_node_index == 1) - { - //_var[1][0] = node_coord[0] + _var[2][0]; - //_var[1][1] = node_coord[1]; - - _var[3][0] = -1.0; - _var[3][1] = 1.0; - } - else if( element_node_index == 2) - { - //_var[1][0] = node_coord[0] + _var[2][0]; - //_var[1][1] = node_coord[1] + _var[2][1]; - - _var[3][0] = -1.0; - _var[3][1] = -1.0; - } - else if( element_node_index == 3) - { - //_var[1][0] = node_coord[0]; - //_var[1][1] = node_coord[1] + _var[2][1]; - - _var[3][0] = 1.0; - _var[3][1] = -1.0; - } - } - - PRISMS::Coordinate<2> min() const - { - PRISMS::Coordinate<2> coord = _var[1]; - - if( _var[3][0] == -1.0) - coord[0] -= _var[2][0]; - if( _var[3][1] == -1.0) - coord[1] -= _var[2][1]; - - return coord; - } - - PRISMS::Coordinate<2> max() const - { - PRISMS::Coordinate<2> coord = _var[1]; - - if( _var[3][0] == 1.0) - coord[0] += _var[2][0]; - if( _var[3][1] == 1.0) - coord[1] += _var[2][1]; - - return coord; - } - - bool is_in_range(const Coordinate &coord) - { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - double e; - - for( int i=0; i<2; i++) - { - e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; - if( e < 0.0 || e >= 1.0) - return false; - - //if( e == 0.0 && std::signbit(e)) - // return false; - } - - - - //std::cout << "e: " ; - //for( int i=0; i<2; i++) - //{ - // e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; - // std::cout << e << " "; - //} - //std::cout << std::endl; - - return true; + //_var[1][0] = node_coord[0] + _var[2][0]; + //_var[1][1] = node_coord[1]; + + _var[3][0] = -1.0; + _var[3][1] = 1.0; } - - // for the following, - // you are expected to KNOW that the coord is_in_range!!! - - double operator()(const Coordinate &coord) + else if (element_node_index == 2) { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - return (*this->_bfunc)(_var); + //_var[1][0] = node_coord[0] + _var[2][0]; + //_var[1][1] = node_coord[1] + _var[2][1]; + + _var[3][0] = -1.0; + _var[3][1] = -1.0; } - - double grad(const Coordinate &coord, size_type di) + else if (element_node_index == 3) { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - return (*this->_bfunc).grad(_var,di); + //_var[1][0] = node_coord[0]; + //_var[1][1] = node_coord[1] + _var[2][1]; + + _var[3][0] = 1.0; + _var[3][1] = -1.0; } - - double hess(const Coordinate &coord, size_type di, size_type dj) + } + + PRISMS::Coordinate<2> + min() const + { + PRISMS::Coordinate<2> coord = _var[1]; + + if (_var[3][0] == -1.0) + coord[0] -= _var[2][0]; + if (_var[3][1] == -1.0) + coord[1] -= _var[2][1]; + + return coord; + } + + PRISMS::Coordinate<2> + max() const + { + PRISMS::Coordinate<2> coord = _var[1]; + + if (_var[3][0] == 1.0) + coord[0] += _var[2][0]; + if (_var[3][1] == 1.0) + coord[1] += _var[2][1]; + + return coord; + } + + bool + is_in_range(const Coordinate &coord) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + double e; + + for (int i = 0; i < 2; i++) { - _var[0][0] = coord[0]; - _var[0][1] = coord[1]; - return (*this->_bfunc).hess(_var,di,dj); + e = _var[3][i] * (_var[0][i] - _var[1][i]) / _var[2][i]; + if (e < 0.0 || e >= 1.0) + return false; + + // if( e == 0.0 && std::signbit(e)) + // return false; } - }; + // std::cout << "e: " ; + // for( int i=0; i<2; i++) + //{ + // e = _var[3][i]*(_var[0][i] - _var[1][i])/_var[2][i]; + // std::cout << e << " "; + // } + // std::cout << std::endl; + + return true; + } + + // for the following, + // you are expected to KNOW that the coord is_in_range!!! + + double + operator()(const Coordinate &coord) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + return (*this->_bfunc)(_var); + } + + double + grad(const Coordinate &coord, size_type di) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + return (*this->_bfunc).grad(_var, di); + } -} + double + hess(const Coordinate &coord, size_type di, size_type dj) + { + _var[0][0] = coord[0]; + _var[0][1] = coord[1]; + return (*this->_bfunc).hess(_var, di, dj); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PBasisSet.hh b/include/IntegrationTools/pfunction/PBasisSet.hh index 54c13495c..e48a849b9 100644 --- a/include/IntegrationTools/pfunction/PBasisSet.hh +++ b/include/IntegrationTools/pfunction/PBasisSet.hh @@ -2,191 +2,213 @@ #ifndef PBasisSet_HH #define PBasisSet_HH -#include -#include -#include -#include - -#include "./PFunction.hh" #include "./PBasisSetBase.hh" +#include "./PFunction.hh" +#include +#include +#include +#include namespace PRISMS { - /// Evaluate basis functions and their derivatives. Store and access the results. - /// - template< class InType, class OutType> - class PBasisSet - { - public: - - - std::string name() const - { - return (*ptr).name(); - } - std::string description() const - { - return (*ptr)._description(); - } - int size() const - { - return (*ptr).size(); - } - - - virtual void resize( int N) - { - (*ptr).resize(N); - } - - virtual int max_size() - { - // default to (essentially) no limit - return (*ptr).max_size(); - } - - - // Output PFunction for individual basis functions: - - - PFunction basis_function(int term) - { - return (*ptr).basis_function(term); - }; - - - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - OutType operator()(int term, const InType &var) - { - return (*ptr)(term,var); - } - OutType grad(int term, const InType &var) - { - return (*ptr).grad(term, var); - } - OutType hess(int term, const InType &var) - { - return (*ptr).hess(term, var); - } - - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - - // By default, evaluate each individual term one-by-one. - // These are virtual so derived classes may implement more efficient methods. - // Returns vector containing results - virtual const std::vector& eval(const InType &var) - { - return (*ptr).eval(var); - } - virtual const std::vector& eval_grad(const InType &var) - { - return (*ptr).eval_grad(var); - } - virtual const std::vector& eval_hess(const InType &var) - { - return (*ptr).eval_hess(var); - } + /// Evaluate basis functions and their derivatives. Store and access the results. + /// + template + class PBasisSet + { + public: + std::string + name() const + { + return (*ptr).name(); + } - // Getters for individual terms - OutType operator()(int term) const - { - return (*ptr)(term); - } - OutType grad(int term) const - { - return (*ptr).grad(term); - } - OutType hess(int term) const - { - return (*ptr).hess(term); - } - - // Getters returning vector containing all terms - const std::vector& operator()() const - { - return (*ptr)(); - } - const std::vector& grad() const - { - return (*ptr).grad(); - } - const std::vector& hess() const - { - return (*ptr).hess(); - } - - // PFunction unique members ------------------------------------------ + std::string + description() const + { + return (*ptr)._description(); + } - PBasisSet &operator=(const PBasisSet &RHS) - { - if(ptr != NULL) - delete ptr; - ptr = RHS.ptr->clone(); - return *this; - } + int + size() const + { + return (*ptr).size(); + } - template PBasisSet &operator=(const T &RHS) - { - RHS.is_derived_from_PBasisSetBase(); + virtual void + resize(int N) + { + (*ptr).resize(N); + } - if(ptr != NULL) - delete ptr; - ptr = RHS.clone(); - return *this; - } + virtual int + max_size() + { + // default to (essentially) no limit + return (*ptr).max_size(); + } - // If you use this, PBasisSet becomes the 'owner' of the function RHS points to - // and it will delete it - PBasisSet &set( PBasisSet *RHS) - { - if(RHS == NULL) - { - std::cout << "Error in PBasisSet::set. RHS == NULL." << std::endl; - exit(1); - } - if(ptr != NULL) - delete ptr; - ptr = RHS; - return *this; - } + // Output PFunction for individual basis functions: - PBasisSet() - { - ptr = NULL; - } + PFunction + basis_function(int term) + { + return (*ptr).basis_function(term); + }; - PBasisSet(const PBasisSet &RHS) - { - if( RHS.ptr != NULL) - ptr = RHS.ptr->clone(); - else - ptr = NULL; - } - - template PBasisSet(const T &RHS) - { - RHS.is_derived_from_PBasisSetBase(); + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + OutType + operator()(int term, const InType &var) + { + return (*ptr)(term, var); + } - ptr = RHS.clone(); - } + OutType + grad(int term, const InType &var) + { + return (*ptr).grad(term, var); + } + + OutType + hess(int term, const InType &var) + { + return (*ptr).hess(term, var); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + + // By default, evaluate each individual term one-by-one. + // These are virtual so derived classes may implement more efficient methods. + // Returns vector containing results + virtual const std::vector & + eval(const InType &var) + { + return (*ptr).eval(var); + } + + virtual const std::vector & + eval_grad(const InType &var) + { + return (*ptr).eval_grad(var); + } + + virtual const std::vector & + eval_hess(const InType &var) + { + return (*ptr).eval_hess(var); + } + + // Getters for individual terms + OutType + operator()(int term) const + { + return (*ptr)(term); + } + + OutType + grad(int term) const + { + return (*ptr).grad(term); + } + + OutType + hess(int term) const + { + return (*ptr).hess(term); + } + + // Getters returning vector containing all terms + const std::vector & + operator()() const + { + return (*ptr)(); + } + + const std::vector & + grad() const + { + return (*ptr).grad(); + } + + const std::vector & + hess() const + { + return (*ptr).hess(); + } - ~PBasisSet() + // PFunction unique members ------------------------------------------ + + PBasisSet & + operator=(const PBasisSet &RHS) + { + if (ptr != NULL) + delete ptr; + ptr = RHS.ptr->clone(); + return *this; + } + + template + PBasisSet & + operator=(const T &RHS) + { + RHS.is_derived_from_PBasisSetBase(); + + if (ptr != NULL) + delete ptr; + ptr = RHS.clone(); + return *this; + } + + // If you use this, PBasisSet becomes the 'owner' of the function RHS points to + // and it will delete it + PBasisSet & + set(PBasisSet *RHS) + { + if (RHS == NULL) { - if(ptr != NULL) - delete ptr; + std::cout << "Error in PBasisSet::set. RHS == NULL." << std::endl; + exit(1); } + if (ptr != NULL) + delete ptr; + ptr = RHS; + return *this; + } + PBasisSet() + { + ptr = NULL; + } - private: - - PBasisSetBase *ptr; - - }; + PBasisSet(const PBasisSet &RHS) + { + if (RHS.ptr != NULL) + ptr = RHS.ptr->clone(); + else + ptr = NULL; + } + + template + PBasisSet(const T &RHS) + { + RHS.is_derived_from_PBasisSetBase(); + + ptr = RHS.clone(); + } + + ~PBasisSet() + { + if (ptr != NULL) + delete ptr; + } -} + private: + PBasisSetBase *ptr; + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PBasisSetBase.hh b/include/IntegrationTools/pfunction/PBasisSetBase.hh index c4ef5c3f6..bf96d42da 100644 --- a/include/IntegrationTools/pfunction/PBasisSetBase.hh +++ b/include/IntegrationTools/pfunction/PBasisSetBase.hh @@ -2,182 +2,213 @@ #ifndef PBasisSetBase_HH #define PBasisSetBase_HH -#include -#include -#include -#include - #include "./PFunction.hh" +#include +#include +#include +#include namespace PRISMS { - /// Evaluate basis functions and their derivatives. Store and access the results. - /// - template< class InType, class OutType> - class PBasisSetBase - { - public: - - typedef typename std::vector::size_type size_type; - - std::string _name; - std::string _description; - - std::vector _val; - std::vector _grad_val; - std::vector _hess_val; - - std::string name() const - { - return _name; - } - std::string description() const - { - return _description; - } - size_type size() const - { - return _val.size(); - } - - PBasisSetBase( size_type N) - { - resize(N); - } - - void is_derived_from_PBasisSetBase() const - { - return; - } - - virtual void resize( size_type N) - { - _val.resize(N); - _grad_val.resize(N); - _hess_val.resize(N); - } - - virtual ~PBasisSetBase() - { - - } - - virtual size_type max_size() const - { - // default to (essentially) no limit - return std::numeric_limits::max(); - } - - virtual PBasisSetBase* clone() const - { - return new PBasisSetBase(*this); - } - - virtual PFunction basis_function(size_type term) const - { - undefined("const PFunction& basis_function(size_type term) const"); - return PFunction(); - } - - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - OutType operator()(size_type term, const InType &var) - { - return _val[term] = eval(term, var); - } - OutType grad(size_type term, const InType &var) - { - return _grad_val[term] = eval_grad(term, var); - } - OutType hess(size_type term, const InType &var) - { - return _hess_val[term] = eval_hess(term, var); - } - - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - - // By default, evaluate each individual term one-by-one. - // These are virtual so derived classes may implement more efficient methods. - // Returns vector containing results - virtual const std::vector& eval(const InType &var) - { - for( size_type i=0; i<_val.size(); i++) - (*this)(i,var); - return _val; - } - virtual const std::vector& eval_grad(const InType &var) - { - for( size_type i=0; i<_val.size(); i++) - (*this).grad(i,var); - return _grad_val; - } - virtual const std::vector& eval_hess(const InType &var) - { - for( size_type i=0; i<_val.size(); i++) - (*this).hess(i,var); - return _hess_val; - } - - // Getters for individual terms - OutType operator()(size_type term) const - { - return _val[term]; - } - OutType grad(size_type term) const - { - return _grad_val[term]; - } - OutType hess(size_type term) const - { - return _hess_val[term]; - } - - // Getters returning vector containing all terms - const std::vector& operator()() const - { - return _val; - } - const std::vector& grad() const - { - return _grad_val; - } - const std::vector& hess() const - { - return _hess_val; - } - - private: - - /// ---------------------------------------------------------- - /// !!! Derived classes must define these functions !!! - /// Usually this is done using PBasisSetWriter - virtual OutType eval(size_type term, const InType &var) - { - undefined("OutType PBasisSetBase::eval(size_type term, const InType &var)"); - return OutType(); - } - virtual OutType eval_grad(size_type term, const InType &var) - { - undefined("OutType PBasisSetBase::eval_grad(size_type term, const InType &var)"); - return OutType(); - } - virtual OutType eval_hess(size_type term, const InType &var) - { - undefined("OutType PBasisSetBase::eval_hess(size_type term, const InType &var)"); - return OutType(); - } - - void undefined(std::string fname) const - { - std::cout << "Error. In PBasisSetBase '" << _name << "'." << std::endl; - std::cout << " The member function '" << fname << "' has not been defined." << std::endl; - exit(1); - } - }; - -} + /// Evaluate basis functions and their derivatives. Store and access the results. + /// + template + class PBasisSetBase + { + public: + typedef typename std::vector::size_type size_type; + + std::string _name; + std::string _description; + + std::vector _val; + std::vector _grad_val; + std::vector _hess_val; + + std::string + name() const + { + return _name; + } + + std::string + description() const + { + return _description; + } + + size_type + size() const + { + return _val.size(); + } + + PBasisSetBase(size_type N) + { + resize(N); + } + + void + is_derived_from_PBasisSetBase() const + { + return; + } + + virtual void + resize(size_type N) + { + _val.resize(N); + _grad_val.resize(N); + _hess_val.resize(N); + } + + virtual ~PBasisSetBase() + {} + + virtual size_type + max_size() const + { + // default to (essentially) no limit + return std::numeric_limits::max(); + } + + virtual PBasisSetBase * + clone() const + { + return new PBasisSetBase(*this); + } + + virtual PFunction + basis_function(size_type term) const + { + undefined("const PFunction& basis_function(size_type term) const"); + return PFunction(); + } + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + OutType + operator()(size_type term, const InType &var) + { + return _val[term] = eval(term, var); + } + + OutType + grad(size_type term, const InType &var) + { + return _grad_val[term] = eval_grad(term, var); + } + + OutType + hess(size_type term, const InType &var) + { + return _hess_val[term] = eval_hess(term, var); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + + // By default, evaluate each individual term one-by-one. + // These are virtual so derived classes may implement more efficient methods. + // Returns vector containing results + virtual const std::vector & + eval(const InType &var) + { + for (size_type i = 0; i < _val.size(); i++) + (*this)(i, var); + return _val; + } + + virtual const std::vector & + eval_grad(const InType &var) + { + for (size_type i = 0; i < _val.size(); i++) + (*this).grad(i, var); + return _grad_val; + } + + virtual const std::vector & + eval_hess(const InType &var) + { + for (size_type i = 0; i < _val.size(); i++) + (*this).hess(i, var); + return _hess_val; + } + + // Getters for individual terms + OutType + operator()(size_type term) const + { + return _val[term]; + } + + OutType + grad(size_type term) const + { + return _grad_val[term]; + } + + OutType + hess(size_type term) const + { + return _hess_val[term]; + } + + // Getters returning vector containing all terms + const std::vector & + operator()() const + { + return _val; + } + + const std::vector & + grad() const + { + return _grad_val; + } + + const std::vector & + hess() const + { + return _hess_val; + } + + private: + /// ---------------------------------------------------------- + /// !!! Derived classes must define these functions !!! + /// Usually this is done using PBasisSetWriter + virtual OutType + eval(size_type term, const InType &var) + { + undefined("OutType PBasisSetBase::eval(size_type term, const InType &var)"); + return OutType(); + } + + virtual OutType + eval_grad(size_type term, const InType &var) + { + undefined("OutType PBasisSetBase::eval_grad(size_type term, const InType &var)"); + return OutType(); + } + + virtual OutType + eval_hess(size_type term, const InType &var) + { + undefined("OutType PBasisSetBase::eval_hess(size_type term, const InType &var)"); + return OutType(); + } + + void + undefined(std::string fname) const + { + std::cout << "Error. In PBasisSetBase '" << _name << "'." << std::endl; + std::cout << " The member function '" << fname << "' has not been defined." + << std::endl; + exit(1); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PFlexFunction.hh b/include/IntegrationTools/pfunction/PFlexFunction.hh index 4da065846..f10fcf229 100644 --- a/include/IntegrationTools/pfunction/PFlexFunction.hh +++ b/include/IntegrationTools/pfunction/PFlexFunction.hh @@ -2,203 +2,219 @@ #ifndef PFlexFunction_HH #define PFlexFunction_HH -#include -#include -#include -#include - #include "./PFuncBase.hh" #include "./PSimpleFunction.hh" +#include +#include +#include +#include namespace PRISMS { - /// A class to create functions consisting of PSimpleFunctions - /// Used to create basis_functions from PSimpleFunctions for f, grad_f, and hess_f - /// - template - class PFlexFunction : public PFuncBase< VarContainer, OutType> - { - PSimpleFunction< VarContainer, OutType> _val; - std::vector< PSimpleFunction< VarContainer, OutType> > _grad_val; - std::vector< std::vector< PSimpleFunction< VarContainer, OutType> > > _hess_val; - - public: - - typedef typename PFuncBase< VarContainer, OutType>::size_type size_type; - - PFlexFunction() - { - - } - - PFlexFunction( const std::string &name, - const std::vector &var_name, - const std::vector &var_description, - const PSimpleFunction< VarContainer, OutType> &simplef, - const std::vector< PSimpleFunction< VarContainer, OutType> > &grad_simplef, - const std::vector< std::vector< PSimpleFunction< VarContainer, OutType> > > &hess_simplef) : - PFuncBase< VarContainer, OutType>(name, var_name, var_description), - _val(simplef), _grad_val(grad_simplef), _hess_val(hess_simplef) - { - check(); - } - - void clear() - { - this->_name = ""; - this->_var_name.clear(); - this->_var_description.clear(); - _val = PSimpleFunction< VarContainer, OutType>(); - _grad_val.clear(); - _hess_val.clear(); - } - - void set( const std::string &name, - const std::vector &var_name, - const std::vector &var_description, - const PSimpleFunction< VarContainer, OutType> &simplef, - const std::vector< PSimpleFunction< VarContainer, OutType> > &grad_simplef, - const std::vector< std::vector< PSimpleFunction< VarContainer, OutType> > > &hess_simplef) - { - this->_name = name; - this->_var_name = var_name; - this->_var_description = var_description; - _val = simplef; - _grad_val = grad_simplef; - _hess_val = hess_simplef; - - check(); - } + /// A class to create functions consisting of PSimpleFunctions + /// Used to create basis_functions from PSimpleFunctions for f, grad_f, and hess_f + /// + template + class PFlexFunction : public PFuncBase + { + PSimpleFunction _val; + std::vector> _grad_val; + std::vector>> _hess_val; + + public: + typedef typename PFuncBase::size_type size_type; + + PFlexFunction() + {} + + PFlexFunction(const std::string &name, + const std::vector &var_name, + const std::vector &var_description, + const PSimpleFunction &simplef, + const std::vector> &grad_simplef, + const std::vector>> + &hess_simplef) + : PFuncBase(name, var_name, var_description) + , _val(simplef) + , _grad_val(grad_simplef) + , _hess_val(hess_simplef) + { + check(); + } - PFlexFunction(const PFlexFunction &RHS ) - { - this->_name = RHS._name; - this->_var_name = RHS._var_name; - this->_var_description = RHS._var_description; - - _val = RHS._val; - _grad_val = RHS._grad_val; - _hess_val = RHS._hess_val; - } + void + clear() + { + this->_name = ""; + this->_var_name.clear(); + this->_var_description.clear(); + _val = PSimpleFunction(); + _grad_val.clear(); + _hess_val.clear(); + } + + void + set(const std::string &name, + const std::vector &var_name, + const std::vector &var_description, + const PSimpleFunction &simplef, + const std::vector> &grad_simplef, + const std::vector>> + &hess_simplef) + { + this->_name = name; + this->_var_name = var_name; + this->_var_description = var_description; + _val = simplef; + _grad_val = grad_simplef; + _hess_val = hess_simplef; - PFlexFunction& operator=(const PFlexFunction &RHS ) - { - this->_name = RHS._name; - this->_var_name = RHS._var_name; - this->_var_description = RHS._var_description; - - _val = RHS._val; - _grad_val = RHS._grad_val; - _hess_val = RHS._hess_val; - - } + check(); + } - ~PFlexFunction() - { - - } + PFlexFunction(const PFlexFunction &RHS) + { + this->_name = RHS._name; + this->_var_name = RHS._var_name; + this->_var_description = RHS._var_description; - PFlexFunction* clone() const - { - return new PFlexFunction(*this); - } + _val = RHS._val; + _grad_val = RHS._grad_val; + _hess_val = RHS._hess_val; + } - PSimpleFunction< VarContainer, OutType> simplefunction() const - { - return _val; - } + PFlexFunction & + operator=(const PFlexFunction &RHS) + { + this->_name = RHS._name; + this->_var_name = RHS._var_name; + this->_var_description = RHS._var_description; - PSimpleFunction< VarContainer, OutType> grad_simplefunction(size_type di) const - { - return _grad_val[di]; - } + _val = RHS._val; + _grad_val = RHS._grad_val; + _hess_val = RHS._hess_val; + } - PSimpleFunction< VarContainer, OutType> hess_simplefunction(size_type di, size_type dj) const - { - return _hess_val[di][dj]; - } + ~PFlexFunction() + {} - OutType operator()(const VarContainer &var) - { - return _val(var); - } + PFlexFunction * + clone() const + { + return new PFlexFunction(*this); + } - OutType grad(const VarContainer &var, size_type di) - { - return _grad_val[di](var); - } + PSimpleFunction + simplefunction() const + { + return _val; + } - OutType hess(const VarContainer &var, size_type di, size_type dj) - { - return _hess_val[di][dj](var); - } - - void eval(const VarContainer &var) - { - (*this)(var); - } + PSimpleFunction + grad_simplefunction(size_type di) const + { + return _grad_val[di]; + } - void eval_grad(const VarContainer &var) - { - for( size_type i=0; i<_grad_val.size(); i++) - _grad_val[i](var); - } + PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + return _hess_val[di][dj]; + } - void eval_hess(const VarContainer &var) - { - for( size_type i=0; i<_hess_val.size(); i++) - for( size_type j=0; j<_hess_val[i].size(); j++) - _hess_val[i][j](var); - } + OutType + operator()(const VarContainer &var) + { + return _val(var); + } + + OutType + grad(const VarContainer &var, size_type di) + { + return _grad_val[di](var); + } + + OutType + hess(const VarContainer &var, size_type di, size_type dj) + { + return _hess_val[di][dj](var); + } + + void + eval(const VarContainer &var) + { + (*this)(var); + } + + void + eval_grad(const VarContainer &var) + { + for (size_type i = 0; i < _grad_val.size(); i++) + _grad_val[i](var); + } + + void + eval_hess(const VarContainer &var) + { + for (size_type i = 0; i < _hess_val.size(); i++) + for (size_type j = 0; j < _hess_val[i].size(); j++) + _hess_val[i][j](var); + } + + OutType + operator()() const + { + return _val(); + } + + OutType + grad(size_type di) const + { + return _grad_val[di](); + } - OutType operator()() const + OutType + hess(size_type di, size_type dj) const + { + return _hess_val[di][dj](); + } + + private: + void + check() + { + size_type n = this->_var_name.size(); + if (this->_var_description.size() != n) { - return _val(); + std::cerr + << "Error in PFlexFunction. _var_name.size() != _var_description.size()." + << std::endl; + exit(1); } - - OutType grad(size_type di) const + if (_grad_val.size() != n) { - return _grad_val[di](); + std::cerr << "Error in PFlexFunction. _var_name.size() != _grad_val.size()." + << std::endl; + exit(1); } - - OutType hess(size_type di, size_type dj) const + if (_hess_val.size() != n) { - return _hess_val[di][dj](); + std::cerr << "Error in PFlexFunction. _var_name.size() != _hess_val.size()." + << std::endl; + exit(1); } - - private: - - void check() + for (size_type i = 0; i < _hess_val.size(); i++) { - size_type n = this->_var_name.size(); - if( this->_var_description.size() != n) + if (_hess_val[i].size() != n) { - std::cerr << "Error in PFlexFunction. _var_name.size() != _var_description.size()." << std::endl; - exit(1); - } - if( _grad_val.size() != n) - { - std::cerr << "Error in PFlexFunction. _var_name.size() != _grad_val.size()." << std::endl; - exit(1); - } - if( _hess_val.size() != n) - { - std::cerr << "Error in PFlexFunction. _var_name.size() != _hess_val.size()." << std::endl; - exit(1); - } - for( size_type i=0; i<_hess_val.size(); i++) - { - if( _hess_val[i].size() != n) - { - std::cerr << "Error in PFlexFunction. _var_name.size() != _hess_val[" << i << "].size()." << std::endl; - exit(1); - } + std::cerr << "Error in PFlexFunction. _var_name.size() != _hess_val[" << i + << "].size()." << std::endl; + exit(1); } } + } + }; - }; - -} - +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PFuncBase.hh b/include/IntegrationTools/pfunction/PFuncBase.hh index 56455ec1e..81618f32c 100644 --- a/include/IntegrationTools/pfunction/PFuncBase.hh +++ b/include/IntegrationTools/pfunction/PFuncBase.hh @@ -2,152 +2,187 @@ #ifndef PFuncBase_HH #define PFuncBase_HH -#include -#include -#include -#include -#include - #include "./PSimpleFunction.hh" +#include +#include +#include +#include +#include namespace PRISMS -{ - /// A Base class for a function, including grad & hess - /// - template< class VarContainer, class OutType> - class PFuncBase - { - protected: - std::string _name; - std::vector _var_name; - std::vector _var_description; - - public: - typedef std::vector::size_type size_type; - - PFuncBase(){} - - PFuncBase( const std::string &name, - const std::vector &var_name, - const std::vector &var_description) : - _name(name), _var_name(var_name), _var_description(var_description) - { - } - - virtual ~PFuncBase(){} - - std::string name() - { - return _name; - } - size_type size() const - { - return _var_name.size(); - } - std::vector var_name() - { - return _var_name; - } - std::string var_name(size_type i) - { - return _var_name[i]; - } - std::vector var_description() - { - return _var_description; - } - std::string var_description(size_type i) - { - return _var_description[i]; - } - - void is_derived_from_PFuncBase() const - { - return; - } - - virtual PFuncBase *clone() const - { - return new PFuncBase(*this); - } - - virtual PSimpleFunction simplefunction() const - { - undefined("PSimpleFunction simplefunction() const"); - return PSimpleFunction(); - } - - virtual PSimpleFunction grad_simplefunction(size_type di) const - { - undefined("PSimpleFunction grad_simplefunction() const"); - return PSimpleFunction(); - } - - virtual PSimpleFunction hess_simplefunction(size_type di, size_type dj) const - { - undefined("PSimpleFunction hess_simplefunction(size_type di, size_type dj) const"); - return PSimpleFunction(); - } - - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - virtual OutType operator()(const VarContainer &var) - { - undefined("OutType operator()(const VarContainer &var)"); - return OutType(); - } - virtual OutType grad(const VarContainer &var, size_type di) - { - undefined("OutType grad(const VarContainer &var, size_type di)"); - return OutType(); - } - virtual OutType hess(const VarContainer &var, size_type di, size_type dj) - { - undefined("OutType hess(const VarContainer &var, size_type di, size_type dj)"); - return OutType(); - } - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - virtual void eval(const VarContainer &var) - { - undefined("void eval_grad( const VarContainer &var)"); - } - virtual void eval_grad(const VarContainer &var) - { - undefined("void eval_grad( const VarContainer &var)"); - } - virtual void eval_hess(const VarContainer &var) - { - undefined("void eval_hess( const VarContainer &var)"); - } - - virtual OutType operator()() const - { - undefined("OutType operator()"); - return OutType(); - } - virtual OutType grad(size_type di) const - { - undefined("OutType grad(size_type di)"); - return OutType(); - } - virtual OutType hess(size_type di, size_type dj) const - { - undefined("OutType hess(size_type di, size_type dj)"); - return OutType(); - } - - private: - void undefined(std::string fname) const - { - std::string msg = "Error in PFuncBase '" + _name + "'.\n" + " The member function '" + fname + "' has not been defined.\n"; - throw std::runtime_error(msg); - } - - }; - -} +{ + /// A Base class for a function, including grad & hess + /// + template + class PFuncBase + { + protected: + std::string _name; + std::vector _var_name; + std::vector _var_description; + + public: + typedef std::vector::size_type size_type; + + PFuncBase() + {} + + PFuncBase(const std::string &name, + const std::vector &var_name, + const std::vector &var_description) + : _name(name) + , _var_name(var_name) + , _var_description(var_description) + {} + + virtual ~PFuncBase() + {} + + std::string + name() + { + return _name; + } + + size_type + size() const + { + return _var_name.size(); + } + + std::vector + var_name() + { + return _var_name; + } + + std::string + var_name(size_type i) + { + return _var_name[i]; + } + + std::vector + var_description() + { + return _var_description; + } + + std::string + var_description(size_type i) + { + return _var_description[i]; + } + + void + is_derived_from_PFuncBase() const + { + return; + } + + virtual PFuncBase * + clone() const + { + return new PFuncBase(*this); + } + + virtual PSimpleFunction + simplefunction() const + { + undefined("PSimpleFunction simplefunction() const"); + return PSimpleFunction(); + } + + virtual PSimpleFunction + grad_simplefunction(size_type di) const + { + undefined("PSimpleFunction grad_simplefunction() const"); + return PSimpleFunction(); + } + + virtual PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + undefined("PSimpleFunction hess_simplefunction(size_type " + "di, size_type dj) const"); + return PSimpleFunction(); + } + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + virtual OutType + operator()(const VarContainer &var) + { + undefined("OutType operator()(const VarContainer &var)"); + return OutType(); + } + + virtual OutType + grad(const VarContainer &var, size_type di) + { + undefined("OutType grad(const VarContainer &var, size_type di)"); + return OutType(); + } + + virtual OutType + hess(const VarContainer &var, size_type di, size_type dj) + { + undefined("OutType hess(const VarContainer &var, size_type di, size_type dj)"); + return OutType(); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + virtual void + eval(const VarContainer &var) + { + undefined("void eval_grad( const VarContainer &var)"); + } + + virtual void + eval_grad(const VarContainer &var) + { + undefined("void eval_grad( const VarContainer &var)"); + } + + virtual void + eval_hess(const VarContainer &var) + { + undefined("void eval_hess( const VarContainer &var)"); + } + + virtual OutType + operator()() const + { + undefined("OutType operator()"); + return OutType(); + } + + virtual OutType + grad(size_type di) const + { + undefined("OutType grad(size_type di)"); + return OutType(); + } + + virtual OutType + hess(size_type di, size_type dj) const + { + undefined("OutType hess(size_type di, size_type dj)"); + return OutType(); + } + + private: + void + undefined(std::string fname) const + { + std::string msg = "Error in PFuncBase '" + _name + "'.\n" + + " The member function '" + fname + "' has not been defined.\n"; + throw std::runtime_error(msg); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PFunction.hh b/include/IntegrationTools/pfunction/PFunction.hh index fbad8cb60..c2238b225 100644 --- a/include/IntegrationTools/pfunction/PFunction.hh +++ b/include/IntegrationTools/pfunction/PFunction.hh @@ -2,204 +2,235 @@ #ifndef PFunction_HH #define PFunction_HH -#include -#include -#include -#include - -#include "./PSimpleFunction.hh" #include "./PFuncBase.hh" +#include "./PSimpleFunction.hh" +#include +#include +#include +#include namespace PRISMS -{ - - /// A class that contains a ptr to a PFuncBase object - /// - like a smart ptr class - /// - same interface as PFuncBase - /// - allows for using PFuncBase objects polymorphically - /// without dereferencing and without worrying about new/delete - /// - VarContainer is either a scalar or something whose elements can be accessed with operator[] - /// - /// example: MyFuncA, MyFuncB, MyFuncC, etc. are defined: - /// template< class VarContainer> - /// MyFuncX : public PFuncBase - /// - /// // Then you can do things like this: - /// - /// MyFuncA > my_func_a; - /// MyFuncB > my_func_b; - /// - /// PFuncBase, double>* my_func_c_ptr; - /// my_func_c_ptr = new MyFuncC, double >(); - /// - /// PFunction, double > f, g, h; - /// - /// f = my_func_a; - /// f = my_func_b; - /// g.set(my_func_c_ptr->clone()); - /// h.set(my_func_c_ptr); - /// double result = f(3.0) + g(4.0) + h(5.0); - /// - /// - No deletions are used in this example. - /// PFunction::set makes PFunction the 'owner' of the MyFuncC object and it will delete it. - /// - template< class VarContainer, class OutType> - class PFunction - { - public: - - std::string name() const - { - return (*ptr).name(); - } - int size() const - { - return (*ptr).size(); - } - std::vector var_name() - { - return (*ptr).var_name(); - } - std::string var_name(int i) - { - return (*ptr).var_name(i); - } - std::vector var_description() - { - return (*ptr).var_description(); - } - std::string var_description(int i) - { - return (*ptr).var_description(i); - } - - PSimpleFunction simplefunction() const - { - return (*ptr).simplefunction(); - } - - PSimpleFunction grad_simplefunction(int di) const - { - return (*ptr).grad_simplefunction(di); - } - - PSimpleFunction hess_simplefunction(int di, int dj) const - { - return (*ptr).hess_simplefunction(di, dj); - } +{ + + /// A class that contains a ptr to a PFuncBase object + /// - like a smart ptr class + /// - same interface as PFuncBase + /// - allows for using PFuncBase objects polymorphically + /// without dereferencing and without worrying about new/delete + /// - VarContainer is either a scalar or something whose elements can be accessed with + /// operator[] + /// + /// example: MyFuncA, MyFuncB, MyFuncC, etc. are defined: + /// template< class VarContainer> + /// MyFuncX : public PFuncBase + /// + /// // Then you can do things like this: + /// + /// MyFuncA > my_func_a; + /// MyFuncB > my_func_b; + /// + /// PFuncBase, double>* my_func_c_ptr; + /// my_func_c_ptr = new MyFuncC, double >(); + /// + /// PFunction, double > f, g, h; + /// + /// f = my_func_a; + /// f = my_func_b; + /// g.set(my_func_c_ptr->clone()); + /// h.set(my_func_c_ptr); + /// double result = f(3.0) + g(4.0) + h(5.0); + /// + /// - No deletions are used in this example. + /// PFunction::set makes PFunction the 'owner' of the MyFuncC object and it will + /// delete it. + /// + template + class PFunction + { + public: + std::string + name() const + { + return (*ptr).name(); + } - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - OutType operator()(const VarContainer &var) - { - return (*ptr)(var); - } - OutType grad(const VarContainer &var, int di) - { - return (*ptr).grad(var, di); - } - OutType hess(const VarContainer &var, int di, int dj) - { - return (*ptr).hess(var, di, dj); - } + int + size() const + { + return (*ptr).size(); + } - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - void eval(const VarContainer &var) - { - return (*ptr).eval(var); - } - void eval_grad(const VarContainer &var) - { - return (*ptr).eval_grad(var); - } - void eval_hess(const VarContainer &var) - { - return (*ptr).eval_hess(var); - } + std::vector + var_name() + { + return (*ptr).var_name(); + } - OutType operator()() const - { - return (*ptr)(); - } - OutType grad(int di) const - { - return (*ptr).grad(di); - } - OutType hess(int di, int dj) const - { - return (*ptr).hess(di, dj); - } + std::string + var_name(int i) + { + return (*ptr).var_name(i); + } + std::vector + var_description() + { + return (*ptr).var_description(); + } - // PFunction unique members ------------------------------------------ + std::string + var_description(int i) + { + return (*ptr).var_description(i); + } - PFunction &operator=(const PFunction &RHS) - { - if(ptr != NULL) - delete ptr; - ptr = RHS.ptr->clone(); - return *this; - } + PSimpleFunction + simplefunction() const + { + return (*ptr).simplefunction(); + } - template PFunction &operator=(const T &RHS) - { - RHS.is_derived_from_PFuncBase(); + PSimpleFunction + grad_simplefunction(int di) const + { + return (*ptr).grad_simplefunction(di); + } - if(ptr != NULL) - delete ptr; - ptr = RHS.clone(); - return *this; - } + PSimpleFunction + hess_simplefunction(int di, int dj) const + { + return (*ptr).hess_simplefunction(di, dj); + } - // If you use this, PFunction becomes the 'owner' of the function RHS points to - // and it will delete it - PFunction &set( PFuncBase *RHS) - { - if(RHS == NULL) - { - std::cout << "Error in PFunction::set. RHS == NULL." << std::endl; - exit(1); - } - if(ptr != NULL) - delete ptr; - ptr = RHS; - return *this; - } + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + OutType + operator()(const VarContainer &var) + { + return (*ptr)(var); + } - PFunction() - { - ptr = NULL; - } + OutType + grad(const VarContainer &var, int di) + { + return (*ptr).grad(var, di); + } - PFunction(const PFunction &RHS) - { - if( RHS.ptr != NULL) - ptr = RHS.ptr->clone(); - else - ptr = NULL; - } - - template PFunction(const T &RHS) - { - RHS.is_derived_from_PFuncBase(); + OutType + hess(const VarContainer &var, int di, int dj) + { + return (*ptr).hess(var, di, dj); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + void + eval(const VarContainer &var) + { + return (*ptr).eval(var); + } - ptr = RHS.clone(); - - } + void + eval_grad(const VarContainer &var) + { + return (*ptr).eval_grad(var); + } + + void + eval_hess(const VarContainer &var) + { + return (*ptr).eval_hess(var); + } + + OutType + operator()() const + { + return (*ptr)(); + } + + OutType + grad(int di) const + { + return (*ptr).grad(di); + } + + OutType + hess(int di, int dj) const + { + return (*ptr).hess(di, dj); + } - ~PFunction() + // PFunction unique members ------------------------------------------ + + PFunction & + operator=(const PFunction &RHS) + { + if (ptr != NULL) + delete ptr; + ptr = RHS.ptr->clone(); + return *this; + } + + template + PFunction & + operator=(const T &RHS) + { + RHS.is_derived_from_PFuncBase(); + + if (ptr != NULL) + delete ptr; + ptr = RHS.clone(); + return *this; + } + + // If you use this, PFunction becomes the 'owner' of the function RHS points to + // and it will delete it + PFunction & + set(PFuncBase *RHS) + { + if (RHS == NULL) { - if(ptr != NULL) - delete ptr; + std::cout << "Error in PFunction::set. RHS == NULL." << std::endl; + exit(1); } + if (ptr != NULL) + delete ptr; + ptr = RHS; + return *this; + } + + PFunction() + { + ptr = NULL; + } - private: - PFuncBase *ptr; + PFunction(const PFunction &RHS) + { + if (RHS.ptr != NULL) + ptr = RHS.ptr->clone(); + else + ptr = NULL; + } + + template + PFunction(const T &RHS) + { + RHS.is_derived_from_PFuncBase(); - }; + ptr = RHS.clone(); + } + + ~PFunction() + { + if (ptr != NULL) + delete ptr; + } -} + private: + PFuncBase *ptr; + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PSeriesFunction.hh b/include/IntegrationTools/pfunction/PSeriesFunction.hh index b8d1eed79..a001ae5fc 100644 --- a/include/IntegrationTools/pfunction/PSeriesFunction.hh +++ b/include/IntegrationTools/pfunction/PSeriesFunction.hh @@ -2,566 +2,617 @@ #ifndef PSeriesFunction_HH #define PSeriesFunction_HH -#include -#include -#include -#include - -#include "./PFuncBase.hh" -#include "./PBasisSet.hh" #include "../datastruc/PNDArray.hh" +#include "./PBasisSet.hh" +#include "./PFuncBase.hh" +#include +#include +#include +#include namespace PRISMS { - /// This class is for a series function that is the tensor product of indepedent bases - /// - /// Example for order-3 tensor: - /// - /// f(x,y,z) = sum_i sum_j sum_k T_ijk * phix_i(x) * phiy_j(y) * phiz_k(z) - /// - /// summations run from i to _dim[i] - /// x, y, z are independent variables (dx/dy = 0, etc.) - /// - /// Template parameters: - /// - InType: input type for BasisSets (must be the same for all BasisSets) - /// - OutType: output type for SeriesFunction and BasisSets - /// - VarContainer: container for input - /// - IndexContainer: container for tensor indices - /// - /// Example template parameters are: , std::vector > + /// This class is for a series function that is the tensor product of indepedent bases + /// + /// Example for order-3 tensor: + /// + /// f(x,y,z) = sum_i sum_j sum_k T_ijk * phix_i(x) * phiy_j(y) * phiz_k(z) + /// + /// summations run from i to _dim[i] + /// x, y, z are independent variables (dx/dy = 0, etc.) + /// + /// Template parameters: + /// - InType: input type for BasisSets (must be the same for all BasisSets) + /// - OutType: output type for SeriesFunction and BasisSets + /// - VarContainer: container for input + /// - IndexContainer: container for tensor indices + /// + /// Example template parameters are: , + /// std::vector > + /// + template + class PSeriesFunction : public PFuncBase + { + protected: + // coefficient tensor + PNDArray _coeff; + + OutType _identity_val; + OutType _zero_val; + + // basis sets + std::vector> _basis_set; + + // evaluated values + OutType _val; + std::vector _grad_val; + std::vector> _hess_val; + + public: + PSeriesFunction(OutType zero, OutType identity) + { + _zero_val = zero; + _identity_val = identity; + } + + PSeriesFunction(OutType zero, + OutType identity, + const std::vector> &basis_set) + { + _zero_val = zero; + _identity_val = identity; + + set(basis_set); + } + + virtual ~PSeriesFunction() + { + clear(); + } + + void + clear() + { + _basis_set.clear(); + + _coeff.clear(); + + _val = _zero_val; + _grad_val.clear(); + _hess_val.clear(); + } + + /// generate the PSeriesFunction, by cloning 'basis_set' and resizing everything to + /// match /// - template< class InType, class OutType, class VarContainer, class IndexContainer> - class PSeriesFunction : public PFuncBase - { - protected: - - // coefficient tensor - PNDArray _coeff; - - OutType _identity_val; - OutType _zero_val; - - // basis sets - std::vector< PBasisSet > _basis_set; - - // evaluated values - OutType _val; - std::vector _grad_val; - std::vector< std::vector< OutType> > _hess_val; - - public: - PSeriesFunction(OutType zero, OutType identity) - { - _zero_val = zero; - _identity_val = identity; - } - - PSeriesFunction( OutType zero, OutType identity, const std::vector > &basis_set) - { - _zero_val = zero; - _identity_val = identity; + void + set(const std::vector> &basis_set) + { + // std::cout << "begin PSeriesFunction::set()" << std::endl; + clear(); - set(basis_set); - } - - virtual ~PSeriesFunction() - { - clear(); - } - - void clear() + std::vector dim(basis_set.size()); + for (int i = 0; i < basis_set.size(); i++) + dim[i] = basis_set[i].size(); + _coeff.resize(dim); + + // std::cout << " resize: " << basis_set.size() << std::endl; + _basis_set.resize(basis_set.size()); + + _grad_val.resize(_coeff.order()); + _hess_val.resize(_coeff.order()); + + for (int i = 0; i < _coeff.order(); i++) { - _basis_set.clear(); - - _coeff.clear(); - - _val = _zero_val; - _grad_val.clear(); - _hess_val.clear(); - + // std::cout << " i: " << i << std::endl; + _hess_val[i].resize(_coeff.order()); + // std::cout << " copy basis set" << std::endl; + _basis_set[i] = basis_set[i]; } - - /// generate the PSeriesFunction, by cloning 'basis_set' and resizing everything to match - /// - void set( const std::vector > &basis_set) + } + + PNDArray & + coeff() + { + return _coeff; + } + + virtual PSeriesFunction * + clone() const + { + return new PSeriesFunction(*this); + } + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + virtual OutType + operator()(const VarContainer &var) + { + // evaluate basis functions needed + eval_basis(var); + + // evaluate basis function products, multiply by _coeff_tensor & sum + return calc_val(); + } + + virtual OutType + grad(const VarContainer &var, int di) + { + // evaluate basis functions needed + for (int i = 0; i < _coeff.order(); i++) { - //std::cout << "begin PSeriesFunction::set()" << std::endl; - clear(); - - std::vector dim(basis_set.size()); - for( int i=0; i& coeff() - { - return _coeff; } - virtual PSeriesFunction *clone() const - { - return new PSeriesFunction(*this); - } + // evaluate basis function products, multiply by _coeff_tensor & sum + return calc_grad_val(di); + } - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - virtual OutType operator()(const VarContainer &var) - { - // evaluate basis functions needed - eval_basis(var); - - // evaluate basis function products, multiply by _coeff_tensor & sum - return calc_val(); - } - - virtual OutType grad(const VarContainer &var, int di) + virtual OutType + hess(const VarContainer &var, int di, int dj) + { + // evaluate basis functions needed + for (int i = 0; i < _coeff.order(); i++) { - // evaluate basis functions needed - for( int i=0; i<_coeff.order(); i++) + if (i == di || i == dj) { - if( i == di) - { - _basis_set[i].eval_grad(var[i]); + if (di == dj) + { + _basis_set[i].eval_hess(var[i]); } - else + else { - _basis_set[i](var[i]); + _basis_set[i].eval_grad(var[i]); } } - - // evaluate basis function products, multiply by _coeff_tensor & sum - return calc_grad_val(di); - } - - virtual OutType hess(const VarContainer &var, int di, int dj) - { - // evaluate basis functions needed - for( int i=0; i<_coeff.order(); i++) + else { - if( i == di || i == dj) - { - if( di == dj) - { - _basis_set[i].eval_hess(var[i]); - } - else - { - _basis_set[i].eval_grad(var[i]); - } - } - else - { - _basis_set[i](var[i]); - } + _basis_set[i](var[i]); } - - // evaluate basis function products, multiply by _coeff_tensor & sum - return calc_hess_val(di,dj); - } - - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - virtual void eval(const VarContainer &var) - { - (*this)(var); - } - virtual void eval_grad(const VarContainer &var) - { - eval_basis(var); - - eval_basis_grad(var); - - for( int i=0; i<_coeff.order(); i++) - (*this).calc_grad_val(i); - } - virtual void eval_hess(const VarContainer &var) - { - eval_basis(var); - - eval_basis_grad(var); - - eval_basis_hess(var); - - - for( int i=0; i<_coeff.order(); i++) - for( int j=0; j<_coeff.order(); j++) - (*this).calc_hess_val(i,j); } - virtual OutType operator()() const - { - return _val; - } - virtual OutType grad(int di) const - { - return _grad_val[di]; - } - virtual OutType hess(int di, int dj) const - { - return _hess_val[di][dj]; - } + // evaluate basis function products, multiply by _coeff_tensor & sum + return calc_hess_val(di, dj); + } + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + virtual void + eval(const VarContainer &var) + { + (*this)(var); + } - // Functions for evaluating basis functions & their derivatives: + virtual void + eval_grad(const VarContainer &var) + { + eval_basis(var); - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value + eval_basis_grad(var); - // use basis index and term index for individual basis function - virtual OutType basis(int bindex, int term, const VarContainer &var) - { - return _basis_set[bindex](term,var[bindex]); - } - virtual OutType basis_grad(int bindex, int term, const VarContainer &var) - { - return _basis_set[bindex].grad(term, var[bindex]); - } - virtual OutType basis_hess(int bindex, int term, const VarContainer &var) - { - return _basis_set[bindex].hess(term, var[bindex]); - } + for (int i = 0; i < _coeff.order(); i++) + (*this).calc_grad_val(i); + } - // or use tensor indices to evaluate basis function product - virtual OutType basis(const IndexContainer &term, const VarContainer &var) - { - // evaluate basis functions needed - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) - tmp *= _basis_set[i](term[i],var[i]); - - return tmp; - } - virtual OutType basis_grad(const IndexContainer &term, const VarContainer &var, int di) + virtual void + eval_hess(const VarContainer &var) + { + eval_basis(var); + + eval_basis_grad(var); + + eval_basis_hess(var); + + for (int i = 0; i < _coeff.order(); i++) + for (int j = 0; j < _coeff.order(); j++) + (*this).calc_hess_val(i, j); + } + + virtual OutType + operator()() const + { + return _val; + } + + virtual OutType + grad(int di) const + { + return _grad_val[di]; + } + + virtual OutType + hess(int di, int dj) const + { + return _hess_val[di][dj]; + } + + // Functions for evaluating basis functions & their derivatives: + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + + // use basis index and term index for individual basis function + virtual OutType + basis(int bindex, int term, const VarContainer &var) + { + return _basis_set[bindex](term, var[bindex]); + } + + virtual OutType + basis_grad(int bindex, int term, const VarContainer &var) + { + return _basis_set[bindex].grad(term, var[bindex]); + } + + virtual OutType + basis_hess(int bindex, int term, const VarContainer &var) + { + return _basis_set[bindex].hess(term, var[bindex]); + } + + // or use tensor indices to evaluate basis function product + virtual OutType + basis(const IndexContainer &term, const VarContainer &var) + { + // evaluate basis functions needed + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) + tmp *= _basis_set[i](term[i], var[i]); + + return tmp; + } + + virtual OutType + basis_grad(const IndexContainer &term, const VarContainer &var, int di) + { + // evaluate basis functions needed + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) { - // evaluate basis functions needed - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) + if (i == di) { - if( i == di) - { - tmp *= _basis_set[i].grad(term[i], var[i]); - } - else - { - tmp *= _basis_set[i](term[i], var[i]); - } + tmp *= _basis_set[i].grad(term[i], var[i]); + } + else + { + tmp *= _basis_set[i](term[i], var[i]); } - - return tmp; } - virtual OutType basis_hess(const IndexContainer &term, const VarContainer &var, int di, int dj) + + return tmp; + } + + virtual OutType + basis_hess(const IndexContainer &term, const VarContainer &var, int di, int dj) + { + // evaluate basis functions needed + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) { - // evaluate basis functions needed - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) + if (i == di || i == dj) { - if( i == di || i == dj) - { - if( di == dj) - { - tmp *= _basis_set[i].hess(term[i], var[i]); - } - else - { - tmp *= _basis_set[i].grad(term[i], var[i]); - } + if (di == dj) + { + tmp *= _basis_set[i].hess(term[i], var[i]); } - else + else { - tmp *= _basis_set[i](term[i], var[i]); + tmp *= _basis_set[i].grad(term[i], var[i]); } } - - return tmp; + else + { + tmp *= _basis_set[i](term[i], var[i]); + } } - // ---------------------------------------------------------- - // Use these functions to evaluate all basis functions, - // then use following methods to access results. + return tmp; + } - virtual void eval_basis(const VarContainer &var) - { - // evaluate basis functions - for( int i=0; i<_coeff.order(); i++) - _basis_set[i].eval(var[i]); - - } - virtual void eval_basis(const VarContainer &var, int i) - { - // evaluate basis functions - _basis_set[i].eval(var[i]); - - } - virtual void eval_basis_grad(const VarContainer &var) - { - // evaluate basis grad functions - for( int i=0; i<_coeff.order(); i++) - _basis_set[i].eval_grad(var[i]); - } - virtual void eval_basis_grad(const VarContainer &var, int i) - { - // evaluate basis grad functions - _basis_set[i].eval_grad(var[i]); - } - virtual void eval_basis_hess(const VarContainer &var) - { - // evaluate basis hess functions - for( int i=0; i<_coeff.order(); i++) - _basis_set[i].eval_hess(var[i]); - } - virtual void eval_basis_hess(const VarContainer &var, int i) - { - // evaluate basis hess functions - _basis_set[i].eval_hess(var[i]); - } + // ---------------------------------------------------------- + // Use these functions to evaluate all basis functions, + // then use following methods to access results. - // use basis index and term index for individual basis function - virtual OutType basis(int bindex, int term) const - { - return _basis_set[bindex](term); - } - virtual OutType basis_grad(int bindex, int term) const - { - return _basis_set[bindex].grad(term); - } - virtual OutType basis_hess(int bindex, int term) const - { - return _basis_set[bindex].hess(term); - } + virtual void + eval_basis(const VarContainer &var) + { + // evaluate basis functions + for (int i = 0; i < _coeff.order(); i++) + _basis_set[i].eval(var[i]); + } - // or use tensor indices to evaluate basis function product - virtual OutType basis(const IndexContainer &term) const - { - // evaluate basis function product - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) - tmp *= _basis_set[i](term[i]); - - return tmp; - } - virtual OutType basis_grad(const IndexContainer &term, int di) const + virtual void + eval_basis(const VarContainer &var, int i) + { + // evaluate basis functions + _basis_set[i].eval(var[i]); + } + + virtual void + eval_basis_grad(const VarContainer &var) + { + // evaluate basis grad functions + for (int i = 0; i < _coeff.order(); i++) + _basis_set[i].eval_grad(var[i]); + } + + virtual void + eval_basis_grad(const VarContainer &var, int i) + { + // evaluate basis grad functions + _basis_set[i].eval_grad(var[i]); + } + + virtual void + eval_basis_hess(const VarContainer &var) + { + // evaluate basis hess functions + for (int i = 0; i < _coeff.order(); i++) + _basis_set[i].eval_hess(var[i]); + } + + virtual void + eval_basis_hess(const VarContainer &var, int i) + { + // evaluate basis hess functions + _basis_set[i].eval_hess(var[i]); + } + + // use basis index and term index for individual basis function + virtual OutType + basis(int bindex, int term) const + { + return _basis_set[bindex](term); + } + + virtual OutType + basis_grad(int bindex, int term) const + { + return _basis_set[bindex].grad(term); + } + + virtual OutType + basis_hess(int bindex, int term) const + { + return _basis_set[bindex].hess(term); + } + + // or use tensor indices to evaluate basis function product + virtual OutType + basis(const IndexContainer &term) const + { + // evaluate basis function product + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) + tmp *= _basis_set[i](term[i]); + + return tmp; + } + + virtual OutType + basis_grad(const IndexContainer &term, int di) const + { + // evaluate basis function product + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) { - // evaluate basis function product - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) + if (i == di) { - if( i == di) - { - tmp *= _basis_set[i].grad(term[i]); - } - else - { - tmp *= _basis_set[i](term[i]); - } + tmp *= _basis_set[i].grad(term[i]); + } + else + { + tmp *= _basis_set[i](term[i]); } - - return tmp; } - virtual OutType basis_hess(const IndexContainer &term, int di, int dj) const + + return tmp; + } + + virtual OutType + basis_hess(const IndexContainer &term, int di, int dj) const + { + // evaluate basis function product + OutType tmp = _identity_val; + for (int i = 0; i < _coeff.order(); i++) { - // evaluate basis function product - OutType tmp = _identity_val; - for( int i=0; i<_coeff.order(); i++) + if (i == di || i == dj) { - if( i == di || i == dj) - { - if( di == dj) - { - tmp *= _basis_set[i].hess(term[i]); - } - else - { - tmp *= _basis_set[i].grad(term[i]); - } + if (di == dj) + { + tmp *= _basis_set[i].hess(term[i]); } - else + else { - tmp *= _basis_set[i](term[i]); + tmp *= _basis_set[i].grad(term[i]); } } - - return tmp; - } - - // ---------------------------------------------------------- - // Read and write routines: - // Write routines only uses 'getters', they assume you have - // already evaluated the necessary basis functions - // - // These are not included in PExtern - - void print_basis(std::ostream &sout) - { - std::vector term = std::vector(_coeff.order(),0); - for( int i=0; i<_coeff.volume(); i++) + else { - _coeff.tensor_indices(i,term); - for( int j=0; j term = std::vector(_coeff.order(), 0); + for (int i = 0; i < _coeff.volume(); i++) + { + _coeff.tensor_indices(i, term); + for (int j = 0; j < term.size(); j++) + sout << term[j] << " "; + sout << basis(term) << "\n"; + } + } + + void + print_basis_grad(std::ostream &sout, int di) + { + std::vector term = std::vector(_coeff.order(), 0); + for (int i = 0; i < _coeff.volume(); i++) { - std::vector term = std::vector(_coeff.order(),0); - for( int i=0; i<_coeff.volume(); i++) - { - _coeff.tensor_indices(i,term); - for( int j=0; j term = std::vector(_coeff.order(), 0); + for (int i = 0; i < _coeff.volume(); i++) { - std::vector term = std::vector(_coeff.order(),0); - for( int i=0; i<_coeff.volume(); i++) - { - _coeff.tensor_indices(i,term); - for( int j=0; j term = std::vector(_coeff.order(), 0); + for (int i = 0; i < _coeff.volume(); i++) { - std::vector term = std::vector(_coeff.order(),0); - for( int i=0; i<_coeff.volume(); i++) - { - _coeff.tensor_indices(i,term); - for( int j=0; j term = std::vector(_coeff.order(), 0); + for (int i = 0; i < _coeff.volume(); i++) { - // no error checking, - // assumes format identical to print_coeff output - std::vector term = std::vector(_coeff.order(),0); - for( int i=0; i<_coeff.volume(); i++) + // ignores tensor indices + for (int j = 0; j < term.size(); j++) { - // ignores tensor indices - for( int j=0; j> term[j]; - } - sin >> _coeff(i); + sin >> term[j]; } + sin >> _coeff(i); } - - - /// - private: - - // ---------------------------------------------------------- - // These assume you've evaluated the necessary _basis, _basis_grad, and _basis_hess functions - // and just need to take products & sum - - // evaluate basis function products, multiply by _coeff_tensor & sum - OutType calc_val() + } + + /// + + private: + // ---------------------------------------------------------- + // These assume you've evaluated the necessary _basis, _basis_grad, and _basis_hess + // functions + // and just need to take products & sum + + // evaluate basis function products, multiply by _coeff_tensor & sum + OutType + calc_val() + { + std::vector tindex(_coeff.order()); + OutType tmp; + _val = _zero_val; + for (int i = 0; i < _coeff.volume(); i++) { - std::vector tindex(_coeff.order()); - OutType tmp; - _val = _zero_val; - for( int i=0; i<_coeff.volume(); i++) + tmp = _coeff(i); + _coeff.tensor_indices(i, tindex); + + for (int j = 0; j < _coeff.order(); j++) { - tmp = _coeff(i); - _coeff.tensor_indices(i, tindex); - - for( int j=0; j<_coeff.order(); j++) - { - tmp *= _basis_set[j](tindex[j]); - } - _val += tmp; + tmp *= _basis_set[j](tindex[j]); } - - return _val; + _val += tmp; } - - // evaluate basis function products, multiply by _coeff_tensor & sum - OutType calc_grad_val(int di) - { - std::vector tindex(_coeff.order()); - OutType tmp; - _grad_val[di] = _zero_val; - for( int i=0; i<_coeff.volume(); i++) + + return _val; + } + + // evaluate basis function products, multiply by _coeff_tensor & sum + OutType + calc_grad_val(int di) + { + std::vector tindex(_coeff.order()); + OutType tmp; + _grad_val[di] = _zero_val; + for (int i = 0; i < _coeff.volume(); i++) + { + tmp = _coeff(i); + _coeff.tensor_indices(i, tindex); + for (int j = 0; j < _coeff.order(); j++) { - tmp = _coeff(i); - _coeff.tensor_indices(i, tindex); - for( int j=0; j<_coeff.order(); j++) + if (j == di) { - if( j == di) - { - tmp *= _basis_set[j].grad(tindex[j]); - } - else - { - tmp *= _basis_set[j](tindex[j]); - } + tmp *= _basis_set[j].grad(tindex[j]); + } + else + { + tmp *= _basis_set[j](tindex[j]); } - - _grad_val[di] += tmp; - } - return _grad_val[di]; + + _grad_val[di] += tmp; } - - // evaluate basis function products, multiply by _coeff_tensor & sum - OutType calc_hess_val(int di, int dj) - { - std::vector tindex(_coeff.order()); - OutType tmp; - _hess_val[di][dj] = _zero_val; - for( int i=0; i<_coeff.volume(); i++) + return _grad_val[di]; + } + + // evaluate basis function products, multiply by _coeff_tensor & sum + OutType + calc_hess_val(int di, int dj) + { + std::vector tindex(_coeff.order()); + OutType tmp; + _hess_val[di][dj] = _zero_val; + for (int i = 0; i < _coeff.volume(); i++) + { + tmp = _coeff(i); + _coeff.tensor_indices(i, tindex); + for (int j = 0; j < _coeff.order(); j++) { - tmp = _coeff(i); - _coeff.tensor_indices(i, tindex); - for( int j=0; j<_coeff.order(); j++) + if (j == di || j == dj) { - if( j == di || j == dj) + if (di == dj) { - if( di == dj) - { - tmp *= _basis_set[j].hess(tindex[j]); - } - else - { - tmp *= _basis_set[j].grad(tindex[j]); - } + tmp *= _basis_set[j].hess(tindex[j]); } - else + else { - tmp *= _basis_set[j](tindex[j]); + tmp *= _basis_set[j].grad(tindex[j]); } } - _hess_val[di][dj] += tmp; + else + { + tmp *= _basis_set[j](tindex[j]); + } } - return _hess_val[di][dj]; + _hess_val[di][dj] += tmp; } - }; - - -} + return _hess_val[di][dj]; + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PSimpleBase.hh b/include/IntegrationTools/pfunction/PSimpleBase.hh index 055547e74..041bd1663 100644 --- a/include/IntegrationTools/pfunction/PSimpleBase.hh +++ b/include/IntegrationTools/pfunction/PSimpleBase.hh @@ -2,73 +2,98 @@ #ifndef PSimpleBase_HH #define PSimpleBase_HH -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include namespace PRISMS { - /// Base classes for functions that can be hard-coded, - /// then shared and used elsewhere - - /// A simple expression evaluator - /// - template< class VarContainer, class OutType> - class PSimpleBase + /// Base classes for functions that can be hard-coded, + /// then shared and used elsewhere + + /// A simple expression evaluator + /// + template + class PSimpleBase + { + protected: + std::string _name; + OutType _val; + + public: + virtual ~PSimpleBase() {}; + + std::string + name() const + { + return _name; + } + + OutType + operator()(const VarContainer &var) + { + return _val = eval(var); + } + + OutType + operator()() const + { + return _val; + } + + void + is_derived_from_PSimpleBase() const + { + return; + } + + virtual PSimpleBase * + clone() const + { + return new PSimpleBase(*this); + } + + virtual std::string + csrc() const + { + undefined("std::string csrc()"); + return std::string(); + } + + virtual std::string + sym() const + { + undefined("std::string sym()"); + return std::string(); + } + + virtual std::string + latex() const + { + undefined("std::string latex()"); + return std::string(); + } + + private: + virtual OutType + eval(const VarContainer &var) const + { + undefined("OutType eval( const VarContainer &var)"); + return OutType(); + } + + void + undefined(std::string fname) const { - protected: - std::string _name; - OutType _val; - - public: - virtual ~PSimpleBase(){}; - - std::string name() const { return _name;} - OutType operator()( const VarContainer &var){ return _val = eval(var);} - OutType operator()() const { return _val;} - - void is_derived_from_PSimpleBase() const - { - return; - } - - virtual PSimpleBase* clone() const - { - return new PSimpleBase(*this); - } - - virtual std::string csrc() const - { - undefined("std::string csrc()"); - return std::string(); - } - - virtual std::string sym() const - { - undefined("std::string sym()"); - return std::string(); - } - - virtual std::string latex() const - { - undefined("std::string latex()"); - return std::string(); - } - - private: - virtual OutType eval( const VarContainer &var) const { undefined("OutType eval( const VarContainer &var)"); return OutType();} - - void undefined(std::string fname) const - { - std::string msg = "Error in PSimpleBase '" + _name + "'.\n" + " The member function '" + fname + "' has not been defined.\n"; - throw std::runtime_error(msg); - } - }; - -} + std::string msg = "Error in PSimpleBase '" + _name + "'.\n" + + " The member function '" + fname + "' has not been defined.\n"; + throw std::runtime_error(msg); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/pfunction/PSimpleFunction.hh b/include/IntegrationTools/pfunction/PSimpleFunction.hh index d85ae32b6..efd93227d 100644 --- a/include/IntegrationTools/pfunction/PSimpleFunction.hh +++ b/include/IntegrationTools/pfunction/PSimpleFunction.hh @@ -2,128 +2,135 @@ #ifndef PSimpleFunction_HH #define PSimpleFunction_HH -#include -#include -#include -#include - #include "./PSimpleBase.hh" +#include +#include +#include +#include namespace PRISMS { - template< class VarContainer, class OutType> - class PSimpleFunction + template + class PSimpleFunction + { + private: + PSimpleBase *ptr; + + public: + std::string + name() const { - private: - PSimpleBase *ptr; - - public: - - std::string name() const - { - return (*ptr).name(); - } - std::string csrc() const - { - return (*ptr).csrc(); - } - std::string sym() const - { - return (*ptr).sym(); - } - std::string latex() const - { - return (*ptr).latex(); - } - + return (*ptr).name(); + } - // ---------------------------------------------------------- - // Use this function if you want to evaluate, - // return and store result - OutType operator()(const VarContainer &var) - { - return (*ptr)(var); - } - - // ---------------------------------------------------------- - // Then use 'get' methods to access results later - void eval(const VarContainer &var) - { - (*ptr)(var); - } - - OutType operator()() const - { - return (*ptr)(); - } + std::string + csrc() const + { + return (*ptr).csrc(); + } - // PFunction unique members ------------------------------------------ + std::string + sym() const + { + return (*ptr).sym(); + } - PSimpleFunction& operator=(const PSimpleFunction &RHS) - { - if(ptr != NULL) - delete ptr; - ptr = RHS.ptr->clone(); - return *this; - } + std::string + latex() const + { + return (*ptr).latex(); + } + + // ---------------------------------------------------------- + // Use this function if you want to evaluate, + // return and store result + OutType + operator()(const VarContainer &var) + { + return (*ptr)(var); + } - template - PSimpleFunction& operator=(const T &RHS) - { - RHS.is_derived_from_PSimpleBase(); + // ---------------------------------------------------------- + // Then use 'get' methods to access results later + void + eval(const VarContainer &var) + { + (*ptr)(var); + } - if(ptr != NULL) - delete ptr; - ptr = RHS.clone(); - return *this; - } + OutType + operator()() const + { + return (*ptr)(); + } - // If you use this, PSimpleFunction becomes the 'owner' of the function RHS points to - // and it will delete it - PSimpleFunction& set( PSimpleBase *RHS) - { - if(RHS == NULL) - { - std::cout << "Error in PSimpleFunction::set. RHS == NULL." << std::endl; - exit(1); - } - if(ptr != NULL) - delete ptr; - ptr = RHS; - return *this; - } - - PSimpleFunction() - { - ptr = NULL; - } + // PFunction unique members ------------------------------------------ - PSimpleFunction(const PSimpleFunction &RHS) + PSimpleFunction & + operator=(const PSimpleFunction &RHS) + { + if (ptr != NULL) + delete ptr; + ptr = RHS.ptr->clone(); + return *this; + } + + template + PSimpleFunction & + operator=(const T &RHS) + { + RHS.is_derived_from_PSimpleBase(); + + if (ptr != NULL) + delete ptr; + ptr = RHS.clone(); + return *this; + } + + // If you use this, PSimpleFunction becomes the 'owner' of the function RHS points to + // and it will delete it + PSimpleFunction & + set(PSimpleBase *RHS) + { + if (RHS == NULL) { - if( RHS.ptr != NULL) - ptr = RHS.ptr->clone(); - else - ptr = NULL; + std::cout << "Error in PSimpleFunction::set. RHS == NULL." << std::endl; + exit(1); } - - template PSimpleFunction(const T &RHS) - { - RHS.is_derived_from_PSimpleBase(); + if (ptr != NULL) + delete ptr; + ptr = RHS; + return *this; + } - ptr = RHS.clone(); - - } + PSimpleFunction() + { + ptr = NULL; + } - ~PSimpleFunction() - { - if(ptr != NULL) - delete ptr; - } + PSimpleFunction(const PSimpleFunction &RHS) + { + if (RHS.ptr != NULL) + ptr = RHS.ptr->clone(); + else + ptr = NULL; + } + + template + PSimpleFunction(const T &RHS) + { + RHS.is_derived_from_PSimpleBase(); - - }; + ptr = RHS.clone(); + } -} + ~PSimpleFunction() + { + if (ptr != NULL) + delete ptr; + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/piecewise/PPieceWiseFuncBase.hh b/include/IntegrationTools/piecewise/PPieceWiseFuncBase.hh index 10c701e9f..bd802ab74 100644 --- a/include/IntegrationTools/piecewise/PPieceWiseFuncBase.hh +++ b/include/IntegrationTools/piecewise/PPieceWiseFuncBase.hh @@ -2,148 +2,172 @@ #ifndef PPieceWiseFuncBase_HH #define PPieceWiseFuncBase_HH -#include -#include -#include - #include "../pfunction/PFuncBase.hh" -#include "./SimplePiece.hh" #include "./PPieceWiseSimpleBase.hh" #include "./Piece.hh" +#include "./SimplePiece.hh" +#include +#include +#include namespace PRISMS -{ - - /// Class to define a PieceWise Function - /// - /// Contains a vector of 'Piece'. Throws a domain_error if it - /// is evaluated outside of the valid domain of any piece. - /// - template< class VarContainer, class OutType> - class PPieceWiseFuncBase : public PFuncBase - { - public: - typedef typename PFuncBase::size_type size_type; - - size_type _curr_piece; - std::vector > _piece; - - PPieceWiseFuncBase() {} - - PPieceWiseFuncBase( const std::vector > &piece) - { - _piece = piece; - } - - bool in_piece( const VarContainer &var) const - { - for( size_type i=0; i<_piece.size(); i++) - { - if( _piece[i].in_piece(var) ) - return true; - } - return false; - } - - size_type piece(const VarContainer &var) - { - for( size_type i=0; i<_piece.size(); i++) - { - if( _piece[i].in_piece(var)) - return _curr_piece = i; - } - - throw std::domain_error("PPieceWiseFuncBase: Not in any piece"); - } - - virtual PPieceWiseFuncBase *clone() const - { - return new PPieceWiseFuncBase(*this); - } - - virtual PSimpleFunction simplefunction() const - { - std::vector > piece; - - for( size_type i=0; i<_piece.size(); i++) - { - piece.push_back( _piece[i].simplepiece() ); - } - - return PSimpleFunction( PPieceWiseSimpleBase(piece) ); - } - - virtual PSimpleFunction grad_simplefunction(size_type di) const - { - std::vector > piece; - - for( size_type i=0; i<_piece.size(); i++) - { - piece.push_back( _piece[i].grad_simplepiece(di)); - } - - return PSimpleFunction( PPieceWiseSimpleBase(piece)); - } - - virtual PSimpleFunction hess_simplefunction(size_type di, size_type dj) const - { - std::vector > piece; - - for( size_type i=0; i<_piece.size(); i++) - { - piece.push_back( _piece[i].hess_simplepiece(di,dj)); - } - - return PSimpleFunction( PPieceWiseSimpleBase(piece) ); - } +{ - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - - virtual OutType operator()(const VarContainer &var) - { - return _piece[piece(var)](var); - } - virtual OutType grad(const VarContainer &var, size_type di) - { - return _piece[piece(var)].grad(var, di); - } - virtual OutType hess(const VarContainer &var, size_type di, size_type dj) - { - return _piece[piece(var)].hess(var, di, dj); - } + /// Class to define a PieceWise Function + /// + /// Contains a vector of 'Piece'. Throws a domain_error if it + /// is evaluated outside of the valid domain of any piece. + /// + template + class PPieceWiseFuncBase : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - virtual void eval(const VarContainer &var) - { - _piece[piece(var)].eval(var); - } - virtual void eval_grad(const VarContainer &var) + size_type _curr_piece; + std::vector> _piece; + + PPieceWiseFuncBase() + {} + + PPieceWiseFuncBase(const std::vector> &piece) + { + _piece = piece; + } + + bool + in_piece(const VarContainer &var) const + { + for (size_type i = 0; i < _piece.size(); i++) { - _piece[piece(var)].eval_grad(var); + if (_piece[i].in_piece(var)) + return true; } - virtual void eval_hess(const VarContainer &var) + return false; + } + + size_type + piece(const VarContainer &var) + { + for (size_type i = 0; i < _piece.size(); i++) { - _piece[piece(var)].eval_hess(var); + if (_piece[i].in_piece(var)) + return _curr_piece = i; } - - /// These don't recheck the domain - virtual OutType operator()() const + + throw std::domain_error("PPieceWiseFuncBase: Not in any piece"); + } + + virtual PPieceWiseFuncBase * + clone() const + { + return new PPieceWiseFuncBase(*this); + } + + virtual PSimpleFunction + simplefunction() const + { + std::vector> piece; + + for (size_type i = 0; i < _piece.size(); i++) { - return _piece[_curr_piece](); + piece.push_back(_piece[i].simplepiece()); } - virtual OutType grad(size_type di) const + + return PSimpleFunction( + PPieceWiseSimpleBase(piece)); + } + + virtual PSimpleFunction + grad_simplefunction(size_type di) const + { + std::vector> piece; + + for (size_type i = 0; i < _piece.size(); i++) { - return _piece[_curr_piece].grad(di); + piece.push_back(_piece[i].grad_simplepiece(di)); } - virtual OutType hess(size_type di, size_type dj) const + + return PSimpleFunction( + PPieceWiseSimpleBase(piece)); + } + + virtual PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + std::vector> piece; + + for (size_type i = 0; i < _piece.size(); i++) { - return _piece[_curr_piece].hess(di, dj); + piece.push_back(_piece[i].hess_simplepiece(di, dj)); } - }; -} + return PSimpleFunction( + PPieceWiseSimpleBase(piece)); + } + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + + virtual OutType + operator()(const VarContainer &var) + { + return _piece[piece(var)](var); + } + + virtual OutType + grad(const VarContainer &var, size_type di) + { + return _piece[piece(var)].grad(var, di); + } + + virtual OutType + hess(const VarContainer &var, size_type di, size_type dj) + { + return _piece[piece(var)].hess(var, di, dj); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + virtual void + eval(const VarContainer &var) + { + _piece[piece(var)].eval(var); + } + + virtual void + eval_grad(const VarContainer &var) + { + _piece[piece(var)].eval_grad(var); + } + + virtual void + eval_hess(const VarContainer &var) + { + _piece[piece(var)].eval_hess(var); + } + + /// These don't recheck the domain + virtual OutType + operator()() const + { + return _piece[_curr_piece](); + } + + virtual OutType + grad(size_type di) const + { + return _piece[_curr_piece].grad(di); + } + + virtual OutType + hess(size_type di, size_type dj) const + { + return _piece[_curr_piece].hess(di, dj); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/piecewise/PPieceWiseSimpleBase.hh b/include/IntegrationTools/piecewise/PPieceWiseSimpleBase.hh index 9cae8ab8d..6672b930d 100644 --- a/include/IntegrationTools/piecewise/PPieceWiseSimpleBase.hh +++ b/include/IntegrationTools/piecewise/PPieceWiseSimpleBase.hh @@ -2,138 +2,143 @@ #ifndef PPieceWiseSimpleBase_HH #define PPieceWiseSimpleBase_HH -#include -#include -#include - #include "../pfunction/PSimpleBase.hh" #include "./SimplePiece.hh" +#include +#include +#include namespace PRISMS { - /// Class to define a PieceWise SimpleFunction - /// - /// Contains a vector of 'SimplePiece'. Throws a domain_error if it - /// is evaluated outside of the valid domain of any piece. - /// - template< class VarContainer, class OutType> - class PPieceWiseSimpleBase : public PSimpleBase + /// Class to define a PieceWise SimpleFunction + /// + /// Contains a vector of 'SimplePiece'. Throws a domain_error if it + /// is evaluated outside of the valid domain of any piece. + /// + template + class PPieceWiseSimpleBase : public PSimpleBase + { + public: + mutable std::vector> _piece; + + PPieceWiseSimpleBase() + {} + + PPieceWiseSimpleBase(const std::vector> &piece) { - public: - - mutable std::vector > _piece; - - PPieceWiseSimpleBase() {} - - PPieceWiseSimpleBase( const std::vector > &piece) - { - _piece = piece; - } - - virtual std::string csrc() const + _piece = piece; + } + + virtual std::string + csrc() const + { + std::string str = ""; + for (int i = 0; i < _piece.size(); i++) { - std::string str = ""; - for( int i=0; i<_piece.size(); i++) + if (i == 0) { - if( i == 0) - { - str += _piece[i].csrc(); - } - else if( i == _piece.size()-1 ) - { - str += "; and " + _piece[i].csrc(); - } - else - { - str += "; " + _piece[i].csrc(); - } + str += _piece[i].csrc(); } - return str; - } - - virtual std::string sym() const - { - std::string str = ""; - for( int i=0; i<_piece.size(); i++) + else if (i == _piece.size() - 1) { - if( i == 0) - { - str += _piece[i].sym(); - } - else if( i == _piece.size()-1 ) - { - str += "; and " + _piece[i].sym(); - } - else - { - str += "; " + _piece[i].sym(); - } + str += "; and " + _piece[i].csrc(); } - return str; - } - - virtual std::string latex() const - { - std::string str = ""; - for( int i=0; i<_piece.size(); i++) + else { - if( i == 0) - { - str += "\\left\\{ \\begin{array}{ll} " + _piece[i].latex(); - } - else if( i == _piece.size()-1 ) - { - str += " \\\\ " + _piece[i].latex(); - } - else - { - str += " \\\\ " + _piece[i].latex(); - } + str += "; " + _piece[i].csrc(); } - - str += " \\end{array} \\right."; - return str; - } - - virtual PPieceWiseSimpleBase* clone() const - { - return new PPieceWiseSimpleBase(*this); } - - bool in_piece( const VarContainer &var) const + return str; + } + + virtual std::string + sym() const + { + std::string str = ""; + for (int i = 0; i < _piece.size(); i++) { - for( int i=0; i<_piece.size(); i++) + if (i == 0) { - if( _piece[i].in_piece(var) ) - return true; + str += _piece[i].sym(); + } + else if (i == _piece.size() - 1) + { + str += "; and " + _piece[i].sym(); + } + else + { + str += "; " + _piece[i].sym(); } - return false; } - - int piece(const VarContainer &var) + return str; + } + + virtual std::string + latex() const + { + std::string str = ""; + for (int i = 0; i < _piece.size(); i++) { - for( int i=0; i<_piece.size(); i++) + if (i == 0) { - if( _piece[i].in_piece(var)) - return i; + str += "\\left\\{ \\begin{array}{ll} " + _piece[i].latex(); } - - throw std::domain_error("PPieceWiseSimpleBase: Not in any piece"); - } - - private: - virtual OutType eval( const VarContainer &var) const - { - for( int i=0; i<_piece.size(); i++) + else if (i == _piece.size() - 1) { - if( _piece[i].in_piece(var)) - return _piece[i](var); + str += " \\\\ " + _piece[i].latex(); } - - throw std::domain_error("PPieceWiseSimpleBase: Not in any piece"); + else + { + str += " \\\\ " + _piece[i].latex(); + } + } + + str += " \\end{array} \\right."; + return str; + } + + virtual PPieceWiseSimpleBase * + clone() const + { + return new PPieceWiseSimpleBase(*this); + } + + bool + in_piece(const VarContainer &var) const + { + for (int i = 0; i < _piece.size(); i++) + { + if (_piece[i].in_piece(var)) + return true; + } + return false; + } + + int + piece(const VarContainer &var) + { + for (int i = 0; i < _piece.size(); i++) + { + if (_piece[i].in_piece(var)) + return i; + } + + throw std::domain_error("PPieceWiseSimpleBase: Not in any piece"); + } + + private: + virtual OutType + eval(const VarContainer &var) const + { + for (int i = 0; i < _piece.size(); i++) + { + if (_piece[i].in_piece(var)) + return _piece[i](var); } - }; -} + throw std::domain_error("PPieceWiseSimpleBase: Not in any piece"); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/piecewise/Piece.hh b/include/IntegrationTools/piecewise/Piece.hh index 9be936e8e..61ed1e72d 100644 --- a/include/IntegrationTools/piecewise/Piece.hh +++ b/include/IntegrationTools/piecewise/Piece.hh @@ -2,147 +2,177 @@ #ifndef Piece_HH #define Piece_HH -#include -#include - #include "../pfunction/PFunction.hh" #include "./SimplePiece.hh" +#include +#include namespace PRISMS -{ - - /// Class to contain a Function and the piece in which it is valid. - /// - /// This can be evaluated in or out of the piece in which it is declared valid - /// - template< class VarContainer, class OutType> - class Piece : public PFuncBase - { - protected: - - PFunction _expr; - mutable std::vector< PSimpleFunction > _condition; - - typedef typename std::vector< PSimpleFunction >::size_type cond_size_type; - - public: - - typedef typename PFuncBase< VarContainer, OutType>::size_type size_type; - - Piece( const PFunction &expr, const std::vector > &condition) - { - _expr = expr; - _condition = condition; - this->_name = _expr.name(); - this->_var_name = _expr.var_name(); - this->_var_description = _expr.var_description(); - } - - bool in_piece( const VarContainer &var) const - { - for( cond_size_type i=0; i<_condition.size(); i++) - { - if( !_condition[i](var) ) - return false; - } - return true; - } - - PFunction expr() const - { - return _expr; - } - - std::vector< PSimpleFunction > condition() const - { - return _condition; - } - - SimplePiece simplepiece() const - { - return SimplePiece(_expr.simplefunction(), _condition); - } - - SimplePiece grad_simplepiece(size_type di) const - { - return SimplePiece(_expr.grad_simplefunction(di), _condition); - } - - SimplePiece hess_simplepiece(size_type di, size_type dj) const - { - return SimplePiece(_expr.hess_simplefunction(di, dj), _condition); - } - - virtual Piece *clone() const - { - return new Piece(*this); - } - - virtual PSimpleFunction simplefunction() const - { - return PSimpleFunction( SimplePiece(_expr.simplefunction(), _condition)); - } - - virtual PSimpleFunction grad_simplefunction(size_type di) const - { - return PSimpleFunction( SimplePiece(_expr.grad_simplefunction(di), _condition)); - } - - virtual PSimpleFunction hess_simplefunction(size_type di, size_type dj) const - { - return PSimpleFunction( SimplePiece(_expr.hess_simplefunction(di,dj), _condition)); - } +{ - // ---------------------------------------------------------- - // Use these functions if you want to evaluate a single value - - /// These will return '_expr' evaluated anywhere. Must check in_piece first. We - /// don't check it here to avoid double checking when evaluating PPieceWiseFuncBase - /// - virtual OutType operator()(const VarContainer &var) - { - return _expr(var); - } - virtual OutType grad(const VarContainer &var, size_type di) - { - return _expr.grad(var, di); - } - virtual OutType hess(const VarContainer &var, size_type di, size_type dj) - { - return _expr.hess(var, di, dj); - } + /// Class to contain a Function and the piece in which it is valid. + /// + /// This can be evaluated in or out of the piece in which it is declared valid + /// + template + class Piece : public PFuncBase + { + protected: + PFunction _expr; + mutable std::vector> _condition; - // ---------------------------------------------------------- - // Use these functions to evaluate several values, then use 'get' methods to access results - virtual void eval(const VarContainer &var) - { - _expr(var); - } - virtual void eval_grad(const VarContainer &var) - { - _expr.eval_grad(var); - } - virtual void eval_hess(const VarContainer &var) - { - _expr.eval_hess(var); - } - - /// These don't recheck the domain - virtual OutType operator()() const - { - return _expr(); - } - virtual OutType grad(size_type di) const - { - return _expr.grad(di); - } - virtual OutType hess(size_type di, size_type dj) const + typedef + typename std::vector>::size_type cond_size_type; + + public: + typedef typename PFuncBase::size_type size_type; + + Piece(const PFunction &expr, + const std::vector> &condition) + { + _expr = expr; + _condition = condition; + this->_name = _expr.name(); + this->_var_name = _expr.var_name(); + this->_var_description = _expr.var_description(); + } + + bool + in_piece(const VarContainer &var) const + { + for (cond_size_type i = 0; i < _condition.size(); i++) { - return _expr.hess(di, dj); + if (!_condition[i](var)) + return false; } - }; + return true; + } + + PFunction + expr() const + { + return _expr; + } -} + std::vector> + condition() const + { + return _condition; + } + + SimplePiece + simplepiece() const + { + return SimplePiece(_expr.simplefunction(), _condition); + } + + SimplePiece + grad_simplepiece(size_type di) const + { + return SimplePiece(_expr.grad_simplefunction(di), + _condition); + } + + SimplePiece + hess_simplepiece(size_type di, size_type dj) const + { + return SimplePiece(_expr.hess_simplefunction(di, dj), + _condition); + } + + virtual Piece * + clone() const + { + return new Piece(*this); + } + + virtual PSimpleFunction + simplefunction() const + { + return PSimpleFunction( + SimplePiece(_expr.simplefunction(), _condition)); + } + + virtual PSimpleFunction + grad_simplefunction(size_type di) const + { + return PSimpleFunction( + SimplePiece(_expr.grad_simplefunction(di), _condition)); + } + + virtual PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const + { + return PSimpleFunction( + SimplePiece(_expr.hess_simplefunction(di, dj), + _condition)); + } + + // ---------------------------------------------------------- + // Use these functions if you want to evaluate a single value + + /// These will return '_expr' evaluated anywhere. Must check in_piece first. We + /// don't check it here to avoid double checking when evaluating PPieceWiseFuncBase + /// + virtual OutType + operator()(const VarContainer &var) + { + return _expr(var); + } + + virtual OutType + grad(const VarContainer &var, size_type di) + { + return _expr.grad(var, di); + } + + virtual OutType + hess(const VarContainer &var, size_type di, size_type dj) + { + return _expr.hess(var, di, dj); + } + + // ---------------------------------------------------------- + // Use these functions to evaluate several values, then use 'get' methods to access + // results + virtual void + eval(const VarContainer &var) + { + _expr(var); + } + + virtual void + eval_grad(const VarContainer &var) + { + _expr.eval_grad(var); + } + + virtual void + eval_hess(const VarContainer &var) + { + _expr.eval_hess(var); + } + + /// These don't recheck the domain + virtual OutType + operator()() const + { + return _expr(); + } + + virtual OutType + grad(size_type di) const + { + return _expr.grad(di); + } + + virtual OutType + hess(size_type di, size_type dj) const + { + return _expr.hess(di, dj); + } + }; +} // namespace PRISMS #endif \ No newline at end of file diff --git a/include/IntegrationTools/piecewise/SimplePiece.hh b/include/IntegrationTools/piecewise/SimplePiece.hh index 82ce7c3be..9c1a20e2e 100644 --- a/include/IntegrationTools/piecewise/SimplePiece.hh +++ b/include/IntegrationTools/piecewise/SimplePiece.hh @@ -2,135 +2,141 @@ #ifndef SimplePiece_HH #define SimplePiece_HH -#include -#include - #include "../pfunction/PSimpleFunction.hh" +#include +#include namespace PRISMS -{ - /// Class to contain a SimpleFunction and the piece in which it is valid. - /// - /// This can be evaluated in or out of the piece in which it is declared valid - /// - template< class VarContainer, class OutType> - class SimplePiece : public PSimpleBase +{ + /// Class to contain a SimpleFunction and the piece in which it is valid. + /// + /// This can be evaluated in or out of the piece in which it is declared valid + /// + template + class SimplePiece : public PSimpleBase + { + protected: + mutable PSimpleFunction _expr; + mutable std::vector> _condition; + typedef + typename std::vector>::size_type size_type; + + public: + SimplePiece(const PSimpleFunction &expr, + const std::vector> &condition) { - protected: - - mutable PSimpleFunction _expr; - mutable std::vector< PSimpleFunction > _condition; - typedef typename std::vector< PSimpleFunction >::size_type size_type; - - public: - - SimplePiece( const PSimpleFunction &expr, const std::vector< PSimpleFunction > &condition) - { - _expr = expr; - _condition = condition; - this->_name = _expr.name(); - } - - virtual std::string csrc() const + _expr = expr; + _condition = condition; + this->_name = _expr.name(); + } + + virtual std::string + csrc() const + { + std::string str = _expr.csrc(); + for (size_type i = 0; i < _condition.size(); i++) { - std::string str = _expr.csrc(); - for( size_type i=0; i<_condition.size(); i++) + if (i == 0) { - if( i == 0) - { - str += " if " + _condition[i].csrc(); - } - else if( i == _condition.size()-1 ) - { - str += " and " + _condition[i].csrc(); - } - else - { - str += ", " + _condition[i].csrc(); - } + str += " if " + _condition[i].csrc(); } - return str; - } - - virtual std::string sym() const - { - std::string str = _expr.sym(); - for( size_type i=0; i<_condition.size(); i++) + else if (i == _condition.size() - 1) { - if( i == 0) - { - str += " if " + _condition[i].sym(); - } - else if( i == _condition.size()-1 ) - { - str += " and " + _condition[i].sym(); - } - else - { - str += ", " + _condition[i].sym(); - } + str += " and " + _condition[i].csrc(); } - return str; - } - - virtual std::string latex() const - { - std::string str = _expr.latex(); - for( size_type i=0; i<_condition.size(); i++) + else { - if( i == 0) - { - str += " & \\mbox{ if } " + _condition[i].latex(); - } - else if( i == _condition.size()-1 ) - { - str += " \\mbox{ and } " + _condition[i].latex(); - } - else - { - str += " \\mbox{, } " + _condition[i].sym(); - } + str += ", " + _condition[i].csrc(); } - return str; - } - - virtual SimplePiece* clone() const - { - return new SimplePiece(*this); } - - bool in_piece( const VarContainer &var) const + return str; + } + + virtual std::string + sym() const + { + std::string str = _expr.sym(); + for (size_type i = 0; i < _condition.size(); i++) { - for( size_type i=0; i<_condition.size(); i++) + if (i == 0) + { + str += " if " + _condition[i].sym(); + } + else if (i == _condition.size() - 1) + { + str += " and " + _condition[i].sym(); + } + else { - if( !_condition[i](var) ) - return false; + str += ", " + _condition[i].sym(); } - return true; } - - PSimpleFunction expr() const + return str; + } + + virtual std::string + latex() const + { + std::string str = _expr.latex(); + for (size_type i = 0; i < _condition.size(); i++) { - return _expr; + if (i == 0) + { + str += " & \\mbox{ if } " + _condition[i].latex(); + } + else if (i == _condition.size() - 1) + { + str += " \\mbox{ and } " + _condition[i].latex(); + } + else + { + str += " \\mbox{, } " + _condition[i].sym(); + } } - - std::vector< PSimpleFunction > condition() const + return str; + } + + virtual SimplePiece * + clone() const + { + return new SimplePiece(*this); + } + + bool + in_piece(const VarContainer &var) const + { + for (size_type i = 0; i < _condition.size(); i++) { - return _condition; - } - - private: - - /// This will return '_expr' evaluated anywhere. Must check in_piece first. We - /// don't check it here to avoid double checking when evaluating PPieceWiseSimpleBase - /// - virtual OutType eval( const VarContainer &var) const - { - return _expr(var); + if (!_condition[i](var)) + return false; } - }; + return true; + } + + PSimpleFunction + expr() const + { + return _expr; + } -} + std::vector> + condition() const + { + return _condition; + } + + private: + /// This will return '_expr' evaluated anywhere. Must check in_piece first. We + /// don't check it here to avoid double checking when evaluating + /// PPieceWiseSimpleBase + /// + virtual OutType + eval(const VarContainer &var) const + { + return _expr(var); + } + }; +} // namespace PRISMS #endif \ No newline at end of file