diff --git a/specific/configuration/flexible_acquisition.cpp b/specific/configuration/flexible_acquisition.cpp
index b16a83972e6a2c533218213d7ec8364a06748a8b..21911ec32085ad7dac30f50f31de8b8ae6f2166c 100644
--- a/specific/configuration/flexible_acquisition.cpp
+++ b/specific/configuration/flexible_acquisition.cpp
@@ -67,21 +67,26 @@ namespace flexible_acquisition
// 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");
+ // TODO: define as supersurface and subsurface
std::string unknownRegion = "subsurface";
gmshFem.userDefinedParameter(unknownRegion, "unknown");
if (unknownRegion == "subsurface")
{
- _model_known[Support::BLK] = Domain("supersurface");
- _model_known[Support::BND] = Domain("supersurface_bnd");
+ _model_known[Support::BLK] = _supersurface[Support::BLK];
+ _model_known[Support::BND] = _supersurface[Support::BND];
- _model_unknown[Support::BLK] = Domain("subsurface");
- _model_unknown[Support::BND] = Domain("subsurface_bnd");
+ _model_unknown[Support::BLK] = _subsurface[Support::BLK];
+ _model_unknown[Support::BND] = _subsurface[Support::BND];
}
else if (unknownRegion == "none")
{
- _model_known[Support::BLK] = Domain("supersurface") | Domain("subsurface");
- _model_known[Support::BND] = Domain("supersurface_bnd") | Domain("subsurface_bnd");
+ _model_known[Support::BLK] = _supersurface[Support::BLK] | _subsurface[Support::BLK];
+ _model_known[Support::BND] = _supersurface[Support::BND] | _subsurface[Support::BND];
}
else {
throw Exception("Invalid unknown region type: " + unknownRegion);
@@ -96,6 +101,10 @@ namespace flexible_acquisition
/*
* 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");
/*
_m_super.resize(model_size());
@@ -436,6 +445,108 @@ namespace flexible_acquisition
*/
}
+ 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");
+ 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);
+ }
+ }
+
std::array<unsigned int,2> Configuration::data_coordinate_to_index(double xs, double xr) const
{
throw Exception("Unimplemented Configuration::data_coordinate_to_index");
diff --git a/specific/configuration/flexible_acquisition.h b/specific/configuration/flexible_acquisition.h
index b8c40af0747f84092c65397cda5bc16ebea0a8cb..608b9e903ccfcd602714a9e29a12b0e78e874de7 100644
--- a/specific/configuration/flexible_acquisition.h
+++ b/specific/configuration/flexible_acquisition.h
@@ -77,10 +77,17 @@ namespace flexible_acquisition
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);