Flexible acqusition (Rectangular domain with arbitrary sources and receivers)

CMakeLists.txt

@@ -185,6 +185,15 @@ endif(NOT GMSH_INC)
 list(APPEND EXTRA_INCS ${GMSH_INC})
 list(APPEND EXTRA_LIBS ${GMSH_LIB})
 
+# Yaml-cpp
+option(ENABLE_YAML_CPP "Enable YAML-CPP library" ON)
+
+if(ENABLE_YAML_CPP)
+    find_package(yaml-cpp REQUIRED)
+    add_definitions(-DENABLE_YAML_CPP)
+endif()
+
+
 # GmshFem
 
 find_library(GMSHFEM_LIB gmshfem)
@@ -227,30 +236,31 @@ include_directories(${EXTRA_INCS})
 file( GLOB LIB_SOURCES contrib/blossom5v2.05/*.cpp contrib/blossom5v2.05/*/*.cpp common/*/*/*.cpp common/*/*.cpp common/*.cpp specific/*/*/*.cpp specific/*/*.cpp specific/*.cpp)
 file( GLOB LIB_HEADERS contrib/blossom5v2.05/*.h contrib/blossom5v2.05/*/*.h common/*/*/*.h common/*/*.h common/*.h specific/*/*/*.h specific/*/*.h specific/*.h)
 add_library( fwi ${LIB_SOURCES} ${LIB_HEADERS})
+target_link_libraries(fwi yaml-cpp)
 
-add_executable( synthetics synthetics.cpp ${EXTRA_INCS})
+add_executable( synthetics synthetics.cpp )
 target_link_libraries( synthetics fwi  ${EXTRA_LIBS})
 #add_executable( time_synthetics_wave time_synthetics_wave.cpp ${EXTRA_INCS})
 #target_link_libraries( time_synthetics_wave fwi  ${EXTRA_LIBS})
-add_executable( directional directional.cpp ${EXTRA_INCS})
+add_executable( directional directional.cpp )
 target_link_libraries( directional fwi  ${EXTRA_LIBS})
-add_executable( directional_field directional_field.cpp ${EXTRA_INCS})
+add_executable( directional_field directional_field.cpp )
 target_link_libraries( directional_field fwi  ${EXTRA_LIBS})
-add_executable( gradient gradient.cpp ${EXTRA_INCS})
+add_executable( gradient gradient.cpp )
 target_link_libraries( gradient fwi  ${EXTRA_LIBS})
 #add_executable( convergence convergence.cpp ${EXTRA_INCS})
 #target_link_libraries( convergence fwi  ${EXTRA_LIBS})
 #add_executable( preconditioner preconditioner.cpp ${EXTRA_INCS})
 #target_link_libraries( preconditioner fwi  ${EXTRA_LIBS})
-add_executable( ip_comparison ip_comparison.cpp ${EXTRA_INCS})
+add_executable( ip_comparison ip_comparison.cpp )
 target_link_libraries( ip_comparison fwi  ${EXTRA_LIBS})
 #add_executable( ob_comparison ob_comparison.cpp ${EXTRA_INCS})
 #target_link_libraries( ob_comparison fwi  ${EXTRA_LIBS})
-add_executable( multiscale multiscale.cpp ${EXTRA_INCS})
+add_executable( multiscale multiscale.cpp )
 target_link_libraries( multiscale fwi  ${EXTRA_LIBS})
-add_executable( inversion inversion.cpp ${EXTRA_INCS})
+add_executable( inversion inversion.cpp )
 target_link_libraries( inversion fwi  ${EXTRA_LIBS})
-add_executable( error error.cpp ${EXTRA_INCS})
+add_executable( error error.cpp)
 target_link_libraries( error fwi  ${EXTRA_LIBS})
 
 

common/configuration.cpp

@@ -53,7 +53,7 @@ bool ConfigurationInterface::pntIsValid(unsigned int e_r) const
 {
   if( !(e_r<_np) )
   {
-    throw Exception("Emitter-receiver "+std::to_string(e_r)+" is out of scope.");
+    throw Exception("Emitter-receiver "+std::to_string(e_r)+" is out of scope. There are only "+std::to_string(_np)+" points.");
   }
   return true;
 }

common/configuration.h

@@ -51,6 +51,12 @@ public:
         gmshFem.userDefinedParameter(remesh, "remesh");
         _remesh = ((bool) remesh);
     };
+
+    ConfigurationInterface() : _name("dummy"), _parametrization(nullptr), _mc(1,0.)
+    {
+        _remesh = true;
+    };
+
     virtual ~ConfigurationInterface() {delete _parametrization;};
 
     virtual void mesh() const;

specific/configuration/flexible_acquisition.cpp

+//Standard Library
+#include <numeric>
+#include <array>
+
+//GmshFEM Library
+#include "GmshFem.h"
+#include "Exception.h"
+#include "Message.h"
+#include "Domain.h"
+#include "Function.h"
+#include "CSVio.h"
+
+//Gmsh Library
+#include "gmsh.h"
+
+//GmshFWI Library
+#include "../../common/model/macro.h"
+#include "../../common/data/element.h"
+#include "../wave/correlation.h"
+#include "flexible_acquisition.h"
+#include "yaml_interface.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 flexible_acquisition
+{
+    
+    /*
+    *  class Configuration
+    */
+    Configuration::Configuration(std::string name, const ParametrizationInterface* const parametrization, const GmshFem& gmshFem) : ConfigurationInterface(name, parametrization, gmshFem)
+    {
+        msg::print << "Initialize flexible configuration" << msg::endl;
+
+        /*
+        * MESH
+        */
+        if(
+          !(
+            gmshFem.userDefinedParameter(_H, "H") &&
+            gmshFem.userDefinedParameter(_L, "L") &&
+            gmshFem.userDefinedParameter(_Hsub, "Hsub") &&
+            gmshFem.userDefinedParameter(_h, "h")
+            )
+        )
+        {
+          throw common::Exception("A geometric parameter could not be found.");
+        }
+
+        // SETUP FROM YAML
+        std::string shotsPath;
+        if (!gmshFem.userDefinedParameter(shotsPath, "shots_config"))
+        {
+          throw common::Exception("Path to a shots_config is needed in flexible_acquisition.");
+        }
+        ShotsConfigurationYAML shotsConfig(shotsPath);
+        _ns = shotsConfig.numShots();
+
+        // Configure points
+        for (const auto &coords : shotsConfig.points())
+        {
+          _er_positions.push_back(coords);
+        }
+
+        _np = _er_positions.size();
+
+
+        mesh();
+
+        
+        // Configure shots
+        for (unsigned s = 0; s < shotsConfig.numShots(); ++s)
+        {
+            _emitter.push_back(shotsConfig.emitters(s));
+            _receiver.push_back(shotsConfig.receivers(s));
+        }
+
+        for (unsigned p = 0; p < _er_positions.size(); ++p) {
+            _point.push_back("emitter_receiver_"+std::to_string(p));
+            _points |= _point[p];
+        }
+
+
+
+        /*
+        * DOMAIN
+        */
+       // TODO: actual data mesh!
+        _data_omega = Domain("supersurface") | Domain("subsurface");
+
+        _subsurface[Support::BLK] = Domain("subsurface");
+        _subsurface[Support::BND] = Domain("subsurface_bnd");
+        _supersurface[Support::BLK] = Domain("supersurface");
+        _supersurface[Support::BND] = Domain("supersurface_bnd");
+
+        std::string unknownRegion = "subsurface";
+        gmshFem.userDefinedParameter(unknownRegion, "unknown");
+        if (unknownRegion == "subsurface")
+        {
+          _model_known[Support::BLK] = _supersurface[Support::BLK];
+          _model_known[Support::BND] = _supersurface[Support::BND];
+
+          _model_unknown[Support::BLK] = _subsurface[Support::BLK];
+          _model_unknown[Support::BND] = _subsurface[Support::BND];
+        }
+        else if (unknownRegion == "none")
+        {
+          _model_known[Support::BLK] = _supersurface[Support::BLK] | _subsurface[Support::BLK];
+          _model_known[Support::BND] = _supersurface[Support::BND] | _subsurface[Support::BND];
+        }
+        else if (unknownRegion == "all")
+        {
+          _model_unknown[Support::BLK] = _supersurface[Support::BLK] | _subsurface[Support::BLK];
+          _model_unknown[Support::BND] = _supersurface[Support::BND] | _subsurface[Support::BND];
+        }
+        else {
+            throw Exception("Invalid unknown region type: " + unknownRegion);
+        }
+
+        _wave_omega[Support::BLK] = _model_known[Support::BLK] | _model_unknown[Support::BLK];
+        _wave_omega[Support::BND] = _model_known[Support::BND] | _model_unknown[Support::BND];
+        
+        // Once Dirichlet BCs are added
+        //_named_domains() ...
+        _named_domains["top_bnd"] = Domain("top_bnd");
+
+        /*
+         * Reference model function
+         */
+        // TODO: remove debugging stuff
+        setupInitialModel(parametrization, gmshFem);
+        gmshfem::post::save(_m0[0], _model_known[Support::BLK], "model_known");
+        gmshfem::post::save(_m0[0], _model_unknown[Support::BLK], "model_unknown");
+
+        
+    }
+
+
+    void Configuration::wave_mesh() const
+    {
+        int p0 = factory::addPoint(0., -_H, 0., _h);
+        int p1 = factory::addPoint(_L, -_H, 0., _h);
+        int p2 = factory::addPoint(_L, -_Hsub, 0., _h);
+        int p3 = factory::addPoint(0., -_Hsub, 0., _h);
+        int p4 = factory::addPoint(_L, 0., 0., _h);
+        int p5 = factory::addPoint(0., 0., 0., _h);
+
+        // All horizontal lines defined left to right, vertical bottom to top
+        int lineBottom = factory::addLine(p0, p1);
+        int lineSub = factory::addLine(p3, p2);
+        int lineTop = factory::addLine(p5, p4);
+        int lineLeftSub = factory::addLine(p0, p3);
+        int lineLeftSuper = factory::addLine(p3, p5);
+        int lineRightSub = factory::addLine(p1, p2);
+        int lineRightSuper = factory::addLine(p2, p4);
+
+        int clSub = factory::addCurveLoop({lineBottom, lineRightSub, -lineSub, -lineLeftSub});
+        int clSuper = factory::addCurveLoop({lineSub, lineRightSuper, -lineTop, -lineLeftSuper});
+        
+        int sSub = factory::addPlaneSurface({clSub});
+        int sSuper = factory::addPlaneSurface({clSuper});
+
+        std::vector<int> p_er;
+        // Emitters and receivers
+        for (auto pos: _er_positions) {
+            p_er.push_back(factory::addPoint(pos.first, pos.second, 0.0, _h));
+        }
+
+        factory::synchronize();
+
+        // Embed the e/r. Check if they're in the subsurface or supersurface
+        for (unsigned i = 0; i < p_er.size(); ++i) {
+            bool isSub = _er_positions[i].second < (-_Hsub);
+            gmodel::mesh::embed(0, {p_er[i]}, 2, isSub ? sSub : sSuper);
+
+        }
+
+        gmodel::addPhysicalGroup(2, {sSub}, 21);
+        gmodel::setPhysicalName(2, 21, "subsurface");
+        gmodel::addPhysicalGroup(2, {sSuper}, 22);
+        gmodel::setPhysicalName(2, 22, "supersurface");
+        
+        gmodel::addPhysicalGroup(1, {lineLeftSub, lineBottom, lineRightSub}, 11);
+        gmodel::setPhysicalName(1, 11, "subsurface_bnd");
+        gmodel::addPhysicalGroup(1, {lineLeftSuper, lineTop, lineRightSuper}, 12);
+        gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+
+        gmodel::addPhysicalGroup(1, {lineTop}, 13);
+        gmodel::setPhysicalName(1, 13, "top_bnd");
+
+        for (unsigned p = 0; p < p_er.size(); ++p) {
+            gmodel::addPhysicalGroup(0, {p_er[p]}, 1000+p);
+            gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p));
+        }
+
+    }
+    void Configuration::data_mesh() const
+    {
+        // TODO
+        /*
+        double DAngle = 2. * M_PI /((double) _np);
+        double Ds_e = DAngle * _emitter_skip * _rer;
+        double Ds_r = DAngle * _receiver_skip * _rer;
+        double Le0 = DAngle * _emitter_offset * _rer;
+        double Lr0 = DAngle * _receiver_offset * _rer;
+        double Le = (_ne-1) * Ds_e;
+        double Lr = (_nr-1) * Ds_r;
+
+        int p0 = factory::addPoint(Le0, Lr0, 1.);
+        int p1 = factory::addPoint(Le0+Le, Lr0, 1.);
+        int p2 = factory::addPoint(Le0+Le, Lr0+Lr, 1.);
+        int p3 = factory::addPoint(Le0, Lr0+Lr, 1.);
+
+        int l0 = factory::addLine(p0, p1);
+        int l1 = factory::addLine(p1, p2);
+        int l2 = factory::addLine(p2, p3);
+        int l3 = factory::addLine(p3, p0);
+
+        int cl = factory::addCurveLoop({l0,l1,l2,l3});
+        int s1 = factory::addPlaneSurface({cl});
+
+        factory::mesh::setTransfiniteCurve(l0,_ne);
+        factory::mesh::setTransfiniteCurve(l1,_nr);
+        factory::mesh::setTransfiniteCurve(l2,_ne);
+        factory::mesh::setTransfiniteCurve(l3,_nr);
+        factory::mesh::setTransfiniteSurface(s1,"Left",{p0,p1,p2,p3});
+
+        factory::mesh::setRecombine(2,s1);
+
+        gmodel::addPhysicalGroup(2, {s1}, 1);
+        gmodel::setPhysicalName(2, 1, "data_omega");
+        
+        */
+    }
+
+    void Configuration::setupInitialModel(const ParametrizationInterface* const parametrization, const GmshFem& gmshFem) {
+        /*
+         * Reference model function
+         */
+        _m_super.resize(model_size());
+        _m_sub.resize(model_size());
+        for (unsigned int c = 0; c < model_size(); c++)
+        {
+            ScalarPiecewiseFunction<std::complex<double>> m0;
+            std::string suffix = "c" + std::to_string(c);
+            double Rem_super, Imm_super, Rem_sub, Imm_sub;
+            if (!(
+                    gmshFem.userDefinedParameter(Rem_super, "Re(m_super" + suffix + ")") && gmshFem.userDefinedParameter(Imm_super, "Im(m_super" + suffix + ")") &&
+                    gmshFem.userDefinedParameter(Rem_sub, "Re(m_sub" + suffix + ")") && gmshFem.userDefinedParameter(Imm_sub, "Im(m_sub" + suffix + ")")))
+            {
+                throw common::Exception("A model component " + suffix + " could not be found.");
+            }
+            _m_super[c] = Rem_super + im * Imm_super;
+            _m_sub[c] = Rem_sub + im * Imm_sub;
+            _mc[c] = _m_sub[c];
+
+            m0.addFunction(_m_super[c], _supersurface[Support::BLK] | _supersurface[Support::BND] | _points);
+
+            std::string m0_type;
+            if (!gmshFem.userDefinedParameter(m0_type, "m0_type" + suffix))
+            {
+                throw Exception("Reference model type could not be found.");
+            }
+
+            if (m0_type == "file")
+            {
+                std::string path = "";
+                if (!gmshFem.userDefinedParameter(path, "m0_path" + suffix))
+                {
+                    throw common::Exception("Path to subsurface data could not be found.");
+                }
+                m0.addFunction(bilinearInterpolation(path), _subsurface[Support::BLK] | _subsurface[Support::BND]);
+            }
+            else if (m0_type == "file.pos")
+            {
+                std::string path = "";
+                if (!gmshFem.userDefinedParameter(path, "m0_path" + suffix))
+                {
+                    throw common::Exception("Path to subsurface data could not be found.");
+                }
+
+                gmsh::merge(path + suffix + ".pos");
+                // Fix for weird tags when loading .pos files.
+                std::vector<int> tags;
+                gmsh::view::getTags(tags);
+                auto view = tags.back();
+                ScalarFunction<std::complex<double>> mpos = probeScalarView<std::complex<double>>(view);
+
+                m0.addFunction(mpos, _subsurface[Support::BLK] | _subsurface[Support::BND]);
+            }
+            /*
+            else if (m0_type == "inverse_linear_squared")
+            {
+                double Rea_0, Ima_0;
+                if (!(
+                        gmshFem.userDefinedParameter(Rea_0, "Re(a0" + suffix + ")") && gmshFem.userDefinedParameter(Ima_0, "Im(a0" + suffix + ")")))
+                {
+                    throw common::Exception("Initial model parameter (a0) could not be found.");
+                }
+
+                double Rea_H, Ima_H;
+                if (!(
+                        gmshFem.userDefinedParameter(Rea_H, "Re(aH" + suffix + ")") && gmshFem.userDefinedParameter(Ima_H, "Im(aH" + suffix + ")")))
+                {
+                    throw common::Exception("Initial model parameter (aH) could not be found.");
+                }
+
+                ScalarFunction<std::complex<double>> num = (Rea_0 + im * Ima_0) - ((Rea_H + im * Ima_H) - (Rea_0 + im * Ima_0)) / H() * y<std::complex<double>>();
+                m0.addFunction(1. / pow(num, 2), _subsurface[Support::BLK] | _subsurface[Support::BND]);
+            }
+            else if (m0_type == "linear")
+            {
+                double Rea_0, Ima_0;
+                if (!(
+                        gmshFem.userDefinedParameter(Rea_0, "Re(a0" + suffix + ")") && gmshFem.userDefinedParameter(Ima_0, "Im(a0" + suffix + ")")))
+                {
+                    throw common::Exception("Initial model parameter (a0) could not be found.");
+                }
+
+                double Rea_H, Ima_H;
+                if (!(
+                        gmshFem.userDefinedParameter(Rea_H, "Re(aH" + suffix + ")") && gmshFem.userDefinedParameter(Ima_H, "Im(aH" + suffix + ")")))
+                {
+                    throw common::Exception("Initial model parameter (aH) could not be found.");
+                }
+
+                ScalarFunction<std::complex<double>> lin = (Rea_0 + im * Ima_0) - ((Rea_H + im * Ima_H) - (Rea_0 + im * Ima_0)) / H() * (y<std::complex<double>>() + _ymin);
+                m0.addFunction(lin, _subsurface[Support::BLK] | _subsurface[Support::BND]);
+            }*/
+            else if (m0_type == "constant")
+            {
+                m0.addFunction(_m_sub[c], _subsurface[Support::BLK] | _subsurface[Support::BND]);
+            }
+            else
+            {
+                throw common::Exception("Initial model type (" + suffix + ") " + m0_type + " is unknown.");
+            }
+            _m0.push_back(m0);
+        }
+    }
+
+    std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
+    {
+        throw Exception("Unimplemented Configuration::data_coordinate_to_index");
+        /*
+        if( (!_receiver_on_emitter) )
+        {
+            throw Exception("Data space when an emitter is not also a receiver is not implemented yet");
+        }
+        std::array<unsigned int,2> index;
+
+        double DAngle = 2. * M_PI /((double) _np);
+        double Ds_e = DAngle * _emitter_skip * _rer;
+        double Ds_r = DAngle * _receiver_skip * _rer;
+        double Le0 = DAngle * _emitter_offset * _rer;
+        double Lr0 = DAngle * _receiver_offset * _rer;
+        index[0] = std::round( (xs-Le0) / Ds_e );
+        index[1] = std::round( (xr-Lr0) / Ds_r );
+
+        return index;
+        */
+    }
+
+    std::array<double,2> Configuration::index_to_data_coordinate(unsigned int s, unsigned int r) const
+    {
+        throw Exception("Unimplemented Configuration::index_to_data_coordinate");
+
+        /*
+        if( (!_receiver_on_emitter) )
+        {
+            throw Exception("Data space when an emitter is not also a receiver is not implemented yet");
+        }
+
+        std::array<double,2> index;
+
+        double DAngle = 2. * M_PI /((double) _np);
+        double Ds_e = DAngle * _emitter_skip * _rer;
+        double Ds_r = DAngle * _receiver_skip * _rer;
+        double Le0 = DAngle * _emitter_offset * _rer;
+        double Lr0 = DAngle * _receiver_offset * _rer;
+        index[0] = Le0 + ((double)s) * Ds_e;
+        index[1] = Lr0 + ((double)r) * Ds_r;
+
+        return index;*/
+    }
+
+    double Configuration::data_area() const
+    {
+        throw Exception("Unimplemented Configuration::data_area");
+        /*
+        double DAngle = 2. * M_PI /((double) _np);
+        double Ds_e = DAngle * _emitter_skip * _rer;
+        double Ds_r = DAngle * _receiver_skip * _rer;
+        double Le = (_ne-1) * Ds_e;
+        double Lr = (_nr-1) * Ds_r;
+
+        return Le * Lr;
+        */
+    };
+
+    double Configuration::datapoint_area() const
+    {
+        throw Exception("Unimplemented Configuration::datapoint_area");
+        /*
+                double DAngle = 2. * M_PI /((double) _np);
+                double Ds_e = DAngle * _emitter_skip * _rer;
+                double Ds_r = DAngle * _receiver_skip * _rer;
+
+                return Ds_e * Ds_r;*/
+    };
+
+    bool Configuration::data_coordinate_isValid(double xs,double xr) const
+    {
+        throw Exception("Unimplemented Configuration::data_coordinate_isValid");
+
+        /*
+        double DAngle = 2. * M_PI /((double) _np);
+        double Ds_e = DAngle * _emitter_skip * _rer;
+        double Ds_r = DAngle * _receiver_skip * _rer;
+        double Le0 = DAngle * _emitter_offset * _rer;
+        double Lr0 = DAngle * _receiver_offset * _rer;
+        double Le = (_ne-1) * Ds_e;
+        double Lr = (_nr-1) * Ds_r;
+
+        double eps = 1e-8;
+        gmsh::option::getNumber("Geometry.MatchMeshTolerance",eps);
+        if( (Le0-eps <= xs) && (xs <= Le0+Le+eps) && (Lr0-eps <= xr) && (xr <= Lr0+Lr+eps) )
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+
+        */
+    }
+
+    template<Physic T_Physic>
+    ModelMonoFunction green0_preconditioner(const Data<T_Physic>& dd, const WaveMultiField<T_Physic>& g, const Configuration* const config)
+    {
+        return autocorrelate<T_Physic>(g,config->emitter_idx())*autocorrelate<T_Physic>(g,config->receiver_idx());
+    }
+
+    template ModelMonoFunction green0_preconditioner<Physic::acoustic>(const Data<Physic::acoustic>&, const WaveMultiField<Physic::acoustic>&, const Configuration* const);
+    template ModelMonoFunction green0_preconditioner<Physic::electromagnetic>(const Data<Physic::electromagnetic>&, const WaveMultiField<Physic::electromagnetic>&, const Configuration* const);
+    template ModelMonoFunction green0_preconditioner<Physic::elastodynamic>(const Data<Physic::elastodynamic>&, const WaveMultiField<Physic::elastodynamic>&, const Configuration* const);
+
+} // namespace circular_acquisition

