Skip to content
Snippets Groups Projects
Select Git revision
  • 1a6edd28a9f95556b376bfbd0b336660d84847fb
  • 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

MeshNodeBoundary.hpp

Blame
  • MeshNodeBoundary.hpp 10.36 KiB
    #ifndef MESH_NODE_BOUNDARY_HPP
    #define MESH_NODE_BOUNDARY_HPP
    
    #include <Array.hpp>
    #include <ItemValue.hpp>
    
    #include <Kokkos_Vector.hpp>
    #include <TinyVector.hpp>
    
    #include <RefNodeList.hpp>
    #include <RefFaceList.hpp>
    
    #include <ConnectivityMatrix.hpp>
    #include <IConnectivity.hpp>
    #include <iostream>
    
    template <size_t dimension>
    class MeshNodeBoundary
    {
     protected:
      Array<const NodeId> m_node_list;
     public:
      MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
      MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
    
      const Array<const NodeId>& nodeList() const
      {
        return m_node_list;
      }
    
      template <typename MeshType>
      MeshNodeBoundary(const MeshType& mesh,
                       const RefFaceList& ref_face_list)
      {
        static_assert(dimension == MeshType::dimension);
        const auto& face_to_cell_matrix
            = mesh.connectivity().faceToCellMatrix();
    
        const Array<const FaceId>& face_list = ref_face_list.faceList();
        Kokkos::parallel_for(face_list.size(), KOKKOS_LAMBDA(const int& l){
            const auto& face_cells = face_to_cell_matrix[face_list[l]];
            if (face_cells.size()>1) {
              std::cerr << "internal faces cannot be used to define mesh boundaries\n";
              std::exit(1);
            }
          });
    
        Kokkos::vector<unsigned int> node_ids;
        // not enough but should reduce significantly the number of resizing
        node_ids.reserve(dimension*face_list.size());
        const auto& face_to_node_matrix
            = mesh.connectivity().faceToNodeMatrix();
    
        for (size_t l=0; l<face_list.size(); ++l) {
          const size_t face_number = face_list[l];
          const auto& face_nodes = face_to_node_matrix[face_number];
    
          for (size_t r=0; r<face_nodes.size(); ++r) {
            node_ids.push_back(face_nodes[r]);
          }
        }
        std::sort(node_ids.begin(), node_ids.end());
        auto last = std::unique(node_ids.begin(), node_ids.end());
        node_ids.resize(std::distance(node_ids.begin(), last));
    
        Array<NodeId> node_list(node_ids.size());
        Kokkos::parallel_for(node_ids.size(), KOKKOS_LAMBDA(const int& r){
            node_list[r] = node_ids[r];
          });
        m_node_list = node_list;
      }
    
      template <typename MeshType>
      MeshNodeBoundary(const MeshType& mesh,
                       const RefNodeList& ref_node_list)
          : m_node_list(ref_node_list.nodeList())
      {
        static_assert(dimension == MeshType::dimension);
      }
    
      MeshNodeBoundary() = default;
      MeshNodeBoundary(const MeshNodeBoundary&) = default;
      MeshNodeBoundary(MeshNodeBoundary&&) = default;
      virtual ~MeshNodeBoundary() = default;
    };
    
    
    template <size_t dimension>
    class MeshFlatNodeBoundary
        : public MeshNodeBoundary<dimension>
    {
     public:
      using Rd = TinyVector<dimension, double>;
    
     private:
      const Rd m_outgoing_normal;
    
      template <typename MeshType>
      inline Rd _getNormal(const MeshType& mesh);
    
      template <typename MeshType>
      inline void _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
                                       const TinyVector<2,double>& xmin,
                                       const TinyVector<2,double>& xmax,
                                       const MeshType& mesh) const;
    
      template <typename MeshType>
      inline Rd _getOutgoingNormal(const MeshType& mesh);
     public:
      const Rd& outgoingNormal() const
      {
        return m_outgoing_normal;
      }
    
      MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
      MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
    
      template <typename MeshType>
      MeshFlatNodeBoundary(const MeshType& mesh,
                           const RefFaceList& ref_face_list)
          : MeshNodeBoundary<dimension>(mesh, ref_face_list),
            m_outgoing_normal(_getOutgoingNormal(mesh))
      {
        ;
      }
    
      template <typename MeshType>
      MeshFlatNodeBoundary(const MeshType& mesh,
                           const RefNodeList& ref_node_list)
          : MeshNodeBoundary<dimension>(mesh, ref_node_list),
            m_outgoing_normal(_getOutgoingNormal(mesh))
      {
        ;
      }
    
      MeshFlatNodeBoundary() = default;
      MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
      MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
      virtual ~MeshFlatNodeBoundary() = default;
    };
    
    template<>
    template <typename MeshType>
    void MeshFlatNodeBoundary<2>::
    _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
                         const TinyVector<2,double>& xmin,
                         const TinyVector<2,double>& xmax,
                         const MeshType& mesh) const
    {
      static_assert(MeshType::dimension == 2);
      using R2 = TinyVector<2,double>;
    
      const R2 origin = 0.5*(xmin+xmax);
      const double length = l2Norm(xmax-xmin);
    
      const NodeValue<const R2>& xr = mesh.xr();
    
      Kokkos::parallel_for(m_node_list.size(), KOKKOS_LAMBDA(const size_t& r) {
          const R2& x = xr[m_node_list[r]];
          if ((x-origin,normal)>1E-13*length) {
            std::cerr << "this FlatBoundary is not flat!\n";
            std::exit(1);
          }
        });
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<1,double>
    MeshFlatNodeBoundary<1>::
    _getNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 1);
      using R = TinyVector<1,double>;
    
      if (m_node_list.size() != 1) {
        std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
        std::exit(1);
      }
    
      return R(1);
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<2,double>
    MeshFlatNodeBoundary<2>::
    _getNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 2);
      using R2 = TinyVector<2,double>;
    
      const NodeValue<const R2>& xr = mesh.xr();
    
      R2 xmin(std::numeric_limits<double>::max(),
              std::numeric_limits<double>::max());
    
      R2 xmax(-std::numeric_limits<double>::max(),
              -std::numeric_limits<double>::max());
    
      for (size_t r=0; r<m_node_list.size(); ++r) {
        const R2& x = xr[m_node_list[r]];
        if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
          xmin = x;
        }
        if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
          xmax = x;
        }
      }
    
      if (xmin == xmax) {
        std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
        std::exit(1);
      }
    
      R2 dx = xmax-xmin;
      dx *= 1./l2Norm(dx);
    
      R2 normal(-dx[1], dx[0]);
    
      this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
    
      return normal;
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<3,double>
    MeshFlatNodeBoundary<3>::
    _getNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 3);
      using R3 = TinyVector<3,double>;
    
    
      R3 xmin(std::numeric_limits<double>::max(),
              std::numeric_limits<double>::max(),
              std::numeric_limits<double>::max());
      R3 ymin = xmin;
      R3 zmin = xmin;;
    
      R3 xmax =-xmin;
      R3 ymax = xmax;
      R3 zmax = xmax;
    
      const NodeValue<const R3>& xr = mesh.xr();
    
      for (size_t r=0; r<m_node_list.size(); ++r) {
        const R3& x = xr[m_node_list[r]];
        if (x[0]<xmin[0]) {
          xmin = x;
        }
        if (x[1]<ymin[1]) {
          ymin = x;
        }
        if (x[2]<zmin[2]) {
          zmin = x;
        }
        if (x[0]>xmax[0]) {
          xmax = x;
        }
        if (x[1]>ymax[1]) {
          ymax = x;
        }
        if (x[2]>zmax[2]) {
          zmax = x;
        }
      }
    
      const R3 u = xmax-xmin;
      const R3 v = ymax-ymin;
      const R3 w = zmax-zmin;
    
      const R3 uv = crossProduct(u,v);
      const double uv_l2 = (uv, uv);
    
      R3 normal = uv;
      double normal_l2 = uv_l2;
    
      const R3 uw = crossProduct(u,w);
      const double uw_l2 = (uw,uw);
    
      if (uw_l2 > uv_l2) {
        normal = uw;
        normal_l2 = uw_l2;
      }
    
      const R3 vw = crossProduct(v,w);
      const double vw_l2 = (vw,vw);
    
      if (vw_l2 > normal_l2) {
        normal = vw;
        normal_l2 = vw_l2;
      }
    
    
      if (normal_l2 ==0) {
        std::cerr << "not able to compute normal!\n";
        std::exit(1);
      }
    
      normal *= 1./sqrt(normal_l2);
    
    #warning Add flatness test
      // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
    
      return normal;
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<1,double>
    MeshFlatNodeBoundary<1>::
    _getOutgoingNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 1);
      using R = TinyVector<1,double>;
    
      const R normal = this->_getNormal(mesh);
    
      const NodeValue<const R>& xr = mesh.xr();
      const auto& cell_to_node_matrix
          = mesh.connectivity().cellToNodeMatrix();
    
      const auto& node_to_cell_matrix
          = mesh.connectivity().nodeToCellMatrix();
    
      const NodeId r0 = m_node_list[0];
      const CellId j0 = node_to_cell_matrix[r0][0];
      const auto& j0_nodes = cell_to_node_matrix[j0];
      double max_height = 0;
      for (size_t r=0; r<j0_nodes.size(); ++r) {
        const double height = (xr[j0_nodes[r]]-xr[r0], normal);
        if (std::abs(height) > std::abs(max_height)) {
          max_height = height;
        }
      }
      if (max_height > 0) {
        return -normal;
      } else {
        return normal;
      }
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<2,double>
    MeshFlatNodeBoundary<2>::
    _getOutgoingNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 2);
      using R2 = TinyVector<2,double>;
    
      const R2 normal = this->_getNormal(mesh);
    
      const NodeValue<const R2>& xr = mesh.xr();
      const auto& cell_to_node_matrix
          = mesh.connectivity().cellToNodeMatrix();
    
      const auto& node_to_cell_matrix
          = mesh.connectivity().nodeToCellMatrix();
    
      const NodeId r0 = m_node_list[0];
      const CellId j0 = node_to_cell_matrix[r0][0];
      const auto& j0_nodes = cell_to_node_matrix[j0];
      double max_height = 0;
      for (size_t r=0; r<j0_nodes.size(); ++r) {
        const double height = (xr[j0_nodes[r]]-xr[r0], normal);
        if (std::abs(height) > std::abs(max_height)) {
          max_height = height;
        }
      }
      if (max_height > 0) {
        return -normal;
      } else {
        return normal;
      }
    }
    
    template <>
    template <typename MeshType>
    inline TinyVector<3,double>
    MeshFlatNodeBoundary<3>::
    _getOutgoingNormal(const MeshType& mesh)
    {
      static_assert(MeshType::dimension == 3);
      using R3 = TinyVector<3,double>;
    
      const R3 normal = this->_getNormal(mesh);
    
      const NodeValue<const R3>& xr = mesh.xr();
      const auto& cell_to_node_matrix
          = mesh.connectivity().cellToNodeMatrix();
    
      const auto& node_to_cell_matrix
          = mesh.connectivity().nodeToCellMatrix();
    
      const NodeId r0 = m_node_list[0];
      const CellId j0 = node_to_cell_matrix[r0][0];
      const auto& j0_nodes = cell_to_node_matrix[j0];
      double max_height = 0;
      for (size_t r=0; r<j0_nodes.size(); ++r) {
        const double height = (xr[j0_nodes[r]]-xr[r0], normal);
        if (std::abs(height) > std::abs(max_height)) {
          max_height = height;
        }
      }
      if (max_height > 0) {
        return -normal;
      } else {
        return normal;
      }
    }
    
    
    #endif // MESH_NODE_BOUNDARY_HPP