Initial commit (reproducing result of paper 1)

common/wave/state.h

+#ifndef H_COMMON_WAVE_STATE
+#define H_COMMON_WAVE_STATE
+
+//FWI Library
+#include "../type.h"
+#include "element.h"
+#include "../configuration.h"
+
+//Forward declaration
+template<Physic T_Physic>
+class WaveUpdater;
+
+/*
+* WaveStateEvaluator
+*/
+template<Physic T_Physic>
+class WaveStateEvaluator
+{
+public:
+    virtual const WaveMultiField<T_Physic>& state(Type type) const = 0;
+    void write(Type type, std::string name) const {state(type).write(name);};
+};
+
+/*
+* WaveState
+*/
+template<Physic T_Physic>
+class WaveState final: public WaveStateEvaluator<T_Physic>
+{
+private:
+    std::array<WaveMultiField<T_Physic>,4> _state;
+    std::array<bool,4> _isUpToDate;
+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>()),
+    }, _isUpToDate{false,false,false,false} {};
+
+    virtual const WaveMultiField<T_Physic>& state(Type type) const;
+
+    friend class WaveUpdater<T_Physic>;
+};
+
+#endif // H_COMMON_WAVE_STATE

common/wave/updater.cpp

+//Standard library
+
+//GmshFem library
+#include "Exception.h"
+
+//GmshFWI Library
+#include "updater.h"
+#include "../data/updater.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+
+
+template<Physic T_Physic>
+void WaveUpdater<T_Physic>::update(Type type, const ModelStateEvaluator& m)
+{
+    if(_ws._isUpToDate[type]){return;}
+
+    std::array<bool,4> dataNeedToBeUpToDate = {false,false,false,false};
+    switch (type)
+    {
+        case Type::F:
+            break;
+        case Type::A:
+            dataNeedToBeUpToDate[Type::A] = true;
+            break;
+        case Type::PF:
+            update(Type::F,m);
+            break;
+        case Type::PA:
+            update(Type::A,m);
+            dataNeedToBeUpToDate[Type::PA] = true;
+            break;
+    }
+
+    unsigned int s = 0;
+    for (auto it = _ws._state[type].begin(); it != _ws._state[type].end(); it++)
+    {
+        *it = _equation->update_wave(type,s,_du->get(dataNeedToBeUpToDate,m),m,_ws);
+        s++;
+    }
+    _ws._isUpToDate[type] = true;
+}
+
+template<Physic T_Physic>
+const WaveStateEvaluator<T_Physic>& WaveUpdater<T_Physic>::get(std::array<bool,4> needToBeUpToDate, const ModelStateEvaluator& m)
+{
+    for (unsigned int t = 0; t < 4; t++)
+    {
+        Type type = ((Type) t);
+        if(needToBeUpToDate[type]){update(type,m);}
+    }
+    return _ws;
+}
+
+template<Physic T_Physic>
+void WaveUpdater<T_Physic>::isObsolete(std::array<bool,4> NoMoreUpToDate)
+{
+    for (unsigned int t = 0; t < 4; t++)
+    {
+        if(NoMoreUpToDate[t]){_ws._isUpToDate[t] = false;}
+    }
+}
+
+/*
+template<Physic T_Physic>
+ModelFunction WaveUpdater<T_Physic>::update_sensitivity(Order order, Support support, const ModelStateEvaluator& ms)
+{
+    switch (order)
+    {
+        case Order::SCD:
+            update(Type::PF,ms);
+            update(Type::PA,ms);
+            //no break;
+        case Order::FST:
+            update(Type::F,ms);
+            update(Type::A,ms);
+    }
+    return _equation->update_sensitivity(order, support, ms, _ws);
+}
+
+template<Physic T_Physic>
+ModelFunction WaveUpdater<T_Physic>::update_preconditioner(const ModelStateEvaluator& ms)
+{
+    return _equation->update_preconditioner(ms);
+}
+*/
+
+template class WaveUpdater<Physic::acoustic>;

common/wave/updater.h

+#ifndef H_COMMON_WAVE_UPDATER
+#define H_COMMON_WAVE_UPDATER
+
+//GmshFWI Library
+#include "equation/equation.h"
+#include "element.h"
+//#include "../configuration.h"
+
+//Forward declaration
+template<Physic T_Physic>
+class DataUpdater;
+
+/*
+* WaveUpdater
+*/
+template<Physic T_Physic>
+class WaveUpdater final
+{
+private:
+    WaveState<T_Physic> _ws;
+    DataUpdater<T_Physic>* const _du;
+    EquationInterface<T_Physic>* const _equation;
+
+    void update(Type type, const ModelStateEvaluator& ms);
+public:
+    WaveUpdater(const ConfigurationInterface* const config,DataUpdater<T_Physic>* du, EquationInterface<T_Physic>* equation) : _ws(config), _du(du), _equation(equation) {};
+
+    const WaveStateEvaluator<T_Physic>& get(std::array<bool,4> needToBeUpToDate, const ModelStateEvaluator& ms);
+
+    void isObsolete(std::array<bool,4> NoMoreUpToDate);
+};
+
+#endif // H_COMMON_WAVE_ELEMENT

convergence.cpp

