Skip to content
Snippets Groups Projects
Select Git revision
  • 6e40cab6b242754dd6c687ed2bce201e4ef113f2
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

PViewDataRemote.h

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    main.cpp 12.26 KiB
    #include <gmsh.h>
    #include <gmshddm/Formulation.h>
    #include <gmshddm/GmshDdm.h>
    #include <gmshddm/MPIInterface.h>
    #include <gmshfem/Message.h>
    #include <gmshfem/AnalyticalFunction.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::analytics;
    using namespace gmshfem::post;
    
    void sphere(const double Rin, const double Rext, const double lc)
    {
      gmsh::model::add("sphere");
      int s_in = gmsh::model::occ::addSphere(0, 0, 0, Rin);
      int s_ext = gmsh::model::occ::addSphere(0, 0, 0, Rext);
    
      std::vector< std::pair< int, int > > ov_phy;
      std::vector< std::vector< std::pair< int, int > > > ovv_phy;
      gmsh::model::occ::cut({{3, s_ext}}, {{3, s_in}}, ov_phy, ovv_phy, 5);
    
      gmsh::model::occ::synchronize();
    
      gmsh::model::addPhysicalGroup(2, {s_ext}, 1);
      gmsh::model::setPhysicalName(2, 1, "gamma_scat");
      gmsh::model::addPhysicalGroup(2, {s_in}, 2);
      gmsh::model::setPhysicalName(2, 2, "gamma_ext");
    
      gmsh::model::addPhysicalGroup(3, {5}, 1);
      gmsh::model::setPhysicalName(3, 1, "omega");
    
      gmsh::option::setNumber("Mesh.MeshSizeMax", lc);
    
      gmsh::model::occ::synchronize();
      gmsh::model::mesh::generate(3);
    }
    
    int main(int argc, char **argv)
    {
      GmshDdm gmshDdm(argc, argv);
    
      const double pi = 3.14159265358979323846264338327950288;
      const std::complex< double > im(0., 1.);
    
      std::string mesh = "";
      gmshDdm.userDefinedParameter(mesh, "mesh");
    
      int nDom = 2;
      gmshDdm.userDefinedParameter(nDom, "nDom");
    
      double lc = 0.1;
      gmshDdm.userDefinedParameter(lc, "lc");
    
      double Rin = 0.5;
      gmshDdm.userDefinedParameter(Rin, "Rin");
    
      double Rext = 1.5;
      gmshDdm.userDefinedParameter(Rin, "Rext");
    
      double k = 1;
      gmshDdm.userDefinedParameter(k, "k");
    
      double theta = pi/3.; // orientation of the incident wave - azimuthal angle [0, 2*pi]
      gmshDdm.userDefinedParameter(theta, "theta");
      double phi = pi/2.; // orientation of the incident wave - polar angle [0, pi]
      gmshDdm.userDefinedParameter(phi, "phi");
    
      bool ABC2 = true; // false: Sommerfeld ABC, true: 2nd order ABC
      gmshDdm.userDefinedParameter(ABC2, "ABC2");
      if(ABC2 && !getMPIRank())
        msg::info << " - use 2nd order exterior absorbing boundary condition" << msg::endl;
    
      int FEMorder = 2;
      gmshDdm.userDefinedParameter(FEMorder, "FEMorder");
      std::string gauss = "Gauss" + std::to_string(2 * FEMorder + 1);
    
      gmsh::option::setNumber("Mesh.ElementOrder", FEMorder);
      gmsh::option::setNumber("Mesh.SecondOrderLinear", 1); // straight-sided elements
    
      std::string Transmission = "Taylor2";
      gmshDdm.userDefinedParameter(Transmission, "TC");
    
      double alpha = -pi / 2.;
      gmshDdm.userDefinedParameter(alpha, "alpha");
    
      // if no mesh is given, create a partitioned mesh of a sphere on rank 0
      if(mesh.empty()) {
        if(!getMPIRank()) {
          float pointsByWl = 2 * pi * FEMorder / k / lc;
          msg::info << " - dofs density by wavelength = " << pointsByWl << "" << msg::endl;
          if(pointsByWl < 6)
            msg::warning << " - less than 6 dofs per wavelength !" << msg::endl;
          // build geometry and mesh
          sphere(Rin, Rext, lc);
          // partition
          gmsh::model::mesh::partition(nDom);
          // save all, as partitioning will create some entities without physical groups
          gmsh::option::setNumber("Mesh.SaveAll", 1);
          // save partitioned mesh in single file for sequential runs
          gmsh::option::setNumber("Mesh.PartitionSplitMeshFiles", 0);
          gmsh::write("sphere.msh");
          // save partitioned mesh in separate file for distributed runs
          gmsh::option::setNumber("Mesh.PartitionSplitMeshFiles", 1);
          gmsh::write("sphere.msh");
          gmsh::clear();
        }
        mesh = "sphere";
      }
      barrier();
    
      // read partitioned mesh
      if(getMPISize() == 1) {
        gmsh::merge(mesh + ".msh");
      }
      else {
        for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
          if (mpi::isItMySubdomain(i)) {
            gmsh::merge(mesh + "_" + std::to_string(i + 1) + ".msh");
          }
        }
      }
      if(gmsh::model::getNumberOfPartitions() != nDom)
        msg::error << "Wrong number of partitions in mesh file " << gmsh::model::getNumberOfPartitions() << " vs. " << nDom << msg::endl;
    
      Domain omegaMono = gmshfem::domain::Domain("omega");
      Domain gammaScatMono = gmshfem::domain::Domain("gamma_scat");
      Domain gammaExtMono = gmshfem::domain::Domain("gamma_ext");
    
      Subdomain omega = Subdomain::buildSubdomainOf(omegaMono);
      Subdomain gammaScat = Subdomain::buildSubdomainOf(gammaScatMono);
      Subdomain gammaExt = Subdomain::buildSubdomainOf(gammaExtMono);
    
      std::vector< std::vector< unsigned int > > topology;
      Interface sigma = Interface::buildInterface(topology);
    
      gmshddm::problem::Formulation< std::complex< double > > formulation("Helmholtz", topology);
    
      SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaScat | gammaExt | sigma, FunctionSpaceTypeForm0::Lagrange);
      InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::Lagrange);
      formulation.addInterfaceField(g);
    
      // Define Analytical solution
      Function< std::complex< double >, Degree::Degree0 > *u_exact = nullptr;
      if(mesh == "sphere") {
        u_exact = new AnalyticalFunction< helmholtz3D::ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta, phi); // -k because we use e^iwt convention
        u_exact->activateMemorization();
      }
    
      // Define incident field
      Function< std::complex< double >, Degree::Degree0 > uinc, dn_uinc;
    
      if(mesh == "sphere") {
        uinc = exp< std::complex< double > >(im * k * ( sin(phi)*cos(theta)*x< std::complex< double > >() + sin(phi)*sin(theta)*y< std::complex< double > >() + cos(phi)*z< std::complex< double > >() ) );
        dn_uinc = AnalyticalFunction< helmholtz3D::dr_ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta, phi);
      }
      else {
        uinc = exp< std::complex< double > >(im * k * z< std::complex< double > >());
        dn_uinc = normal< std::complex< double > >() * vector< std::complex< double > > (0, 0, im * k * uinc);
      }
    
      ScalarFunction< std::complex< double > > Si = 0., dS = 0.;
      if(Transmission == "Taylor0") {
        if(!getMPIRank())
          msg::info << " - Use Sommerfeld transmission condition with rotation angle = " << alpha << "" << msg::endl;
        Si = im * k * exp(im * alpha / 2.);
        dS = 0.;
      }
      else if(Transmission == "Taylor2") {
        if(!getMPIRank())
          msg::info << " - Use second order transmission condition with rotation angle = " << alpha << "" << msg::endl;
        Si = im * k * cos(alpha / 2.);
        dS = im * exp(-im * alpha / 2.) / (2.0 * k);
      }
    
      for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
        // Volume weak form
        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);
    
        // Sommerfeld BC
        formulation(i).integral(im * k * dof(u(i)), tf(u(i)), gammaExt(i), gauss);
        if (ABC2) {
          formulation(i).integral( ( 1. / Rext ) * dof(u(i)), tf(u(i)), gammaExt(i), gauss);
          formulation(i).integral( ( -im / (2. * k) + 1. / (2. * k * k *Rext) ) * grad(dof(u(i))), grad(tf(u(i))), gammaExt(i), gauss);
        }
    
        // Incident field
        formulation(i).integral( formulation.physicalSource(dn_uinc) , tf(u(i)), gammaScat(i), gauss);
    
        // Artificial source
        for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
          const unsigned int j = topology[i][jj];
          if(!sigma(i, j).isEmpty()) {
            formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss);
            // Transmission conditions terms
            formulation(i).integral(Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss);
            formulation(i).integral(-dS * grad(dof(u(i))), grad(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];
          if(!sigma(i, j).isEmpty()) {
            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);
            // Transmission conditions terms
            formulation(i, j).integral(-2. * Si * u(i), tf(g(i, j)), sigma(i, j), gauss);
            formulation(i, j).integral(2. * dS * grad(u(i)), grad(tf(g(i, j))), sigma(i, j), gauss);
          }
        }
      }
    
      // Solve DDM
      double tolerence = 1e-6;
      int maxIter = 1000;
      formulation.pre();
      formulation.solve("gmres", tolerence, maxIter, true);
    
      if(u_exact) {
        // Compute analytic L2-error
        for(auto i = 0; i < nDom; ++i) {
          if(gmshddm::mpi::isItMySubdomain(i)) {
            std::complex< double > denLocal = integrate(pow(abs(*u_exact), 2), omega(i), gauss);
            std::complex< double > numLocal = integrate(pow(abs(*u_exact - u(i)), 2), omega(i), gauss);
            double num = numLocal.real();
            double den = denLocal.real();
            double local_err = sqrt(num / den);
            msg::info << "Analytic-ddm L2 error (%) = " << 100.*local_err << " on subdomain " << i << msg::endl;
          }
        }
      }
    
      // save solutions
      gmsh::option::setNumber("Mesh.Binary", 1);
      gmsh::option::setNumber("PostProcessing.Binary", 1);
      gmsh::option::setNumber("PostProcessing.SaveMesh", 0);
    
      for(auto i = 0; i < nDom; ++i) {
        if(gmshddm::mpi::isItMySubdomain(i)) {
          int tag = save(u(i), omega(i), "u", "", "", true, 0, 0., i + 1);
          gmsh::view::write(tag, "u_" + std::to_string(i + 1) + ".msh");
    
          // also save a cut
          for(int step = 0; step < 2; step++) {
            gmsh::view::option::setNumber(tag, "TimeStep", step);
            gmsh::view::option::setNumber(tag, "AdaptVisualizationGrid", 1);
            gmsh::plugin::setNumber("CutPlane", "View", gmsh::view::getIndex(tag));
            gmsh::plugin::setNumber("CutPlane", "A", 0);
            gmsh::plugin::setNumber("CutPlane", "B", 0);
            gmsh::plugin::setNumber("CutPlane", "C", 1);
            gmsh::plugin::setNumber("CutPlane", "D", 0);
            gmsh::plugin::setNumber("CutPlane", "RecurLevel", FEMorder);
            gmsh::plugin::setNumber("CutPlane", "TargetError", -1); // 1e-4 to get decent AMR
            int tag2 = gmsh::plugin::run("CutPlane");
            std::string name2 = std::string("u_cut_step") + std::to_string(step) + "_" + std::to_string(i + 1);
            gmsh::view::option::setString(tag2, "Name", name2);
            gmsh::view::write(tag2, name2 + ".pos");
          }
        }
      }
    
      bool ComputeMono = false; // monodomain solution
      gmshDdm.userDefinedParameter(ComputeMono, "ComputeMono");
      if (ComputeMono) {
        // monodomain solution
        gmshfem::field::Field< std::complex< double >, Form::Form0 > uMono("uMono", omegaMono | gammaScatMono | gammaExtMono, FunctionSpaceTypeForm0::Lagrange);
        gmshfem::problem::Formulation< std::complex< double > > monodomain("Helmholtz");
    
        monodomain.integral(grad(dof(uMono)), grad(tf(uMono)), omegaMono, gauss);
        monodomain.integral(-k * k * dof(uMono), tf(uMono), omegaMono, gauss);
    
        // Sommerfeld BC
        monodomain.integral(im * k * dof(uMono), tf(uMono), gammaExtMono, gauss);
        if (ABC2) { // 2nd order ABC
          monodomain.integral( ( 1. / Rext ) * dof(uMono), tf(uMono), gammaExtMono, gauss);
          monodomain.integral( ( -im / (2. * k) + 1. / (2. * k * k *Rext) ) * grad(dof(uMono)), grad(tf(uMono)), gammaExtMono, gauss);
        }
    
        // Incident field
        monodomain.integral(dn_uinc, tf(uMono), gammaScatMono, gauss);
        // solve
        monodomain.pre();
        monodomain.assemble();
        monodomain.solve();
    
        // Compute analytic L2-error
        double numM = 0., denM = 0.;
        double mono_err; // L2 local error
        for (unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
          if(gmshddm::mpi::isItMySubdomain(i)) {
            std::complex< double > denMono = integrate(pow(abs(uMono), 2), omega(i), gauss);
            std::complex< double > numMono = integrate(pow(abs(uMono - u(i)), 2), omega(i), gauss);
            numM += numMono.real();
            denM += denMono.real();
            mono_err = sqrt(numM / denM);
            msg::info << "Monodomain-ddm L2 error (%) = " << 100.*mono_err << " on subdomain " << i << msg::endl;
          }
        }
    
      } // end Monodomain
    
      return 0;
    }