Review paper 1 + emetters/receivers mistake fixing + beginning of presentation 2

common/configuration.cpp

 //Standard library
-//GmshFem library
-/*
-#include "Exception.h"
-#include "Message.h"
-#include "Formulation.h"
-*/
-//FWI Library
+//GmshFEM Library
+//GmshFWI Library
+
 #include "configuration.h"
-/*
-using namespace gmshfem::field;
-using namespace gmshfem::function;
-using namespace gmshfem::problem;
-using namespace gmshfem::equation;
-*/
+
 using namespace gmshfem;
 using namespace gmshfem::common;
 using namespace gmshfem::domain;
@@ -61,34 +52,34 @@ const std::vector<unsigned int>& ConfigurationInterface::receiver_idx(unsigned i
 }
 unsigned int ConfigurationInterface::isEmitter(unsigned int point) const
 {
-  std::vector<unsigned int> idx = emitter_idx();
+  std::vector<unsigned int> idx = all_emitter_idx();
   return (std::find(idx.begin(), idx.end(), point) != idx.end());
 }
 unsigned int ConfigurationInterface::isReceiver(unsigned int point) const
 {
-  std::vector<unsigned int> idx = receiver_idx();
+  std::vector<unsigned int> idx = all_receiver_idx();
   return (std::find(idx.begin(), idx.end(), point) != idx.end());
 }
 Domain ConfigurationInterface::emitter(unsigned int shot, unsigned int em) const
 {
   emIsValid(shot,em);
-  return _point[emitter_idx(shot)[em]];
+  return _point[_emitter[shot][em]];
 }
 Domain ConfigurationInterface::receiver(unsigned int shot, unsigned int rec) const
 {
   recIsValid(shot,rec);
-  return _point[receiver_idx(shot)[rec]];
+  return _point[_receiver[shot][rec]];
 }
 Domain ConfigurationInterface::point(unsigned int idx) const
 {
   pntIsValid(idx);
   return _point[idx];
 }
