Element.h

#pragma once
 
#include <Caribou/config.h>
#include <Caribou/traits.h>
#include <Caribou/macros.h>
#include <Eigen/Dense>
#include <vector>
 
#include <iostream>
 
namespace caribou::geometry {
 
template<typename T>
struct traits;
 
template<class T, class Enable = void>
struct element_has_boundaries : std::false_type {};
 
template< class T>
constexpr bool element_has_boundaries_v = element_has_boundaries<T>::value;
 
template<typename Derived, typename ScalarType = FLOATING_POINT_TYPE>
struct Element {
    // Types
    using Scalar = ScalarType;
 
    template <INTEGER_TYPE Dim>
    using Vector = Eigen::Matrix<Scalar, Dim, 1>;
 
    template <INTEGER_TYPE Rows, INTEGER_TYPE Cols, int Options = Eigen::ColMajor>
    using Matrix = Eigen::Matrix<Scalar, Rows, Cols, Options>;
 
    template <INTEGER_TYPE Rows, INTEGER_TYPE Cols, int Options = Eigen::ColMajor>
    using MatrixI = Eigen::Matrix<Scalar, Rows, Cols, Options>;
 
    static constexpr auto CanonicalDimension = traits<Derived>::CanonicalDimension;
    static constexpr auto Dimension = traits<Derived>::Dimension;
    static constexpr auto NumberOfNodesAtCompileTime = traits<Derived>::NumberOfNodesAtCompileTime;
    static constexpr auto NumberOfGaussNodesAtCompileTime = traits<Derived>::NumberOfGaussNodesAtCompileTime;
    using LocalCoordinates = Vector<CanonicalDimension>;
    using WorldCoordinates = Vector<Dimension>;
 
    struct GaussNode {
        LocalCoordinates position;
        Scalar weight;
    };
 
 
    // Functions
    [[nodiscard]]
    inline auto number_of_nodes() const -> UNSIGNED_INTEGER_TYPE {return self().get_number_of_nodes();}
 
    [[nodiscard]]
    inline auto number_of_gauss_nodes() const -> UNSIGNED_INTEGER_TYPE {return self().get_number_of_gauss_nodes();}
 
    inline auto node(const UNSIGNED_INTEGER_TYPE & index) const {return self().get_node(index);};
 
    inline auto nodes() const {return self().get_nodes();};
 
    inline auto gauss_node(const UNSIGNED_INTEGER_TYPE & index) const -> const GaussNode & {return self().gauss_nodes()[index];};
 
    inline auto gauss_nodes() const -> const std::vector<GaussNode> & {return self().get_gauss_nodes();};
 
    inline auto number_of_boundary_elements() const {
        if constexpr (element_has_boundaries_v<Derived>) {
            return self().get_number_of_boundary_elements();
        } else {
            return 0;
        }
    }
 
    inline auto boundary_elements_node_indices() const -> const auto & {
        static_assert(element_has_boundaries_v<Derived>, "This element type has no boundary elements defined.");
        return self().get_boundary_elements_nodes();
    }
 
    inline auto boundary_element(const UNSIGNED_INTEGER_TYPE & boundary_id) const {
        static_assert(element_has_boundaries_v<Derived>, "This element type has no boundary elements defined.");
        using BoundaryElement = typename traits<Derived>::BoundaryElementType;
        const auto & node_indices = boundary_elements_node_indices()[boundary_id];
        Matrix<BoundaryElement::NumberOfNodesAtCompileTime, Dimension> nodes;
        for (Eigen::Index boundary_node_id = 0; boundary_node_id < nodes.rows(); ++boundary_node_id) {
            nodes.row(boundary_node_id) = self().node(node_indices[boundary_node_id]);
        }
        return BoundaryElement(nodes);
    }
 
    inline auto  L(const LocalCoordinates & xi) const -> Vector<NumberOfNodesAtCompileTime> {return self().get_L(xi);};
 
    inline auto dL(const LocalCoordinates & xi) const -> Matrix<NumberOfNodesAtCompileTime, CanonicalDimension> {return self().get_dL(xi);};
 
