BarycentricContainer.h

#pragma once
 
#include <memory>
#include <vector>
#include <unordered_map>
 
#include <Caribou/config.h>
#include <Caribou/macros.h>
#include <Caribou/constants.h>
#include <Caribou/Topology/Domain.h>
#include <Caribou/Topology/HashGrid.h>
 
namespace caribou::topology {
 
 template <typename Domain>
class BarycentricContainer {
public:
    using Mesh = typename Domain::MeshType;
    using ContainerElement = typename Domain::ElementType;
 
    static constexpr INTEGER_TYPE Dimension = ContainerElement::Dimension;
    using ElementIndex = INTEGER_TYPE;
    using LocalCoordinates = typename ContainerElement::LocalCoordinates;
    using WorldCoordinates = typename ContainerElement::WorldCoordinates;
    using HashGridT = HashGrid<ContainerElement>;
 
    struct BarycentricPoint {
        BarycentricPoint() : element_index(-1), local_coordinates(LocalCoordinates::Zero()) {}
        BarycentricPoint(const ElementIndex & index, const LocalCoordinates & coordinates)
        : element_index(index), local_coordinates(coordinates) {}
 
        ElementIndex element_index;
        LocalCoordinates local_coordinates;
    };
 
    BarycentricContainer() = delete;
 
    template <typename Derived>
    BarycentricContainer(const Domain * container_domain, const Eigen::MatrixBase<Derived> & embedded_points)
    : BarycentricContainer(container_domain) {
        // Set the embedded points
        set_embedded_points(embedded_points);
    }
 
    explicit BarycentricContainer(const Domain * container_domain)
        : p_container_domain(container_domain) {
        if (container_domain->number_of_elements() == 0) {
            throw std::runtime_error("Trying to create a barycentric container from an empty domain.");
        }
 
        // Get the mean size of the elements
        WorldCoordinates H_mean = WorldCoordinates::Zero();
        for (UNSIGNED_INTEGER_TYPE element_id = 0; element_id < container_domain->number_of_elements(); ++element_id) {
            const ContainerElement element = container_domain->element(element_id);
            const auto element_nodes = element.nodes();
            const auto min = element_nodes.colwise().minCoeff();
            const auto max = element_nodes.colwise().maxCoeff();
            H_mean += (max - min);
        }
        H_mean /= container_domain->number_of_elements();
 
        // Create the Hash grid
        p_hash_grid = std::make_unique<HashGridT>(H_mean.maxCoeff(), container_domain->number_of_elements());
 
        for (UNSIGNED_INTEGER_TYPE element_id = 0; element_id < container_domain->number_of_elements(); ++element_id) {
            p_hash_grid->add(container_domain->element(element_id), element_id);
        }
    }
 
    auto barycentric_point(const WorldCoordinates & p) const -> BarycentricPoint {
        // List of candidate elements that could contain the point
        const auto candidate_element_indices = p_hash_grid->get(p);
 
 
        for (const auto & element_index : candidate_element_indices) {
            const ContainerElement e = p_container_domain->element(element_index);
            const LocalCoordinates local_coordinates = e.local_coordinates(p);
            if (e.contains_local(local_coordinates)) {
                // Found one, let's return it
                return {static_cast<ElementIndex>(element_index), local_coordinates};
            }
        }
 
        // No elements are containing the queried point
        return {-1, LocalCoordinates::Zero()};
    }
 
    auto closest_elements(const WorldCoordinates & p) const -> std::vector<BarycentricPoint> {
        // List of candidate elements that could contain the point
        const auto candidate_element_indices = p_hash_grid->get(p);
 
        std::vector<BarycentricPoint> closest_elements;
        closest_elements.reserve(candidate_element_indices.size());
        for (const auto & element_index : candidate_element_indices) {
            const ContainerElement e = p_container_domain->element(element_index);
            const LocalCoordinates local_coordinates = e.local_coordinates(p);
            closest_elements.emplace_back(static_cast<ElementIndex>(element_index), local_coordinates);
        }
 
        return closest_elements;
    }
 
