Auxilary field final

common/wave/element.cpp

+//GmshFWI Library
+#include "element.h"
+
+/*
+* WaveMultiField
+*/
+template<Physic T_Physic>
+void WaveMultiField<T_Physic>::initializeField(const WaveField<T_Physic>& field)
+{
+    _field = field;
+    for (auto it = _dofValues.begin(); it != _dofValues.end(); it++)
+    {
+        it->resize(_field.numberOfDofs());
+    }
+    _isInitialized = true;
+}
+
+template<Physic T_Physic>
+const WaveField<T_Physic>& WaveMultiField<T_Physic>::operator[](const unsigned int index) const
+{
+  _field.setAllVector(_dofValues[index], gmshfem::dofs::RawOrder::Hash);
+  return _field;
+}
+
+template<Physic T_Physic>
+void WaveMultiField<T_Physic>::setValues(const unsigned int index, gmshfem::algebra::Vector< std::complex< double > > &vector)
+{
+  _dofValues[index] = std::move(vector);
+}
+
+template<Physic T_Physic>
+void WaveMultiField<T_Physic>::write(std::string name) const
+{
+    for (unsigned int s = 0; s < size(); s++)
+    {
+        (*this)[s].write(name+"s"+std::to_string(s));
+    }
+}
+
+/*
+* WaveAuxilaryField
+*/
+template<Physic T_Physic>
+void WaveAuxilaryField<T_Physic>::assignValues(const std::vector< std::complex<double> > &values)
+{
+    unsigned int i = 0;
+    for(auto it = _wmf->_field.begin(); it != _wmf->_field.end(); it++)
+    {
+        if(it->first->type() == gmshfem::dofs::Type::Unknown)
+        {
+            _wmf->_dofValues[_index][i++] = values[it->first->numDof() - 1];
+        }
+    }
+}
+
+template class WaveMultiField<Physic::acoustic>;
+template class WaveAuxilaryField<Physic::acoustic>;

common/wave/element.h

@@ -6,70 +6,34 @@
 #include "../../specific/wave/element.h"
 #include "../configuration.h"
 
+//Forward declaration
+template<Physic T_Physic> class WaveAuxilaryField;
+
 /*
 * WaveMultiField
 */
 template<Physic T_Physic>
 class WaveMultiField
 {
-public:
+private:
     mutable WaveField<T_Physic> _field;
     bool _isInitialized;
     std::vector< gmshfem::algebra::Vector< std::complex< double > > > _dofValues;
+public:
     WaveMultiField(unsigned int n) : _field(), _isInitialized(false), _dofValues(n) {};
 
-    void setField(const WaveField<T_Physic>& field)
-    {
-        _field = field;
-        for (auto it = _dofValues.begin(); it != _dofValues.end(); it++)
-        {
-            it->resize(_field.numberOfDofs());
-        }
-        _isInitialized = true;
-    };
-
+    void initializeField(const WaveField<T_Physic>& field);
     bool isInitialized() const {return _isInitialized;};
 
-    unsigned int size() const
-    {
-      return _dofValues.size();
-    }
+    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;
-    }
+    const WaveField<T_Physic> &operator[](const unsigned int index) const;
 
-    /*
-    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);
 
-    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();
-    }
+    void write(std::string name) const;
 
-    std::vector< gmshfem::algebra::Vector< std::complex< double > > >::const_iterator end() const
-    {
-      return _dofValues.end();
-    }
-*/
-    void write(std::string name) const
-    {
-        for (unsigned int s = 0; s < size(); s++)
-        {
-            (*this)[s].write(name+"s"+std::to_string(s));
-        }
-    };
+    friend class WaveAuxilaryField<T_Physic>;
 };
 
 /*
@@ -85,17 +49,8 @@ public:
     WaveAuxilaryField(std::string name, gmshfem::domain::Domain domain, std::string gmodel, const wave::Discretization<T_Physic>& w_discret): WaveField<T_Physic>(name,domain,gmodel,w_discret), _wmf(nullptr), _index(0) {};
     void setField(WaveMultiField<T_Physic>* wmf){_wmf=wmf;};
     void setIndex(unsigned int index){_index=index;};
-    virtual void assignValues(const std::vector< std::complex<double> > &values) override
-    {
-        unsigned int i = 0;
-        for(auto it = _wmf->_field.begin(); it != _wmf->_field.end(); it++)
-        {
-            if(it->first->type() == gmshfem::dofs::Type::Unknown)
-            {
-                _wmf->_dofValues[_index][i++] = values[it->first->numDof() - 1];
-            }
-        }
-    };
+    virtual void assignValues(const std::vector< std::complex<double> > &values) override;
 };
 
+
 #endif // H_COMMON_WAVE_ELEMENT

common/wave/equation/equation.cpp

@@ -104,7 +104,7 @@ DifferentialEquationInterface<T_Physic>::updateWithSystem(WaveMultiField<T_Physi
       _formulation.removeTerms();
     }
 
-    if(!output.isInitialized()){output.setField(_v);};
+    if(!output.isInitialized()){output.initializeField(_v);};
     _v.setField(&output);
 
     for (unsigned int s = 0; s < _config->ns(); s++)
@@ -178,31 +178,33 @@ 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();
+    }
+
+    if(!_g.isInitialized()){_g.initializeField(_v);};
+    _v.setField(&_g);
+
+    for (auto it = p_idx.begin(); it != p_idx.end(); it++)
+    {
+        _v.setIndex(*it);
+        setGreenRHS(*it);
+        _formulation.assemble();
+        _formulation.removeTerms();
+
+        _formulation.solve(_systemIsFactorized);
+        _systemIsFactorized=true;
+        _formulation.setRHSToZero();
+    }
 }
 
 template<Physic T_Physic>