    inline auto center() const -> WorldCoordinates {return self().get_center();}
 
    inline auto world_coordinates(const LocalCoordinates & coordinates) const -> WorldCoordinates {
        return WorldCoordinates(self().interpolate(coordinates, self().nodes()));
    }
 
    inline auto local_coordinates(const WorldCoordinates & coordinates) const -> LocalCoordinates {
        return local_coordinates(coordinates, LocalCoordinates::Constant(0), 1e-5, 5);
    }
 
    inline auto local_coordinates(
        const WorldCoordinates & coordinates,
        const LocalCoordinates & starting_point,
        const FLOATING_POINT_TYPE & residual_tolerance,
        const UNSIGNED_INTEGER_TYPE & maximum_number_of_iterations) const -> LocalCoordinates {
 
        using namespace Eigen;
 
        LocalCoordinates xi = starting_point;
        UNSIGNED_INTEGER_TYPE iteration = 0;
        FLOATING_POINT_TYPE squared_threshold = residual_tolerance*residual_tolerance;
 
        // Initial residual
        WorldCoordinates residual = coordinates - world_coordinates(xi);
        FLOATING_POINT_TYPE r_norm_2 = residual.squaredNorm();
        FLOATING_POINT_TYPE r0_norm_2 = r_norm_2;
 
        if (r_norm_2 < squared_threshold) {
            // The initial guess is good enough
            return xi;
        }
 
        // Start the iterations
        do {
            Matrix<Dimension, CanonicalDimension> J = jacobian(xi);
 
            LocalCoordinates dxi;
            if constexpr (Dimension == CanonicalDimension) {
                dxi.noalias() = J.inverse() * residual;
            } else {
                dxi.noalias() = (J.transpose()*J).inverse() * (J.transpose()*residual);
            }
 
            xi.noalias() = (xi + dxi).eval();
            residual.noalias() = coordinates - world_coordinates(xi);
            r_norm_2 = residual.squaredNorm();
 
            ++iteration;
        }
        while (iteration < maximum_number_of_iterations and r_norm_2 >= r0_norm_2*squared_threshold);
 
        return xi;
    }
 
    inline auto contains_local(const LocalCoordinates & xi, const FLOATING_POINT_TYPE & eps = 1e-10) const -> bool {
        return self().get_contains_local(xi, eps);
    }
 
    template <typename MatrixType>
    inline auto interpolate (const LocalCoordinates & coordinates,
                             const Eigen::MatrixBase<MatrixType> & values) const {
        static_assert(Eigen::MatrixBase<MatrixType>::RowsAtCompileTime == NumberOfNodesAtCompileTime,
                      "The matrix containing the values at each nodes must have one node-value per row.");
        constexpr auto NbCols = Eigen::MatrixBase<MatrixType>::ColsAtCompileTime;
        using MatrixScalar = typename Eigen::MatrixBase<MatrixType>::Scalar;
        const auto result = ((values.array().colwise() * self().L(coordinates).array().template cast<typename Eigen::MatrixBase<MatrixType>::Scalar>()).matrix().colwise().sum().transpose()).eval();
        if constexpr (NbCols == 1) {
            return static_cast<MatrixScalar>(result[0]);
        } else {
            return result.template cast<MatrixScalar>();
        }
    }
 
    inline auto jacobian (const LocalCoordinates & coordinates) const -> Matrix<Dimension, CanonicalDimension>
    {
        const auto shape_derivatives = self().dL(coordinates);
        return self().nodes().transpose() * shape_derivatives;
    }
private:
    auto self() -> Derived& { return *static_cast<Derived*>(this); }
    auto self() const -> const Derived& { return *static_cast<const Derived*>(this); }
};
 
 
template <class T>
using element_boundary_type_t = typename traits<T>::BoundaryElementType;
 
template<class T>
struct element_has_boundaries<T, CLASS_REQUIRES(
    caribou::internal::is_detected_v<element_boundary_type_t, T>
)> : std::true_type {};
 
 
}  /// namespace caribou::geometry