+//GmshFEM Library
+#include "GmshFem.h"
+#include "Post.h"
+#include "CSVio.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::post;
+using namespace gmshfem::function;
+
+//FWI Library
+#include "specific/physic.h"
+#include "specific/configuration/newConfiguration.h"
+#include "specific/wave/equation/newEquation.h"
+#include "specific/data/objective/newObjective.h"
+#include "common/statefunctional.h"
+
+static const std::complex< double > im = std::complex< double >(0., 1.);
+
+enum class Interval: unsigned int { Linear=0, Log=1 };
+Interval to_interval(const GmshFem& gmshFem)
+{
+  std::string interval_buffer;
+  if(gmshFem.userDefinedParameter(interval_buffer, "interval"))
+  {
+    if(interval_buffer=="linear"){return Interval::Linear;}
+    else if(interval_buffer=="log"){return Interval::Log;}
+  }
+  return Interval::Linear;
+}
+
+template <Physic T_Physic>
+int convergence(const GmshFem& gmshFem)
+{
+    std::string name = "noname";
+    gmshFem.userDefinedParameter(name, "name");
+
+    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+
+    wave::Discretization<T_Physic> w_discret(gmshFem);
+
+    double frequency = 0.;
+    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    {
+        throw Exception("Frequency could not be found.");
+    }
+
+    EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
+
+    std::string filename = "noname";
+    gmshFem.userDefinedParameter(filename, "data");
+    Data<T_Physic> d0(filename,configuration);
+
+    ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);
+
+    StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,equation,objective);
+    FunctionalInterface* const functional = statefunctional;
+
+    double Redm,Imdm;
+    double eps0,epsN;
+    unsigned int N;
+    if(!(
+        gmshFem.userDefinedParameter(Redm, "Re(dm)") &&
+        gmshFem.userDefinedParameter(Imdm, "Im(dm)") &&
+        gmshFem.userDefinedParameter(eps0, "eps0") &&
+        gmshFem.userDefinedParameter(epsN, "epsN") &&
+        gmshFem.userDefinedParameter(N, "N")
+    ) )
+    {
+        throw common::Exception("A directional parameter could not be found.");
+    }
+
+    Interval interval = to_interval(gmshFem);
+
+    ScalarPiecewiseFunction< std::complex< double > > dm;
+    dm.addFunction(Redm+im*Imdm,configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
+
+    functional->setModelPerturbation(dm);
+
+    CSVio output(name+"_convergence", ';', common::OpeningMode::NewFile);
+
+    ModelFunction m = configuration->m0();
+
+    double j = functional->performance( m );
+    msg::print << "j: "<< j << msg::endl;
+    
+    double dj_v = statefunctional->directional(Order::FST, Support::BLK, dm);
+    msg::print << "dj_v: "<< dj_v << msg::endl;
+
+    double dj_s = statefunctional->directional(Order::FST, Support::BND, dm);
+    msg::print << "dj_s: "<< dj_s << msg::endl;
+    msg::print << "jd: "<< dj_s+dj_v << msg::endl;
+
+    double d2j_s = 0.;
+    msg::print << "d2j_s: "<< d2j_s << msg::endl;
+    double d2j_v = 0.;
+    msg::print << "d2j_v: "<< d2j_v << msg::endl;
+    output <<j << dj_s << dj_v << dj_s+dj_v << d2j_s << d2j_v << d2j_s+d2j_v << csv::endl;
+
+    double base =  epsN / eps0;
+    double span =  epsN - eps0;
+    for (unsigned int n = 0; n <= N; n++)
+    {
+        double eps;
+        switch (interval)
+        {
+          case Interval::Log:
+            eps  = eps0 * std::pow(base, ((double) n) / ((double) N));
+            break;
+          case Interval::Linear:
+            eps = eps0 + ((double)n) / ((double)N) * span;
+            break;
+        }
+        msg::print << "--- Directional derivative validation convergence point "<< n << ": eps = " << eps << " ---" << msg::endl;
+
+        double jp =  functional->performance(m + eps*dm);
+        msg::print << "jp: "<< jp << msg::endl;
+        double jm =  functional->performance(m - eps*dm);
+        msg::print << "jm: "<< jm << msg::endl;
+        msg::print << "jFD: "<< (jp - jm) / 2. / eps << msg::endl;
+
+        //msg::print << "j2FD: "<< (jp - 2.*j + jm) / eps /eps << msg::endl;
+
+        output << eps << (jp - jm) / 2. / eps << std::abs(dj_s+dj_v - (jp - jm) / 2. / eps) << std::abs(dj_v - (jp - jm) / 2. / eps) << (jp - 2.*j + jm) / eps / eps  << std::abs(d2j_s+d2j_v - (jp - 2.*j + jm) / eps /eps) << csv::endl;
+        msg::print << msg::endl;
+        msg::print << msg::endl;
+    }
+    output.close();
+
+    delete statefunctional;
+    delete configuration;
+    delete equation;
+    delete objective;
+
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+  Physic T_Physic = to_physic(gmshFem);
+  switch (T_Physic)
+  {
+    case Physic::acoustic: default:
+      return convergence<Physic::acoustic>(gmshFem);
+  }
+}

convergence_plot.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 19 14:22:17 2019
+
+@author: xavier
+"""
+import numpy as np
+from enum import Enum
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+#H = 4.0;
+#L = 5.0;
+#Value for 2 width figure
+H = 3.0;
+L = 3.20;
+
+class Config:
+    weakly_inclusion=1
+    weakly_background=2
+    highly_inclusion=3
+    highly_background=4
+
+config = Config.weakly_background;
+path = '../Paper1/data/convergence/'
+filename = {
+    Config.weakly_inclusion: 'weakly_inclusion',
+    Config.weakly_background: 'weakly_background',
+    Config.highly_inclusion: 'highly_inclusion',
+    Config.highly_background: 'highly_background',
+}[config]
+
+ref = np.loadtxt(path+filename+'_convergence.csv',delimiter=';',skiprows=0,max_rows=1);
+x = np.loadtxt(path+filename+'_convergence.csv',delimiter=';',skiprows=1);
+
+
+#Cost function
+f0 = plt.figure(figsize=(L,H),tight_layout=False);
+plt.subplots_adjust(top=0.93,right=0.95,bottom=0.2,left=0.25)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ax = plt.gca().invert_xaxis()
+plt.ylim(1e-9,2.)
+plt.loglog(x[:,0],x[:,2],color='k',marker='+',markersize=5,linewidth=0)
+plt.hlines(ref[1],x[0,0],x[-1,2],colors='tab:blue', linestyles='dashed')
+plt.xlabel(r'$\epsilon$')
+plt.grid()
+plt.show()
+
+f1 = plt.figure(figsize=(L,H),tight_layout=False);
+plt.subplots_adjust(top=0.93,right=0.95,bottom=0.2,left=0.25)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ax = plt.gca().invert_xaxis()
+plt.ylim(1e-9,2.)
+plt.loglog(x[:,0],x[:,3],color='k',marker='+',markersize=5,linewidth=0)
+plt.hlines(ref[1],x[0,0],x[-1,0],colors='tab:blue', linestyles='dashed')
+plt.xlabel(r'$\epsilon$')
+plt.grid()
+plt.show()

directional.cpp

+//GmshFem Library
+#include "GmshFem.h"
+#include "Post.h"
+#include "CSVio.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::post;
+using namespace gmshfem::function;
+
+//FWI Library
+#include "specific/physic.h"
+#include "specific/configuration/newConfiguration.h"
+#include "specific/wave/equation/newEquation.h"
+#include "specific/data/objective/newObjective.h"
+#include "common/statefunctional.h"
+
+static const std::complex< double > im = std::complex< double >(0., 1.);
+
+enum class Interval: unsigned int { Linear=0, Log=1 };
+Interval to_interval(const GmshFem& gmshFem)
+{
+  std::string interval_buffer;
+  if(gmshFem.userDefinedParameter(interval_buffer, "interval"))
+  {
+    if(interval_buffer=="linear"){return Interval::Linear;}
+    else if(interval_buffer=="log"){return Interval::Log;}
+  }
+  return Interval::Linear;
+}
+
+template <Physic T_Physic>
+int directional(const GmshFem& gmshFem)
+{
+    std::string name = "noname";
+    gmshFem.userDefinedParameter(name, "name");
+
+    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+
+    wave::Discretization<T_Physic> w_discret(gmshFem);
+
+    double frequency = 0.;
+    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    {
+        throw Exception("Frequency could not be found.");
+    }
+
+    EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
+
+    std::string filename = "noname";
+    gmshFem.userDefinedParameter(filename, "data");
+    Data<T_Physic> d0(filename,configuration);
+
+    ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);
+
+    StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,equation,objective);
+    FunctionalInterface* const functional = statefunctional;
+
+    double Remu0,Immu0,RemuN,ImmuN,Redm,Imdm,eps;
+    unsigned int N;
+    if(!(
+        gmshFem.userDefinedParameter(Remu0, "Re(mu)") &&
+        gmshFem.userDefinedParameter(Immu0, "Im(mu)") &&
+        gmshFem.userDefinedParameter(RemuN, "Re(muN)") &&
+        gmshFem.userDefinedParameter(ImmuN, "Im(muN)") &&
+        gmshFem.userDefinedParameter(Redm, "Re(dm)") &&
+        gmshFem.userDefinedParameter(Imdm, "Im(dm)") &&
+        gmshFem.userDefinedParameter(N, "N") &&
+        gmshFem.userDefinedParameter(eps, "eps")
+        ) )
+    {
+        throw Exception("A directional parameter could not be found.");
+    }
+
+
+    Interval interval = to_interval(gmshFem);
+
+    std::complex<double> mu0 = Remu0 + im * Immu0;
+    std::complex<double> muN = RemuN + im * ImmuN;
+    double span = std::abs(muN - mu0);
+    std::complex<double> dir = (muN - mu0) / span;
+
+    ScalarPiecewiseFunction< std::complex< double > > dm;
+    dm.addFunction(Redm+im*Imdm,configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
+
+    functional->setModelPerturbation(dm);
+
+    CSVio output(name+"_directional", ';', common::OpeningMode::NewFile);
+    for (unsigned int n = 0; n <= N; n++)
+    {
+        msg::print << "--- Directional derivatives validation point "<< n << " ---" << msg::endl;
+        double step=0.;
+        switch (interval)
+        {
+          case Interval::Log:
+            step = std::pow(span, ((double)n) / ((double)N) ) - 1.;
+            break;
+          case Interval::Linear:
+            step = ((double)n) / ((double)N) * span;
+            break;
+        }
+        std::complex<double> mn = mu0 + step*dir;
+        msg::print << "mn = " << mn << msg::endl;
+        ScalarPiecewiseFunction<std::complex<double>> m;
+        m.addFunction(mn,configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
+
+
+        double j =  functional->performance(m);
+        msg::print << "j: "<< j << msg::endl;
+
+        double dj_s = statefunctional->directional(Order::FST, Support::BND, dm);
+        double dj_v = statefunctional->directional(Order::FST, Support::BLK, dm);
+
+        /*
+        double jfd = 0.;
+        double d2j_s = 0.;
+        double d2j_v = 0.;
+        double j2fd = 0.;
+        */
+
+        double d2j_s = statefunctional->directional(Order::SCD, Support::BND, dm);
+        double d2j_v = statefunctional->directional(Order::SCD, Support::BLK, dm);
+
+        double jp =  functional->performance(m + eps*dm);
+        //msg::print << "jp: "<< jp << msg::endl;
+        double jm =  functional->performance(m - eps*dm);
+        //msg::print << "jm: "<< jm << msg::endl;
+        double jfd = (jp - jm) / 2. / eps;
+        msg::print << "dj_s: "<< dj_s << msg::endl;
+        msg::print << "dj_v: "<< dj_v << msg::endl;
+        msg::print << "jd: "<< dj_s+dj_v << msg::endl;
+        msg::print << "jFD: "<< jfd << msg::endl;
+
+        double jp1 =  functional->performance(m + 2.*eps*dm);
+        //msg::print << "jp1: "<< jp1 << msg::endl;
+        double jm1 =  functional->performance(m - 2.*eps*dm);
+        //msg::print << "jm1: "<< jm1 << msg::endl;
+
+        double j2fd = (jm1 -2.*jm +2.*j -2.*jp + jp1) / eps / eps / 2.;
+        msg::print << "d2j_s: "<< d2j_s << msg::endl;
+        msg::print << "d2j_v: "<< d2j_v << msg::endl;
+        msg::print << "d2j: "<< d2j_s+d2j_v << msg::endl;
+        msg::print << "j2FD: "<< j2fd << msg::endl;
+
+        /*
+        */
+
+        output << mn  << j << dj_s << dj_v << dj_s+dj_v << jfd << dj_s+dj_v - jfd  << d2j_s << d2j_v << d2j_s+d2j_v << j2fd  << d2j_s+d2j_v - j2fd << csv::endl;
+        msg::print << msg::endl;
+        msg::print << msg::endl;
+    }
+    output.close();
+
+    delete statefunctional;
+    delete configuration;
+    delete equation;
+    delete objective;
+
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+  Physic T_Physic = to_physic(gmshFem);
+  switch (T_Physic)
+  {
+    case Physic::acoustic: default:
+      return directional<Physic::acoustic>(gmshFem);
+  }
+}

directional_plot.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 19 14:22:17 2019
+
+@author: xavier
+"""
+import numpy as np
+from enum import Enum
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+#H = 4.0;
+#L = 5.0;
+#Value for 2 width figure
+H = 3.0;
+L = 3.20;
+
+class Config:
+    weakly_inclusion=1
+    weakly_background=2
+    highly_inclusion=3
+    highly_background=4
+
+config = Config.weakly_background;
+path = '../Paper1/data/directional/'
+filename = {
+    Config.weakly_inclusion: 'weakly_inclusion',
+    Config.weakly_background: 'weakly_background',
+    Config.highly_inclusion: 'highly_inclusion',
+    Config.highly_background: 'highly_background',
+}[config]
+
+xc = {
+    Config.weakly_inclusion: 1,
+    Config.weakly_background: 0,
+    Config.highly_inclusion: 0,
+    Config.highly_background: 0,
+}[config]
+
+ax = {
+    Config.weakly_inclusion: -1,
+    Config.weakly_background: 1,
+    Config.highly_inclusion: 1,
+    Config.highly_background: 1,
+}[config]
+
+var = {
+    Config.weakly_inclusion: '\\alpha_c',
+    Config.weakly_background: 's^2_0',
+    Config.highly_inclusion: 's^2_c',
+    Config.highly_background: 's^2_0',
+}[config]
+
+x = np.loadtxt(path+filename+'_directional.csv',delimiter=';');
+
+
+#Cost function
+f0 = plt.figure(figsize=(L,H),tight_layout=False);
+plt.subplots_adjust(top=0.93,right=0.95,bottom=0.2,left=0.25)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+plt.plot(ax*x[:,xc],x[:,2],color='k',marker='o',markersize=3,linewidth=0)
+if config == Config.weakly_inclusion:
+    plt.xscale('log')
+plt.xlabel(r'$'+var+'$')
+plt.ylabel(r'$\mathcal{J}('+var+')$')
+plt.grid()
+plt.show()
+
+#First derivative
+f1 = plt.figure(figsize=(L,H))
+plt.subplots_adjust(top=0.93,right=0.95,bottom=0.2,left=0.25)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+#Finite difference
+plt.plot(ax*x[:,xc],x[:,6],color='tab:green',marker='x',markersize=5,linewidth=0,label='FD')
+
+
+#Adjoint state method
+if config == Config.highly_inclusion:
+    #Volume cylinder
+    plt.plot(ax*x[:,xc],x[:,5],color='tab:brown',marker='o',markersize=3,linewidth=0,label='AS')
+elif config == Config.weakly_inclusion:
+    #Surface cylinder
+    plt.plot(ax*x[:,xc],x[:,5],color='tab:brown',marker='o',markersize=3,linewidth=0,label='AS',fillstyle='none')
+    plt.xscale('log')
+elif config == Config.highly_background or config == Config.weakly_background:
+    #Volume background
+    plt.plot(ax*x[:,xc],x[:,4],color='tab:blue',marker='o',markersize=3,linewidth=0,label='ASv')
+    #Surface background
+    plt.plot(ax*x[:,xc],x[:,3],color='tab:blue',marker='o',fillstyle='none',markersize=3,linewidth=0,label='ASs')
+
+plt.grid()
+#plt.legend(loc=3,ncol=1);
+plt.xlabel(r'$'+var+'$')
+plt.ylabel(r'$D_{'+var+'} \mathcal{J}('+var+')$')
+plt.show();
+
+#Second derivative
+f2 = plt.figure(figsize=(L,H));
+plt.subplots_adjust(top=0.93,right=0.95,bottom=0.2,left=0.25)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+#Finite difference
+plt.plot(ax*x[:,xc],x[:,11],color='tab:green',marker='x',markersize=5,linewidth=0,label='FD')
+
+#Adjoint state method
+if config == Config.highly_inclusion:
+    #Volume cylinder
+    plt.plot(ax*x[:,xc],x[:,9],color='tab:brown',marker=r'$\bullet$',markersize=3,linewidth=0,label='AS')
+elif config == Config.weakly_inclusion:
+    #Surface cylinder
+    plt.plot(ax*x[:,xc],x[:,8],color='tab:brown',marker=r'$\circ$',markersize=3,linewidth=0,label='AS')
+    plt.xscale('log')
+elif config == Config.highly_background or config == Config.weakly_background:
+    #Volume background
+    plt.plot(ax*x[:,xc],x[:,9],color='tab:blue',marker=r'$\blacksquare$',markersize=3,linewidth=0,label='ASv')
+    #Surface background
+    plt.plot(ax*x[:,xc],x[:,8],color='tab:blue',marker=r'$\square$',markersize=3,linewidth=0,label='ASs')
+
+plt.grid()
+#plt.legend(loc=3,ncol=1);
+plt.xlabel(r'$'+var+'$')
+plt.ylabel(r'$D^2_{'+var+'} \mathcal{J}('+var+')$')
+plt.show();

