Elements.hpp

Go to the documentation of this file.

#pragma once
// Elements.hpp -- Finite element definitions for the Geom module.
//
// Architecture
// ------------
// Element metadata is split across three layers:
//
//   ElemEnum.hpp        ElemType / ParamSpace enumerations
//   ElementTraits.hpp   Per-element compile-time data (nodes, faces, coords,
//                       elevation spans, bisection topology, VTK mapping, ...)
//   Elements.hpp        (this file) Runtime API consumed by the rest of the
//                       codebase: thin wrappers that dispatch to ElementTraits
//                       via DispatchElementType(), the Element struct, shape
//                       function dispatch, and geometry helpers.
//
// Shape functions are auto-generated by tools/gen_shape_functions/ and live
// in per-element headers under Geom/Elements/<ElemName>.hpp.  Each header
// provides a full specialization of ShapeFuncImpl<ElemType> with static
// Diff0..Diff3 methods.  ShapeFunc_DiNj() resolves the runtime ElemType to
// the correct specialization via DispatchElementType().
//
// Remaining runtime-only logic
// ----------------------------
// GetO2ElemBisectVariant() -- chooses the shortest-diagonal bisection variant
// for Tet10 (3 variants) and Pyramid14 (2 variants) based on physical
// coordinates.  This cannot be made constexpr because it depends on runtime
// geometry.
 
#include <type_traits>
#include <cstdint>
#include <array>
#include <tuple>
#include "DNDS/HardEigen.hpp"
 
#include "DNDS/Defines.hpp"
#include "Geometric.hpp"
#include "BaseFunction.hpp"
#include "ElemEnum.hpp"
#include "ElementTraits.hpp"
 
// Auto-generated per-element shape function headers
#include "Geom/Elements/Line2.hpp"
#include "Geom/Elements/Line3.hpp"
#include "Geom/Elements/Tri3.hpp"
#include "Geom/Elements/Tri6.hpp"
#include "Geom/Elements/Quad4.hpp"
#include "Geom/Elements/Quad9.hpp"
#include "Geom/Elements/Tet4.hpp"
#include "Geom/Elements/Tet10.hpp"
#include "Geom/Elements/Hex8.hpp"
#include "Geom/Elements/Hex27.hpp"
#include "Geom/Elements/Prism6.hpp"
#include "Geom/Elements/Prism18.hpp"
#include "Geom/Elements/Pyramid5.hpp"
#include "Geom/Elements/Pyramid14.hpp"
 
namespace DNDS::Geom::Elem
{
 
    // ================================================================
    // Thin wrappers that delegate to ElementTraits via DispatchElementType.
    // These preserve the old free-function signatures for backward compat.
    // ================================================================
 
