Skip to content
Snippets Groups Projects
Select Git revision
  • 972ab2f158f193ecdc14247c01ea434e6da66e6e
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • 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
  • 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

MeshData.hpp

Blame
  • Stephane Del Pino's avatar
    Stéphane Del Pino authored
    Separate the notion of Nlr (Eucclhyd normals) of the calculation of
    Cjr in 3D. It remain to implement its calculation in 2D to finish this
    part.
    
    This should allow a trivial implementation of the Eucclhyd scheme and
    simplify coding off pressure boundary conditions.
    972ab2f1
    History
    MeshData.hpp 9.91 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;
      std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
      std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
      std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
      std::shared_ptr<CellValue<const Rd>> m_xj;
      std::shared_ptr<CellValue<const double>> m_Vj;
      std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr;
    
      PUGS_INLINE
      void
      _computeIsobarycenter()
      {   // 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 = std::make_shared<CellValue<const Rd>>(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 = std::make_shared<CellValue<const Rd>>(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 = std::make_shared<CellValue<const double>>(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) {
          throw NotImplementedError("Nlr are not yet computed in 2d");
        } 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 = std::make_shared<NodeValuePerFace<const Rd>>(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 = std::make_shared<NodeValuePerCell<const Rd>>(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 = std::make_shared<NodeValuePerCell<const Rd>>(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 = std::make_shared<NodeValuePerCell<const Rd>>(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 = std::make_shared<NodeValuePerCell<const double>>(ljr);
    
        } else {
          NodeValuePerCell<double> ljr(m_mesh.connectivity());
          parallel_for(
            Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
          m_ljr = std::make_shared<NodeValuePerCell<const double>>(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 = std::make_shared<NodeValuePerCell<const Rd>>(njr);
        }
      }
    
      void
      _checkCellVolume()
      {
        auto Vj = this->Vj();
    
        bool is_valid = [&] {
          for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
            if (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
      NodeValuePerFace<const Rd>
      Nlr()
      {
        if (not m_Nlr) {
          this->_computeNlr();
        }
        return *m_Nlr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const Rd>
      Cjr()
      {
        if (not m_Cjr) {
          this->_computeCjr();
        }
        return *m_Cjr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const double>
      ljr()
      {
        if (not m_ljr) {
          this->_compute_ljr();
        }
        return *m_ljr;
      }
    
      PUGS_INLINE
      NodeValuePerCell<const Rd>
      njr()
      {
        if (not m_njr) {
          this->_compute_njr();
        }
        return *m_njr;
      }
    
      PUGS_INLINE
      CellValue<const Rd>
      xj()
      {
        if (not m_xj) {
          this->_computeIsobarycenter();
        }
        return *m_xj;
      }
    
      PUGS_INLINE
      CellValue<const double>
      Vj()
      {
        if (not m_Vj) {
          this->_computeCellVolume();
          this->_checkCellVolume();
        }
        return *m_Vj;
      }
    
      void
      updateAllData()
      {
        ;
      }
    
      MeshData(const MeshType& mesh) : m_mesh(mesh) {}
    
      MeshData() = delete;
    
      MeshData(const MeshData&) = default;
      MeshData(MeshData&&)      = default;
    
      ~MeshData()
      {
        ;
      }
    };
    
    #endif   // MESH_DATA_HPP