diff --git a/common/configuration.h b/common/configuration.h
index 8818d6b341c14aebcc1dee9fb91e28cac8b91a17..fb346bbbb019bee519b447c80d5c18677987d12d 100644
--- a/common/configuration.h
+++ b/common/configuration.h
@@ -6,6 +6,7 @@
#include "Domain.h"
//Standard Library
#include <array>
+#include <map>
//GmshFWI Library
#include "model/element.h"
#include "model/parametrization.h"
@@ -21,6 +22,7 @@ protected:
std::array<gmshfem::domain::Domain,2> _wave_omega;
std::array<gmshfem::domain::Domain,2> _model_known;
std::array<gmshfem::domain::Domain,2> _model_unknown;
+ std::map<std::string, gmshfem::domain::Domain> _named_domains;
gmshfem::domain::Domain _data_omega;
unsigned int _np;
@@ -74,6 +76,7 @@ public:
gmshfem::domain::Domain model_evaldomain() const {return _model_unknown[Support::BLK] | _model_unknown[Support::BND];};
gmshfem::domain::Domain wave_evaldomain() const {return _wave_omega[Support::BLK] | _wave_omega[Support::BND] | _points;};
gmshfem::domain::Domain data_evaldomain() const {return _data_omega;};
+ const std::map<std::string, gmshfem::domain::Domain>& named_domains() const {return _named_domains;}
gmshfem::domain::Domain emitter(unsigned int shot, unsigned int em) const;
gmshfem::domain::Domain receiver(unsigned int shot, unsigned int rec) const;
diff --git a/specific/configuration/surface_acquisition.cpp b/specific/configuration/surface_acquisition.cpp
index 0a81f151044eba1020d48131d49c5578c27c307e..37ed89dcaf96ef81a2eb274393bdaccbc744a1a4 100644
--- a/specific/configuration/surface_acquisition.cpp
+++ b/specific/configuration/surface_acquisition.cpp
@@ -113,6 +113,8 @@ namespace surface_acquisition
_supersurface[Support::BLK] = Domain(2, 22);
_supersurface[Support::BND] = Domain(1, 12);
+ _named_domains["top_bnd"] = Domain("top_bnd");
+
_wave_omega[Support::BLK] = _subsurface[Support::BLK] | _supersurface[Support::BLK];
_wave_omega[Support::BND] = _subsurface[Support::BND] | _supersurface[Support::BND];
@@ -341,6 +343,7 @@ namespace surface_acquisition
lymin_reverse[n] = -lymin[n];
}
+ // Nomenclature: line from one y to another, left or right. All lines go up
int lHyminl = factory::addLine(pHl,pyminl);
int lymin0l = factory::addLine(pyminl,p0l);
int l0ymaxl = factory::addLine(p0l,pymaxl);
@@ -430,6 +433,10 @@ namespace surface_acquisition
gmodel::addPhysicalGroup(1, lsuper, 12);
gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+ // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+ gmodel::addPhysicalGroup(1, lymax, 13);
+ gmodel::setPhysicalName(1, 13, "top_bnd");
+
gmodel::addPhysicalGroup(2, {s_sub}, 21);
gmodel::setPhysicalName(2, 21, "subsurface_vol");
gmodel::addPhysicalGroup(2, {s_supermin,s_supermax}, 22);
@@ -507,6 +514,10 @@ namespace surface_acquisition
gmodel::addPhysicalGroup(1, lsuper, 12);
gmodel::setPhysicalName(1, 12, "supersurface_bnd");
+ // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+ gmodel::addPhysicalGroup(1, lymax, 13);
+ gmodel::setPhysicalName(1, 13, "top_bnd");
+
gmodel::addPhysicalGroup(2, {s_sub}, 21);
gmodel::setPhysicalName(2, 21, "subsurface_vol");
gmodel::addPhysicalGroup(2, {s_supermax}, 22);
diff --git a/specific/wave/equation/simplewave.cpp b/specific/wave/equation/simplewave.cpp
index 499781c96d0d767370a5f4d323e75ba678e11e92..b2fe3c972cdd03307245449e0c39aa6cfe79019c 100644
--- a/specific/wave/equation/simplewave.cpp
+++ b/specific/wave/equation/simplewave.cpp
@@ -72,6 +72,15 @@ namespace simplewave
{
throw Exception("Impossible to compute boundary terms with real impedance: sensitivities are wrong.");
}
+
+ if (gmshFem.userDefinedParameter(_topBC, "top_bc")) {
+ if (_topBC != "Absorbing" && _topBC != "Dirichlet" && _topBC != "Neumann") {
+ throw Exception("Invalid value of parameter top_bc: " + _topBC);
+ }
+ }
+ else {
+ msg::warning << "SimpleWave without specified top_bc defaults to Absorbing." << msg::endl;
+ }
}
bool Equation::modelIsObsolete()
@@ -104,16 +113,35 @@ 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())
+ throw Exception("Asked for specific BCs on top but this entity isn't defined.");
+ }
+
+ if (_topBC == "Dirichlet")
+ {
+ _v.addConstraint(top->second, 0.0);
+ }
+
_formulation.integral(-1.*grad(dof(_v)), grad(tf(_v)), _config->wave_omega(Support::BLK),integrationType(2*_w_degree));
_formulation.integral(_pulsation * _pulsation * m0 * dof(_v), tf(_v), _config->model_known(Support::BLK), integrationType(_integrationDegreeBlk));
+
+ auto absorbingPartKnown = _config->model_known(Support::BND);
+ auto absorbingPartUnknown = _config->model_unknown(Support::BND);
+ if (_topBC == "Neumann") {
+ absorbingPartKnown &= ~(top->second);
+ absorbingPartUnknown &= ~(top->second);
+ }
if(_realImpedance)
{
- _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
else
{
- _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
if(!_config->model_unknown(Support::BLK).isEmpty())
{
@@ -123,11 +151,11 @@ namespace simplewave
{
if(_realImpedance)
{
- _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) * dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
else
{
- _formulation.integral(-1. * im * _pulsation * sqrt(m) * dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+ _formulation.integral(-1. * im * _pulsation * sqrt(m) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
}
}
diff --git a/specific/wave/equation/simplewave.h b/specific/wave/equation/simplewave.h
index 0acee2dd67c44b19b7f41fa51d4fa0e51c731e78..2cbe35d648b5d907c4ebe90f62096b9ffa384909 100644
--- a/specific/wave/equation/simplewave.h
+++ b/specific/wave/equation/simplewave.h
@@ -21,6 +21,7 @@ namespace simplewave
bool _realImpedance;
double _const_preconditioner;
using EquationInterface<Physic::acoustic>::_boundary;
+ std::string _topBC;
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 = "");