main.cpp

#include "gmsh.h"

#include <gmshddm/Formulation.h>
#include <gmshddm/GmshDdm.h>
#include <gmshddm/MPIInterface.h>
#include <gmshfem/Message.h>

using namespace gmshddm;
using namespace gmshddm::common;
using namespace gmshddm::domain;
using namespace gmshddm::problem;
using namespace gmshddm::field;
using namespace gmshddm::mpi;

using namespace gmshfem;
using namespace gmshfem::domain;
using namespace gmshfem::equation;
using namespace gmshfem::problem;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::post;

#include <cstring>
#include <fstream>
#include <sstream>

static std::string freqToStr(double f) {
  char buffer[100];
  sprintf(buffer, "%.2f", f);
  return std::string(buffer);
}

struct sourcePosition{
  double x;
  double y;
  sourcePosition(double x, double y) : x(x), y(y) {}
};

std::vector<sourcePosition> readSources(const std::string& path) {
    std::ifstream inFile(path);
    std::vector<sourcePosition> positions;

    if (!inFile.is_open()) {
        throw std::runtime_error("Unable to open file: " + path);
    }

    std::string line;
    int lineNumber = 0;
    while (std::getline(inFile, line)) {
        lineNumber++;
        std::istringstream lineStream(line);
        double x, y;
        char delimiter;

        if (lineStream >> x >> delimiter >> y && delimiter == ';') {
            positions.push_back({x, y});
        } else {
            throw std::runtime_error("Invalid line format in file: " + path + " at line " + std::to_string(lineNumber));
        }
    }

    inFile.close();
    return positions;
}


/**
 * Returns ID of the domains containing the source
*/
std::vector<int> mesh(double L, double H, double x0, double y0, const std::vector<sourcePosition>& sourcePositions, unsigned nX = 2, unsigned nY = 2, double lc = 0.1) {
  gmsh::model::add("freeSpace");
  // Disable warnings for large meshes
  gmsh::option::setNumber("General.ExpertMode", 1);

  namespace gmodel = gmsh::model;

  for(unsigned i = 0; i < sourcePositions.size(); ++i) {
      double x_src = sourcePositions[i].x;
      double y_src = sourcePositions[i].y;
        if(x_src < x0 || y_src < y0 || x_src > (x0 + L) || y_src > (y0 + H)) {
            throw gmshfem::common::Exception("Invalid source location.");
        }
  }

  std::vector<std::vector<int>> points(nX + 1);
  std::vector<std::vector<int>> edgesH(nX);
  std::vector<std::vector<int>> edgesV(nX+1);
  std::vector<std::vector<int>> curvedLoops(nX);
  // Indices from 0 to N (N+1 points)
  for (unsigned i = 0; i <= nX; ++i)
    points[i].resize(nY + 1);

  const double dx = L / nX;
  const double dy = H / nY;

  // Create all points
  for (unsigned i = 0; i <= nX; ++i) {
    for (unsigned j = 0; j <= nY; ++j) {
      points[i][j] = gmodel::geo::addPoint(x0 + i * dx, y0 + j * dy, 0.0, lc);
    }
  }

  // Create all horizontal edges
  for (unsigned i = 0; i < nX; ++i) {
    edgesH[i].resize(nY+1);
    for (unsigned j = 0; j <= nY; ++j) {
      edgesH[i][j] = gmodel::geo::addLine(points.at(i).at(j), points.at(i+1).at(j));
    }
  }

  // Create all vertical edges
  for (unsigned i = 0; i <= nX; ++i) {
    edgesV[i].resize(nY);
    for (unsigned j = 0; j < nY; ++j) {
      edgesV[i][j] = gmodel::geo::addLine(points.at(i).at(j), points.at(i).at(j+1));
    }
  }

  // Create all surfaces
  for (unsigned i = 0; i < nX; ++i) {
    curvedLoops[i].resize(nY);
    for (unsigned j = 0; j < nY; ++j) {
      curvedLoops[i][j] = gmodel::geo::addCurveLoop({edgesH.at(i).at(j), edgesV.at(i+1).at(j), -edgesH.at(i).at(j+1), -edgesV.at(i).at(j)});
      gmodel::geo::addPlaneSurface({curvedLoops[i][j]}, curvedLoops[i][j]); // Same indices
    }
  }

  // Create sources. TODO: optimize
  std::vector<int> sources;
  std::vector<int> sourcesDomains;
  for (const auto& pos : sourcePositions) {
    int source = gmodel::geo::addPoint(pos.x, pos.y, 0, lc);
    int i_source = int(pos.x / dx);
    int j_source = int(pos.y / dy);
    gmodel::geo::synchronize();
    gmodel::mesh::embed(0, {source}, 2, curvedLoops[i_source][j_source]);
    sources.push_back(source);
    sourcesDomains.push_back(i_source + j_source * nX);
  }

  gmodel::geo::synchronize();

  /* Define physical entities */

  

  // Domains
  std::vector<int> omegaAll;
  for (unsigned i = 0; i < nX; ++i) {
    for (unsigned j = 0; j < nY; ++j) {
      int omegaIJ = gmodel::addPhysicalGroup(2, {curvedLoops[i][j]}, -1);
      omegaAll.push_back(omegaIJ);
      gmodel::setPhysicalName(2, omegaIJ, "omega_" + std::to_string(i) + "_" + std::to_string(j));
    }
  }

  int omegaAllIdx = gmodel::addPhysicalGroup(2, omegaAll, -1);
  gmodel::setPhysicalName(2, omegaAllIdx, "omega");

  // Outer boundary
  std::vector<int> gammaOut;
  for (unsigned i = 0; i < nX; ++i)
  {
    gammaOut.push_back(edgesH[i][0]);
    gammaOut.push_back(edgesH[i][nY]);
  }

  for (unsigned j = 0; j < nY; ++j)
  {
    gammaOut.push_back(edgesV[0][j]);
    gammaOut.push_back(edgesV[nX][j]);
  }

  int gammaOutIdx = gmodel::addPhysicalGroup(1, gammaOut, -1);
  gmodel::setPhysicalName(1, gammaOutIdx, "gammaOut");

  // All eges
  // H-edges
  for (unsigned i = 0; i < nX; ++i) {
    for (unsigned j = 0; j <= nY; ++j) {
      int HIJ = gmodel::addPhysicalGroup(1, {edgesH[i][j]}, -1);
      gmodel::setPhysicalName(1, HIJ, "edge_h_" + std::to_string(i) + "_" + std::to_string(j));
    }
  }

  // V-edges
  for (unsigned i = 0; i <= nX; ++i) {
    for (unsigned j = 0; j < nY; ++j) {
      int VIJ = gmodel::addPhysicalGroup(1, {edgesV[i][j]}, -1);
      gmodel::setPhysicalName(1, VIJ, "edge_v_" + std::to_string(i) + "_" + std::to_string(j));
    }
  }

  // Points
  for (unsigned i = 0; i <= nX; ++i) {
    for (unsigned j = 0; j <= nY; ++j) {
      int VIJ = gmodel::addPhysicalGroup(1, {points[i][j]}, -1);
      gmodel::setPhysicalName(1, VIJ, "point_" + std::to_string(i) + "_" + std::to_string(j));
    }
  }

  // Src
  for (unsigned srcId = 0; srcId < sources.size(); ++srcId) {
    int srcEntity = gmodel::addPhysicalGroup(0, {sources[srcId]});
    gmodel::setPhysicalName(0, srcEntity, "source_" + std::to_string(srcId));

  }

  gmodel::geo::synchronize();
  gmodel::mesh::generate();

  return sourcesDomains;

}

