Pyramid5.hpp

Go to the documentation of this file.

#pragma once
// Auto-generated by tools/gen_shape_functions -- DO NOT EDIT
// Element: Pyramid5
// Regenerate: /usr/bin/python3 -m tools.gen_shape_functions.generate
 
#include "DNDS/Defines.hpp"
#include "Geom/Geometric.hpp"
#include "Geom/ElemEnum.hpp"
#include "Geom/ElementTraitsBase.hpp"
 
namespace DNDS::Geom::Elem
{
 
    // Forward declaration (primary template is in ElementTraitsBase.hpp)
    template <ElemType> struct ShapeFuncImpl;
 
    // <GEN_SHAPE_FUNCS_BEGIN>
    template <>
    struct ShapeFuncImpl<Pyramid5>
    {
        template <class TPoint, class TArray>
        static inline void Diff0(const TPoint &p, TArray &&v)
        {
            t_real xi = p[0];
            t_real et = p[1];
            t_real zt = p[2];
            // singularity guard
            if (std::abs(1 - zt) < 1e-15)
                zt -= DNDS::signP(1 - zt) * 1e-15;
            const t_real _t0 = zt - 1;
            const t_real _t1 = _t0 + et;
            const t_real _t2 = _t0 + xi;
            const t_real _t3 = (1.0 / (4*zt - 4));
            const t_real _t4 = 1 - zt;
            const t_real _t5 = _t4 + xi;
            const t_real _t6 = ((0.25))/_t0;
            const t_real _t7 = _t4 + et;
            v(0, 0) = -_t1*_t2*_t3;
            v(0, 1) = _t1*_t5*_t6;
            v(0, 2) = -_t3*_t5*_t7;
            v(0, 3) = _t2*_t6*_t7;
            v(0, 4) = zt;
        }
 
        template <class TPoint, class TArray>
        static inline void Diff1(const TPoint &p, TArray &&v)
        {
            t_real xi = p[0];
            t_real et = p[1];
            t_real zt = p[2];
            // singularity guard
            if (std::abs(1 - zt) < 1e-15)
                zt -= DNDS::signP(1 - zt) * 1e-15;
            const t_real _t0 = zt - 1;
            const t_real _t1 = _t0 + et;
            const t_real _t2 = ((0.25))/_t0;
            const t_real _t3 = 1 - zt;
            const t_real _t4 = _t3 + et;
            const t_real _t5 = _t0 + xi;
            const t_real _t6 = _t3 + xi;
            const t_real _t7 = 2*zt;
            const t_real _t8 = et*xi;
            const t_real _t9 = ((zt) * (zt));
            const t_real _t10 = -_t7 + _t9 + 1;
            const t_real _t11 = ((0.25))/_t10;
            const t_real _t12 = _t11*(_t7 + _t8 - _t9 - 1);
            const t_real _t13 = _t11*(-_t10 - _t8);
            v(0, 0) = -_t1*_t2;
            v(0, 1) = _t1*_t2;
            v(0, 2) = -_t2*_t4;
            v(0, 3) = _t2*_t4;
            v(1, 0) = -_t2*_t5;
            v(1, 1) = _t2*_t6;
            v(1, 2) = -_t2*_t6;
            v(1, 3) = _t2*_t5;
            v(2, 0) = _t12;
            v(2, 1) = _t13;
            v(2, 2) = _t12;
            v(2, 3) = _t13;
            v(2, 4) = 1;
        }
 