-std::vector<unsigned int> ConfigurationInterface::emitter_idx(bool includeReceivers) const
+std::vector<unsigned int> ConfigurationInterface::all_emitter_idx(bool includeReceivers) const
 {
   std::vector<unsigned int> em_idx;
   std::vector<unsigned int> rec_idx;
-  if(!includeReceivers){rec_idx = receiver_idx(true);}
+  if(!includeReceivers){rec_idx = all_receiver_idx(true);}
   for (unsigned int s = 0; s < ns(); s++)
   {
     for (unsigned int e = 0; e < ne(s); e++)
@@ -106,11 +97,11 @@ std::vector<unsigned int> ConfigurationInterface::emitter_idx(bool includeReceiv
   }
   return em_idx;
 }
-std::vector<unsigned int> ConfigurationInterface::receiver_idx(bool includeEmitters) const
+std::vector<unsigned int> ConfigurationInterface::all_receiver_idx(bool includeEmitters) const
 {
   std::vector<unsigned int> rec_idx;
   std::vector<unsigned int> em_idx;
-  if(!includeEmitters){em_idx = emitter_idx(true);}
+  if(!includeEmitters){em_idx = all_emitter_idx(true);}
   for (unsigned int s = 0; s < ns(); s++)
   {
     for (unsigned int r = 0; r < nr(s); r++)

common/configuration.h

 #ifndef H_CONFIGURATION_INTERFACE
 #define H_CONFIGURATION_INTERFACE
 
-//GmshFEM library
+//GmshFEM Library
 #include "GmshFem.h"
 #include "Domain.h"
 //Standard Library
 #include <array>
-/*
-//Standard Library
-#include <string>
-#include "Function.h"
-#include "Exception.h"
-//#include "Message.h"
-#include "FieldForm0.h"
-#include "discretization.h"
-#include "../specific/physics.h"
-*/
-//GmshFWI
+//GmshFWI Library
 #include "model/element.h"
 #include "support.h"
 
@@ -63,6 +53,7 @@ public:
 
     virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const = 0;
     virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const = 0;
+    virtual bool data_coordinate_isValid(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];};
@@ -78,8 +69,8 @@ public:
     gmshfem::domain::Domain point(unsigned int idx) const;
     const std::vector<unsigned int>& emitter_idx(unsigned int shot) const;
     const std::vector<unsigned int>& receiver_idx(unsigned int shot) const;
-    std::vector<unsigned int> emitter_idx(bool includeReceivers = false) const;
-    std::vector<unsigned int> receiver_idx(bool includeEmitters = false) const;
+    std::vector<unsigned int> all_emitter_idx(bool includeReceivers = false) const;
+    std::vector<unsigned int> all_receiver_idx(bool includeEmitters = false) const;
     unsigned int isEmitter(unsigned int point) const;
     unsigned int isReceiver(unsigned int point) const;
     unsigned int np() const {return _np;};

common/data/element.cpp

@@ -100,6 +100,8 @@ PointData<T_Physic> Data<T_Physic>::value(unsigned int s,unsigned int r) const
 
 template<Physic T_Physic>
 PointData<T_Physic> Data<T_Physic>::value(double xs, double xr) const
+{
+    if(_config->data_coordinate_isValid(xs,xr))
     {
         std::array<unsigned int,2> index = _config->data_coordinate_to_index(xs,xr);
         unsigned int s = index[0];
@@ -107,6 +109,12 @@ PointData<T_Physic> Data<T_Physic>::value(double xs, double xr) const
         isValid(s,r);
         return _value[s][r];
     }
+    else
+    {
+        msg::warning << "Data coordinate: ("<< std::to_string(xs) << "," << std::to_string(xr) << ")" << " is not valid." << msg::endl;
+        return std::complex<double>(0.);
+    }
+}
 
 template<Physic T_Physic>
 void Data<T_Physic>::value(unsigned int s,unsigned int r,const PointData<T_Physic>& value)

common/wave/equation/equation.cpp

@@ -169,14 +169,14 @@ void DifferentialEquationInterface<T_Physic>::updateGreen(const ModelStateEvalua
 template<Physic T_Physic>
 void DifferentialEquationInterface<T_Physic>::updateGreenEmitter(const ModelStateEvaluator& m)
 {
-    updateGreen(m,_config->emitter_idx(!_grAreUpToDate));
+    updateGreen(m,_config->all_emitter_idx(!_grAreUpToDate));
     _geAreUpToDate=true;
 }
 
 template<Physic T_Physic>
 void DifferentialEquationInterface<T_Physic>::updateGreenReceiver(const ModelStateEvaluator& m)
 {
-    updateGreen(m,_config->receiver_idx(!_geAreUpToDate));
+    updateGreen(m,_config->all_receiver_idx(!_geAreUpToDate));
     _grAreUpToDate=true;
 }
 

common/wave/state.cpp

@@ -19,4 +19,18 @@ const WaveMultiField<T_Physic>& WaveState<T_Physic>::state(Type type) const
     }
 }
 
+template<Physic T_Physic>
+const WaveField<T_Physic>& WaveState<T_Physic>::state(Type type, unsigned int index) const
+{
+    if(!(index < _state[((unsigned int) type)].size()))
+    {
+        throw Exception("Wave field index is bigger than vector size.");
+    }
+    if(!_isUpToDate[((unsigned int) type)])
+    {
+        throw Exception(type_to_string(type) + " wave is not up to date while required.");
+    }
+    return (_state[((unsigned int) type)])[index];
+}
+
 template class WaveState<Physic::acoustic>;

common/wave/state.h

@@ -18,6 +18,7 @@ class WaveStateEvaluator
 {
 public:
     virtual const WaveMultiField<T_Physic>& state(Type type) const = 0;
+    virtual const WaveField<T_Physic>& state(Type type, unsigned int index) const = 0;
     void write(Type type, std::string name) const {state(type).write(name);};
 };
 
@@ -40,6 +41,7 @@ public:
     }, _isUpToDate{false,false,false,false} {};
 
     virtual const WaveMultiField<T_Physic>& state(Type type) const;
+    virtual const WaveField<T_Physic>& state(Type type, unsigned int index) const;
 
     friend class WaveUpdater<T_Physic>;
 };

directional.cpp

@@ -29,6 +29,18 @@ Interval to_interval(const GmshFem& gmshFem)
   return Interval::Linear;
 }
 
