diff --git a/common/sensitivity/state.h b/common/sensitivity/state.h
index ddf4a1bc1c23bcc2ef859379c96de97b0645fcfc..7d0221700dfc0434b5d6e96c72c6ad430476d832 100644
--- a/common/sensitivity/state.h
+++ b/common/sensitivity/state.h
@@ -1,11 +1,13 @@
#ifndef H_COMMON_SENSITIVITY_STATE
#define H_COMMON_SENSITIVITY_STATE
+//Standerd Library
+#include <array>
+
//GmshFWI Library
#include "element.h"
#include "../enum.h"
-#include <array>
/*
* SensitivityStateEvaluator
diff --git a/common/wave/element.h b/common/wave/element.h
index 4613bae011df52d134ce5b8559ac57f4f418adea..e69260daa5b271a2ff1bb71016ac42eaea033d65 100644
--- a/common/wave/element.h
+++ b/common/wave/element.h
@@ -13,48 +13,88 @@ template<Physic T_Physic>
class WaveMultiField
{
public:
- mutable WaveField<T_Physic> *_field;
+ mutable WaveField<T_Physic> _field;
+ bool _isInitialized;
std::vector< gmshfem::algebra::Vector< std::complex< double > > > _dofValues;
- WaveMultiField(unsigned int n,const WaveField<T_Physic>& other) : _field(), _dofValues(n) {};
-
+ 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;
+ };
+
+ bool isInitialized() const {return _isInitialized;};
+
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;
+ _field.setAllVector(_dofValues[index], gmshfem::dofs::RawOrder::Hash);
+ return _field;
}
-
- void setDofs(const std::vector< gmshfem::dofs::RawDof > &dofs)
+
+ /*
+ 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++));
-// }
+ for (unsigned int s = 0; s < size(); s++)
+ {
+ (*this)[s].write(name+"s"+std::to_string(s));
+ }
+ };
+};
+
+/*
+* WaveAuxilaryField
+*/
+template<Physic T_Physic>
+class WaveAuxilaryField : public WaveField<T_Physic>
+{
+private:
+ WaveMultiField<T_Physic>* _wmf;
+ unsigned int _index;
+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];
+ }
+ }
};
};
diff --git a/common/wave/equation/equation.cpp b/common/wave/equation/equation.cpp
index ac10edc4441fb05f474cfd9dcab9ad3e3d8a4b35..dd52426c0d171c2d86c09ea361728501c8f5cda4 100644
--- a/common/wave/equation/equation.cpp
+++ b/common/wave/equation/equation.cpp
@@ -68,22 +68,28 @@ std::complex<double> EquationInterface<T_Physic>::directional(const Sensitivity&
*/
template<Physic T_Physic>
DifferentialEquationInterface<T_Physic>::DifferentialEquationInterface(const ConfigurationInterface* const config, const wave::Discretization<T_Physic>& w_discret, const GmshFem& gmshFem, std::string suffix)
-: EquationInterface<T_Physic>(config, gmshFem, suffix), _w_discret(w_discret), _v("_v",_config->wave_evaldomain(),_config->wave_gmodel(),_w_discret),_g(_config->np(),_v),/* _dofg(_config->np(),gmshfem::algebra::Vector<std::complex<double>>()),*/_geAreUpToDate(false),_grAreUpToDate(false), _formulation("wave_formulation"),_systemIsUpToDate(false),_systemIsFactorized(false)
+: EquationInterface<T_Physic>(config, gmshFem, suffix), _w_discret(w_discret), _v("_v",_config->wave_evaldomain(),_config->wave_gmodel(),_w_discret),_g(_config->np()),/* _dofg(_config->np(),gmshfem::algebra::Vector<std::complex<double>>()),*/_geAreUpToDate(false),_grAreUpToDate(false), _formulation("wave_formulation"),_systemIsUpToDate(false),_systemIsFactorized(false)
{
+ /*
unsigned int useGreen_buffer = 0;
gmshFem.userDefinedParameter(useGreen_buffer,"useGreen");
_useGreen = ((bool) useGreen_buffer);
+ */
}
template<Physic T_Physic>
-const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::update_wave(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& d, const ModelStateEvaluator& m, const WaveStateEvaluator<T_Physic>& w)
+void DifferentialEquationInterface<T_Physic>::update_wave(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws)
{
+ /*
if(_useGreen){return updateWithGreen(type,s,d,m,w);}
else {return updateWithSystem(type,s,d,m,w);}
+ */
+ updateWithSystem(output,type,ds,ms,ws);
}
template<Physic T_Physic>
-const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::updateWithSystem(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& d, const ModelStateEvaluator& m, const WaveStateEvaluator<T_Physic>& w)
+void
+DifferentialEquationInterface<T_Physic>::updateWithSystem(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws)
{
if(!_systemIsUpToDate)
{
@@ -91,22 +97,27 @@ const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::updateWithSy
_formulation.initSystem();
_systemIsFactorized=false;
- setLHS(m);
+ setLHS(ms);
_formulation.pre();
_formulation.assemble();
_systemIsUpToDate=true;
_formulation.removeTerms();
}
- setRHS(type,s,d,m,w);
- _formulation.assemble();
- _formulation.removeTerms();
+ if(!output.isInitialized()){output.setField(_v);};
+ _v.setField(&output);
- _formulation.solve(_systemIsFactorized);
- _systemIsFactorized=true;
- _formulation.setRHSToZero();
+ for (unsigned int s = 0; s < _config->ns(); s++)
+ {
+ _v.setIndex(s);
+ setRHS(type,s,ds,ms,ws);
+ _formulation.assemble();
+ _formulation.removeTerms();
- return _v;
+ _formulation.solve(_systemIsFactorized);
+ _systemIsFactorized=true;
+ _formulation.setRHSToZero();
+ }
}
/*
@@ -136,7 +147,7 @@ multiplyAddAssignementEmitters(algebra::Vector<std::complex<double>>& dofv, cons
}
}
*/
-
+/*
template<Physic T_Physic>
const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::updateWithGreen(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& d, const ModelStateEvaluator& m, const WaveStateEvaluator<T_Physic>& w)
{
@@ -146,28 +157,24 @@ const WaveField<T_Physic>& DifferentialEquationInterface<T_Physic>::updateWithGr
case Type::F:
{
if(!_geAreUpToDate){updateGreenEmitter(m);}
- /*
- algebra::Vector<std::complex<double>> _dofv(_v.size());
- multiplyAddAssignementEmitters(_dofv,_dofg,_config,s);
- _v.setVector(_dofv);
- */
+ //algebra::Vector<std::complex<double>> _dofv(_v.size());
+ //multiplyAddAssignementEmitters(_dofv,_dofg,_config,s);
+ //_v.setVector(_dofv);
return _v;
}
case Type::A:
{
if(!_grAreUpToDate){updateGreenReceiver(m);}
- /*
- algebra::Vector<std::complex<double>> _dofv(_v.size());
- multiplyAddAssignementReceivers(_dofv,_dofg,d.state(Type::A),_config,s);
- _v.setVector(_dofv);
- */
+ //algebra::Vector<std::complex<double>> _dofv(_v.size());
+ //multiplyAddAssignementReceivers(_dofv,_dofg,d.state(Type::A),_config,s);
+ //_v.setVector(_dofv);
return _v;
}
case Type::PF: case Type::PA: default:
return updateWithSystem(type,s,d,m,w);
}
}
-
+*/
template<Physic T_Physic>
void DifferentialEquationInterface<T_Physic>::updateGreen(const ModelStateEvaluator& m, const std::vector<unsigned int>& p_idx)
{
@@ -276,10 +283,11 @@ ParametrizedEquation<T_Physic>::ParametrizedEquation(const ConfigurationInterfac
//Wave related
template<Physic T_Physic>
-const WaveField<T_Physic>& ParametrizedEquation<T_Physic>::update_wave(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws)
+void
+ParametrizedEquation<T_Physic>::update_wave(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws)
{
update_model(ms);
- return _equation->update_wave(type,s,ds,_msc,ws);
+ _equation->update_wave(output,type,ds,_msc,ws);
}
//Model/sensitivity related
diff --git a/common/wave/equation/equation.h b/common/wave/equation/equation.h
index b49264f9b1760154310e670ba6261fb946c6d62a..90ac03c5ad3299808712df818c773a2b7d7756f4 100644
--- a/common/wave/equation/equation.h
+++ b/common/wave/equation/equation.h
@@ -34,7 +34,7 @@ public:
virtual bool compatible(const ParametrizationInterface* const parametrization) const = 0;
//Wave related
- virtual const WaveField<T_Physic>& update_wave(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws) = 0;
+ virtual void update_wave(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws) = 0;
//Model/sensitivity related
virtual Sensitivity update_sensitivity(Order order, Support support, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws) = 0;
@@ -55,9 +55,9 @@ protected:
using EquationInterface<T_Physic>::_config;
const wave::Discretization<T_Physic> _w_discret;
- WaveField<T_Physic> _v;
+ WaveAuxilaryField<T_Physic> _v;
- bool _useGreen;
+ //bool _useGreen;
WaveMultiField<T_Physic> _g;
//std::vector< gmshfem::algebra::Vector<std::complex<double>> > _dofg;
bool _geAreUpToDate;
@@ -70,12 +70,12 @@ public:
DifferentialEquationInterface(const ConfigurationInterface* const config, const wave::Discretization<T_Physic>& w_discret, const gmshfem::common::GmshFem& gmshFem,std::string suffix = "");
//Wave related
- virtual const WaveField<T_Physic>& update_wave(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& d, const ModelStateEvaluator& m, const WaveStateEvaluator<T_Physic>& w);
+ virtual void update_wave(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
bool modelIsObsolete() override;//Returns true, if sensitivity[DIAG] depends on model
protected:
- const WaveField<T_Physic>& updateWithSystem(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
- const WaveField<T_Physic>& updateWithGreen(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
+ void updateWithSystem(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
+ //const WaveField<T_Physic>& updateWithGreen(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
virtual void setLHS(const ModelStateEvaluator& m) = 0;
virtual void setRHS(Type type,unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws) = 0;
@@ -124,7 +124,7 @@ public:
virtual bool compatible(const ParametrizationInterface* const parametrization) const {return false;};
//Wave related
- virtual const WaveField<T_Physic>& update_wave(Type type, unsigned int s, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
+ virtual void update_wave(WaveMultiField<T_Physic>& output, Type type, const DataStateEvaluator<T_Physic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<T_Physic>& ws);
//Model/sensitivity related
void update_model(const ModelStateEvaluator& ms);
diff --git a/common/wave/state.h b/common/wave/state.h
index 7cebb3a5d813d4bc480acd7da8fbcc1aaee885a1..8a9506c93798eb0d13aa49dad4e277c1281c7d42 100644
--- a/common/wave/state.h
+++ b/common/wave/state.h
@@ -37,10 +37,10 @@ private:
public:
WaveState(const ConfigurationInterface* const config):
_state{
- WaveMultiField<T_Physic>(config->ns(),WaveField<T_Physic>()),
- WaveMultiField<T_Physic>(config->ns(),WaveField<T_Physic>()),
- WaveMultiField<T_Physic>(config->ns(),WaveField<T_Physic>()),
- WaveMultiField<T_Physic>(config->ns(),WaveField<T_Physic>())
+ WaveMultiField<T_Physic>(config->ns()),
+ WaveMultiField<T_Physic>(config->ns()),
+ WaveMultiField<T_Physic>(config->ns()),
+ WaveMultiField<T_Physic>(config->ns())
}, _isUpToDate{false,false,false,false} {};
virtual const WaveMultiField<T_Physic>& state(Type type) const;
diff --git a/common/wave/updater.cpp b/common/wave/updater.cpp
index 2eaacdb85c5561a6e8a477ea1a941a6c79d89710..ef3abeeb0cf0cf20bf68c815f0de29be5613be94 100644
--- a/common/wave/updater.cpp
+++ b/common/wave/updater.cpp
@@ -36,27 +36,12 @@ void WaveUpdater<T_Physic>::update(Type type, const ModelStateEvaluator& m)
break;
}
- unsigned int s = 0;
- for (auto it = _ws._state[type].begin(); it != _ws._state[type].end(); it++)
+ if( (_du==nullptr) && (type != Type::F) )
{
- if( (_du==nullptr) && (type != Type::F) )
- {
- throw Exception("Wave updater can only update forward state when no data updater is provided (_du = nullptr)");
- }
- 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++;
+ throw Exception("Wave updater can only update forward state when no data updater is provided (_du = nullptr)");
}
+
+ _equation->update_wave(_ws._state[type],type,_du->get(dataNeedToBeUpToDate,m),m,_ws);
_ws._isUpToDate[type] = true;
}