Initial model setup (flexible acquisition)

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");

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);