diff --git a/common/configuration.h b/common/configuration.h index fb346bbbb019bee519b447c80d5c18677987d12d..fe5d0986a3ada8aa6b95c6927c1a5272380127e7 100644 --- a/common/configuration.h +++ b/common/configuration.h @@ -51,6 +51,12 @@ public: gmshFem.userDefinedParameter(remesh, "remesh"); _remesh = ((bool) remesh); }; + + ConfigurationInterface() : _name("dummy"), _parametrization(nullptr), _mc(1,0.) + { + _remesh = true; + }; + virtual ~ConfigurationInterface() {delete _parametrization;}; virtual void mesh() const; diff --git a/common/data/element.h b/common/data/element.h index d200c18c42b8b31107877d8f86ff0d3aaa3a902f..0089970b201531a573cbf936b5463a1336c8af47 100644 --- a/common/data/element.h +++ b/common/data/element.h @@ -182,4 +182,39 @@ Data<T_Physic> arg(const Data<T_Physic>& d) return out; }; +/** + * @brief Returns a vector ([f][s]) that tells for each shot what coefficient minimizes the L² misfit for that shot. + * If d0 is the computed value for a unit shot and d the reference data, we look for a complexe number "alpha" such that + * alpha * d0 is the best approximation to d. + * It is given by alpha = <d0, d> / <d0, d0> with < , > the complex inner product +*/ +template<Physic T_Physic> +std::vector<std::vector<std::complex<double>>> leastSquareShotIntensity(const Data<T_Physic>& d0, const Data<T_Physic>& d) +{ + areCompatible(d0, d); + std::vector<std::vector<std::complex<double>>> alpha; + alpha.resize(d.nf()); + for (auto &shots: alpha) + shots.resize(d.ns()); + + + for (unsigned int f = 0; f < d.nf(); f++) + { + for (unsigned int s = 0; s < d.ns(); s++) + { + std::complex<double> num = 0.; + std::complex<double> den = 0.; + + for (unsigned int r = 0; r < d.nr(s); r++) + { + num += conj<T_Physic>(d0.value(f,s,r)) * d.value(f,s,r); + den += conj<T_Physic>(d0.value(f,s,r)) * d0.value(f,s,r); + } + + alpha[f][s] = num / den; + } + } + return alpha; +}; + #endif // H_COMMON_DATA_ELEMENT diff --git a/tests/test.cpp b/tests/test.cpp index 90c0a351a8915bd26f239058458196d16bbf037b..942cb3ce202bf4cf6ce5788f2a26fccf6e8b9810 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -7,6 +7,10 @@ #include <gmshfem/Formulation.h> #include <gmshfem/GmshFem.h> #include "../common/wave/element.h" +#include "../common/data/element.h" +#include "../specific/data/element.h" + +#include "dummyAcquisition.h" #define INIT_GMSHFEM() \ namespace geo = gmsh::model::geo; \ @@ -102,4 +106,26 @@ BOOST_AUTO_TEST_CASE(dummy_test) BOOST_CHECK(1 + 1 == 2); } +BOOST_AUTO_TEST_CASE(LeastSquareShotIntensityFitAcoustic) +{ + std::complex<double> im(0, 1); + // A Data needs a configuration with some shots + DummyConfigurationOneSource config(2, {1, 2}); // Check perfect fit if one rec, least square fit for 2 + std::vector<double> freqs{0}; + Data<Physic::acoustic> target(freqs, &config); + Data<Physic::acoustic> unitShot(freqs, &config); + target.value(0, 0, 0, im); + unitShot.value(0, 0, 0, 1.); // Expected coef is "i" + + target.value(0, 1, 0, 1. + im); + target.value(0, 1, 1, 1. - im); + + unitShot.value(0, 1, 0, -2); + unitShot.value(0, 1, 1, -2); // Expected coef is -0.5 + + auto alphas = leastSquareShotIntensity(unitShot, target); + BOOST_CHECK_SMALL(std::norm(alphas[0][0] - im), 1e-4); + BOOST_CHECK_SMALL(std::norm(alphas[0][1] - (-0.5)), 1e-4); +} + BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file