specific/configuration/flexible_acquisition.h

+#ifndef H_CONFIGURATION_FLEXIBLE_ACQUISITION
+#define H_CONFIGURATION_FLEXIBLE_ACQUISITION
+
+//GmshFEM Library
+#include "GmshFem.h"
+//#include "Function.h"
+
+//GmshFWI Library
+#include "../../common/configuration.h"
+#include "../../common/wave/element.h"
+#include "inclusion/inclusion.h"
+
+//Forward declaration
+template<Physic T_Physic>
+class Data;
+
+namespace flexible_acquisition
+{
+
+    class Configuration final : public ConfigurationInterface
+    {
+    private:
+    
+        /**
+         * Emetters and receivers are stored together, and for each shot, one point is set as emitter, 
+         * and some others (potentially including the emitter) are enabled as receivers
+        */
+        std::vector<std::pair<double, double>> _er_positions;
+
+        struct ShotConfig {
+            unsigned int emetterIdx;
+            std::vector<unsigned int> enabledReceivers;
+        };
+
+        std::vector<ShotConfig> _shotsConfigs;
+
+        /* Domain is a rectangle. "x" goes from 0 to L and "y" from -H to 0. Known part is y between 0 and -Hsub*/
+        double _H;
+        double _L;
+        double _Hsub;
+
+        // Mesh resolution
+        double _h;
+
+        // TODO
+        std::vector<unsigned int> _emitter_idx;
+        std::vector<unsigned int> _receiver_idx;
+
+        bool _areFilled;
+
+        // Subdomains: all in the parent class
+    /*
+        double _rer;
+        double _rext;
+
+        unsigned int _ne;
+        unsigned int _nr;
+        unsigned int _emitter_offset;
+        unsigned int _receiver_offset;
+        unsigned int _emitter_skip;
+        unsigned int _receiver_skip;
+        bool _receiver_on_emitter;
+
+        unsigned int _ni;
+        
+        std::vector<const InclusionInterface*> _inclusion_geo;
+
+
+        UnknownRegion _unknown_region;
+        std::array<gmshfem::domain::Domain,2> _background;
+        std::array<gmshfem::domain::Domain,2> _inclusions;
+        std::vector<std::array<gmshfem::domain::Domain,2>> _inclusion;
+
+
+
+        std::vector<std::complex<double>> _mb;
+        std::vector<std::vector<std::complex<double>>> _mi;
+*/
+
+        std::array<gmshfem::domain::Domain, 2> _subsurface;
+        std::array<gmshfem::domain::Domain, 2> _supersurface;
+
+        std::vector<std::complex<double>> _m_super;
+        std::vector<std::complex<double>> _m_sub;
+
+        virtual void wave_mesh() const;
+        virtual void data_mesh() const;
+
+
+        void setupInitialModel(const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem);
+    public:
+        Configuration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem);
+
+        const std::vector<unsigned int>& emitter_idx() const {return _emitter_idx;};
+        const std::vector<unsigned int>& receiver_idx() const {return _receiver_idx;};
+
+        bool areFilled() const {return _areFilled;};
+        unsigned int ner() const {return _np;};
+
+//        virtual void mesh() const;
+
+        virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+        virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
+        virtual bool data_coordinate_isValid(double xs,double xr) const;
+
+        virtual double area() const {return _H * _L;};
+        virtual double data_area() const override;
+        virtual double datapoint_area() const;
+        virtual double array() const {return 2. * (_L + _H);};
+        virtual double depth() const {return _H;};
+
+        virtual std::string wave_gmodel() const {return _name;};
+        virtual std::string model_gmodel() const {return _name;};
+        virtual std::string data_gmodel() const {return _name;};
+    };
+
+    template<Physic T_Physic>
+    ModelMonoFunction green0_preconditioner(const Data<T_Physic>& dd, const WaveMultiField<T_Physic>& g, const Configuration* const config);
+} // namespace soil
+
+#endif // H_CONFIGURATION_FLEXIBLE_ACQUISITION

