Select Git revision
inclusion.cpp
inclusion.cpp 19.53 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 "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(std::string name, const GmshFem& gmshFem) : ConfigurationInterface(name, 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 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;
}
void Configuration::mesh() const
{
msg::print << "Generate meshes" << msg::endl;
gmsh::option::setNumber("General.Terminal", 1);
std::string name = _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");
}
}
factory::removeAllDuplicates();
}
void Configuration::data_mesh() const
{
if(nxe() != 0)
{
double ors = (_H - _He) / 2.;
double orr = (_H - _Hr) / 2.;
//Vertical array
int p00 = factory::addPoint(ors, orr, 0.);
int pHe0 = factory::addPoint(ors+_He,orr, 0.);
int pHeHr = factory::addPoint(ors+_He, orr+_Hr, 0.);
int p0Hr = factory::addPoint(ors, orr+_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");
factory::removeAllDuplicates();
}
if(nye() != 0)
{
double ors = _H + (_L - _Le) / 2.;
double orr = _H + (_L - _Lr) / 2.; //Horizontal array
int pHeHr = factory::addPoint(ors, orr, 0.);
int pLe0 = factory::addPoint(ors+_Le, orr, 0.);
int pLeLr = factory::addPoint(ors+_Le, orr+_Lr, 0.);
int p0Lr = factory::addPoint(ors+0., orr+_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");
}
}
std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
{
std::array<unsigned int,2> index;
if(xs<_H && xr < _H)
{
double xe0 = (_H-_He)/2.;
double stepe = _He/(nxe()-1);
double xr0 = (_H-_Hr)/2.;
double stepr = _Hr/(nxr()-1);
index[0] = std::round((xs-xe0)/stepe);
index[1] = std::round((xr-xr0)/stepr);
}
else if(xs>=_H && xr>=_H)
{
double ye0 = _H + (_L-_Le)/2.;
double stepe = _Le/(nye()-1);
double yr0 = _H + (_L-_Lr)/2.;
double stepr = _Lr/(nyr()-1);
index[0] = std::round((xs-ye0)/stepe);
index[1] = std::round((xr-yr0)/stepr);
}
else
{
throw Exception("Invalid data coordinate.");
}
return index;
}
std::array<double,2> Configuration::index_to_data_coordinate(unsigned int s, unsigned int r) const
{
std::array<double,2> index;
if(s < nxe() && r < nxr())
{
double xe0 = (_H-_He)/2.;
double stepe = _He/(nxe()-1);
double xr0 = (_H-_Hr)/2.;
double stepr = _Hr/(nxr()-1);
index[0] = xe0 + s*stepe;
index[1] = xr0 + r*stepr;
}
else if( s>=nxe() && r>=nxr())
{
double xe0 = _H + (_L-_Le)/2.; double stepe = _Le/(nye()-1);
double xr0 = _H + (_L-_Lr)/2.;
double stepr = _Lr/(nyr()-1);
index[0] = xe0 + (s-nxe())*stepe;
index[1] = xr0 + (r-nxr())*stepr;
}
else
{
throw Exception("Invalid index.");
}
return index;
}
template<Physic T_Physic>
ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
{
ModelFunction termX;
if( !config->receiversAreEmittersX() )
{
termX = autocorrelate<T_Physic>(g,config->emitter_idx_X())*autocorrelate<T_Physic>(g,config->receiver_idx_X());
}
else
{
termX = pow(autocorrelate<T_Physic>(g,config->emitter_idx_X()),2);
}
ModelFunction termY;
if( !config->receiversAreEmittersY() )
{
termY = autocorrelate<T_Physic>(g,config->emitter_idx_Y())*autocorrelate<T_Physic>(g,config->receiver_idx_Y());
}
else
{
termY = pow(autocorrelate<T_Physic>(g,config->emitter_idx_Y()),2);
}
return termX + termY;
}
template ModelFunction green0_preconditioner<Physic::acoustic>(const WaveMultiField<Physic::acoustic>&, const Configuration* const);
} // namespace soil