int main(int argc, char **argv)
{
  GmshDdm gmshDdm(argc, argv);
  constexpr std::complex<double> im(0, 1);
  constexpr double pi = 3.14159265359;

  unsigned rank = gmshddm::mpi::getMPIRank();

  int order = 2;
  gmshDdm.userDefinedParameter(order, "order");
  double f = 2;
  gmshDdm.userDefinedParameter(f, "f");

  double lc = (1. / f) / 10;


  gmshDdm.userDefinedParameter(lc, "lc");

  std::string gauss = "Gauss10";
  gmshDdm.userDefinedParameter(gauss, "gauss");

  // Output parameters
  bool mustSolve{true};
  gmshDdm.userDefinedParameter(mustSolve, "noSolve");


  auto nX = 2;
  auto nY = 2;
  gmshDdm.userDefinedParameter(nX, "nX");
  gmshDdm.userDefinedParameter(nY, "nY");

  double L{1}, H{1};
  gmshDdm.userDefinedParameter(L, "L");
  gmshDdm.userDefinedParameter(H, "H");

  

  bool save_matrix = false;
  gmshDdm.userDefinedParameter(save_matrix, "save_iteration_matrix");

  // TODO : Weird behavior for negative y-axis -> Unsupported x0, y0
  std::vector<sourcePosition> sources;
  std::string sourcePath = "../srcs.csv";
  gmshDdm.userDefinedParameter(sourcePath, "sourcePath");
  sources = readSources(sourcePath);
  auto d_source = mesh(L, H, 0.0, 0, sources, nX, nY, lc);


  // Configure wavenumber. If no model is given, assume homogenuous velocity of 1
  // (i.e. k = 2 * pi * f)
  // If a model is given, it is assumed to be a squared slowness.
  gmshfem::function::Function<std::complex<double>, gmshfem::Degree::Degree0> k;

  std::string model_path = "";
  std::string model_name;
  gmshDdm.userDefinedParameter(model_path, "model_path");
  if (model_path != "") {
    model_name = model_path.substr(0, model_path.size() - 4);
    gmsh::merge(model_path);
    std::vector<int> tags;
    gmsh::view::getTags(tags);

    auto view = tags.back();
    ScalarFunction<std::complex<double>> model = probeScalarView<std::complex<double>>(view);
    // k = omega * slowness
    k = 2 * pi * f * sqrt(model);
  }
  else {
    // Unit velocity
    k = 2 * pi * f;
    model_name = "free_space";
  }
  
  gmshDdm.userDefinedParameter(model_name, "model_name");
  msg::info << "Model name is " << model_name << gmshfem::msg::endl;
  

  auto nDom = nX * nY;
  // Define domain
  Subdomain omega(nDom);
  Subdomain gammaOut(nDom);
  Interface sigma(nDom);

  #define OMEGA(i,j) "omega_" + std::to_string((i)) + "_" + std::to_string((j))
  #define POINT(i,j) "point_" + std::to_string((i)) + "_" + std::to_string((j))
  #define EDGE_H(i,j) "edge_h_" + std::to_string((i)) + "_" + std::to_string((j))
  #define EDGE_V(i,j) "edge_v_" + std::to_string((i)) + "_" + std::to_string((j))

  for(int i = 0; i < nX; ++i) {
    for (int j = 0; j < nY; ++j) {
      int d = i + nX * j;
      omega(d) = Domain(OMEGA(i,j));

      // Left border / interface
      if(i == 0) {
        gammaOut(d) |= Domain(EDGE_V(i, j));
      }
      else {
        sigma(d, d - 1) =  Domain(EDGE_V(i, j));
      }
      // Right border / interface
      if(i == nX - 1) {
        gammaOut(d) |= Domain(EDGE_V(i + 1, j));
      }
      else {
        sigma(d, d + 1) =  Domain(EDGE_V(i + 1, j));
      }
      // Bottom border / interface
      if(j == 0) {
        gammaOut(d) |= Domain(EDGE_H(i, j));
      }
      else {
        sigma(d, d - nX) =  Domain(EDGE_H(i, j));
      }
      // Top border / interface
      if(j == nY - 1) {
        gammaOut(d) |= Domain(EDGE_H(i, j + 1));
      }
      else {
        sigma(d, d + nX) =  Domain(EDGE_H(i, j + 1));
      }

      
    }
  }

  // Define topology
  std::vector< std::vector< unsigned int > > topology(nDom);
  for(unsigned int i = 0; i < nX; ++i) {
    for (unsigned int j = 0; j < nY; ++j) {
        unsigned int d = i + nX * j;
        if (i != 0) topology[d].push_back(d-1);
        if (i != nX-1) topology[d].push_back(d+1);
        if (j != 0) topology[d].push_back(d-nX);
        if (j != nY-1) topology[d].push_back(d+nX);
    }
  }


  // Create DDM formulation
  gmshddm::problem::Formulation< std::complex< double > > formulation("Helmholtz", topology);

  // Create fields
  SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaOut, FunctionSpaceTypeForm0::HierarchicalH1, order);
  InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order);
  // Tell to the formulation that g is field that have to be exchanged between subdomains
  formulation.addInterfaceField(g);

  // Add terms common to all RHS
  for(unsigned int i = 0; i < nDom; ++i) {
        // VOLUME TERMS
        formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss);
        formulation(i).integral(-k * k * dof(u(i)), tf(u(i)), omega(i), gauss);
        formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), gammaOut(i), gauss);

        // Artificial source
        for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
        const unsigned int j = topology[i][jj];
        formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), sigma(i, j), gauss);
        formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss);
        }

        // INTERFACE TERMS
        for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
        const unsigned int j = topology[i][jj];
        formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss);
        formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss);
        formulation(i, j).integral(2. * im * k * u(i), tf(g(i, j)), sigma(i, j), gauss);
        }
    }

  // Add all sources
  formulation.setShotNumber(sources.size());
  for(unsigned isrc = 0; isrc < sources.size(); ++isrc) {
    formulation.numberedPhysicalSourceTerm(isrc, formulation(d_source[isrc]).integral(formulation.physicalSource(-1), tf(u(d_source[isrc])), Domain("source_" + std::to_string(isrc)), gauss));
  }

      formulation.pre(false);


  for(unsigned isrc = 0; isrc < sources.size(); ++isrc) {
    formulation.disableAllShots();
    formulation.enableShot(isrc);

    // Solve the DDM formulation
    formulation.assemble();

    if(mustSolve) {
        int iterMax = 1000;
        std::string solver = "gmres";
        gmshDdm.userDefinedParameter(solver, "solver");
        gmshDdm.userDefinedParameter(iterMax, "iter");
        formulation.solve(solver, 1e-6, iterMax);
         save(u.buildUnifiedField(), u.getLocalDomain(), "u_rhs_" + std::to_string(isrc) + "_rank_" + std::to_string(rank));
    }


    std::string suffix = model_name + "_" + freqToStr(f) + "_" + std::to_string(nX) + "_" + std::to_string(nY);
    std::string mat_name{"matrix_" + suffix}, vec_name{"rhs_" + std::to_string(isrc) + "_" + suffix};
    if(save_matrix) {
        // Export all RHS but only one matrix
        if(isrc == 0) {
        auto mat = formulation.computeMatrix();
        msg::info << "nnz = " << mat.numberOfNonZero() << ", symmetric = " << mat.isSymmetric() << ", hermitian = " << mat.isHermitian() << msg::endl;
        mat.save(mat_name);
        }

        formulation.getRHS().save(vec_name);
    }
  }


  return 0;
}