Avoid field copy

common/sensitivity/state.h

@@ -5,6 +5,8 @@
 #include "element.h"
 #include "../enum.h"
 
+#include <array>
+
 /*
 * SensitivityStateEvaluator
 */

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) {};
     
-    void write(std::string name) const
+    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)
     {
-        unsigned int s = 0;
-        for (auto it = this->begin(); it != this->end(); it++)
+      _dofValues[index] = std::move(vector);
+    }
+    
+    std::vector< gmshfem::algebra::Vector< std::complex< double > > >::const_iterator begin() const
     {
-            it->write(name+"s"+std::to_string(s++));
+      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++));
+//        }
     };
 };
 

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>

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;

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)
     {

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;

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 >(){};