    DNDS_DEVICE_CALLABLE constexpr ParamSpace ElemType_to_ParamSpace(const ElemType t)
    {
        return DispatchElementType(t, [](auto tr) -> ParamSpace
                                   { return decltype(tr)::paramSpace; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr ElemType GetFaceType(ElemType t_v, t_index iFace)
    {
        return DispatchElementType(t_v, [iFace](auto tr) -> ElemType
                                   { return decltype(tr)::GetFaceType(iFace); });
    }
 
    /// Get the number of edges for a 3D element. Returns 0 for 1D/2D elements.
    DNDS_DEVICE_CALLABLE constexpr t_index GetNumEdges(ElemType t_v)
    {
        return DispatchElementType(t_v, [](auto tr) -> t_index
                                   {
            using T = decltype(tr);
            if constexpr (T::dim >= 3)
                return T::numEdges;
            else
                return 0; });
    }
 
    /// Get the element type of edge iEdge for a 3D element. Returns UnknownElem for 1D/2D.
    DNDS_DEVICE_CALLABLE constexpr ElemType GetEdgeType(ElemType t_v, t_index iEdge)
    {
        return DispatchElementType(t_v, [iEdge](auto tr) -> ElemType
                                   {
            using T = decltype(tr);
            if constexpr (T::dim >= 3)
                return T::GetEdgeType(iEdge);
            else
                return UnknownElem; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr t_real ParamSpaceVol(ParamSpace ps)
    {
        return Elem::ParamSpaceVolume(ps);
    }
 
    DNDS_DEVICE_CALLABLE constexpr int GetElemElevation_O1O2_NumNode(ElemType t)
    {
        return DispatchElementType(t, [](auto tr) -> int
                                   { return decltype(tr)::numElevNodes; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr ElemType GetElemElevation_O1O2_ElevatedType(ElemType t)
    {
        return DispatchElementType(t, [](auto tr) -> ElemType
                                   { return decltype(tr)::elevatedType; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr ElemType GetElemElevation_O1O2_NodeSpanType(ElemType t, t_index ine)
    {
        return DispatchElementType(t, [ine](auto tr) -> ElemType
                                   {
            using T = decltype(tr);
            if constexpr (T::numElevNodes > 0)
                return T::elevNodeSpanTypes[ine];
            else
                return UnknownElem; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr ElemType GetElemO1(ParamSpace ps)
    {
        return ParamSpaceO1Elem(ps);
    }
 
    DNDS_DEVICE_CALLABLE constexpr t_index GetO2ElemBisectNum(ElemType t)
    {
        return DispatchElementType(t, [](auto tr) -> t_index
                                   {
            using T = decltype(tr);
            if constexpr (T::order == 2)
                return T::numBisect;
            else
                return 0; });
    }
 
    DNDS_DEVICE_CALLABLE constexpr ElemType GetO2ElemBisectElem(ElemType t, t_index i)
    {
        return DispatchElementType(t, [i](auto tr) -> ElemType
                                   {
            using T = decltype(tr);
            if constexpr (T::order == 2)
                return T::GetBisectElemType(i);
            else
                return UnknownElem; });
    }
 
    DNDS_DEVICE_CALLABLE static constexpr real _iipow(real x, int y)
    {
        if (y == 0)
            return 1.;
        if (y == 1)
            return x;
        if (y > 1)
            return x * _iipow(x, y - 1);
        if (y == -1)
            return 1. / x;
        if (y < -1)
            return 1. / x * _iipow(x, y + 1);
        return UnInitReal;
    }
 
    /// Primary template for per-element shape function implementations.
    /// Full specializations are auto-generated in Geom/Elements/<Elem>.hpp.
    template <ElemType>
    struct ShapeFuncImpl;
 
    /**
     * @brief Dispatch to per-element shape function code via ShapeFuncImpl<T>.
     *
     * Each ShapeFuncImpl specialization provides static Diff0..Diff3 methods.
     * DispatchElementType (from ElementTraits.hpp) resolves the runtime ElemType
     * to the correct compile-time specialization.
     */
    template <int diffOrder, class TPoint, class TArray>
    DNDS_DEVICE_CALLABLE void ShapeFunc_DiNj(ElemType t, const TPoint &p, TArray &&v)
    {
        static_assert(diffOrder == 0 || diffOrder == 1 || diffOrder == 2 || diffOrder == 3);
 
        DispatchElementType(t, [&](auto traits)
                            {
            constexpr ElemType eType = decltype(traits)::elemType;
            if constexpr (diffOrder == 0)
                ShapeFuncImpl<eType>::Diff0(p, std::forward<TArray>(v));
            else if constexpr (diffOrder == 1)
                ShapeFuncImpl<eType>::Diff1(p, std::forward<TArray>(v));
            else if constexpr (diffOrder == 2)
                ShapeFuncImpl<eType>::Diff2(p, std::forward<TArray>(v));
            else if constexpr (diffOrder == 3)
                ShapeFuncImpl<eType>::Diff3(p, std::forward<TArray>(v)); });
    }
 
    using tNj = Eigen::RowVector<t_real, Eigen::Dynamic>;
    using tD1Nj = Eigen::Matrix<t_real, 3, Eigen::Dynamic>;
    using tD01Nj = Eigen::Matrix<t_real, 4, Eigen::Dynamic>;
    using tDiNj = Eigen::Matrix<t_real, Eigen::Dynamic, Eigen::Dynamic>;
 
    struct Element
    {
        ElemType type = UnknownElem;
 
        DNDS_DEVICE_TRIVIAL_COPY_DEFINE(Element, Element)
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr ParamSpace GetParamSpace() const
        {
            return ElemType_to_ParamSpace(type);
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetDim() const
        {
            return DispatchElementType(type, [](auto tr) -> t_index
                                       { return decltype(tr)::dim; });
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetOrder() const
        {
            return DispatchElementType(type, [](auto tr) -> t_index
                                       { return decltype(tr)::order; });
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetNumVertices() const
        {
            return DispatchElementType(type, [](auto tr) -> t_index
                                       { return decltype(tr)::numVertices; });
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetNumNodes() const
        {
            return DispatchElementType(type, [](auto tr) -> t_index
                                       { return decltype(tr)::numNodes; });
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetNumFaces() const
        {
            return DispatchElementType(type, [](auto tr) -> t_index
                                       { return decltype(tr)::numFaces; });
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetNumElev_O1O2() const
        {
            return GetElemElevation_O1O2_NumNode(type);
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] t_index GetO2NumBisect() const
        {
            return GetO2ElemBisectNum(type);
        }
 
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr Element ObtainFace(t_index iFace) const
        {
            DNDS_assert(iFace < this->GetNumFaces());
            return Element{GetFaceType(type, iFace)};
        }
 
        /// Number of edges (3D elements only; returns 0 for 1D/2D).
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr t_index GetNumEdges() const
        {
            return Geom::Elem::GetNumEdges(type);
        }
 
        /// Edge sub-element type for edge iEdge (Line2 or Line3).
        DNDS_DEVICE_CALLABLE [[nodiscard]] constexpr Element ObtainEdge(t_index iEdge) const
        {
            DNDS_assert(iEdge < this->GetNumEdges());
            return Element{Geom::Elem::GetEdgeType(type, iEdge)};
        }
 
        /// Extract the global node indices of edge iEdge from a cell's node list.
        /// @warning Assumes edgeNodesOut has sufficient size.
        template <class TIn, class TOut>
        void ExtractEdgeNodes(t_index iEdge, const TIn &nodes, TOut &edgeNodesOut)
        {
            DNDS_assert(iEdge < this->GetNumEdges());
            DispatchElementType(type, [&](auto tr)
                                {
                using T = decltype(tr);
                if constexpr (T::dim >= 3)
                {
                    for (t_index i = 0; i < ObtainEdge(iEdge).GetNumNodes(); i++)
                        edgeNodesOut[i] = nodes[T::edgeNodes[iEdge][i]];
                } });
        }
 
        /**
         * @warning assuming Out has correct size
         */
        template <class TIn, class TOut>
        void ExtractFaceNodes(t_index iFace, const TIn &nodes, TOut &faceNodesOut)
        {
            DNDS_assert(iFace < this->GetNumFaces());
            DispatchElementType(type, [&](auto tr)
                                {
                using T = decltype(tr);
                if constexpr (T::numFaces > 0)
                {
                    for (t_index i = 0; i < ObtainFace(iFace).GetNumNodes(); i++)
                        faceNodesOut[i] = nodes[T::faceNodes[iFace][i]];
                } });
        }
 
        [[nodiscard]] constexpr Element ObtainElevNodeSpan(t_index iNodeElev) const
        {
            DNDS_assert(iNodeElev < this->GetNumElev_O1O2());
            return Element{GetElemElevation_O1O2_NodeSpanType(type, iNodeElev)};
        }
 
        [[nodiscard]] constexpr Element ObtainElevatedElem() const
        {
            return Element{GetElemElevation_O1O2_ElevatedType(type)};
        }
 
        [[nodiscard]] constexpr Element ObtainO1Elem() const
        {
            return Element{GetElemO1(GetParamSpace())};
        }
 
        [[nodiscard]] Element ObtainO2BisectElem(t_index iSubElem) const
        {
            DNDS_assert(iSubElem < this->GetO2NumBisect());
            return Element{GetO2ElemBisectElem(type, iSubElem)};
        }
 
        template <class TIn, class TOut>
        void ExtractElevNodeSpanNodes(t_index iNodeElev, const TIn &nodes, TOut &spanNodes)
        {
            DNDS_assert(iNodeElev < this->GetNumElev_O1O2());
            DispatchElementType(type, [&](auto tr)
                                {
                using T = decltype(tr);
                if constexpr (T::numElevNodes > 0)
                {
                    auto spanElem = Element{T::elevNodeSpanTypes[iNodeElev]};
                    for (t_index i = 0; i < spanElem.GetNumNodes(); i++)
                        spanNodes[i] = nodes[T::elevSpans[iNodeElev][i]];
                } });
        }
 
        template <class TIn, class TOut>
        void ExtractO2BisectElemNodes(t_index iSubElem, t_index iVariant, const TIn &nodes, TOut &subNodes)
        {
            DNDS_assert(iSubElem < this->GetO2NumBisect());
            DispatchElementType(type, [&](auto tr)
                                {
                using T = decltype(tr);
                if constexpr (T::order == 2)
                {
                    DNDS_assert(iVariant < T::numBisectVariants);
                    t_index offset = iSubElem + T::numBisect * iVariant;
                    auto subElem = Element{T::GetBisectElemType(iSubElem)};
                    for (t_index i = 0; i < subElem.GetNumNodes(); i++)
                        subNodes[i] = nodes[T::bisectElements[offset][i]];
                } });
        }
 
        /**
         * @warning Nj resized within
         */
        DNDS_DEVICE_CALLABLE void GetNj(const tPoint &pParam, tNj &Nj) const
        {
            Nj.setZero(1, this->GetNumNodes());
            ShapeFunc_DiNj<0>(type, pParam, Nj);
        }
 
        /**
         * @warning D1Nj resized within
         */
        DNDS_DEVICE_CALLABLE void GetD1Nj(const tPoint &pParam, tD1Nj &D1Nj) const
        {
            D1Nj.setZero(3, this->GetNumNodes());
            ShapeFunc_DiNj<1>(type, pParam, D1Nj);
        }
 
        /**
         * @warning DiNj resized within
         */
        void GetD01Nj(const tPoint &pParam, tD01Nj &D01Nj) const
        {
            tNj Nj;
            tD1Nj D1Nj;
            this->GetNj(pParam, Nj);
            this->GetD1Nj(pParam, D1Nj);
 
            D01Nj.resize(4, this->GetNumNodes());
            D01Nj(0, EigenAll) = Nj;
            D01Nj({1, 2, 3}, EigenAll) = D1Nj;
        }
        /**
         * @warning DiNj resized within
         */
        void GetDiNj(const tPoint &pParam, tDiNj &DiNj, int maxOrder) const
        {
            if (this->GetDim() == 2)
            {
                int nDiffCur = Base::ndiffSizC2D.at(maxOrder);
                DiNj.setZero(nDiffCur, this->GetNumNodes());
                for (; maxOrder >= 0; maxOrder--)
                {
                    if (maxOrder == 3)
                        ShapeFunc_DiNj<3>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC2D.at(2), Base::ndiffSizC2D.at(3) - 1),
                                EigenAll));
                    if (maxOrder == 2)
                        ShapeFunc_DiNj<2>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC2D.at(1), Base::ndiffSizC2D.at(2) - 1),
                                EigenAll));
                    if (maxOrder == 1)
                        ShapeFunc_DiNj<1>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC2D.at(0), Base::ndiffSizC2D.at(1) - 1),
                                EigenAll));
                    if (maxOrder == 0)
                        ShapeFunc_DiNj<0>(
                            type, pParam,
                            DiNj(
                                0,
                                EigenAll));
                }
            }
            else if (this->GetDim() == 3)
            {
                int nDiffCur = Base::ndiffSizC.at(maxOrder);
                DiNj.setZero(nDiffCur, this->GetNumNodes());
                for (; maxOrder >= 0; maxOrder--)
                {
                    if (maxOrder == 3)
                        ShapeFunc_DiNj<3>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC.at(2), Base::ndiffSizC.at(3) - 1),
                                EigenAll));
                    if (maxOrder == 2)
                        ShapeFunc_DiNj<2>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC.at(1), Base::ndiffSizC.at(2) - 1),
                                EigenAll));
                    if (maxOrder == 1)
                        ShapeFunc_DiNj<1>(
                            type, pParam,
                            DiNj(
                                Eigen::seq(Base::ndiffSizC.at(0), Base::ndiffSizC.at(1) - 1),
                                EigenAll));
                    if (maxOrder == 0)
                        ShapeFunc_DiNj<0>(
                            type, pParam,
                            DiNj(
                                0,
                                EigenAll));
                }
            }
            else if (this->GetDim() == 1)
            {
                int nDiffCur = maxOrder + 1;
                DiNj.setZero(nDiffCur, this->GetNumNodes());
                for (; maxOrder >= 0; maxOrder--)
                {
                    if (maxOrder == 3)
                        ShapeFunc_DiNj<3>(
                            type, pParam,
                            DiNj(
                                3,
                                EigenAll));
                    if (maxOrder == 2)
                        ShapeFunc_DiNj<2>(
                            type, pParam,
                            DiNj(
                                2,
                                EigenAll));
                    if (maxOrder == 1)
                        ShapeFunc_DiNj<1>(
                            type, pParam,
                            DiNj(
                                1,
                                EigenAll));
                    if (maxOrder == 0)
                        ShapeFunc_DiNj<0>(
                            type, pParam,
                            DiNj(
                                0,
                                EigenAll));
                }
            }
        }
    };
 
    Eigen::Matrix<t_real, 3, Eigen::Dynamic> GetStandardCoord(ElemType t);
}
 
namespace DNDS::Geom::Elem
{
    inline bool cellsAreFaceConnected(
        const std::vector<DNDS::index> &nodes_A,
        const std::vector<DNDS::index> &nodes_B,
        Element eA,
        Element eB)
    {
        DNDS_assert(nodes_A.size() >= eA.GetNumNodes());
        DNDS_assert(nodes_B.size() >= eB.GetNumVertices());
        std::vector<DNDS::index> nodes_B_Vert{nodes_B.begin(), nodes_B.begin() + eB.GetNumVertices()};
        std::sort(nodes_B_Vert.begin(), nodes_B_Vert.end());
        for (int iF = 0; iF < eA.GetNumFaces(); iF++)
        {
            auto eF = eA.ObtainFace(iF);
            std::vector<DNDS::index> fNodes(eF.GetNumNodes());
            eA.ExtractFaceNodes(iF, nodes_A, fNodes);
            std::sort(fNodes.begin(), fNodes.begin() + eF.GetNumVertices());
            if (std::includes(
                    nodes_B_Vert.begin(), nodes_B_Vert.end(),
                    fNodes.begin(), fNodes.begin() + eF.GetNumVertices()))
                return true;
        }
        return false;
    }
 
    inline bool cellsAreFaceConnected(
        const std::vector<DNDS::index> &verts_A,
        const std::vector<DNDS::index> &nodes_B)
    {
 
        return false;
    }
 
    template <class tCoordsIn>
    tJacobi ShapeJacobianCoordD01Nj(const tCoordsIn &cs, Eigen::Ref<const tD01Nj> DiNj)
    {
        return cs * DiNj({1, 2, 3}, EigenAll).transpose();
    }
 
    template <class tCoordsIn>
    tPoint PPhysicsCoordD01Nj(const tCoordsIn &cs, Eigen::Ref<const tD01Nj> DiNj)
    {
        return cs * DiNj(0, EigenAll).transpose();
    }
 
    template <int dim>
    real JacobiDetFace(Eigen::Ref<const tJacobi> J)
    {
        static_assert(dim == 2 || dim == 3);
        if constexpr (dim == 2)
            return J(EigenAll, 0).stableNorm();
        else
            return J(EigenAll, 0).cross(J(EigenAll, 1)).stableNorm();
    }
 
    // TODO: add a integration-based counterpart
    inline tPoint GetElemNodeMajorSpan(const tSmallCoords &coords)
    {
        tPoint c = coords.rowwise().mean();
        tSmallCoords coordsC = coords.colwise() - c;
        tPoint ve0 = coordsC(EigenAll, 1) - coordsC(EigenAll, 0);
        tJacobi inertia = coordsC * coordsC.transpose();
        real cond01 = HardEigen::Eigen3x3RealSymEigenDecompositionGetCond01(inertia); // ratio of 2 largest eigenvalues
        if (cond01 < 1 + 1e-6)
            inertia += ve0 * ve0.transpose() * 1e-4; // first edge gets priority
        auto decRet = HardEigen::Eigen3x3RealSymEigenDecompositionNormalized(inertia);
        coordsC = decRet.transpose() * coordsC;
        tPoint ret = coordsC.rowwise().maxCoeff() - coordsC.rowwise().minCoeff();
        std::sort(ret.begin(), ret.end(), std::greater_equal<real>());
        return ret;
    }
 
    inline t_index GetO2ElemBisectVariant(Element e, const tSmallCoords &coords)
    {
        DNDS_assert(e.GetOrder() == 2);
        //! if serendipity ?
        switch (e.type)
        {
        case Tet10:
        {
            Eigen::Vector<real, 3> lengths{
                (coords(EigenAll, 5 - 1) - coords(EigenAll, 10 - 1)).norm(),
                (coords(EigenAll, 6 - 1) - coords(EigenAll, 8 - 1)).norm(),
                (coords(EigenAll, 7 - 1) - coords(EigenAll, 9 - 1)).norm(),
            };
            Eigen::Index iMin{-1};
            real minLen = lengths.minCoeff(&iMin);
            (void)minLen;
            return t_index(iMin);
        }
        break;
        case Pyramid14:
        {
            Eigen::Vector<real, 2> lengths{
                (coords(EigenAll, 10 - 1) - coords(EigenAll, 12 - 1)).norm(),
                (coords(EigenAll, 11 - 1) - coords(EigenAll, 13 - 1)).norm(),
            };
            Eigen::Index iMin{-1};
            real minLen = lengths.minCoeff(&iMin);
            (void)minLen;
            return t_index(iMin);
        }
        break;
        default:
            return 0;
        }
    }
 
    template <class TIn>
    std::pair<int, std::vector<index>> ToVTKVertsAndData(Element e, const TIn &vin)
    {
        return DispatchElementType(e.type, [&](auto tr) -> std::pair<int, std::vector<index>>
                                   {
            using T = decltype(tr);
            constexpr int nOut = static_cast<int>(
                std::tuple_size<decltype(T::vtkNodeOrder)>::value);
            std::vector<index> v(nOut);
            for (int i = 0; i < nOut; i++)
                v[i] = vin[T::vtkNodeOrder[i]];
            return {T::vtkCellType, std::move(v)}; });
    }
}