EBGeometry_DCEL_Reconcile.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2026 Robert Marskar <robert.marskar@sintef.no>
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef EBGEOMETRY_DCEL_RECONCILE_HPP
#define EBGEOMETRY_DCEL_RECONCILE_HPP
 
// Std includes
#include <cstddef>
#include <cstdint>
 
// Our includes
#include "EBGeometry_DCEL_Mesh.hpp"
#include "EBGeometry_DCEL_Vertex.hpp"
#include "EBGeometry_GPUTransform.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_PODVector.hpp"
#include "EBGeometry_Pool.hpp"
 
namespace EBGeometry {
 
namespace DCEL {
 
namespace ReconcileDetail {
 
struct Corner
{
  uint32_t m_face;
 
  uint32_t m_first;
 
  uint32_t m_second;
};
 
} // namespace ReconcileDetail
 
class ReconcileSchedule
{
public:
  template <class T>
  EBGEOMETRY_HOST
  ReconcileSchedule(Pool& a_pool, const MeshT<T>& a_mesh);
 
private:
  template <class T>
  friend void
  reconcile(Pool& a_pool, const MeshT<T>& a_mesh, const ReconcileSchedule& a_schedule);
 
  PODVector<uint32_t> m_cornerBegin;
 
  PODVector<ReconcileDetail::Corner> m_corners;
 
  PODVector<uint32_t> m_status;
 
  uint64_t m_rootId = 0;
 
  uint64_t m_vertexOffset = 0;
 
  uint64_t m_edgeOffset = 0;
 
  uint64_t m_faceOffset = 0;
 
  uint32_t m_numVertices = 0;
 
  uint32_t m_numEdges = 0;
 
  uint32_t m_numFaces = 0;
};
 
template <class T>
EBGEOMETRY_HOST
inline void
reconcile(Pool& a_pool, const MeshT<T>& a_mesh, const ReconcileSchedule& a_schedule);
 
namespace ReconcileDetail {
 
template <class T>
struct Faces
{
  mutable MeshT<T> m_mesh;
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  operator()(const std::size_t a_i) const
  {
    m_mesh.getFace(static_cast<uint32_t>(a_i)).reconcile(m_mesh);
  }
};
 
template <class T>
struct Edges
{
  mutable MeshT<T> m_mesh;
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  operator()(const std::size_t a_i) const noexcept
  {
    m_mesh.getEdge(static_cast<uint32_t>(a_i)).reconcile(m_mesh);
  }
};
 
template <class T>
struct Vertices
{
  mutable MeshT<T> m_mesh;
 
  void* m_base;
 
  PODVector<uint32_t> m_cornerBegin;
 
  PODVector<Corner> m_corners;
 
  PODVector<uint32_t> m_status;
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  operator()(const std::size_t a_i) const noexcept
  {
    const auto     v     = static_cast<uint32_t>(a_i);
    const uint32_t begin = m_cornerBegin.at(m_base, v);
    const uint32_t end   = m_cornerBegin.at(m_base, v + 1);
 
    VertexT<T>&    vertex = m_mesh.getVertex(v);
    const Vec3T<T> x0     = vertex.getPosition();
 
    // A vertex no face uses gets a zero normal, as in MeshT::reconcile().
    Vec3T<T> normal = Vec3T<T>::zeros();
 
    for (uint32_t c = begin; c < end; c++) {
      const Corner&   corner = m_corners.at(m_base, c);
      const FaceT<T>& face   = m_mesh.getFace(corner.m_face);
 
      // A zero-area face contributes nothing, as in VertexT::computeVertexNormalAngleWeighted().
      if (face.getNormal().length2() == T(0)) {
        continue;
      }
 
      const Vec3T<T>& x1 = m_mesh.getVertex(corner.m_first).getPosition();
      const Vec3T<T>& x2 = m_mesh.getVertex(corner.m_second).getPosition();
 
      // Every thread that finds a coincidence writes the same value, so the order does not matter.
      if (x0 == x1 || x0 == x2 || x1 == x2) {
        m_status.at(m_base, 0) = 1;
 
        continue;
      }
 
      normal += VertexDetail::angleWeightedTerm(x0, x1, x2, face.getNormal());
    }
 
    vertex.setNormal(normal);
 
    if (end > begin) {
      vertex.normalizeNormalVector();
    }
  }
};
 
} // namespace ReconcileDetail
 
} // namespace DCEL
 
} // namespace EBGeometry
 
#include "EBGeometry_DCEL_ReconcileImplem.hpp"
 
#endif