patchkernel: add a method to find all the closest cells to a point

src/patchkernel/patch_skd_tree.cpp

@@ -729,6 +729,47 @@ void SkdNode::updatePointClosestCell(const std::array<double, 3> &point, bool in
     }
 }
 
+/*!
+* Given the specified point find if, among the cells contained in the
+* bounding box associated to the node, there are cells closer to the
+* one received in input.
+*
+* In output all the cells at minimal distance from the given point are
+* provided.
+*
+* \param point is the point
+* \param interiorCellsOnly if set to true, only interior cells will be considered,
+* it will be possible to consider non-interior cells only if the tree has been
+* instantiated with non-interior cells support enabled
+* \param[in,out] closestIds are the indices of the closest cells, on output it will
+* be updated if the specified cell is closer than the current closest cells or
+* it is at the same distance.
+* \param[in,out] closestDistance is the distance of the current closest cells,
+* on output it will be updated if closer cells are found
+*/
+void SkdNode::updatePointClosestCells(const std::array<double, 3> &point, bool interiorCellsOnly,
+                                     std::vector<long> &closestId, double *closestDistance) const
+{
+    if (getCellCount() == 0) {
+        return;
+    }
+
+    const PatchKernel &patch = m_patchInfo->getPatch();
+    const PiercedVector<Cell> &cells = patch.getCells();
+    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
+
+    for (std::size_t n = m_cellRangeBegin; n < m_cellRangeEnd; n++) {
+        std::size_t cellRawId = cellRawIds[n];
+        const Cell &cell = cells.rawAt(cellRawId);
+        if (interiorCellsOnly && !cell.isInterior()) {
+            continue;
+        }
+
+
+        updatePointClosestCells(point, cell, closestId, closestDistance);
+    }
+}
+
 /*!
 * Given the specified point find if the given cell is closer than the current
 * one.
@@ -850,6 +891,98 @@ void SkdNode::updatePointClosestCell(const std::array<double, 3> &point, const C
     }
 }
 
+/*!
+* Given the specified point find if the given cell is closer than the current
+* one.
+*
+* If the specified cell has the same distance of the closest cells, this cell
+* is added to the closest cells container
+*
+* \param point is the point
+* \param cell is the cell
+* \param[in,out] closestIds are the indices of the closest cells, on output it will
+* be updated if the specified cell is closer than the current closest cells or
+* it is at the same distance.
+* \param[in,out] closestDistance is the distance of the closest cells, on output
+* it will be updated if the specified cell is closer than the current closest
+* cells
+*/
+void SkdNode::updatePointClosestCells(const std::array<double, 3> &point, const Cell &cell,
+        std::vector<long> &closestIds, double *closestDistance) const
+{
+    const PatchKernel &patch = m_patchInfo->getPatch();
+
+    // Cell id
+    long cellId = cell.getId();
+
+    // Project the point on the cell
+    int nCellVertices = cell.getVertexCount();
+    BITPIT_CREATE_WORKSPACE(cellVertexCoordinates, std::array<double BITPIT_COMMA 3>, nCellVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
+    patch.getElementVertexCoordinates(cell, cellVertexCoordinates);
+
+    double cellDistance;
+    std::array<double, 3> cellProjection;
+    cell.evalPointProjection(point, cellVertexCoordinates, &cellProjection, &cellDistance);
+
+    // Detect if the specified cell is closer than the current closest cells
+    const int DISTANCE_CLOSER  = - 1;
+    const int DISTANCE_EQUAL   =   0;
+    const int DISTANCE_FARTHER =   1;
+
+    int distanceFlag = DISTANCE_FARTHER;
+    if (utils::DoubleFloatingEqual()(cellDistance, *closestDistance, patch.getTol(), patch.getTol())) {
+        distanceFlag = DISTANCE_EQUAL;
+    } else if (cellDistance < *closestDistance) {
+        distanceFlag = DISTANCE_CLOSER;
+    }
+
+    // Consider the case where no closest cell is defined
+    //
+    // Even if the id of the closest cell is null, we may have an
+    // estimated of the closest cell distance. We need to update
+    // the closest cell information only if the current cells is
+    // closer than the estimate.
+    if (closestIds[0] == Cell::NULL_ID) {
+        if (distanceFlag != DISTANCE_FARTHER) {
+            closestIds[0]       = cellId;
+            *closestDistance = cellDistance;
+        }
+
+        return;
+    }
+
+    // Update closest cells information accordingly to the distance flag:
+    //  - if the specified cell is closer the specified cell will become the
+    //    closest cell
+    //  - if the two cells have the same distance, the specified cell will be
+    //    added to the closest cells container.
+    //  - if the specified cell is farther there is nothing to do;
+    switch (distanceFlag) {
+
+    case DISTANCE_CLOSER:
+    {
+        closestIds.resize(1);
+        closestIds[0]       = cellId;
+        *closestDistance = cellDistance;
+
+        break;
+    }
+
+    case DISTANCE_EQUAL:
+    {
+        closestIds.push_back(cellId);
+
+        break;
+    }
+
+    default:
+    {
+        // Nothing to do
+    }
+
+    }
+}
+
 #if BITPIT_ENABLE_MPI
 /*!
 * \class SkdGlobalCellDistance

src/patchkernel/patch_skd_tree.hpp

@@ -127,6 +127,7 @@ friend class PatchSkdTree;
 
     void findPointClosestCell(const std::array<double, 3> &point, bool interiorCellsOnly, long *closestId, double *closestDistance) const;
     void updatePointClosestCell(const std::array<double, 3> &point, bool interiorCellsOnly, long *closestId, double *closestDistance) const;
+    void updatePointClosestCells(const std::array<double, 3> &point, bool interiorCellsOnly, std::vector<long> &closestId, double *closestDistance) const;
 
 private:
     const SkdPatchInfo *m_patchInfo;
@@ -139,6 +140,7 @@ friend class PatchSkdTree;
     void initializeBoundingBox();
 
     void updatePointClosestCell(const std::array<double, 3> &point, const Cell &cell, long *closestId, double *closestDistance) const;
+    void updatePointClosestCells(const std::array<double, 3> &point, const Cell &cell, std::vector<long> &closestId, double *closestDistance) const;
 
 };