Skip to content
Snippets Groups Projects
Select Git revision
  • b4e98b6b8315d0832efa69fe907661aa70f46304
  • master default protected
  • NewDistributionGmshFWI
  • parametrizationSimpleWave
  • tuto_obstacle
  • everything
  • cleanup_configuuration_mesh
  • fix
  • source_estimation
  • unique_ptr
  • SobolevDirectionalFilter
  • OT
  • newPhysics
  • SimultaneousFrequency
  • SobolevDistance
  • BonesImaging
  • MultiParameter
  • UpdateAntho
  • v2.0
  • v1.0
20 results

flexible_acquisition.cpp

Blame
    • flexible_acquisition.cpp 17.06 KiB 
    //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