Forward interpolation operator

common/configuration.h

@@ -59,6 +59,8 @@ public:
     virtual std::string model_gmodel() const = 0;
     virtual std::string data_gmodel() const = 0;
 
+    virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const = 0;
+
     gmshfem::domain::Domain wave_omega(Support support) const {return _wave_omega[support];};
     gmshfem::domain::Domain model_unknown(Support support) const {return _model_unknown[support];};
     gmshfem::domain::Domain model_known(Support support) const {return _model_known[support];};

common/data/element.cpp

@@ -70,6 +70,8 @@ template<Physic T_Physic> void Data<T_Physic>::write(std::string filename,std::s
 {
     /* msg::print << "Write " + filename << msg::endl; */
     CSVio file(filename,';',OpeningMode::NewFile);
+    if(type=="csv")
+    {
         for (auto s = _value.begin(); s != _value.end(); s++)
         {
             for (auto r = s->begin(); r != s->end(); r++)
@@ -78,6 +80,14 @@ template<Physic T_Physic> void Data<T_Physic>::write(std::string filename,std::s
             }
             file << csv::endl;
         }
+    }
+    else
+    {
+        data::Discretization<T_Physic> d_discret;
+        DataField<T_Physic> d(filename,_config->data_evaldomain(),_config->data_gmodel(),d_discret);
+        d.forward_interpolate(*this);
+        d.write(filename,type,"");
+    }
     file.close();
 }
 
@@ -88,6 +98,16 @@ PointData<T_Physic> Data<T_Physic>::value(unsigned int s,unsigned int r) const
     return _value[s][r];
 }
 
+template<Physic T_Physic>
+PointData<T_Physic> Data<T_Physic>::value(double xs, double xr) const
+{
+    std::array<unsigned int,2> index = _config->data_coordinate_to_index(xs,xr);
+    unsigned int s = index[0];
+    unsigned int r = index[1];
+    isValid(s,r);
+    return _value[s][r];
+}
+
 template<Physic T_Physic>
 void Data<T_Physic>::value(unsigned int s,unsigned int r,const PointData<T_Physic>& value)
 {

common/data/element.h

@@ -24,6 +24,7 @@ public:
     void read(std::string filename);
 
     PointData<T_Physic> value(unsigned int s,unsigned int r) const;
+    PointData<T_Physic> value(double  xs, double xr) const;
     void value(unsigned int s,unsigned int r,const PointData<T_Physic>& value);
     void value(const WaveMultiField<T_Physic>& w);
     unsigned int ns() const {return _value.size();};

input/ZP2018Marmousi2D/data_1x1sides.txt

-#Data
-data0 = VO2009I_1x1sides_data0
-data1 = VO2009I_1x1sides_data1
-data2 = VO2009I_1x1sides_data2
-data3 = VO2009I_1x1sides_data3
-data4 = VO2009I_1x1sides_data4
-data5 = VO2009I_1x1sides_data5

input/ZP2018Marmousi2D/synthetics.txt

@@ -6,4 +6,4 @@ path = ../input/ZP2018Marmousi2D/marmousi
 #Discretization
 h = 0.024
 #Output
-write_fields = 1
+write_fields = 0

specific/configuration/inclusion.cpp

@@ -509,6 +509,11 @@ namespace inclusion
         }
     }
 
+    std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
+    {
+        throw Exception("Data coordinate to index not yet implemented for inclusion.");
+    }
+
     template<Physic T_Physic>
     ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
     {

specific/configuration/inclusion.h

@@ -85,6 +85,9 @@ namespace inclusion
         unsigned int nyr() const {return _nyr;};
 
         virtual void mesh() const;
+
+        virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+
         virtual double area() const {return _L * _H; };
         virtual double array() const {return std::max(_He,_Le);};
         virtual double depth() const {return std::max(_H,_L)/2.;};

specific/configuration/soil.cpp

@@ -415,6 +415,16 @@ namespace soil
         gmodel::setPhysicalName(2, 1, "data_omega");
     }
 