+enum class SpatialDistribution: unsigned int { Constant=0, HeavisideY0=1 };
+SpatialDistribution to_spatialdistribution(const GmshFem& gmshFem)
+{
+  std::string buffer;
+  if(gmshFem.userDefinedParameter(buffer, "spatialdistribution"))
+  {
+    if(buffer=="constant"){return SpatialDistribution::Constant;}
+    else if(buffer=="heavisidey0"){return SpatialDistribution::HeavisideY0;}
+  }
+  return SpatialDistribution::Constant;
+}
+
 template <Physic T_Physic>
 int directional(const GmshFem& gmshFem)
 {
@@ -73,6 +85,28 @@ int directional(const GmshFem& gmshFem)
         throw Exception("A directional parameter could not be found.");
     }
 
+    double Remjp = 0.;
+    double Immjp = 0.;
+    double yjp = 0.;
+    std::complex<double> mjp = 0.;
+    SpatialDistribution spatialdistribution = to_spatialdistribution(gmshFem);
+    switch (spatialdistribution)
+    {
+        case SpatialDistribution::HeavisideY0:
+            if(!(
+                gmshFem.userDefinedParameter(Remjp, "Re(mjump)") &&
+                gmshFem.userDefinedParameter(Immjp, "Im(mjump)") &&
+                gmshFem.userDefinedParameter(yjp, "yjump")
+                )
+            )
+            {
+                throw Exception("A spatial distribution parameter could not be found.");
+            }
+            mjp = Remjp + im * Immjp;
+            break;
+        case SpatialDistribution::Constant:
+            break;
+    }
 
     Interval interval = to_interval(gmshFem);
 
@@ -106,8 +140,17 @@ int directional(const GmshFem& gmshFem)
         std::complex<double> mn = mu0 + step*dir;
         msg::print << "mn = " << mn << msg::endl;
         ScalarPiecewiseFunction<std::complex<double>> m;