        template <class TPoint, class TArray>
        static inline void Diff2(const TPoint &p, TArray &&v)
        {
            t_real xi = p[0];
            t_real et = p[1];
            t_real zt = p[2];
            // singularity guard
            if (std::abs(1 - zt) < 1e-15)
                zt -= DNDS::signP(1 - zt) * 1e-15;
            const t_real _t0 = ((zt) * (zt));
            const t_real _t1 = ((zt) * (zt) * (zt));
            const t_real _t2 = et*xi;
            const t_real _t3 = -_t2/(-6*_t0 + 2*_t1 + 6*zt - 2);
            const t_real _t4 = ((0.5))*_t2/(-3*_t0 + _t1 + 3*zt - 1);
            const t_real _t5 = -1/(4*zt - 4);
            const t_real _t6 = zt - 1;
            const t_real _t7 = ((0.25))/_t6;
            const t_real _t8 = ((0.25))/((_t6) * (_t6));
            const t_real _t9 = _t8*xi;
            const t_real _t10 = -_t9;
            const t_real _t11 = _t8*et;
            const t_real _t12 = -_t11;
            v(2, 0) = _t3;
            v(2, 1) = _t4;
            v(2, 2) = _t3;
            v(2, 3) = _t4;
            v(3, 0) = _t5;
            v(3, 1) = _t7;
            v(3, 2) = _t5;
            v(3, 3) = _t7;
            v(4, 0) = _t9;
            v(4, 1) = _t10;
            v(4, 2) = _t9;
            v(4, 3) = _t10;
            v(5, 0) = _t11;
            v(5, 1) = _t12;
            v(5, 2) = _t11;
            v(5, 3) = _t12;
        }
 
        template <class TPoint, class TArray>
        static inline void Diff3(const TPoint &p, TArray &&v)
        {
            t_real xi = p[0];
            t_real et = p[1];
            t_real zt = p[2];
            // singularity guard
            if (std::abs(1 - zt) < 1e-15)
                zt -= DNDS::signP(1 - zt) * 1e-15;
            const t_real _t0 = ((zt) * (zt) * (zt) * (zt));
            const t_real _t1 = ((zt) * (zt));
            const t_real _t2 = ((zt) * (zt) * (zt));
            const t_real _t3 = et*xi;
            const t_real _t4 = ((1.5))*_t3/(_t0 + 6*_t1 - 4*_t2 - 4*zt + 1);
            const t_real _t5 = -3*_t3/(2*_t0 + 12*_t1 - 8*_t2 - 8*zt + 2);
            const t_real _t6 = zt - 1;
            const t_real _t7 = ((0.5))/((_t6) * (_t6) * (_t6));
            const t_real _t8 = _t7*xi;
            const t_real _t9 = -_t8;
            const t_real _t10 = _t7*et;
            const t_real _t11 = -_t10;
            const t_real _t12 = ((0.25))/((_t6) * (_t6));
            const t_real _t13 = -_t12;
            v(2, 0) = _t4;
            v(2, 1) = _t5;
            v(2, 2) = _t4;
            v(2, 3) = _t5;
            v(6, 0) = _t9;
            v(6, 1) = _t8;
            v(6, 2) = _t9;
            v(6, 3) = _t8;
            v(7, 0) = _t11;
            v(7, 1) = _t10;
            v(7, 2) = _t11;
            v(7, 3) = _t10;
            v(9, 0) = _t12;
            v(9, 1) = _t13;
            v(9, 2) = _t12;
            v(9, 3) = _t13;
        }
    };
    // <GEN_SHAPE_FUNCS_END>
 
 
    /**
     * @brief Element traits for 5-node linear pyramid (Pyramid5)
     *
     * Pyramid5 is a 3D pyramidal element with a quadrilateral base and triangular sides,
     * commonly used for:
     * - Transition between hexahedral and tetrahedral meshes
     * - Octree-based mesh refinement
     * - Adaptive meshing with hanging nodes
     *
     * Geometry:
     * - Reference pyramid: square base at z=0, apex at z=1
     * - 5 nodes: 4 at base corners, 1 at apex
     * - Parametric space volume = 4/3
     *
     * Faces:
     * - 1 quadrilateral base face
     * - 4 triangular side faces
     */
    template <>
    struct ElementTraits<Pyramid5>
    {
        // ============================================================
        // Core Element Identification
        // ============================================================
 