    template <typename Derived1, typename Derived2>
    void interpolate(const Eigen::MatrixBase<Derived1> & container_field_values,
                           Eigen::MatrixBase<Derived2> & embedded_field_values) const {
 
        // Make sure we have one field value per node of the container domain's mesh.
        if (static_cast<unsigned>(container_field_values.rows()) != p_container_domain->mesh().number_of_nodes()) {
            std::stringstream ss;
            ss << "The number of rows of the input matrix must be the same as the number of nodes in the container domain. ";
            ss << "The number of rows is " << container_field_values.rows() << " and the number of nodes is ";
            ss << p_container_domain->mesh().number_of_nodes();
            throw std::runtime_error(ss.str());
        }
 
        // Make sure we can write one field value per node of the embedded mesh.
        if (static_cast<unsigned>(embedded_field_values.rows()) != p_barycentric_points.size()) {
            throw std::runtime_error("The number of rows of the output matrix must be the same as the number of embedded nodes.");
        }
 
        if (container_field_values.cols() != embedded_field_values.cols()) {
            throw std::runtime_error("The number of columns of the input matrix (container field values) must be the "
                                     "same as the number of columns of the output matrix (embedded field values).");
        }
 
        // Interpolate field values
        const auto number_of_embedded_points = static_cast<signed >(p_barycentric_points.size());
        for (Eigen::Index node_id = 0; node_id < number_of_embedded_points; ++node_id) {
            const auto & element_index = p_barycentric_points[node_id].element_index;
            const auto & local_coordinates = p_barycentric_points[node_id].local_coordinates;
 
            if (element_index < 0) {
                continue;
            }
 
            const auto & element_node_indices = p_container_domain->element_indices(element_index);
            const auto element = ContainerElement(p_container_domain->mesh().positions(element_node_indices));
 
            // Get the field values at the nodes of the containing element
            Eigen::Matrix<typename Eigen::MatrixBase<Derived2>::Scalar,
                          ContainerElement::NumberOfNodesAtCompileTime,
                          Eigen::MatrixBase<Derived2>::ColsAtCompileTime>
                values (element.number_of_nodes(), container_field_values.cols());
 
            for (UNSIGNED_INTEGER_TYPE i = 0; i < element.number_of_nodes(); ++i) {
                values.row(i) = container_field_values.row(element_node_indices[i]);
            }
 
            // Interpolate the field value at the local position
            const auto interpolated_value = element.interpolate(local_coordinates, values);
 
            // Write back the interpolated value in the output embedded field
            if constexpr (Eigen::MatrixBase<Derived2>::ColsAtCompileTime == 1) {
                embedded_field_values[node_id] = interpolated_value;
            } else {
                embedded_field_values.row(node_id) = interpolated_value;
            }
        }
    }
 
    auto barycentric_points () const -> const std::vector<BarycentricPoint> & {
        return p_barycentric_points;
    }
 
    [[nodiscard]]
    auto outside_nodes () const -> const std::vector<UNSIGNED_INTEGER_TYPE> & {
        return p_outside_nodes;
    }
 
private:
    template <typename Derived>
    void set_embedded_points(const Eigen::MatrixBase<Derived> & embedded_points) {
        static_assert(Eigen::MatrixBase<Derived>::ColsAtCompileTime == Dimension or Eigen::MatrixBase<Derived>::ColsAtCompileTime == Eigen::Dynamic);
 
        if (Eigen::MatrixBase<Derived>::ColsAtCompileTime == Eigen::Dynamic and embedded_points.cols() != Dimension) {
            throw std::runtime_error("Trying to get the barycentric coordinates of " +
                                     std::to_string(embedded_points.cols()) + "D points from a " +
                                     std::to_string(Dimension) + "D domain");
        }
 
        p_outside_nodes.resize(0);
        p_outside_nodes.reserve(std::floor(embedded_points.rows() / 10.)); // Reserve 10% of the mesh size for outside nodes
 
        p_barycentric_points.resize(0);
        p_barycentric_points.reserve(embedded_points.rows());
 
        const auto number_of_embedded_points = embedded_points.rows();
        for (Eigen::Index node_id = 0; node_id < number_of_embedded_points; ++node_id) {
            const BarycentricPoint bp = barycentric_point(embedded_points.row(node_id).template cast<typename WorldCoordinates::Scalar>());
            if (bp.element_index < 0) {
                p_outside_nodes.emplace_back(node_id);
            }
 
            p_barycentric_points.emplace_back(bp);
        }
    }
 
 
    const Domain * p_container_domain;
    std::vector<BarycentricPoint> p_barycentric_points;
    std::vector<UNSIGNED_INTEGER_TYPE> p_outside_nodes;
    std::unique_ptr<HashGridT> p_hash_grid;
};
 
} // namespace caribou::topology