+        switch (spatialdistribution)
+        {
+            case SpatialDistribution::HeavisideY0:
+                m.addFunction(
+                    mn + (2. * heaviside( y<std::complex<double>>() - yjp ) -1. ) * mjp,
+                    configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
+                break;
+            case SpatialDistribution::Constant:
                 m.addFunction(mn,configuration->model_unknown(Support::BLK) | configuration->model_unknown(Support::BND));
-
+                break;
+        }
 
         double j =  functional->performance(m);
         msg::print << "j: "<< j << msg::endl;

gradient.cpp

@@ -13,7 +13,7 @@ using namespace gmshfem::common;
 #include "specific/configuration/newConfiguration.h"
 #include "specific/wave/equation/newEquation.h"
 #include "specific/data/objective/newObjective.h"
-#include "specific/model/innerproduct/newInnerProduct.h"
+#include "specific/model/innerproduct/sobolev.h"
 #include "common/statefunctional.h"
 
 static const std::complex< double > im = std::complex< double >(0., 1.);

input/MP2008LinearVelocity2D/common.txt

@@ -3,24 +3,31 @@ name = MP2008LinearVelocity2D
 configuration = soil
 #
 xoffset = 0
-Ler = 4.8
-nr = 301
+Ler = 4.356
+nr = 364
 ne = 1
-ymax = 2.4
-L = 7.6
-H = 2.4
+ReceiverOnEmitter=0
+ymax = 2.904
+L = 9.192
+H = 2.904
 #
 Re(v_soil) = 0.
 Im(v_soil) = 0.
 #
+#
+#
 #Equation
 physic = acoustic
 equation = helmholtz
 #
 #Discretization
+h = 0.012
 #Wave
 wave_IntegrationType = Gauss
 wave_FunctionSpaceType = HierarchicalH1
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceDegree = 1
 #
 n_freq = 3
 #Frequency0
@@ -28,9 +35,9 @@ frequency0 = 4.
 wave_FunctionSpaceDegree0 = 1
 #
 #Frequency1
-frequency1 = 7.
-wave_FunctionSpaceDegree1 = 2
+frequency1 = 6.
+wave_FunctionSpaceDegree1 = 1
 #
 #Frequency2
-frequency2 = 13.
-wave_FunctionSpaceDegree2 = 3
+frequency2 = 9.
+wave_FunctionSpaceDegree2 = 2

input/MP2008LinearVelocity2D/multiscale.txt

@@ -11,12 +11,10 @@ v0N = 2.0
 Nv = 0
 alpha0 = 0.35
 alphaN = 0.45
-Na = 20
+Na = 30
 #
 #Objective
 objective = l2distance
-data_scaled = 1
-data_norm = 301
 #
 #Data
 data0 = MP2008LinearVelocity2D_data0
@@ -24,6 +22,5 @@ data1 = MP2008LinearVelocity2D_data1
 data2 = MP2008LinearVelocity2D_data2
 #
 #Discretization
-h = 0.032
 #Output
 write_fields=0

input/MP2008LinearVelocity2D/synthetics.txt

@@ -3,10 +3,9 @@ unknown = none
 m0_type = velocity_depth_noair
 Re(v_0) = 2.0
 Im(v_0) = 0.
-Re(v_H) = 2.96
+Re(v_H) = 3.1616
 Im(v_H) = 0.
 #
-#Discretization
-h = 0.016
 #Output
-write_fields = 1
+write_data_fields = 0
+write_wave_fields = 1

input/paper1/common.txt

 #Equation
 physic = acoustic
 equation = helmholtz
-useGreen = 1
 #
 #Configuration
 configuration = inclusion
@@ -23,6 +22,7 @@ H = 25.0
 wave_FunctionSpaceType = HierarchicalH1
 wave_IntegrationType = Gauss
 wave_FunctionSpaceDegree = 5
+model_IntegrationType = Gauss
 #
 #Frequency
 frequency0 = 1.

input/paper1/gradient_highly.txt

@@ -25,5 +25,4 @@ complex = 0
 #Discretization
 h = 0.25
 model_FunctionSpaceType = HierarchicalH1
-model_IntegrationType = Gauss
 model_FunctionSpaceDegree = 5

input/paper1/synthetic_highly.txt

@@ -13,4 +13,7 @@ Re(mc) = 1.2
 Im(mc) = 0.
 #
 #Discretization
+model_FunctionSpaceDegree = 1
 h = 0.2
+#Output
+write_wave_fields = 0

input/review_paper1_cross2/cluster_inversion.sh

+#!/bin/bash
+#
+#SBATCH --job-name=RevAG2020
+#SBATCH --output=RevAG2020.log
+
+#SBATCH --ntasks=2
+#SBATCH --time=180:00
+#SBATCH --mem-per-cpu=2000
+
+./synthetics ../input/review_paper1/common.txt ../input/review_paper1/synthetic_highly.txt -verbose 1
+./directional ../input/review_paper1/common.txt ../input/review_paper1/directional_highly_inclusion.txt -verbose 1
+./directional ../input/review_paper1/common.txt ../input/review_paper1/directional_highly_background.txt -verbose 1
+./convergence ../input/review_paper1/common.txt ../input/review_paper1/convergence_highly_background.txt -verbose 1
+./gradient ../input/review_paper1/common.txt ../input/review_paper1/gradient_highly.txt -verbose 1
+
+./synthetics ../input/review_paper1/common.txt ../input/review_paper1/synthetic_weakly.txt -verbose 1
+./directional ../input/review_paper1/common.txt ../input/review_paper1/directional_weakly_inclusion.txt -verbose 1
+./directional ../input/review_paper1/common.txt ../input/review_paper1/directional_weakly_background.txt -verbose 1
+./convergence ../input/review_paper1/common.txt ../input/review_paper1/convergence_weakly_background.txt -verbose 1

input/review_paper1_cross2/common.txt

+prename = RevAG2020_
+#Equation
+physic = acoustic
+equation = helmholtz
+useGreen = 1
+#
+#Configuration
+configuration = inclusion
+#
+receiver_on_emitter = 0
+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
+frequency0 = 1.

input/review_paper1_cross2/convergence_highly_background.txt

+name = highly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cross2
+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
+data0 = RevAG2020_synthetic_highly_data0
+#
+#Discretization
+h = 0.25

input/review_paper1_cross2/convergence_weakly_background.txt

+name = weakly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 5.0
+inclusion = cross2
+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
+data0 = RevAG2020_synthetic_weakly_data0
+#
+#Discretization
+h = 0.25

input/review_paper1_cross2/directional_highly_background.txt

+name = highly_background
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cross2
+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
+data0 = RevAG2020_synthetic_highly_data0
+#
+#Discretization
+h = 0.25

input/review_paper1_cross2/directional_highly_inclusion.txt

+name = highly_inclusion
+#Objective
+objective = l2distance
+#Configuration
+xe = 5.0
+xr = 20.0
+inclusion = cross2
+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
+data0 = RevAG2020_synthetic_highly_data0
+#
+#Discretization
+h = 0.25