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

main.cpp

Blame
  • MeshNodeBoundary.hpp 15.13 KiB
    #ifndef MESH_NODE_BOUNDARY_HPP
    #define MESH_NODE_BOUNDARY_HPP
    
    #include <Kokkos_Vector.hpp>
    #include <algebra/TinyVector.hpp>
    #include <mesh/ConnectivityMatrix.hpp>
    #include <mesh/IConnectivity.hpp>
    #include <mesh/ItemValue.hpp>
    #include <mesh/RefItemList.hpp>
    #include <utils/Array.hpp>
    #include <utils/Exceptions.hpp>
    #include <utils/Messenger.hpp>
    
    template <size_t Dimension>
    class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
    {
     protected:
      Array<const NodeId> m_node_list;
      std::string m_boundary_name;
    
     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)
        : m_boundary_name(ref_face_list.refId().tagName())
      {
        static_assert(Dimension == MeshType::Dimension);
        const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
    
        const Array<const FaceId>& face_list = ref_face_list.list();
        parallel_for(
          face_list.size(), PUGS_LAMBDA(int l) {
            const auto& face_cells = face_to_cell_matrix[face_list[l]];
            if (face_cells.size() > 1) {
              std::ostringstream ost;
              ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
                  << ": inner faces cannot be used to define mesh boundaries";
              throw NormalError(ost.str());
            }
          });
    
        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 FaceId 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());
        parallel_for(
          node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
        m_node_list = node_list;
      }
    
      template <typename MeshType>
      MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
        : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
    
      {
        static_assert(Dimension == MeshType::Dimension);
      }
    
      MeshNodeBoundary()          = default;
      virtual ~MeshNodeBoundary() = default;
    
     protected:
      MeshNodeBoundary(const MeshNodeBoundary&) = default;
      MeshNodeBoundary(MeshNodeBoundary&&)      = default;
    };
    
    template <size_t Dimension>
    class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
    {
     public:
      using Rd = TinyVector<Dimension, double>;
    
     private:
      const Rd m_outgoing_normal;
    
      template <typename MeshType>
      PUGS_INLINE Rd _getNormal(const MeshType& mesh);
    
      template <typename MeshType>
      PUGS_INLINE void
      _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
                           const TinyVector<Dimension, double>& origin,
                           const double length,
                           const MeshType& mesh) const
      {
        static_assert(MeshType::Dimension == Dimension);
    
        const NodeValue<const Rd>& xr = mesh.xr();
    
        parallel_for(
          this->m_node_list.size(), PUGS_LAMBDA(const size_t r) {
            const Rd& x = xr[this->m_node_list[r]];
            if (dot(x - origin, normal) > 1E-13 * length) {
              std::ostringstream ost;
              ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
                  << ": boundary is not flat!";
              throw NormalError(ost.str());
            }
          });
      }
    
      template <typename MeshType>
      PUGS_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>
    PUGS_INLINE TinyVector<1, double>
    MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
    {
      static_assert(MeshType::Dimension == 1);
      using R = TinyVector<1, double>;
    
      const size_t number_of_bc_nodes = [&]() {
        size_t number_of_bc_nodes = 0;
        auto node_is_owned        = mesh.connectivity().nodeIsOwned();
        for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
          number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
        }
        return parallel::allReduceMax(number_of_bc_nodes);
      }();
    
      if (number_of_bc_nodes != 1) {
        std::ostringstream ost;
        ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
            << ": node boundaries in 1D require to have exactly 1 node";
        throw NormalError(ost.str());
      }
    
      return R{1};
    }
    
    template <>
    template <typename MeshType>
    PUGS_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());
    
      auto update_xmin = [](const R2& x, R2& xmin) {
        if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
          xmin = x;
        }
      };
    
      auto update_xmax = [](const R2& x, R2& xmax) {
        if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
          xmax = x;
        }
      };
    
      for (size_t r = 0; r < m_node_list.size(); ++r) {
        const R2& x = xr[m_node_list[r]];
        update_xmin(x, xmin);
        update_xmax(x, xmax);
      }
    
      if (parallel::size() > 1) {
        Array<R2> xmin_array = parallel::allGather(xmin);
        Array<R2> xmax_array = parallel::allGather(xmax);
        for (size_t i = 0; i < xmin_array.size(); ++i) {
          update_xmin(xmin_array[i], xmin);
        }
        for (size_t i = 0; i < xmax_array.size(); ++i) {
          update_xmax(xmax_array[i], xmax);
        }
      }
    
      if (xmin == xmax) {
        std::ostringstream ost;
        ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
            << ": unable to compute normal";
        throw NormalError(ost.str());
      }
    
      R2 dx = xmax - xmin;
      dx *= 1. / l2Norm(dx);
    
      const R2 normal(-dx[1], dx[0]);
    
      this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
    
      return normal;
    }
    
    template <>
    template <typename MeshType>
    PUGS_INLINE TinyVector<3, double>
    MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
    {
      static_assert(MeshType::Dimension == 3);
      using R3 = TinyVector<3, double>;
    
      auto update_xmin = [](const R3& x, R3& xmin) {
        // XMIN: X.xmin X.ymax X.zmax
        if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
            ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
          xmin = x;
        }
      };
    
      auto update_xmax = [](const R3& x, R3& xmax) {
        // XMAX: X.xmax X.ymin X.zmin
        if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
            ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
          xmax = x;
        }
      };
    
      auto update_ymin = [](const R3& x, R3& ymin) {
        // YMIN: X.ymin X.zmax X.xmin
        if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
            ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
          ymin = x;
        }
      };
    
      auto update_ymax = [](const R3& x, R3& ymax) {
        // YMAX: X.ymax X.zmin X.xmax
        if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
            ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
          ymax = x;
        }
      };
    
      auto update_zmin = [](const R3& x, R3& zmin) {
        // ZMIN: X.zmin X.xmin X.ymin
        if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
            ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
          zmin = x;
        }
      };
    
      auto update_zmax = [](const R3& x, R3& zmax) {
        // ZMAX: X.zmax X.xmax X.ymax
        if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
            ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
          zmax = x;
        }
      };
    
      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]];
        update_xmin(x, xmin);
        update_xmax(x, xmax);
        update_ymin(x, ymin);
        update_ymax(x, ymax);
        update_zmin(x, zmin);
        update_zmax(x, zmax);
      }
    
      if (parallel::size() > 0) {
        Array<const R3> xmin_array = parallel::allGather(xmin);
        Array<const R3> xmax_array = parallel::allGather(xmax);
        Array<const R3> ymin_array = parallel::allGather(ymin);
        Array<const R3> ymax_array = parallel::allGather(ymax);
        Array<const R3> zmin_array = parallel::allGather(zmin);
        Array<const R3> zmax_array = parallel::allGather(zmax);
    
        for (size_t i = 0; i < xmin_array.size(); ++i) {
          update_xmin(xmin_array[i], xmin);
        }
        for (size_t i = 0; i < ymin_array.size(); ++i) {
          update_ymin(ymin_array[i], ymin);
        }
        for (size_t i = 0; i < zmin_array.size(); ++i) {
          update_zmin(zmin_array[i], zmin);
        }
        for (size_t i = 0; i < xmax_array.size(); ++i) {
          update_xmax(xmax_array[i], xmax);
        }
        for (size_t i = 0; i < ymax_array.size(); ++i) {
          update_ymax(ymax_array[i], ymax);
        }
        for (size_t i = 0; i < zmax_array.size(); ++i) {
          update_zmax(zmax_array[i], zmax);
        }
      }
    
      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 = dot(uv, uv);
    
      R3 normal        = uv;
      double normal_l2 = uv_l2;
    
      const R3 uw        = crossProduct(u, w);
      const double uw_l2 = dot(uw, uw);
    
      if (uw_l2 > uv_l2) {
        normal    = uw;
        normal_l2 = uw_l2;
      }
    
      const R3 vw        = crossProduct(v, w);
      const double vw_l2 = dot(vw, vw);
    
      if (vw_l2 > normal_l2) {
        normal    = vw;
        normal_l2 = vw_l2;
      }
    
      if (normal_l2 == 0) {
        std::ostringstream ost;
        ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
            << ": unable to compute normal";
        throw NormalError(ost.str());
      }
    
      const double length = sqrt(normal_l2);
    
      normal *= 1. / length;
    
      this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
    
      return normal;
    }
    
    template <>
    template <typename MeshType>
    PUGS_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);
    
      double max_height = 0;
    
      if (m_node_list.size() > 0) {
        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];
    
        for (size_t r = 0; r < j0_nodes.size(); ++r) {
          const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
          if (std::abs(height) > std::abs(max_height)) {
            max_height = height;
          }
        }
      }
    
      Array<double> max_height_array = parallel::allGather(max_height);
      for (size_t i = 0; i < max_height_array.size(); ++i) {
        const double height = max_height_array[i];
        if (std::abs(height) > std::abs(max_height)) {
          max_height = height;
        }
      }
    
      if (max_height > 0) {
        return -normal;
      } else {
        return normal;
      }
    }
    
    template <>
    template <typename MeshType>
    PUGS_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);
    
      double max_height = 0;
    
      if (m_node_list.size() > 0) {
        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];
        for (size_t r = 0; r < j0_nodes.size(); ++r) {
          const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
          if (std::abs(height) > std::abs(max_height)) {
            max_height = height;
          }
        }
      }
    
      Array<double> max_height_array = parallel::allGather(max_height);
      for (size_t i = 0; i < max_height_array.size(); ++i) {
        const double height = max_height_array[i];
        if (std::abs(height) > std::abs(max_height)) {
          max_height = height;
        }
      }
    
      if (max_height > 0) {
        return -normal;
      } else {
        return normal;
      }
    }
    
    template <>
    template <typename MeshType>
    PUGS_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);
    
      double max_height = 0;
    
      if (m_node_list.size() > 0) {
        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];
    
        for (size_t r = 0; r < j0_nodes.size(); ++r) {
          const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
          if (std::abs(height) > std::abs(max_height)) {
            max_height = height;
          }
        }
      }
    
      Array<double> max_height_array = parallel::allGather(max_height);
      for (size_t i = 0; i < max_height_array.size(); ++i) {
        const double height = max_height_array[i];
        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