gradient.cpp

+//GmshFem Library
+#include "GmshFem.h"
+//#include "Function.h"
+//#include "Post.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+//using namespace gmshfem::post;
+//using namespace gmshfem::function;
+
+//FWI Library
+#include "specific/physic.h"
+#include "specific/configuration/newConfiguration.h"
+#include "specific/wave/equation/newEquation.h"
+#include "specific/data/objective/newObjective.h"
+#include "specific/model/innerproduct/newInnerProduct.h"
+#include "common/statefunctional.h"
+
+static const std::complex< double > im = std::complex< double >(0., 1.);
+
+enum class Interval: unsigned int { Linear=0, Log=1 };
+Interval to_interval(const GmshFem& gmshFem)
+{
+  std::string interval_buffer;
+  if(gmshFem.userDefinedParameter(interval_buffer, "interval"))
+  {
+    if(interval_buffer=="linear"){return Interval::Linear;}
+    else if(interval_buffer=="log"){return Interval::Log;}
+  }
+  return Interval::Linear;
+}
+
+template <Physic T_Physic>
+int gradient(const GmshFem& gmshFem)
+{
+    std::string name = "noname";
+    gmshFem.userDefinedParameter(name, "name");
+
+    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+
+
+    double frequency = 0.;
+    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    {
+        throw Exception("Frequency could not be found.");
+    }
+
+    wave::Discretization<T_Physic> w_discret(gmshFem);
+    EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
+
+    std::string filename = "noname";
+    gmshFem.userDefinedParameter(filename, "data");
+    Data<T_Physic> d0(filename,configuration);
+
+    ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);
+
+    model::Discretization m_discret(gmshFem);
+    sobolev::InnerProduct innerproduct(configuration,m_discret,gmshFem);
+
+    StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,&innerproduct,equation,objective);
+    FunctionalInterface* const functional = statefunctional;
+    functional->setModel(configuration->m0());
+
+    double scale0,scaleN;
+    unsigned int N;
+    if(!(
+    gmshFem.userDefinedParameter(scale0, "scale0") &&
+    gmshFem.userDefinedParameter(scaleN, "scaleN") &&
+    gmshFem.userDefinedParameter(N, "N")
+    ) )
+    {
+        throw common::Exception("A directional parameter could not be found.");
+    }
+    Interval interval = to_interval(gmshFem);
+
+    double scale = 0.;
+
+    double base =  scaleN / scale0;
+    double span =  scaleN - scale0;
+    for (unsigned int n = 0; n <= N; n++)
+    {
+        if(N==0){scale = scale0;}
+        else
+        {
+            switch (interval)
+            {
+                case Interval::Log:
+                scale  = scale0 * std::pow(base, ((double) n) / ((double) N));
+                break;
+                case Interval::Linear:
+                scale = scale0 + ((double)n) / ((double)N) * span;
+                break;
+            }
+        }
+        msg::print << "--- Gradient analysis # "<< n << ": scale = "<< scale << " ---" << msg::endl;
+
+        innerproduct.alpha1( (scale / (2 * M_PI)) * (scale / (2 * M_PI)) );
+        innerproduct.bulkWeight(1.0);
+        innerproduct.boundaryWeight(0.0);
+        statefunctional->innerproductIsObsolete();
+        functional->gradient().write(name+"_blk_gradientn"+std::to_string(n));
+
+        innerproduct.bulkWeight(0.0);
+        innerproduct.boundaryWeight(1.0);
+        statefunctional->innerproductIsObsolete();
+        functional->gradient().write(name+"_bnd_gradientn"+std::to_string(n));
+    }
+
+    delete statefunctional;
+    delete configuration;
+    delete equation;
+    delete objective;
+
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+  Physic T_Physic = to_physic(gmshFem);
+  switch (T_Physic)
+  {
+    case Physic::acoustic: default:
+      return gradient<Physic::acoustic>(gmshFem);
+  }
+}

