Skip to content
Snippets Groups Projects
Select Git revision
  • 5cbfa289546515953ec10d5cb085ca883a347e28
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

CommandLine.h

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    main.cpp 12.98 KiB
    #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;
    }