diff --git a/common/sensitivity/state.h b/common/sensitivity/state.h index 8b691bed788cf91246d119281d284680ad6336ac..ddf4a1bc1c23bcc2ef859379c96de97b0645fcfc 100644 --- a/common/sensitivity/state.h +++ b/common/sensitivity/state.h @@ -5,6 +5,8 @@ #include "element.h" #include "../enum.h" +#include <array> + /* * SensitivityStateEvaluator */ diff --git a/common/wave/element.h b/common/wave/element.h index 6d138ebf04cd81b939950f5337c2f415b845c7b3..4613bae011df52d134ce5b8559ac57f4f418adea 100644 --- a/common/wave/element.h +++ b/common/wave/element.h @@ -10,18 +10,51 @@ * WaveMultiField */ template<Physic T_Physic> -class WaveMultiField final : public std::vector<WaveField<T_Physic>> +class WaveMultiField { public: - WaveMultiField(unsigned int n,const WaveField<T_Physic>& other) : std::vector<WaveField<T_Physic>>(n,other) {}; + mutable WaveField<T_Physic> *_field; + std::vector< gmshfem::algebra::Vector< std::complex< double > > > _dofValues; + WaveMultiField(unsigned int n,const WaveField<T_Physic>& other) : _field(), _dofValues(n) {}; + + unsigned int size() const + { + return _dofValues.size(); + } + + const WaveField<T_Physic> &operator[](const unsigned int index) const + { + _field->setAllVector(_dofValues[index], gmshfem::dofs::RawOrder::Hash); + return *_field; + } + + void setDofs(const std::vector< gmshfem::dofs::RawDof > &dofs) + { + _field->setAllDofVector(dofs, false); + } + + void setValues(const unsigned int index, gmshfem::algebra::Vector< std::complex< double > > &vector) + { + _dofValues[index] = std::move(vector); + } + + std::vector< gmshfem::algebra::Vector< std::complex< double > > >::const_iterator begin() const + { + return _dofValues.begin(); + } + + std::vector< gmshfem::algebra::Vector< std::complex< double > > >::const_iterator end() const + { + return _dofValues.end(); + } void write(std::string name) const { - unsigned int s = 0; - for (auto it = this->begin(); it != this->end(); it++) - { - it->write(name+"s"+std::to_string(s++)); - } +// unsigned int s = 0; +// for (auto it = this->begin(); it != this->end(); it++) +// { +// it->write(name+"s"+std::to_string(s++)); +// } }; }; diff --git a/common/wave/equation/equation.cpp b/common/wave/equation/equation.cpp index 894c982ce2df39be448fea59a843c8b9b313fc12..ac10edc4441fb05f474cfd9dcab9ad3e3d8a4b35 100644 --- a/common/wave/equation/equation.cpp +++ b/common/wave/equation/equation.cpp @@ -171,31 +171,31 @@ const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::updateWithGr template<Physic T_Physic> void DifferentialEquationInterface<T_Physic>::updateGreen(const ModelStateEvaluator& m, const std::vector<unsigned int>& p_idx) { - if(!_systemIsUpToDate) - { - _formulation.removeSystem(); - _formulation.initSystem(); - - _systemIsFactorized=false; - setLHS(m); - _formulation.pre(); - _formulation.assemble(); - _systemIsUpToDate=true; - _formulation.removeTerms(); - } - - for (auto it = p_idx.begin(); it != p_idx.end(); it++) - { - setGreenRHS(*it); - _formulation.assemble(); - _formulation.removeTerms(); - - _formulation.solve(_systemIsFactorized); - _systemIsFactorized=true; - _formulation.setRHSToZero(); - _g[*it] = _v; - //_v.getVector(_dofg[*it]); - } +// if(!_systemIsUpToDate) +// { +// _formulation.removeSystem(); +// _formulation.initSystem(); +// +// _systemIsFactorized=false; +// setLHS(m); +// _formulation.pre(); +// _formulation.assemble(); +// _systemIsUpToDate=true; +// _formulation.removeTerms(); +// } +// +// for (auto it = p_idx.begin(); it != p_idx.end(); it++) +// { +// setGreenRHS(*it); +// _formulation.assemble(); +// _formulation.removeTerms(); +// +// _formulation.solve(_systemIsFactorized); +// _systemIsFactorized=true; +// _formulation.setRHSToZero(); +// _g[*it] = _v; +// //_v.getVector(_dofg[*it]); +// } } template<Physic T_Physic> diff --git a/common/wave/updater.cpp b/common/wave/updater.cpp index 4853d2f4d3abee395d55c9bed084df824312d1c5..2eaacdb85c5561a6e8a477ea1a941a6c79d89710 100644 --- a/common/wave/updater.cpp +++ b/common/wave/updater.cpp @@ -43,7 +43,18 @@ void WaveUpdater<T_Physic>::update(Type type, const ModelStateEvaluator& m) { throw Exception("Wave updater can only update forward state when no data updater is provided (_du = nullptr)"); } - *it = _equation->update_wave(type,s,_du->get(dataNeedToBeUpToDate,m),m,_ws); + WaveField<T_Physic> &wfield = const_cast< WaveField<T_Physic> & >(_equation->update_wave(type,s,_du->get(dataNeedToBeUpToDate,m),m,_ws)); + if(s == 0) { + _ws._state[type]._field = &wfield; +// std::vector< gmshfem::dofs::RawDof > dofs; +// wfield.getAllDofVector(dofs); +// _ws._state[type].setDofs(dofs); + } + + gmshfem::algebra::Vector< std::complex< double > > vector; + wfield.getAllVector(vector, gmshfem::dofs::RawOrder::Hash); + _ws._state[type].setValues(s, vector); + s++; } _ws._isUpToDate[type] = true; diff --git a/specific/model/innerproduct/sobolev.cpp b/specific/model/innerproduct/sobolev.cpp index 18cb9d13586da136474a23e50f347eabcd7b13f0..e131c292ddd10762838728c0108ccc10d22d67f6 100644 --- a/specific/model/innerproduct/sobolev.cpp +++ b/specific/model/innerproduct/sobolev.cpp @@ -17,8 +17,39 @@ static const std::complex< double > im = std::complex< double >(0., 1.); */ namespace sobolev { - template<class M> - InnerProduct::InnerProduct(const ConfigurationInterface* const config, const M& m_discret, const GmshFem& gmshFem, std::string suffix) : DifferentialInnerProductInterface(config,m_discret,gmshFem,suffix), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _p0(_m_size,1.), _weight(_m_size,1.), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(1.) + InnerProduct::InnerProduct(const ConfigurationInterface* const config, const model::Discretization& m_discret, const GmshFem& gmshFem, std::string suffix) : DifferentialInnerProductInterface(config,m_discret,gmshFem,suffix), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _p0(_m_size,1.), _weight(_m_size,1.), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(1.) + { + double scale = 0.; + gmshFem.userDefinedParameter(scale, "scale"); + gmshFem.userDefinedParameter(scale, "scale"+suffix); + _alpha1 = (scale / 2. / M_PI) * (scale / 2. / M_PI); + + unsigned int diag_preconditioner = 0; + gmshFem.userDefinedParameter(diag_preconditioner, "innerproduct_diag_preconditioner"); + _diag_preconditioner = ((bool) diag_preconditioner); + for (unsigned int c = 0; c < _m_size; c++) + { + gmshFem.userDefinedParameter(_p0[c], "innerproduct_preconditioner_ref"); + gmshFem.userDefinedParameter(_p0[c], "innerproduct_preconditioner_ref"+suffix); + gmshFem.userDefinedParameter(_p0[c], "innerproduct_preconditioner_ref"+std::to_string(c)); + gmshFem.userDefinedParameter(_p0[c], "innerproduct_preconditioner_ref"+std::to_string(c)+suffix); + + gmshFem.userDefinedParameter(_weight[c], "innerproduct_weight"); + gmshFem.userDefinedParameter(_weight[c], "innerproduct_weight"+suffix); + gmshFem.userDefinedParameter(_weight[c], "innerproduct_weight"+std::to_string(c)); + gmshFem.userDefinedParameter(_weight[c], "innerproduct_weight"+std::to_string(c)+suffix); + _weight[c] *= _weight[c]; + msg::print << "inner_product_weight" << c << " = " << _weight[c] << msg::endl; + } + + + gmshFem.userDefinedParameter(_bulkWeight, "innerproduct_weight_blk"); + gmshFem.userDefinedParameter(_bulkWeight, "innerproduct_weight_blk"+suffix); + gmshFem.userDefinedParameter(_boundaryWeight, "innerproduct_weight_bnd"); + gmshFem.userDefinedParameter(_boundaryWeight, "innerproduct_weight_bnd"+suffix); + }; + + InnerProduct::InnerProduct(const ConfigurationInterface* const config, const ModelField& m_discret, const GmshFem& gmshFem, std::string suffix) : DifferentialInnerProductInterface(config,m_discret,gmshFem,suffix), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _p0(_m_size,1.), _weight(_m_size,1.), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(1.) { double scale = 0.; gmshFem.userDefinedParameter(scale, "scale"); @@ -49,8 +80,6 @@ namespace sobolev gmshFem.userDefinedParameter(_boundaryWeight, "innerproduct_weight_bnd"); gmshFem.userDefinedParameter(_boundaryWeight, "innerproduct_weight_bnd"+suffix); }; - template InnerProduct::InnerProduct<model::Discretization>(const ConfigurationInterface* const, const model::Discretization&, const gmshfem::common::GmshFem&, std::string); - template InnerProduct::InnerProduct<ModelField>(const ConfigurationInterface* const, const ModelField&, const gmshfem::common::GmshFem&, std::string); void InnerProduct::link(ModelUpdater* const mu, SensitivityUpdater* const su) { diff --git a/specific/model/innerproduct/sobolev.h b/specific/model/innerproduct/sobolev.h index 9df22be20e550d8a4104f31b5bf7d1e1af29b078..f0a1e93717ae982dfdf7a4a31df20ebed2b1f7a6 100644 --- a/specific/model/innerproduct/sobolev.h +++ b/specific/model/innerproduct/sobolev.h @@ -25,8 +25,8 @@ namespace sobolev double _bulkWeight; double _boundaryWeight; public: - template<class M> - InnerProduct(const ConfigurationInterface* const config, const M& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix=""); + InnerProduct(const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix=""); + InnerProduct(const ConfigurationInterface* const config, const ModelField& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix=""); virtual void link(ModelUpdater* const mu, SensitivityUpdater* const mpu) override; virtual void unlink() override; virtual void modelIsObsolete() override; diff --git a/specific/wave/element.h b/specific/wave/element.h index 041d017513aca59ccdbb17c2820ffcc78090d075..2343e05b420ba034c997d58e10112e6096d98b5d 100644 --- a/specific/wave/element.h +++ b/specific/wave/element.h @@ -18,7 +18,7 @@ class WaveField final /* acoustic */ template<> -class WaveField<Physic::acoustic> final : public gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 > +class WaveField<Physic::acoustic> : public gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 > { public: WaveField<Physic::acoustic>(): gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(){};