input/paper1/common.txt

+#Equation
+physic = acoustic
+equation = helmholtz
+useGreen = 1
+#
+#Configuration
+configuration = inclusion
+#
+ye = 0.0
+yr = 0.0
+nxe = 11
+nxr = 11
+nye = 0
+nyr = 0
+He = 15.0
+Hr = 15.0
+Le = 0.0
+Lr = 0.0
+L = 25.0
+H = 25.0
+#
+#Discretization
+wave_FunctionSpaceType = HierarchicalH1
+wave_IntegrationType = Gauss
+wave_FunctionSpaceDegree = 5
+#
+#Frequency
+frequency = 1.

input/paper1/convergence_highly_background.txt

+name = highly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cylinder
+filled = 1
+r = 2.5
+#
+unknown = background
+Re(mb) = 0.93
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = 0.
+#
+#Directional
+Re(dm) = 1.
+Im(dm) = 0.
+interval = log
+eps0 = 1e-10
+epsN = 1e-1
+N = 9
+#
+#Data
+data = filled_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/convergence_weakly_background.txt

+name = weakly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 5.0
+inclusion = cylinder
+filled = 0
+r = 2.5
+#
+unknown = background
+Re(mb) = 0.93
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = -100.
+#
+#Directional
+Re(dm) = 1.
+Im(dm) = 0.
+interval = log
+eps0 = 1e-10
+epsN = 1e-1
+N = 9
+#
+#Data
+data = empty_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/directional_highly_background.txt

+name = highly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cylinder
+filled = 1
+r = 2.5
+#
+unknown = background
+Re(mb) = 0.
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = 0.
+#
+#Directional
+Re(mu) = 0.9
+Im(mu) = 0.
+Re(muN) = 1.1
+Im(muN) = 0.
+Re(dm) = 1.
+Im(dm) = 0.
+N = 20
+eps = 1e-5
+#
+#Data
+data = filled_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/directional_highly_inclusion.txt

+name = highly_inclusion
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cylinder
+filled = 1
+r = 2.5
+#
+unknown = inclusion
+Re(mb) = 1.
+Im(mb) = 0.
+Re(mc) = 0.
+Im(mc) = 0.
+#
+#Directional
+Re(mu) = 1.0
+Im(mu) = 0.
+Re(muN) = 1.5
+Im(muN) = 0.
+Re(dm) = 1.
+Im(dm) = 0.
+N = 20
+eps = 1e-5
+#
+#Data
+data = filled_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/directional_weakly_background.txt

+name = weakly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 5.0
+inclusion = cylinder
+filled = 0
+r = 2.5
+#
+unknown = background
+Re(mb) = 0.
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = -100.
+#
+#Directional
+Re(mu) = 0.9
+Im(mu) = 0.
+Re(muN) = 1.1
+Im(muN) = 0.
+Re(dm) = 1.
+Im(dm) = 0.
+N = 20
+eps = 1e-5
+#
+#Data
+data = empty_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/directional_weakly_inclusion.txt

+name = weakly_inclusion
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 5.0
+inclusion = cylinder
+filled = 0
+r = 2.5
+#
+unknown = inclusion
+Re(mb) = 1.
+Im(mb) = 0.
+Re(mc) = 0.
+Im(mc) = 0.
+#
+#Directional
+Re(mu) = 1.2
+Im(mu) = -10.
+Re(muN) = 1.2
+Im(muN) = -100001.
+Re(dm) = 0.
+Im(dm) = 1.
+interval = log
+N = 20
+eps = 1e-5
+#
+#Data
+data = empty_inclusion_data
+#
+#Discretization
+h = 0.25

input/paper1/gradient_highly.txt

+name = gradient_highly
+#Objective
+objective = l2distance
+#Configuration
+inclusion = none
+xe = 5.0
+xr = 20.0
+#
+unknown = background
+Re(mb) = 1.0
+Im(mb) = 0.
+#
+#Scale
+scale0 = 1
+scaleN = 100
+interval = log
+N = 2
+#
+#Data
+data = filled_inclusion_data
+#
+#Discretization
+h = 0.25
+model_FunctionSpaceType = HierarchicalH1
+model_IntegrationType = Gauss
+model_FunctionSpaceDegree = 5

input/paper1/synthetic_highly.txt

+name = filled_inclusion
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cylinder
+filled = 1
+r = 2.5
+#
+unknown = none
+Re(mb) = 1.
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = 0.
+#
+#Discretization
+h = 0.2

input/paper1/synthetic_weakly.txt

+name = empty_inclusion
+#Configuration
+xe = 5.0
+xr = 5.0
+inclusion = cylinder
+filled = 0
+r = 2.5
+#
+unknown = none
+Re(mb) = 1.
+Im(mb) = 0.
+Re(mc) = 1.2
+Im(mc) = -100.
+#
+#Discretization
+h = 0.2

specific/configuration/green0_preconditioner.h

+#ifndef H_GREEN0_PRECONDITIONER
+#define H_GREEN0_PRECONDITIONER
+
+//Standard Library
+
+//GmshFEM Library
+#include "Exception.h"
+
+//GmshFWI Library
+#include "soil.h"
+#include "inclusion.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+
+template<Physic T_Physic>
+ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const ConfigurationInterface* const config)
+{
+    if (const soil::Configuration* const soil_config = dynamic_cast<const soil::Configuration* const>(config))
+    {
+        return soil::green0_preconditioner(g,soil_config);
+    }
+    else if (const inclusion::Configuration* const inclusion_config = dynamic_cast<const inclusion::Configuration* const>(config))
+    {
+        return inclusion::green0_preconditioner(g,inclusion_config);
+    }
+    else
+    {
+        throw Exception("Dynamic cast to specific configuration has failed.");
+        return 0.;
+    }
+}
+
+#endif // H_GREEN0_PRECONDITIONER

specific/configuration/inclusion.cpp

