Skip to content
Snippets Groups Projects
Commit 5d810ba5 authored by Boris Martin's avatar Boris Martin
Browse files

cleanup

parent 12a8ff79
No related branches found
No related tags found
1 merge request!6Draft: "Flexible acquisition", a configuration with arbitrary sources and receivers read from a YAML file.
...@@ -116,9 +116,6 @@ namespace flexible_acquisition ...@@ -116,9 +116,6 @@ namespace flexible_acquisition
_wave_omega[Support::BLK] = _model_known[Support::BLK] | _model_unknown[Support::BLK]; _wave_omega[Support::BLK] = _model_known[Support::BLK] | _model_unknown[Support::BLK];
_wave_omega[Support::BND] = _model_known[Support::BND] | _model_unknown[Support::BND]; _wave_omega[Support::BND] = _model_known[Support::BND] | _model_unknown[Support::BND];
// Once Dirichlet BCs are added // Once Dirichlet BCs are added
//_named_domains() ... //_named_domains() ...
...@@ -130,181 +127,7 @@ namespace flexible_acquisition ...@@ -130,181 +127,7 @@ namespace flexible_acquisition
gmshfem::post::save(_m0[0], _model_known[Support::BLK], "model_known"); gmshfem::post::save(_m0[0], _model_known[Support::BLK], "model_known");
gmshfem::post::save(_m0[0], _model_unknown[Support::BLK], "model_unknown"); gmshfem::post::save(_m0[0], _model_unknown[Support::BLK], "model_unknown");
/*
_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");
ScalarFunction<std::complex<double>> mpos = probeScalarView<std::complex<double>>(c);
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);
}
*/
// MODEL ON SUBSURFACE, SURFACE
/*
_mb.resize(model_size());
_mi.resize(model_size());
for (unsigned int c = 0; c < model_size(); c++)
{
std::string suffix = "c"+std::to_string(c);
double Remb, Immb;
if
(!
(
gmshFem.userDefinedParameter(Remb, "Re(mb"+suffix+")") &&
gmshFem.userDefinedParameter(Immb, "Im(mb"+suffix+")")
)
)
{
throw common::Exception("Background (reference) model parameter (component "+suffix+" ) could not be found.");
}
else
{
_mb[c] = Remb + im * Immb;
_mc[c] = _mb[c];
}
ScalarPiecewiseFunction< std::complex< double > > m0;
std::string m0_type = "constant";
if(!gmshFem.userDefinedParameter(m0_type, "m0_type"+suffix))
{
msg::warning << "Reference model type could not be found. Constant background is used (default)." << msg::endl;
}
if(m0_type=="file")
{
std::string path = "";
if(!gmshFem.userDefinedParameter(path, "m0_path"+suffix))
{
throw common::Exception("Path to circular data could not be found.");
}
m0.addFunction(bilinearInterpolation(path),_background[Support::BLK] | _background[Support::BND] );
}
else if(m0_type=="file.pos")
{
std::string path = "";
if(!gmshFem.userDefinedParameter(path, "m0_path"+suffix))
{
throw common::Exception("Path to model file could not be found.");
}
gmsh::merge(path+suffix+".pos");
ScalarFunction<std::complex<double>> mpos = probeScalarView<std::complex<double>>(c);
m0.addFunction(mpos,_background[Support::BLK] | _background[Support::BND]);
}
else if(m0_type=="constant")
{
m0.addFunction(_mb[c],_background[Support::BLK] | _background[Support::BND]);
}
else
{
throw common::Exception("Initial model type ("+ suffix + ") " + m0_type +" is unknown.");
}
_mi[c].resize(_ni);
for (unsigned int i = 0; i < _ni; i++)
{
double Remi=0., Immi=0.;
if
(
!(gmshFem.userDefinedParameter(Remi, "Re(mi"+std::to_string(i)+suffix+")") &&
gmshFem.userDefinedParameter(Immi, "Im(mi"+std::to_string(i)+suffix+")"))
)
{
throw Exception("Inclusion model parameter could not be found.");
}
_mi[c][i] = Remi + im * Immi;
m0.addFunction(_mi[c][i],_inclusion[i][Support::BLK] | _inclusion[i][Support::BND]);
}
_m0.push_back(m0);
}*/
} }
...@@ -347,9 +170,6 @@ namespace flexible_acquisition ...@@ -347,9 +170,6 @@ namespace flexible_acquisition
} }
//gmodel::addPhysicalGroup(2, {}, 1000+p);
//gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p));
gmodel::addPhysicalGroup(2, {sSub}, 21); gmodel::addPhysicalGroup(2, {sSub}, 21);
gmodel::setPhysicalName(2, 21, "subsurface"); gmodel::setPhysicalName(2, 21, "subsurface");
gmodel::addPhysicalGroup(2, {sSuper}, 22); gmodel::addPhysicalGroup(2, {sSuper}, 22);
...@@ -364,71 +184,7 @@ namespace flexible_acquisition ...@@ -364,71 +184,7 @@ namespace flexible_acquisition
gmodel::addPhysicalGroup(0, {p_er[p]}, 1000+p); gmodel::addPhysicalGroup(0, {p_er[p]}, 1000+p);
gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p)); gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p));
} }
/*
//Vertical emitters
std::vector<int> pnt(_np);
double DAngle = 2. * M_PI /((double) _np);
for (unsigned int p = 0; p < _np; p++)
{
double angle = ((double)p) * DAngle;
pnt[p] = factory::addPoint(_rer*std::cos(angle), _rer*std::sin(angle), 0., _h);
}
int sb;
std::vector<int> si;
std::vector<std::vector<int>> li(_ni);
std::vector<int> li_tot;
if(_ni!=0)
{
std::vector<int> cli(_ni,0);
if(_areFilled){si.resize(_ni);}
for (unsigned int i = 0; i < _ni; i++)
{
li[i] = _inclusion_geo[i]->addInclusion();
li_tot.insert(li_tot.end(), li[i].begin(), li[i].end());
cli[i] = factory::addCurveLoop(li[i]);
if(_areFilled)
{
si[i] = factory::addPlaneSurface({cli[i]});
}
}
std::vector<int> clbcli = cli;
clbcli.insert(clbcli.begin(), clb);
sb = factory::addPlaneSurface(clbcli);
}
else
{
sb = factory::addPlaneSurface({clb});
}
factory::synchronize();
gmodel::mesh::embed(0, pnt, 2, sb);
for (unsigned int p = 0; p <_np; p++)
{
gmodel::addPhysicalGroup(0, {pnt[p]}, 1000+p);
gmodel::setPhysicalName(0, 1000+p, "emitter_receiver_"+std::to_string(p));
}
gmodel::addPhysicalGroup(1, {lb1,lb2,lb3,lb4}, 12);
gmodel::setPhysicalName(1, 12, "background_bnd");
gmodel::addPhysicalGroup(2, {sb}, 21);
gmodel::setPhysicalName(2, 21, "background_vol");
for (unsigned int i = 0; i < _ni; i++)
{
if(_areFilled)
{
gmodel::addPhysicalGroup(2, {si[i]}, 1000+i);
gmodel::setPhysicalName(2, 1000+i, "inclusion_vol"+std::to_string(i));
}
else
{
gmodel::addPhysicalGroup(1, li[i], 2000+i);
gmodel::setPhysicalName(1, 2000+i, "inclusion_bnd"+std::to_string(i));
}
}*/
} }
void Configuration::data_mesh() const void Configuration::data_mesh() const
{ {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment