Domain.h

#pragma once
 
#include <Caribou/config.h>
#include <Caribou/macros.h>
#include <Caribou/constants.h>
#include <Caribou/Topology/BaseDomain.h>
#include <Caribou/Topology/BarycentricContainer.h>
#include <Caribou/Geometry/Element.h>
 
#include <vector>
#include <array>
 
namespace caribou::topology {
 
    template <typename Element>
    class DomainStorage {};
 
    template <typename Mesh, typename Element, typename NodeIndex = UNSIGNED_INTEGER_TYPE>
    class Domain final : public BaseDomain, private DomainStorage<Element> {
        friend Mesh;
    public:
        static constexpr INTEGER_TYPE Dimension = geometry::traits<Element>::Dimension;
        using MeshType = Mesh;
        using ElementType = Element;
        using NodeIndexType = NodeIndex;
 
        using ElementsIndices = Eigen::Matrix<NodeIndex, Eigen::Dynamic, geometry::traits<Element>::NumberOfNodesAtCompileTime, Eigen::RowMajor>;
 
        using ElementIndices = Eigen::Matrix<NodeIndex, geometry::traits<Element>::NumberOfNodesAtCompileTime, 1>;
 
        Domain() = delete;
 
        Domain(const Domain & other) noexcept
        : DomainStorage<Element>(other), p_mesh(other.p_mesh), p_buffer(other.p_buffer), p_elements(other.p_elements) {}
 
//        /*! Move constructor */
//        Domain(Domain && other) noexcept {
//            swap(*this, other);
//        }
 
        auto operator=(Domain other) noexcept -> Domain & {
            swap(*this, other);
            return *this;
        }
 
        ~Domain() final = default;
 
        [[nodiscard]]
        auto canonical_dimension() const -> UNSIGNED_INTEGER_TYPE final {
            return geometry::traits<Element>::CanonicalDimension;
        }
 
        [[nodiscard]]
        auto number_of_nodes_per_elements() const -> UNSIGNED_INTEGER_TYPE final;
 
        [[nodiscard]]
        auto number_of_elements() const -> UNSIGNED_INTEGER_TYPE final;
 
         template<typename EigenMatrix>
        inline auto element(const UNSIGNED_INTEGER_TYPE & element_id, const Eigen::DenseBase<EigenMatrix> & positions) const -> Element;
 
        inline auto element(const UNSIGNED_INTEGER_TYPE & element_id) const -> Element;
 
        inline auto element_indices(const UNSIGNED_INTEGER_TYPE & index) const;
 
        inline auto mesh() const -> const Mesh & {
            return *p_mesh;
        }
 
        template <typename Derived>
        inline auto embed(const Eigen::MatrixBase<Derived> & points) const -> BarycentricContainer<Domain> {
            return {this, points};
        }
 
    protected:
 
        friend void swap(Domain & first, Domain& second) noexcept
        {
            // enable ADL
            using std::swap;
            swap(first.p_buffer, second.p_buffer);
            swap(first.p_elements, second.p_elements);
        }
 
        Domain(const Mesh * mesh, const ElementsIndices & elements)
            : p_mesh(mesh), p_buffer(elements), p_elements(p_buffer.data(), p_buffer.rows(), p_buffer.cols(), {p_buffer.cols(), 1}) {}
 
        Domain(const Mesh * mesh, ElementsIndices & elements)
            : p_mesh(mesh), p_buffer(elements), p_elements(p_buffer.data(), p_buffer.rows(), p_buffer.cols(), {p_buffer.cols(), 1}) {}
 
        Domain(const Mesh * mesh, const ElementsIndices * elements)
            : p_mesh(mesh), p_buffer(), p_elements(elements->data(), elements->rows(), elements->cols(), {elements->cols(), 1}) {}
 
        explicit Domain(const Mesh * mesh, const NodeIndex * data, const Eigen::Index & number_of_elements, const Eigen::Index & number_of_nodes_per_elements)
            : p_mesh(mesh), p_buffer(), p_elements(data, number_of_elements, number_of_nodes_per_elements, {number_of_nodes_per_elements, 1}) {}
 
        explicit Domain(const Mesh * mesh, const NodeIndex * data, const Eigen::Index & number_of_elements, const Eigen::Index & number_of_nodes_per_elements, Eigen::Index outer_stride, Eigen::Index inner_stride)
            : p_mesh(mesh), p_buffer(), p_elements(data, number_of_elements, number_of_nodes_per_elements, {outer_stride, inner_stride}) {}
 
        const Mesh * p_mesh;
 
        ElementsIndices p_buffer;
 
        Eigen::Map<const ElementsIndices, Eigen::Unaligned, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> p_elements;
    };
 
    // -------------------------------------------------------------------------------------
    //                                    IMPLEMENTATION
    // -------------------------------------------------------------------------------------
    // Implementation of methods that can be specialized later for an explicit Element type
    // -------------------------------------------------------------------------------------
 
    template <typename Mesh, typename Element, typename NodeIndex>
    auto Domain<Mesh, Element, NodeIndex>::number_of_nodes_per_elements() const -> UNSIGNED_INTEGER_TYPE {
        return p_elements.cols();
    }
 
    template <typename Mesh, typename Element, typename NodeIndex>
    auto Domain<Mesh, Element, NodeIndex>::number_of_elements() const -> UNSIGNED_INTEGER_TYPE {
        return p_elements.rows();
    }
 
    template <typename Mesh, typename Element, typename NodeIndex>
    inline auto Domain<Mesh, Element, NodeIndex>::element_indices(const UNSIGNED_INTEGER_TYPE & index) const {
        caribou_assert(index < number_of_elements() and "Trying to get an element that does not exists.");
 
        return Eigen::Map<const Eigen::Matrix<NodeIndex, geometry::traits<Element>::NumberOfNodesAtCompileTime, 1>, Eigen::Unaligned, Eigen::Stride<1, Eigen::Dynamic>>(
            p_elements.row(index).data(), {1, p_elements.innerStride()}
        );
    }
 
    template <typename Mesh, typename Element, typename NodeIndex>
    template<typename EigenMatrix>
    inline auto Domain<Mesh, Element, NodeIndex>::element(const UNSIGNED_INTEGER_TYPE & element_id, const Eigen::DenseBase<EigenMatrix> & positions) const -> Element {
        caribou_assert(element_id < number_of_elements() &&
                       "Trying to get the element #"+std::to_string(element_id) + ", but the domain only has " +
                       std::to_string(number_of_elements()) + " elements."
        );
 
        using NodeMatrix = typename geometry::Element<Element>::template Matrix<geometry::traits<Element>::NumberOfNodesAtCompileTime, Dimension>;
        NodeMatrix node_positions;
        if constexpr (geometry::traits<Element>::NumberOfNodesAtCompileTime == caribou::Dynamic) {
            node_positions.resize(number_of_nodes_per_elements(), Dimension);
        }
 
        const auto node_indices = element_indices(element_id);
        for (std::size_t i = 0; i < node_indices.size(); ++i) {
            node_positions.row(i) = positions.row(node_indices[i]);
        }
 
        return Element(node_positions);
    }
 
    template <typename Mesh, typename Element, typename NodeIndex>
    inline auto Domain<Mesh, Element, NodeIndex>::element(const UNSIGNED_INTEGER_TYPE & element_id) const -> Element {
        caribou_assert(element_id < number_of_elements() &&
                       ("Trying to get the element #"+std::to_string(element_id) + ", but the domain only has " +
                       std::to_string(number_of_elements()) + " elements.").c_str()
        );
 
        return Element(mesh().positions(element_indices(element_id)));
    }
}