+//Standard Library
+#include <numeric>
+//GmshFem Library
+#include "GmshFem.h"
+#include "Exception.h"
+#include "Message.h"
+#include "Domain.h"
+#include "Function.h"
+//#include "Formulation.h"
+#include "CSVio.h"
+//Gmsh Library
+#include "gmsh.h"
+//Fwi library
+#include "../../common/wave/correlation.h"
+#include "inclusion.h"
+
+namespace gmodel = gmsh::model;
+namespace factory = gmsh::model::geo;
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::domain;
+using namespace gmshfem::function;
+
+static const std::complex< double > im = std::complex< double >(0., 1.);
+
+namespace inclusion
+{
+    Inclusion to_inclusion(std::string inclusion)
+    {
+        if(inclusion=="none"){return Inclusion::None;}
+        else if(inclusion=="cylinder"){return Inclusion::Cylinder;}
+        else
+        {
+            throw common::Exception("Unknown inclusion type:"+inclusion+".");
+            return Inclusion::None;
+        }
+    }
+    UnknownRegion to_unknownregion(const GmshFem& gmshFem)
+    {
+        std::string str;
+        if(! gmshFem.userDefinedParameter(str, "unknown") )
+        {
+        throw common::Exception("The unknown region to consider could not be found.");
+        }
+        if(str=="all") {return UnknownRegion::All;}
+        else if(str=="inclusion") {return UnknownRegion::Inclusion;}
+        else if(str=="background") {return UnknownRegion::Background;}
+        else if(str=="none") {return UnknownRegion::None;}
+        else
+        {
+        throw common::Exception("The unknown region " + str + " is not handled.");
+        }
+        return UnknownRegion::None;
+    }
+    /*
+    *  class Configuration
+    */
+    Configuration::Configuration(const GmshFem& gmshFem) : ConfigurationInterface(gmshFem)
+    {
+        msg::print << "Initialize inclusion configuration." << msg::endl;
+
+        /*
+        * MESH
+        */
+        std::string inclusion;
+        if(
+          !(
+            gmshFem.userDefinedParameter(_xe, "xe") &&
+            gmshFem.userDefinedParameter(_xr, "xr") &&
+            gmshFem.userDefinedParameter(_ye, "ye") &&
+            gmshFem.userDefinedParameter(_yr, "yr") &&
+            gmshFem.userDefinedParameter(_He, "He") &&
+            gmshFem.userDefinedParameter(_Hr, "Hr") &&
+            gmshFem.userDefinedParameter(_Le, "Le") &&
+            gmshFem.userDefinedParameter(_Lr, "Lr") &&
+            gmshFem.userDefinedParameter(_nxe, "nxe") &&
+            gmshFem.userDefinedParameter(_nxr, "nxr") &&
+            gmshFem.userDefinedParameter(_nye, "nye") &&
+            gmshFem.userDefinedParameter(_nyr, "nyr") &&
+            gmshFem.userDefinedParameter(_L, "L") &&
+            gmshFem.userDefinedParameter(_H, "H") &&
+            gmshFem.userDefinedParameter(inclusion, "inclusion") &&
+            gmshFem.userDefinedParameter(_h, "h")
+            )
+        )
+        {
+          throw common::Exception("A geometric parameter could not be found.");
+        }
+        _inclusion_type = to_inclusion(inclusion);
+        if(_inclusion_type != Inclusion::None)
+        {
+          if (!gmshFem.userDefinedParameter(_r, "r"))
+          {
+            throw common::Exception("Inclusion size could not be found.");
+          }
+        }
+        unsigned int filled = 0;
+        gmshFem.userDefinedParameter(filled, "filled");
+        _isFilled = ((bool) filled);
+
+        mesh();
+
+        /*
+        * DOMAIN
+        */
+        if(nxe() != 0){_data_omega |= Domain(2,1);}
+        if(nye() != 0){_data_omega |= Domain(2,2);}
+
+        _background[Support::BLK] = Domain(2,21);
+        _background[Support::BND] = Domain(1,12);
+
+        if(_inclusion_type != Inclusion::None)
+        {
+            if(_isFilled){_inclusion[Support::BLK] = Domain(2,22);}
+            else{_inclusion[Support::BND] = Domain(1,13);}
+        }
+
+        _wave_omega[Support::BLK] = _background[Support::BLK] | _inclusion[Support::BLK];
+        _wave_omega[Support::BND] = _background[Support::BND] | _inclusion[Support::BND];
+        _unknown_region = to_unknownregion(gmshFem);
+        switch (_unknown_region)
+        {
+          case UnknownRegion::Inclusion:
+            _model_known = _background;
+            _model_unknown = _inclusion;
+            break;
+          case UnknownRegion::Background:
+            _model_known = _inclusion;
+            _model_unknown = _background;
+            break;
+          case UnknownRegion::All:
+            _model_unknown[Support::BLK] = _background[Support::BLK] | _inclusion[Support::BLK];
+            _model_unknown[Support::BND] = _background[Support::BND] | _inclusion[Support::BND];
+            break;
+          case UnknownRegion::None: default:
+            _model_known[Support::BLK] = _background[Support::BLK] | _inclusion[Support::BLK];
+            _model_known[Support::BND] = _background[Support::BND] | _inclusion[Support::BND];
+            break;
+        }
+
+        _ns = _nxe+_nye;
+        _np = _nxe+_nye;
+        //Emitters X
+        for (unsigned int e = 0; e < _nxe; e++)
+        {
+          _point.push_back(Domain(0,100+e));
+        }
+        _emitter_idx_X.resize(_nxe);
+        std::iota(_emitter_idx_X.begin(),_emitter_idx_X.end(),0);
+
+        //Emitters Y
+        for (unsigned int e = 0; e < _nye; e++)
+        {
+          _point.push_back(Domain(0,300+e));
+        }
+        _emitter_idx_Y.resize(_nye);
+        std::iota(_emitter_idx_Y.begin(),_emitter_idx_Y.end(),_nxe);
+
+        //Receivers X
+        if(!receiversAreEmittersX())
+        {
+          _np+=_nxr;
+          for (unsigned int r = 0; r < _nxr; r++)
+          {
+            _point.push_back(Domain(0,200+r));
+          }
+          _receiver_idx_X.resize(_nxr);
+          std::iota(_receiver_idx_X.begin(),_receiver_idx_X.end(),_nxe+_nye);
+        }
+        else
+        {
+            _receiver_idx_X = _emitter_idx_X;
+        }
+
+        //Receivers Y
+        if(!receiversAreEmittersY())
+        {
+          _np+=_nyr;
+          for (unsigned int r = 0; r < _nyr; r++)
+          {
+            _point.push_back(Domain(0,400+r));
+          }
+          _receiver_idx_Y.resize(_nyr);
+          std::iota(_receiver_idx_Y.begin(),_receiver_idx_Y.end(),_nxe+_nye+_nxr);
+        }
+        else
+        {
+            _receiver_idx_Y = _emitter_idx_Y;
+        }
+
+        for (unsigned int p = 0; p < _np; p++)
+        {
+          _points |= _point[p];
+        }
+
+        for (unsigned int s = 0; s < _nxe; s++)
+        {
+          _emitter.push_back({s});
+          std::vector<unsigned int> rec(_nxr,0);
+          if(receiversAreEmittersX())
+          {
+            std::iota(rec.begin(),rec.end(),0);
+            rec.erase(rec.begin()+s);
+            _receiver.push_back(rec);
+          }
+          else
+          {
+            std::iota(rec.begin(),rec.end(),_nxe+_nye);
+            _receiver.push_back(rec);
+          }
+        }
+        for (unsigned int s = _nxe; s < _nxe+_nye; s++)
+        {
+          _emitter.push_back({s});
+          std::vector<unsigned int> rec(_nyr,0);
+          if(receiversAreEmittersY())
+          {
+            std::iota(rec.begin(),rec.end(),_nxe);
+            rec.erase(rec.begin()+s-_nxe);
+            _receiver.push_back(rec);
+          }
+          else
+          {
+            std::iota(rec.begin(),rec.end(),_nxe+_nye+_nxr);
+            _receiver.push_back(rec);
+          }
+        }
+
+        /*
+        * Reference model function
+        */
+        double Remb, Immb;
+        if
+        (
+          !(
+            gmshFem.userDefinedParameter(Remb, "Re(mb)") &&
+            gmshFem.userDefinedParameter(Immb, "Im(mb)")
+          )
+        )
+        {
+          throw common::Exception("Background model parameter could not be found.");
+        }
+        _mb = Remb + im * Immb;
+
+        double Remc=0., Immc=0.;
+        if
+        (
+          !(gmshFem.userDefinedParameter(Remc, "Re(mc)") &&
+            gmshFem.userDefinedParameter(Immc, "Im(mc)")
+          )
+          &&
+          (_inclusion_type!=Inclusion::None)
+        )
+        {
+          throw common::Exception("Inclusion model parameter could not be found.");
+        }
+        _mc = Remc + im * Immc;
+
+        ScalarPiecewiseFunction< std::complex< double > > m0;
+        m0.addFunction(_mb,_background[Support::BLK] | _background[Support::BND] | _points);
+        m0.addFunction(_mc,_inclusion[Support::BLK] | _inclusion[Support::BND]);
+        _m0 = m0.copy();
+    }
+
+    void Configuration::mesh() const
+    {
+        msg::print << "Generate meshes" << msg::endl;
+        gmsh::option::setNumber("General.Terminal", 1);
+
+        std::string name = "inclusion";
+        switch (_inclusion_type)
+        {
+            case Inclusion::Cylinder: name += "_cyl"; break;
+            case Inclusion::None: name += "_none"; break;
+            default: break;
+        }
+        if(_isFilled && (_inclusion_type != Inclusion::None)){name += "_filled";}
+        if((!_isFilled) && (_inclusion_type != Inclusion::None)){name += "_empty";}
+
+        gmodel::add("inclusion");
+        if(!_remesh)
+        {
+          gmsh::open(name+".msh");
+          return;
+        }
+
+        wave_mesh();
+        data_mesh();
+
+        factory::synchronize();
+        gmodel::mesh::generate();
+
+        gmsh::write(name+".msh");
+    }
+    void Configuration::wave_mesh() const
+    {
+        int pb1 = factory::addPoint(0., 0., 0., _h);
+        int pb2 = factory::addPoint(_L, 0., 0., _h);
+        int pb3 = factory::addPoint(_L, _H, 0., _h);
+        int pb4 = factory::addPoint(0., _H, 0., _h);
+
+        int lbh1 = factory::addLine(pb1, pb2);
+        int lbv1 = factory::addLine(pb2, pb3);
+        int lbh2 = factory::addLine(pb3, pb4);
+        int lbv2 = factory::addLine(pb4, pb1);
+
+        int clb = factory::addCurveLoop({lbh1,lbv1,lbh2,lbv2});
+
+        //Vertical emitters
+        std::vector<int> px(_nxe);
+        double He0 = (_H - _He)/2.;
+        double deltaHe = _He/(((double) _nxe)-1.);
+        for (unsigned int e = 0; e < _nxe; e++)
+        {
+            double Hee = He0 + ((double)e)*deltaHe;
+            px[e] = factory::addPoint(_xe, Hee, 0., _h);
+        }
+
+        if(!receiversAreEmittersX())
+        {
+        //Vertical receivers
+        px.resize(_nxe+_nxr);
+        double Hr0 = (_H - _Hr)/2.;
+        double deltaHr = _Hr/(((double) _nxr)-1.);
+        for (unsigned int r = 0; r < _nxr; r++)
+        {
+            double Hrr = Hr0 + ((double)r)*deltaHr;
+            px[_nxe+r] = factory::addPoint(_xr, Hrr, 0., _h);
+        }
+        }
+        //Horizontal emitters
+        std::vector<int> py(_nye);
+        double Le0 = (_L - _Le)/2.;
+        double deltaLe = _Le/(((double) _nye)-1.);
+        for (unsigned int e = 0; e < _nye; e++)
+        {
+            double Lee = Le0 + ((double)e)*deltaLe;
+            py[e] = factory::addPoint(Lee, _ye, 0., _h);
+        }
+        if(!receiversAreEmittersY())
+        {
+        //Horizontal receivers
+        py.resize(_nye+_nyr);
+        double Lr0 = (_L - _Lr)/2.;
+        double deltaLr = _Lr/(((double) _nyr)-1.);
+        for (unsigned int r = 0; r < _nyr; r++)
+        {
+            double Lrr = Lr0 + ((double)r)*deltaLr;
+            py[_nye+r] = factory::addPoint(Lrr,_yr, 0., _h);
+        }
+        }
+
+        int si=0;
+        int sb=0;
+        std::vector<int> lci{};
+        switch (_inclusion_type)
+        {
+        case Inclusion::Cylinder:
+        {
+            int pc0 = factory::addPoint(_L/2, _H/2, 0., _h);
+            int pc1 = factory::addPoint(_L/2+_r, _H/2, 0., _h);
+            int pc2 = factory::addPoint(_L/2, _H/2+_r, 0., _h);
+            int pc3 = factory::addPoint(_L/2-_r, _H/2, 0., _h);
+            int pc4 = factory::addPoint(_L/2, _H/2-_r, 0., _h);
+
+            int lc1 = factory::addCircleArc(pc1, pc0, pc2);
+            int lc2 = factory::addCircleArc(pc2, pc0, pc3);
+            int lc3 = factory::addCircleArc(pc3, pc0, pc4);
+            int lc4 = factory::addCircleArc(pc4, pc0, pc1);
+
+            lci = {lc1,lc2,lc3,lc4};
+
+            int cli = factory::addCurveLoop(lci);
+
+            if(_isFilled){si = factory::addPlaneSurface({cli});}
+            sb = factory::addPlaneSurface({clb,cli});
+            break;
+        }
+        case Inclusion::None:
+        {
+            sb = factory::addPlaneSurface({clb});
+            break;
+        }
+        }
+
+        factory::synchronize();
+        gmodel::mesh::embed(0, px, 2, sb);
+        gmodel::mesh::embed(0, py, 2, sb);
+
+        if(receiversAreEmittersX())
+        {
+            for (unsigned int e = 0; e <_nxe; e++)
+            {
+                gmodel::addPhysicalGroup(0, {px[e]}, 100+e);
+                gmodel::setPhysicalName(0, 100+e, "emitter_receiver_x_"+std::to_string(e));
+            }
+        }
+        else
+        {
+        for (unsigned int e = 0; e <_nxe; e++)
+        {
+            gmodel::addPhysicalGroup(0, {px[e]}, 100+e);
+            gmodel::setPhysicalName(0, 100+e,"emitter_x_"+std::to_string(e));
+        }
+        for (unsigned int r = 0; r < _nxr; r++)
+        {
+            gmodel::addPhysicalGroup(0, {px[_nxe+r]}, 200+r);
+            gmodel::setPhysicalName(0, 200+r, "receiver_x_"+std::to_string(r) );
+        }
+        }
+        if(receiversAreEmittersY())
+        {
+            for (unsigned int e = 0; e <_nye; e++)
+            {
+                gmodel::addPhysicalGroup(0, {py[e]}, 300+e);
+                gmodel::setPhysicalName(0, 300+e,"emitter_receiver_y_"+std::to_string(e));
+            }
+        }
+        else
+        {
+            for (unsigned int e = 0; e <_nye; e++)
+            {
+                gmodel::addPhysicalGroup(0, {py[e]}, 300+e);
+                gmodel::setPhysicalName(0, 300+e,"emitter_y_"+std::to_string(e));
+            }
+            for (unsigned int r = 0; r < _nyr; r++)
+            {
+                gmodel::addPhysicalGroup(0, {py[_nye+r]}, 400+r);
+                gmodel::setPhysicalName(0, 400+r, "receiver_y_"+std::to_string(r) );
+            }
+        }
+
+        gmodel::addPhysicalGroup(1, {lbh1,lbv1,lbh2,lbv2}, 12);
+        gmodel::setPhysicalName(1, 12, "background_bnd");
+        gmodel::addPhysicalGroup(2, {sb}, 21);
+        gmodel::setPhysicalName(2, 21, "background_vol");
+        if(_inclusion_type != Inclusion::None)
+        {
+            if(_isFilled)
+            {
+                gmodel::addPhysicalGroup(2, {si}, 22);
+                gmodel::setPhysicalName(2, 22, "inclusion_vol");
+            }
+            else
+            {
+                gmodel::addPhysicalGroup(1, lci, 13);
+                gmodel::setPhysicalName(1, 13, "inclusion_bnd");
+            }
+        }
+    }
+    void Configuration::data_mesh() const
+    {
+        if(nxe() != 0)
+        {
+            //Vertical array
+            int p00 = factory::addPoint(0., 0., 0.);
+            int pHe0 = factory::addPoint(_He, 0., 0.);
+            int pHeHr = factory::addPoint(_He, _Hr, 0.);
+            int p0Hr = factory::addPoint(0., _Hr, 0.);
+
+            int ly0 = factory::addLine(p00, pHe0);
+            int lxHe = factory::addLine(pHe0, pHeHr);
+            int lyHr = factory::addLine(pHeHr, p0Hr);
+            int lx0 = factory::addLine(p0Hr, p00);
+
+            int cl = factory::addCurveLoop({ly0,lxHe,lyHr,lx0});
+            int s1 = factory::addPlaneSurface({cl});
+
+            factory::mesh::setTransfiniteCurve(ly0,_nxe);
+            factory::mesh::setTransfiniteCurve(lxHe,_nxr);
+            factory::mesh::setTransfiniteCurve(lyHr,_nxe);
+            factory::mesh::setTransfiniteCurve(lx0,_nxr);
+
+            factory::mesh::setTransfiniteSurface(s1,"Left",{p00,pHe0,pHeHr,p0Hr});
+
+            factory::mesh::setRecombine(2,s1);
+            gmodel::addPhysicalGroup(2, {s1}, 1);
+            gmodel::setPhysicalName(2, 1, "data_omega1");
+        }
+
+        if(nye() != 0)
+        {
+            //Horizontal array
+            int pHeHr = factory::addPoint(_He, _Hr, 0.);
+            int pLe0 = factory::addPoint(_He+_Le, _Hr+0., 0.);
+            int pLeLr = factory::addPoint(_He+_Le, _Hr+_Lr, 0.);
+            int p0Lr = factory::addPoint(_He+0., _Hr+_Lr, 0.);
+
+            int lyHr2 = factory::addLine(pHeHr, pLe0);
+            int lxLe = factory::addLine(pLe0, pLeLr);
+            int lyLr = factory::addLine(pLeLr, p0Lr);
+            int lxHe2 = factory::addLine(p0Lr, pHeHr);
+
+            int cl2 = factory::addCurveLoop({lyHr2,lxLe,lyLr,lxHe2});
+            int s2 = factory::addPlaneSurface({cl2});
+
+            factory::mesh::setTransfiniteCurve(lyHr2,_nxe);
+            factory::mesh::setTransfiniteCurve(lxLe,_nxr);
+            factory::mesh::setTransfiniteCurve(lyLr,_nxe);
+            factory::mesh::setTransfiniteCurve(lxHe2,_nxr);
+
+            factory::mesh::setTransfiniteSurface(s2,"Left",{pHeHr,pLe0,pLeLr,p0Lr});
+
+            factory::mesh::setRecombine(2,s2);
+            gmodel::addPhysicalGroup(2, {s2}, 2);
+            gmodel::setPhysicalName(2, 2, "data_omega2");
+        }
+    }
+
+    template<Physic T_Physic>
+    ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
+    {
+        return 1.;
+    }
+
+    template ModelFunction green0_preconditioner<Physic::acoustic>(const WaveMultiField<Physic::acoustic>&, const Configuration* const);
+
+} // namespace soil