Skip to content
Snippets Groups Projects
Select Git revision
  • 8a6f831845f7062cf990d45f8dfa32be76f8bd8d
  • develop default protected
  • feature/gmsh-reader
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

AcousticSolver.hpp

Blame
  • MeshData.hpp 11.78 KiB
    #ifndef MESH_DATA_HPP
    #define MESH_DATA_HPP
    
    #include <algebra/TinyVector.hpp>
    #include <mesh/IMeshData.hpp>
    #include <mesh/ItemValue.hpp>
    #include <mesh/SubItemValuePerItem.hpp>
    #include <utils/Exceptions.hpp>
    #include <utils/Messenger.hpp>
    #include <utils/PugsUtils.hpp>
    
    #include <map>
    
    template <size_t Dimension>
    class Connectivity;
    
    template <typename ConnectivityType>
    class Mesh;
    
    template <size_t Dimension>
    class MeshData : public IMeshData
    {
     public:
      static_assert(Dimension > 0, "dimension must be strictly positive");
      static_assert((Dimension <= 3), "only 1d, 2d and 3d are implemented");
    
      using MeshType = Mesh<Connectivity<Dimension>>;
    
      using Rd = TinyVector<Dimension>;
    
      static constexpr double inv_Dimension = 1. / Dimension;
    
     private:
      const MeshType& m_mesh;
      NodeValuePerCell<const Rd> m_Cjr;
      NodeValuePerCell<const double> m_ljr;
      NodeValuePerCell<const Rd> m_njr;
      FaceValue<const Rd> m_xl;
      CellValue<const Rd> m_xj;
      CellValue<const double> m_Vj;
      FaceValue<const double> m_ll;
      NodeValuePerFace<const Rd> m_Nlr;
    
      PUGS_INLINE
      void
      _compute_ll()
      {
        if constexpr (Dimension == 1) {
          static_assert(Dimension != 1, "ll does not make sense in 1d");
        } else {
          const auto& nlr = this->Nlr();
    
          FaceValue<double> ll{m_mesh.connectivity()};
    
          const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
          parallel_for(
            m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
              const auto& face_nodes = face_to_node_matrix[face_id];
    
              double lenght = 0;
              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
                lenght += l2Norm(nlr(face_id, i_node));
              }
    
              ll[face_id] = lenght;
            });
    
          m_ll = ll;
        }
      }
    
      PUGS_INLINE void
      _computeFaceIsobarycenter()
      {   // Computes vertices isobarycenter
        if constexpr (Dimension == 1) {
          static_assert(Dimension != 1, "xl does not make sense in 1d");
        } else {
          const NodeValue<const Rd>& xr = m_mesh.xr();
    
          const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
          FaceValue<Rd> xl(m_mesh.connectivity());
          parallel_for(
            m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
              Rd X                   = zero;
              const auto& face_nodes = face_to_node_matrix[face_id];
              for (size_t R = 0; R < face_nodes.size(); ++R) {
                X += xr[face_nodes[R]];
              }
              xl[face_id] = 1. / face_nodes.size() * X;
            });
          m_xl = xl;
        }
      }
    
      void
      _computeCellIsobarycenter()
      {   // Computes vertices isobarycenter
        if constexpr (Dimension == 1) {
          const NodeValue<const Rd>& xr = m_mesh.xr();
    
          const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
    
          CellValue<Rd> xj(m_mesh.connectivity());
          parallel_for(
            m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
              const auto& cell_nodes = cell_to_node_matrix[j];
              xj[j]                  = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
            });
          m_xj = xj;
        } else {
          const NodeValue<const Rd>& xr = m_mesh.xr();
    
          const CellValue<const double>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes();
    
          const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
          CellValue<Rd> xj(m_mesh.connectivity());
          parallel_for(
            m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
              Rd X                   = zero;
              const auto& cell_nodes = cell_to_node_matrix[j];
              for (size_t R = 0; R < cell_nodes.size(); ++R) {
                X += xr[cell_nodes[R]];
              }
              xj[j] = inv_cell_nb_nodes[j] * X;
            });
          m_xj = xj;
        }
      }
    
      PUGS_INLINE
      void
      _computeCellVolume()
      {
        const NodeValue<const Rd>& xr   = m_mesh.xr();
        const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
    
        auto Cjr = this->Cjr();
    
        CellValue<double> Vj(m_mesh.connectivity());
        parallel_for(
          m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
            double sum_cjr_xr      = 0;
            const auto& cell_nodes = cell_to_node_matrix[j];
    
            for (size_t R = 0; R < cell_nodes.size(); ++R) {
              sum_cjr_xr += (xr[cell_nodes[R]], Cjr(j, R));
            }
            Vj[j] = inv_Dimension * sum_cjr_xr;
          });
        m_Vj = Vj;
      }
    
      PUGS_INLINE
      void
      _computeNlr()
      {
        if constexpr (Dimension == 1) {
          static_assert(Dimension != 1, "Nlr does not make sense in 1d");
        } else if constexpr (Dimension == 2) {
          const NodeValue<const Rd>& xr = m_mesh.xr();
    
          NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
          const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
    
          parallel_for(
            m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId l) {
              const auto& face_nodes = face_to_node_matrix[l];
    
              const Rd xr0 = xr[face_nodes[0]];
              const Rd xr1 = xr[face_nodes[1]];
              const Rd dx  = xr1 - xr0;
    
              const Rd Nr = 0.5 * Rd{dx[1], -dx[0]};
    
              Nlr(l, 0) = Nr;
              Nlr(l, 1) = Nr;
            });
          m_Nlr = Nlr;
        } else {
          const NodeValue<const Rd>& xr = m_mesh.xr();
    
          NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
          const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
    
          parallel_for(
            m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId l) {
              const auto& face_nodes = face_to_node_matrix[l];
              const size_t nb_nodes  = face_nodes.size();
              std::vector<Rd> dxr(nb_nodes);
              for (size_t r = 0; r < nb_nodes; ++r) {
                dxr[r] = xr[face_nodes[(r + 1) % nb_nodes]] - xr[face_nodes[(r + nb_nodes - 1) % nb_nodes]];
              }
              const double inv_12_nb_nodes = 1. / (12. * nb_nodes);
              for (size_t r = 0; r < nb_nodes; ++r) {
                Rd Nr            = zero;
                const Rd two_dxr = 2 * dxr[r];
                for (size_t s = 0; s < nb_nodes; ++s) {
                  Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes[s]]);
                }
                Nr *= inv_12_nb_nodes;
                Nr -= (1. / 6.) * crossProduct(dxr[r], xr[face_nodes[r]]);
                Nlr(l, r) = Nr;
              }
            });
          m_Nlr = Nlr;
        }
      }
    
      PUGS_INLINE
      void
      _computeCjr()
      {
        if constexpr (Dimension == 1) {
          NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
          parallel_for(
            m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
              Cjr(j, 0) = -1;
              Cjr(j, 1) = 1;
            });
          m_Cjr = Cjr;
        } else if constexpr (Dimension == 2) {
          const NodeValue<const Rd>& xr   = m_mesh.xr();
          const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
    
          NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
          parallel_for(
            m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
              const auto& cell_nodes = cell_to_node_matrix[j];
              for (size_t R = 0; R < cell_nodes.size(); ++R) {
                int Rp1         = (R + 1) % cell_nodes.size();
                int Rm1         = (R + cell_nodes.size() - 1) % cell_nodes.size();
                Rd half_xrp_xrm = 0.5 * (xr[cell_nodes[Rp1]] - xr[cell_nodes[Rm1]]);
                Cjr(j, R)       = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
              }
            });
          m_Cjr = Cjr;
        } else if (Dimension == 3) {
          auto Nlr = this->Nlr();
    
          const auto& face_to_node_matrix   = m_mesh.connectivity().faceToNodeMatrix();
          const auto& cell_to_node_matrix   = m_mesh.connectivity().cellToNodeMatrix();
          const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
          const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
    
          NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
          parallel_for(
            Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Cjr[jr] = zero; });
    
          parallel_for(
            m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
              const auto& cell_nodes = cell_to_node_matrix[j];
    
              const auto& cell_faces       = cell_to_face_matrix[j];
              const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
    
              for (size_t L = 0; L < cell_faces.size(); ++L) {
                const FaceId& l        = cell_faces[L];
                const auto& face_nodes = face_to_node_matrix[l];
    
                auto local_node_number_in_cell = [&](NodeId node_number) {
                  for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
                    if (node_number == cell_nodes[i_node]) {
                      return i_node;
                    }
                  }
                  return std::numeric_limits<size_t>::max();
                };
    
                if (face_is_reversed[L]) {
                  for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
                    const size_t R = local_node_number_in_cell(face_nodes[rl]);
                    Cjr(j, R) -= Nlr(l, rl);
                  }
                } else {
                  for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
                    const size_t R = local_node_number_in_cell(face_nodes[rl]);
                    Cjr(j, R) += Nlr(l, rl);
                  }
                }
              }
            });
    
          m_Cjr = Cjr;
        }
      }
    
      PUGS_INLINE
      void
      _compute_ljr()
      {
        auto Cjr = this->Cjr();
        if constexpr (Dimension == 1) {
          NodeValuePerCell<double> ljr(m_mesh.connectivity());
          parallel_for(
            ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
          m_ljr = ljr;
    
        } else {
          NodeValuePerCell<double> ljr(m_mesh.connectivity());
          parallel_for(
            Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
          m_ljr = ljr;
        }
      }
    
      PUGS_INLINE
      void
      _compute_njr()
      {
        auto Cjr = this->Cjr();
        if constexpr (Dimension == 1) {
          // in 1d njr=Cjr (here performs a shallow copy)
          m_njr = m_Cjr;
        } else {
          auto ljr = this->ljr();
    
          NodeValuePerCell<Rd> njr(m_mesh.connectivity());
          parallel_for(
            Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; });
          m_njr = njr;
        }
      }
    
      void
      _checkCellVolume()
      {
        Assert(m_Vj.isBuilt());
    
        bool is_valid = [&] {
          for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
            if (m_Vj[j] <= 0) {
              return false;
            }
          }
          return true;
        }();
    
        if (not parallel::allReduceAnd(is_valid)) {
          throw NormalError("mesh contains cells of non-positive volume");
        }
      }
    
     public:
      PUGS_INLINE
      const MeshType&
      mesh() const
      {
        return m_mesh;
      }
    
      PUGS_INLINE
      FaceValue<const double>
      ll()
      {
        if (not m_ll.isBuilt()) {
          this->_compute_ll();
        }
        return m_ll;
      }
    
      PUGS_INLINE
      NodeValuePerFace<const Rd>
      Nlr()
      {
        if (not m_Nlr.isBuilt()) {
          this->_computeNlr();
        }
        return m_Nlr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const Rd>
      Cjr()
      {
        if (not m_Cjr.isBuilt()) {
          this->_computeCjr();
        }
        return m_Cjr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const double>
      ljr()
      {
        if (not m_ljr.isBuilt()) {
          this->_compute_ljr();
        }
        return m_ljr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const Rd>
      njr()
      {
        if (not m_njr.isBuilt()) {
          this->_compute_njr();
        }
        return m_njr;
      }
    
      PUGS_INLINE
      FaceValue<const Rd>
      xl()
      {
        if (not m_xl.isBuilt()) {
          this->_computeFaceIsobarycenter();
        }
        return m_xl;
      }
    
      PUGS_INLINE
      CellValue<const Rd>
      xj()
      {
        if (not m_xj.isBuilt()) {
          this->_computeCellIsobarycenter();
        }
        return m_xj;
      }
    
      PUGS_INLINE
      CellValue<const double>
      Vj()
      {
        if (not m_Vj.isBuilt()) {
          this->_computeCellVolume();
          this->_checkCellVolume();
        }
        return m_Vj;
      }
    
     private:
      // MeshData **must** be constructed through MeshDataManager
      friend class MeshDataManager;
      MeshData(const MeshType& mesh) : m_mesh(mesh) {}
    
     public:
      MeshData() = delete;
    
      MeshData(const MeshData&) = delete;
      MeshData(MeshData&&)      = delete;
    
      ~MeshData() {}
    };
    
    #endif   // MESH_DATA_HPP