+    std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
+    {
+        std::array<unsigned int,2> index;
+        double x0 = (_L-array())/2.;
+        double step = array()/(ner()-1);
+        index[0] = std::round((xs-x0)/step);
+        index[1] = std::round((xr-x0)/step);
+        return index;
+    }
+
     template<Physic T_Physic>
     ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
     {

specific/configuration/soil.h

@@ -44,6 +44,9 @@ namespace soil
     Configuration(std::string name, const gmshfem::common::GmshFem& gmshFem);
     virtual void mesh() const;
     unsigned int ner() const {return _ner;};
+
+    virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+
     double H() const {return _H;};
     virtual double area() const {return _L * _H;};
     virtual double array() const {return _Ler;};

specific/data/discretization.h

@@ -23,6 +23,8 @@ namespace data
       unsigned int _functionSpaceDegree;
       std::string _integrationType;
   public:
+      Discretization(gmshfem::field::FunctionSpaceTypeForm0 functionSpaceType, unsigned int functionSpaceDegree, std::string integrationType) : _functionSpaceType(functionSpaceType), _functionSpaceDegree(functionSpaceDegree), _integrationType(integrationType) {};
+      Discretization() : Discretization(gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1,1,"") {};
       Discretization(const gmshfem::common::GmshFem& gmshFem);
       gmshfem::field::FunctionSpaceTypeForm0 functionSpaceType() const {return _functionSpaceType;};
       std::string integrationType(unsigned int integrandFieldNumber) const

specific/data/element.cpp

+//GmshFEM
+#include "Exception.h"
+#include "Formulation.h"
+
 //GmshFWI Library
 #include "element.h"
+#include "../../common/data/element.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::equation;
 
 /* acoustic */
 template<>
@@ -13,3 +22,31 @@ PointData<Physic::acoustic> conj(const PointData<Physic::acoustic>& d)
 {
   return std::conj(d);
 }
+
+/*
+* DataField
+*/
+void DataField<Physic::acoustic>::forward_interpolate(const Data<Physic::acoustic>& data)
+{
+    gmshfem::problem::Formulation< std::complex< double > > formulation("init_dof");
+    formulation.galerkin(dof(*this),tf(*this),this->domain(),"None");
+    formulation.pre();
+
+    if(this->getFunctionSpace()->type() != gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1)
+    {
+        throw Exception("Finite element interpolation is possible only for HierarchicalH1 function space type");
+    }
+
+    switch (this->getFunctionSpace()->order())
+    {
+        case 1:
+            for(auto it = this->begin(); it != this->end(); it++)
+            {
+                const dofs::Dof< scalar::Precision< std::complex<double> > > *dof = it->first;
+                it->second = data.value(dof->x(),dof->y());
+            }
+            break;
+        default:
+            throw Exception("Finite element interpolation is possible only for fields of degree 1");
+    }
+}

specific/data/element.h

 #ifndef H_SPECIFIC_DATA_ELEMENT
 #define H_SPECIFIC_DATA_ELEMENT
 
-//GmshFem Library
+//GmshFEM Library
 #include "MathObject.h"
 #include "FieldInterface.h"
 #include "Function.h"
+#include "Post.h"
 
-//FWI Library
+//GmshFWI Library
 #include "discretization.h"
 
 /*
@@ -32,6 +33,10 @@ template<Physic T_Physic>
 PointData<T_Physic> conj(const PointData<T_Physic>& d);
 
 
+//Forward declaration
+template<Physic T_Physic>
+class Data;
+
 /*
 * DataField
 */
@@ -44,8 +49,12 @@ template<>
 class DataField<Physic::acoustic> final : public gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >
 {
 public:
-  DataField<Physic::acoustic>(std::string prefix,const gmshfem::domain::Domain &domain, const data::Discretization<Physic::acoustic>& d_param);
+  DataField<Physic::acoustic>(std::string name,const gmshfem::domain::Domain &domain, std::string gmodel, const data::Discretization<Physic::acoustic>& d_discret) : gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(name, domain, d_discret.functionSpaceType(), d_discret.functionSpaceDegree(), gmodel) {};
   DataField<Physic::acoustic>(const gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >& other) : gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(other) {};
+
+  void write(std::string name, std::string type="pos", std::string path="") const {gmshfem::post::save(*this,this->domain(),name,type,path);}
+
+  void forward_interpolate(const Data<Physic::acoustic>& data);
 };
 
 /*

synthetics.cpp

@@ -46,6 +46,7 @@ int synthetics(const GmshFem& gmshFem)
 
         std::array<bool,4> needUpToDate = {true,false,false,false};
         ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"+suffix);
+        ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"+suffix,"pos");
         
         unsigned int write_fields = 1;
         gmshFem.userDefinedParameter(write_fields, "write_fields");