        static constexpr ElemType elemType = Pyramid5;
        static constexpr int dim = 3;
        static constexpr int order = 1;
        static constexpr int numVertices = 5;
        static constexpr int numNodes = 5;
        static constexpr int numFaces = 5;
        static constexpr ParamSpace paramSpace = PyramidSpace;
        static constexpr t_real paramSpaceVol = 4.0 / 3.0;
 
        // ============================================================
        // Geometry Definition
        // ============================================================
 
        /**
         * @brief Standard coordinates of nodes in parametric space
         *
         * Reference pyramid: square base [-1,1]x[-1,1] at z=0, apex at (0,0,1)
         * Nodes 0-3: base corners (counter-clockwise)
         * Node 4: apex
         */
        static constexpr std::array<t_real, 3 * 5> standardCoords = {
            -1, -1, 0,   // Node 0: base corner
             1, -1, 0,   // Node 1: base corner
             1,  1, 0,   // Node 2: base corner
            -1,  1, 0,   // Node 3: base corner
             0,  0, 1};  // Node 4: apex
 
        // ============================================================
        // Face/Edge Definitions
        // ============================================================
 
        /**
         * @brief Get the element type of a face
         * @param iFace Face index (0 is base, 1-4 are sides)
         * @return Quad4 for base, Tri3 for sides
         */
        static constexpr ElemType GetFaceType(t_index iFace)
        {
            return iFace < 1 ? Quad4 : Tri3;
        }
 
        /**
         * @brief Node indices for each face
         *
         * Face 0: nodes 0-3-2-1 (base, quadrilateral)
         * Faces 1-4: triangular side faces
         */
        static constexpr std::array<std::array<t_index, 10>, 5> faceNodes = {{
            {0, 3, 2, 1},    // Face 0: base (quad)
            {0, 1, 4},       // Face 1: side triangle
            {1, 2, 4},       // Face 2: side triangle
            {2, 3, 4},       // Face 3: side triangle
            {3, 0, 4}}};     // Face 4: side triangle
 
        // ============================================================
        // Order Elevation (P-Refinement)
        // ============================================================
 
        /**
         * @brief Element type after order elevation (O1 -> O2)
         * Pyramid5 elevates to Pyramid14 (14-node quadratic pyramid)
         */
        static constexpr ElemType elevatedType = Pyramid14;
 
        /// @brief Number of additional nodes created during elevation
        static constexpr int numElevNodes = 9;
 
        /**
         * @brief Elevation spans define new node connections
         *
         * Elevation creates 9 new nodes:
         *   Spans 0-3: base edge midpoints (4 edges)
         *   Spans 4-7: lateral edge midpoints (4 edges)
         *   Span 8: base face center
         */
        static constexpr std::array<tElevSpan, 9> elevSpans = {{
            {0, 1}, {1, 2}, {2, 3}, {3, 0},     // Base edges
            {0, 4}, {1, 4}, {2, 4}, {3, 4},     // Lateral edges
            {0, 3, 2, 1}}};                     // Base face center
 
        /// @brief Element type of each elevation span
        static constexpr std::array<ElemType, 9> elevNodeSpanTypes = {
            Line2, Line2, Line2, Line2,    // Base edges
            Line2, Line2, Line2, Line2,    // Lateral edges
            Quad4};                        // Base face center
 
        // ============================================================
        // VTK/Visualization Support
        // ============================================================
 
        /// @brief VTK cell type identifier (14 = VTK_PYRAMID)
        static constexpr int vtkCellType = 14;
 
        /**
         * @brief VTK node ordering map
         *
         * VTK uses the same ordering as DNDS for Pyramid5:
         *   VTK nodes 0-3 = base corners 0-3
         *   VTK node 4 = apex
         */
        static constexpr std::array<int, 5> vtkNodeOrder = {0, 1, 2, 3, 4};
    };
 
} // namespace DNDS::Geom::Elem