specific/configuration/newConfiguration.h

@@ -16,6 +16,7 @@
 #include "circular_acquisition.h"
 #include "line_acquisition.h"
 #include "rectangular_acquisition.h"
+#include "flexible_acquisition.h"
 
 std::unique_ptr<ConfigurationInterface> newConfiguration(std::string name, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem)
 {
@@ -30,6 +31,7 @@ std::unique_ptr<ConfigurationInterface> newConfiguration(std::string name, const
     else if(configuration=="circular_acquisition"){return std::make_unique< circular_acquisition::Configuration>(name, parametrization, gmshFem);}
     else if(configuration=="line_acquisition"){return std::make_unique< line_acquisition::Configuration>(name, parametrization, gmshFem);}
     else if(configuration=="rectangular_acquisition"){return std::make_unique< rectangular_acquisition::Configuration>(name, parametrization, gmshFem);}
+    else if(configuration=="flexible_acquisition"){return std::make_unique< flexible_acquisition::Configuration>(name, parametrization, gmshFem);}
     else
     {
         throw gmshfem::common::Exception("Configuration" + configuration + " is not valid.");

specific/configuration/yaml_interface.cpp

+#include "GmshFem.h"
+#include "yaml_interface.h"
+
+ShotsConfigurationYAML::ShotsConfigurationYAML(const std::string &path)
+{
+#ifndef ENABLE_YAML_CPP
+    throw gmshfem::common::Exception("Cannot use Flexible Acquisition without YAML-CPP.");
+#else
+    YAML::Node config = YAML::LoadFile(path);
+
+    for (const auto &coord : config["coordinates"])
+    {
+        double x = coord[0].as<double>();
+        double y = coord[1].as<double>();
+        _points.emplace_back(x, y);
+    }
+
+    for (const auto &shot : config["shots"])
+    {
+        std::vector<unsigned> emitters, receivers;
+        for (const auto &emitter : shot["emitters"])
+        {
+            emitters.push_back(emitter.as<int>());
+        }
+        for (const auto &receiver : shot["receivers"])
+        {
+            receivers.push_back(receiver.as<int>());
+        }
+
+        _emitters.push_back(std::move(emitters));
+        _receivers.push_back(std::move(receivers));
+    }
+#endif
+}
\ No newline at end of file

specific/configuration/yaml_interface.h

+#ifndef H_YAML_INTERFACE_SHOTS
+#define H_YAML_INTERFACE_SHOTS
+
+#include <vector>
+#include <string>
+#include <utility>
+#include "yaml-cpp/yaml.h"
+// TODO: refactor (abstract + impl)
+
+class ShotsConfigurationInterface //(abstract)
+{
+    public:
+        using Coordinates = std::pair<double, double>;
+
+        /**
+         * Describe the number of shots
+         */
+        virtual unsigned numShots() const = 0;
+
+        /**
+         * Describe the number of points
+         */
+        virtual unsigned numPoints() const = 0;
+
+        /**
+         * emitters(s) is the list of enabled emitters in shot "s".
+         * Array must be of size numShots()
+         */
+        virtual const std::vector<unsigned> &emitters(unsigned shot) const = 0;
+
+        /**
+         * receivers(s) is the list of enabled receivers in shot "s".
+         * Array must be of size numShots()
+         */
+        virtual const std::vector<unsigned> &receivers(unsigned shot) const = 0;
+
+        /**
+         * List of all points
+         * Array must be of size numPoints()
+         */
+        virtual const std::vector<Coordinates> &points() const = 0;
+};
+
+class ShotsConfigurationYAML : public ShotsConfigurationInterface {
+    private:
+        std::vector<std::vector<unsigned>> _emitters, _receivers;
+        std::vector<ShotsConfigurationInterface::Coordinates> _points;
+
+    public:
+
+        ShotsConfigurationYAML(const std::string& path);
+
+        virtual unsigned numShots() const override {
+            return _emitters.size();
+        }
+
+        virtual unsigned numPoints() const override {
+            return _points.size();
+        }
+
+        virtual const std::vector<unsigned> &emitters(unsigned shot) const override {
+            return _emitters.at(shot);
+        }
+
+        virtual const std::vector<unsigned> &receivers(unsigned shot) const override {
+            return _receivers.at(shot);
+        }
+
+        virtual const std::vector<Coordinates> &points() const override {
+            return _points;
+        }
+
+        
+
+};
+
+#endif
\ No newline at end of file

specific/wave/equation/simplewave.cpp

@@ -81,6 +81,12 @@ namespace simplewave
         else {
             msg::warning << "SimpleWave without specified top_bc defaults to Absorbing." << msg::endl;
         }
+
+        _useSecondOrderABC = false;
+        int intArg;
+        gmshFem.userDefinedParameter(intArg, "useSecondOrderABC");
+        _useSecondOrderABC = intArg != 0;
+        msg::warning << "Using Second Order ABC : " << _useSecondOrderABC << msg::endl;
     }
 
     bool Equation::modelIsObsolete()
@@ -114,6 +120,8 @@ namespace simplewave
     void Equation::setLHS(const ModelStateEvaluator& ms)
     {
         auto top = _config->named_domains().find("top_bnd");
+        
+
         if (_topBC == "Dirichlet" || _topBC == "Neumann")
         {
             if (top == _config->named_domains().end())
@@ -138,10 +146,13 @@ namespace simplewave
         if(_realImpedance)
         {
             _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+            if (_useSecondOrderABC) _formulation.integral(1. * im / (2.0 * _pulsation * sqrt(real(m0))) * grad(dof(_v)), grad(tf(_v)), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+
         }
         else
         {
             _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+           if (_useSecondOrderABC)   _formulation.integral(1. * im / (2.0 * _pulsation * sqrt(m0)) * grad(dof(_v)), grad(tf(_v)), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
         }
         if(!_config->model_unknown(Support::BLK).isEmpty())
         {
@@ -152,10 +163,12 @@ namespace simplewave
             if(_realImpedance)
             {
                 _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) *  dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+                if (_useSecondOrderABC) _formulation.integral(1. * im / (2. * _pulsation * sqrt(real(m))) *  grad(dof(_v)), grad(tf(_v)), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
             }
             else
             {
                  _formulation.integral(-1. * im * _pulsation * sqrt(m) *  dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+                if (_useSecondOrderABC) _formulation.integral(1. * im / (2. * _pulsation * sqrt(m)) *  grad(dof(_v)), grad(tf(_v)), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
             }
         }
     }

specific/wave/equation/simplewave.h

@@ -22,6 +22,7 @@ namespace simplewave
         double _const_preconditioner;
         using EquationInterface<Physic::acoustic>::_boundary;
         std::string _topBC;
+        bool _useSecondOrderABC;
     public:
         Equation(unsigned int f_idx,double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");
 

tests/test.cpp

@@ -7,6 +7,10 @@
 #include <gmshfem/Formulation.h>
 #include <gmshfem/GmshFem.h>
 #include "../common/wave/element.h"
+#include "../common/data/element.h"
+#include "../specific/data/element.h"
+
+#include "dummyAcquisition.h"
 
 #define INIT_GMSHFEM()                  \
     namespace geo = gmsh::model::geo;   \