meshentities.h

/******************************************************************************
 *   Copyright (C) 2006-2024 by the GIMLi development team                    *
 *   Carsten Rücker carsten@resistivity.net                                   *
 *                                                                            *
 *   Licensed under the Apache License, Version 2.0 (the "License");          *
 *   you may not use this file except in compliance with the License.         *
 *   You may obtain a copy of the License at                                  *
 *                                                                            *
 *       http://www.apache.org/licenses/LICENSE-2.0                           *
 *                                                                            *
 *   Unless required by applicable law or agreed to in writing, software      *
 *   distributed under the License is distributed on an "AS IS" BASIS,        *
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
 *   See the License for the specific language governing permissions and      *
 *   limitations under the License.                                           *
 *                                                                            *
 ******************************************************************************/
 
#ifndef _GIMLI_MESHENTITIES__H
#define _GIMLI_MESHENTITIES__H
 
#include "gimli.h"
 
#include "vector.h"
#include "matrix.h"
#include "elementmatrix.h"
#include "baseentity.h"
 
#include "polynomial.h"
 
#include <vector>
#include <set>
 
namespace GIMLI{
 
template < class Typname > bool lesserId(const Typname * a, const Typname * b){
    return (a->id() < b->id());
}
 
DLLEXPORT Boundary * findBoundary(const Node & n1);
DLLEXPORT Boundary * findBoundary(const Node & n1, const Node & n2);
DLLEXPORT Boundary * findBoundary(const Node & n1, const Node & n2, const Node & n3);
 
DLLEXPORT Boundary * findBoundary(const std::vector < Node * > & n);

DLLEXPORT std::set < Boundary * >  
findBoundaries(const std::vector < Node * > & n);

DLLEXPORT Boundary * findCommonBoundary(const Cell & c1, const Cell & c2);
 
DLLEXPORT Cell * findCommonCell(const std::vector < Node * > & n, bool warn = true);

class CollectNodeFunctor{
public:
    CollectNodeFunctor(std::set< Node * > & c):c_(& c){}
 
    template < class E > void operator()(E * e){
        for (Index i = 0; i < e->nodeCount(); i++) c_->insert(& e->node(i));
    }
    std::set< Node * > * c_;
};

template < class ContainerOfMeshEntities >
std::set< Node * > commonNodes(const ContainerOfMeshEntities & c){
    std::set< Node * > commons;
 
    CollectNodeFunctor gc(commons);
    for_each(c.begin(), c.end(), gc);
    return commons;
}
 
class DLLEXPORT RegionMarker : public RVector3{
public:
    RegionMarker(const RVector3 & pos, int marker, double area=0.0,
                 bool hole=false)
    : RVector3(pos), marker_(marker), area_(area), isHole_(hole){}
 
    ~RegionMarker(){}
 
    inline void setMarker(SIndex marker) {marker_ = marker;}
    inline int marker() const {return marker_;}
 
    inline void setArea(double area) {area_ = area;}
    inline double area() const {return area_;}
 
    inline void setPos(const Pos & pos) {copy_(pos);}
 
    bool isHole() const { return isHole_; }
    inline void setHole(bool h) { isHole_ = h; }
 
protected:
    int marker_;
    double area_;
    bool isHole_;
};
 
class DLLEXPORT MeshEntity : public BaseEntity {
public:

    MeshEntity();

    virtual ~MeshEntity();

    virtual uint dim() const { return 0; }

    virtual uint rtti() const { return MESH_MESHENTITY_RTTI; }

    virtual uint parentType() const { return MESH_MESHENTITY_RTTI; }
 
    virtual void setNodes(const std::vector < Node * > & nodes);
 
    const std::vector< Node * > & nodes() const { return nodeVector_; }
 
    inline Node & node(uint i) {
        ASSERT_RANGE(i, 0, nodeCount()); return *nodeVector_[i];
    }
    
    inline Node & node(uint i) const {
        ASSERT_RANGE(i, 0, nodeCount()); return *nodeVector_[i];
    }
 
    inline uint nodeCount() const { return nodeVector_.size(); }
 
    inline Shape & shape() { return *shape_; }
 
    inline Shape & shape() const { return *shape_; }
 
    inline Shape * pShape() { return shape_; }

    RVector3 rst(uint i) const;

    RVector3 center() const;

    double size() const;
 
    virtual double attribute() const { return -1.0; }

    IndexArray ids() const;
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;

    virtual RVector N(const RVector3 & rst) const;

    virtual void N(const RVector3 & rst, RVector & n) const;

    virtual RVector dNdL(const RVector3 & rst, uint i) const;

    virtual RMatrix dNdL(const RVector3 & rst) const;

    double pot(const RVector3 & p, const RVector & u) const;

    RVector3 vec(const RVector3 & p, const R3Vector & v) const;

    RVector3 grad(const RVector3 & p, const RVector & u) const;
 
    friend std::ostream & operator << (std::ostream & str, const MeshEntity & c);
 
//    void setUxCache(const ElementMatrix < double >  & mat) const { uxCache_ = mat; }
 
 //   const ElementMatrix < double > & uxCache() const { return uxCache_; }
 
    void setUxCache(const RMatrix & mat) const { uxCache_ = mat; }
 
    const RMatrix & uxCache() const { return uxCache_; }
 
    ElementMatrix < double > & uCache(){ return uCache_; }
 
    ElementMatrix < double > & gradUCache(){ return gradUCache_; }

    void changed();
 
    void addSecondaryNode(Node * n);
 
    void delSecondaryNode(Node * n);
 
    const std::vector < Node * > & secondaryNodes() const;

    const std::vector < Node * > allNodes() const;
 
    Index allNodeCount() const;

    virtual bool enforcePositiveDirection();
 
protected:
    void fillShape_();
 
    virtual void registerNodes_();
    virtual void deRegisterNodes_();
 
    virtual void registerSecNode_(Node *n);
    virtual void deRegisterSecNode_(Node *n);
    
 
    Shape * shape_;
 
    std::vector < Node * > nodeVector_;
    std::vector < Node * > secondaryNodes_;

    mutable ElementMatrix < double > uCache_;
    mutable ElementMatrix < double > gradUCache_;
 
    mutable RMatrix uxCache_;
 
protected:
    MeshEntity(const MeshEntity & ent){
        THROW_TO_IMPL
    }

    MeshEntity & operator = (const MeshEntity & ent){
         if (this != & ent){
             THROW_TO_IMPL
         } return *this;
    }
};


class DLLEXPORT Cell : public MeshEntity {
public:
    DLLEXPORT friend std::ostream & operator << (std::ostream & str, const Cell & c);

    Cell();

    Cell(const std::vector < Node * > & nodes);

    virtual ~Cell();

    bool operator==(const Cell & cell){
        return &cell == this;
    }
 
    virtual uint rtti() const { return MESH_CELL_RTTI; }
    virtual uint parentType() const { return MESH_CELL_RTTI; }
    virtual uint neighborCellCount() const { return 0; }
    inline uint boundaryCount() const { return neighborCellCount(); }
 
    void cleanNeighborInfos();
 
    Cell * neighborCell(const RVector & sf);

    inline Cell * neighborCell(uint i){ return neighborCells_[i]; }

    virtual void findNeighborCell(uint i);
 
    inline double attribute() const { return attribute_; }
 
    inline void setAttribute(double attr) { attribute_ = attr; }

    Node * oppositeTo(const Boundary & bound);

    Boundary * boundaryTo(const RVector & sf);

    Boundary * boundary(Index i);

    virtual std::vector < Node * > boundaryNodes(Index i) const{
        CERR_TO_IMPL
        std::cout << rtti() << std::endl;
        std::vector < Node * > n;
        return n;
    }
 
protected:
    virtual void registerNodes_();
    virtual void deRegisterNodes_();
    virtual void registerSecNode_(Node *n);
    virtual void deRegisterSecNode_(Node *n);
    
    std::vector < Cell * > neighborCells_;
    double attribute_;
 
protected:
    Cell(const Cell & cell){
        std::cerr << "cell(const cell & cell)" << std::endl;
        THROW_TO_IMPL
    }

    Cell & operator = (const Cell & cell){
        if (this != &cell) {
            THROW_TO_IMPL
            std::cerr << "cell=cell" << std::endl;
        }
        return *this;
    }
 
};
 
class DLLEXPORT Boundary : public MeshEntity{
public:
    Boundary();
    Boundary(const std::vector < Node * > & nodes);
    virtual ~Boundary();
 
    virtual uint rtti() const { return MESH_BOUNDARY_RTTI; }
    virtual uint parentType() const { return MESH_BOUNDARY_RTTI; }

    virtual RVector3 rst(uint i) const;

    inline const Cell & leftCell() const { return *leftCell_; }
    inline Cell * leftCell() { return leftCell_; }
 
    inline const Cell & rightCell() const { return *rightCell_; }
    inline Cell * rightCell() { return rightCell_; }
 
    inline void setLeftCell(Cell * cell) { leftCell_  = cell; }
    inline void setRightCell(Cell * cell) { rightCell_ = cell; }
 
    friend std::ostream & operator << (std::ostream & str, const Boundary & e);

    virtual RVector3 norm() const;

    virtual RVector3 norm(const Cell & cell) const;

    bool normShowsOutside(const Cell & cell) const;

    void swapNorm(bool withNeighbours=true);

    bool outside() const { return (leftCell_ != 0) && (rightCell_ == 0); }
 
 
protected:
    virtual void registerNodes_();
    virtual void deRegisterNodes_();
    virtual void registerSecNode_(Node *n);
    virtual void deRegisterSecNode_(Node *n);
    
 
    Cell *leftCell_;
    Cell *rightCell_;
 
protected:
    Boundary(const Boundary & bound){
        std::cerr << "Boundary(const Boundary & bound)" << std::endl;
        THROW_TO_IMPL
    }

    Boundary & operator = (const Boundary & boundary){
        if (this != &boundary) {
            std::cerr << "Assignment for boundaries not yet supported." << std::endl;
            THROW_TO_IMPL
        }
        return *this;
    }
 
};
 
class DLLEXPORT NodeBoundary : public Boundary{
public:
    NodeBoundary(Node & n1);
 
    NodeBoundary(const std::vector < Node * > & nodes);
 
    virtual ~NodeBoundary();
 
    virtual uint dim() const { return 1; }
 
    virtual uint rtti() const { return MESH_BOUNDARY_NODE_RTTI; }
 
    void setNodes(Node & n1);
 
    virtual double size() const { return 1.0; }
 
    friend std::ostream & operator << (std::ostream & str, const NodeBoundary & e);

    virtual RVector3 norm() const;
 
protected:
};
 
class DLLEXPORT Edge : public Boundary{
public:
    Edge(Node & n1, Node & n2);
 
    Edge(const std::vector < Node * > & nodes);
 
    virtual ~Edge();
 
    virtual uint dim() const { return 1; }
 
    virtual uint rtti() const { return MESH_EDGE_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    void setNodes(Node & n1, Node & n2);

    int swap();
 
//     /*! See ref MeshEntity::N() */
//     virtual RVector N(const RVector3 & L) const;
 
//     void shapeFunctionsL(const RVector3 & L, RVector & funct) const;
 
    friend std::ostream & operator << (std::ostream & str, const Edge & e);
 
protected:
};
 
class DLLEXPORT Edge3 : public Edge{
public:
 
    Edge3(const std::vector < Node * > & nodes);
 
    virtual ~Edge3();
 
    virtual uint rtti() const { return MESH_EDGE3_RTTI; }

    virtual RVector3 rst(uint i) const;
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
//     /*! See MeshEntity::N() */
//     virtual RVector N(const RVector3 & L) const;
    /*
    void shapeFunctionsL(const RVector3 & L, RVector & funct) const;*/
 
// friend std::ostream & operator << (std::ostream & str, const Edge & e);
 
protected:
};
 
 
class DLLEXPORT TriangleFace : public Boundary{
public:
    TriangleFace(Node & n1, Node & n2, Node & n3);
 
    TriangleFace(const std::vector < Node * > & nodes);
 
    virtual ~TriangleFace();
 
    virtual uint dim() const { return 2; }
 
    virtual uint rtti() const { return MESH_TRIANGLEFACE_RTTI; }
 
    void setNodes(Node & n1, Node & n2, Node & n3);
 
    friend std::ostream & operator << (std::ostream & str, const TriangleFace & e);

    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
 
     // /*! Don't call this class directly */
    // TriangleFace(const TriangleFace & bound){
        // std::cerr << "TriangleFace(const Boundary & bound)" << std::endl;
    // }
 
    // /*! Don't call this class directly */
    // TriangleFace & operator = (const TriangleFace & boundary){
        // if (this != &boundary) {
            // std::cerr << "TriangleFace boundary=boundary" << std::endl;
        // }
        // return *this;
    // }
};
 
class DLLEXPORT Triangle6Face : public TriangleFace{
public:
    Triangle6Face(const std::vector < Node * > & nodes);
 
    ~Triangle6Face();
 
    virtual uint rtti() const { return MESH_TRIANGLEFACE6_RTTI; }

    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;

    virtual RVector3 rst(uint i) const;
 
protected:
};
 
class DLLEXPORT QuadrangleFace : public Boundary{
public:
    QuadrangleFace(Node & n1, Node & n2, Node & n3, Node & n4);
 
    QuadrangleFace(const std::vector < Node * > & nodes);
 
    virtual ~QuadrangleFace();
 
    virtual uint dim() const { return 2; }
 
    virtual uint rtti() const { return MESH_QUADRANGLEFACE_RTTI; }

    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    void setNodes(Node & n1, Node & n2, Node & n3, Node & n4);
 
    friend std::ostream & operator << (std::ostream & str, const TriangleFace & e);
 
protected:
 
    QuadrangleFace(const QuadrangleFace & quad){
        std::cerr << "QuadrangleFace(const QuadrangleFace & quad)" << std::endl;
    }
};
 
class DLLEXPORT Quadrangle8Face : public QuadrangleFace{
public:
 
    Quadrangle8Face(const std::vector < Node * > & nodes);
 
    virtual ~Quadrangle8Face();
 
    virtual uint rtti() const { return MESH_QUADRANGLEFACE8_RTTI; }

    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;

    virtual RVector3 rst(uint i) const;
 
protected:
 
private:
 
};
 
class DLLEXPORT PolygonFace : public Boundary {
public:
    typedef Pos HoleMarker;
    typedef std::vector < Pos >  HoleMarkerList;
 
    PolygonFace(const std::vector < Node * > & nodes);
 
    ~PolygonFace();
 
    virtual uint dim() const { return 3; }
 
    virtual uint rtti() const { return MESH_POLYGON_FACE_RTTI; }

    void insertNode(Node * node, double tol=TOLERANCE);

    void addSubface(const std::vector < Node * > & nodes, bool isHole=false);
 
    Index subfaceCount() const {return this->subfaces_.size();}
 
    const std::vector < Node * > & subface(Index i) const;

    void addHoleMarker(const RVector3 & pos);
 
    void delHoleMarker(const RVector3 & pos);

    const HoleMarkerList & holeMarkers() const;

    HoleMarkerList & holeMarkers(){ return holeMarker_;}
 
protected:
    std::vector < std::vector < Node * > > subfaces_;
    HoleMarkerList holeMarker_;
private:
 
};
 
class DLLEXPORT EdgeCell : public Cell {
public:
    EdgeCell(Node & n1, Node & n2);
 
    EdgeCell(const std::vector < Node * > & nodes);
 
    virtual ~EdgeCell();
 
    virtual uint dim() const { return 1; }
 
    virtual uint rtti() const { return MESH_EDGE_CELL_RTTI; }
 
    virtual uint neighborCellCount() const { return 2; }
 
//     virtual void findNeighborCell(uint id);
 
    void setNodes(Node & n1, Node & n2);
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
    friend std::ostream & operator << (std::ostream & str, const EdgeCell & t);
 
protected:
};

class DLLEXPORT Edge3Cell : public EdgeCell {
public:
    Edge3Cell(const std::vector < Node * > & nodes);
 
    virtual ~Edge3Cell();
 
    virtual uint rtti() const { return MESH_EDGE3_CELL_RTTI; }
 
    virtual RVector3 rst(uint i) const;
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};


class DLLEXPORT Triangle : public Cell {
public:
    Triangle(Node & n1, Node & n2, Node & n3);
 
    Triangle(const std::vector < Node * > & nodes);
 
    virtual ~Triangle();
 
    virtual uint dim() const { return 2; }
 
    virtual uint rtti() const { return MESH_TRIANGLE_RTTI; }
 
    virtual uint neighborCellCount() const { return 3; }
 
//     virtual void findNeighborCell(uint i);
 
    void setNodes(Node & n1, Node & n2, Node & n3);
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    friend std::ostream & operator << (std::ostream & str, const Triangle & t);

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};


class DLLEXPORT Triangle6 : public Triangle {
public:
    Triangle6(const std::vector < Node * > & nodes);
 
    virtual ~Triangle6();
 
    virtual uint rtti() const { return MESH_TRIANGLE6_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};


class DLLEXPORT Quadrangle : public Cell {
public:
    Quadrangle(Node & n1, Node & n2, Node & n3, Node & n4);
 
    Quadrangle(const std::vector < Node * > & nodes);
 
    virtual ~Quadrangle();
 
    virtual uint dim() const { return 2; }
 
    virtual uint rtti() const { return MESH_QUADRANGLE_RTTI; }
 
    void setNodes(Node & n1, Node & n2, Node & n3, Node & n4);
 
    virtual uint neighborCellCount() const { return 4; }
 
//     virtual void findNeighborCell(uint i);
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    friend std::ostream & operator << (std::ostream & str, const Quadrangle & t);
 
    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};


class DLLEXPORT Quadrangle8 : public Quadrangle {
public:
    Quadrangle8(const std::vector < Node * > & nodes);
 
    virtual ~Quadrangle8();
 
    virtual uint rtti() const { return MESH_QUADRANGLE8_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};


static const uint8 TetrahedronFacesID[4][3] = {
    {1, 2, 3},
    {2, 0, 3},
    {0, 1, 3},
    {0, 2, 1}
};


class DLLEXPORT Tetrahedron : public Cell {
public:
    Tetrahedron(Node & n1, Node & n2, Node & n3, Node & n4);
 
    Tetrahedron(const std::vector < Node * > & nodes);
 
    virtual ~Tetrahedron();
 
    virtual uint dim() const { return 3; }
 
    virtual uint rtti() const { return MESH_TETRAHEDRON_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    void setNodes(Node & n1, Node & n2, Node & n3, Node & n4);
 
    virtual uint neighborCellCount() const { return 4; }
 
//     virtual void findNeighborCell(uint i);
 
    friend std::ostream & operator << (std::ostream & str, const Tetrahedron & t);

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:

    Tetrahedron(const Tetrahedron& cell){ std::cerr << "Tetrahedron cell(const cell & cell)" << std::endl; }

    Tetrahedron & operator = (const Tetrahedron & cell){
        if (this != &cell) {
            std::cerr << "Tetrahedron cell=cell" << std::endl;
        }
        return *this;
    }
 
};
 
//*! VTK,Flaherty,Gimli count: 1-2-3-4, 5(1-2), 6(2-3), 7(3-1), 8(1-4), 9(2-4), 10(3-4)* //
static const uint8 Tet10NodeSplit[10][2] = {
    {0,0},{1,1},{2,2},{3,3},
    {0,1},{1,2},{2,0},
    {0,3},{1,3},{2,3}
};
 
//*! Zienkiewicz count: 1-2-3-4, 5(1-2), 6(1-3), 7(1-4), 8(2-3), 9(3-4), 10(4-2)* //
static const uint8 Tet10NodeSplitZienk[10][2] = {
    {0,0},{1,1},{2,2},{3,3},
    {0,1},{0,2},{0,3},
    {1,2},{2,3},{3,1}
};
 
class DLLEXPORT Tetrahedron10 : public Tetrahedron {
public:
    Tetrahedron10(const std::vector < Node * > & nodes);
 
    virtual ~Tetrahedron10();
 
    virtual uint rtti() const { return MESH_TETRAHEDRON10_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};


 
static const uint8 HexahedronFacesID[6][4] = {
    {1, 2, 6, 5},
    {2, 3, 7, 6},
    {3, 0, 4, 7},
    {0, 1, 5, 4},
    {4, 5, 6, 7},
    {0, 3, 2, 1}
};
 
class DLLEXPORT Hexahedron: public Cell {
public:
    Hexahedron(const std::vector < Node * > & nodes);
 
    virtual ~Hexahedron();
 
    virtual uint dim() const { return 3; }
 
    virtual uint rtti() const { return MESH_HEXAHEDRON_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    virtual uint neighborCellCount() const { return 6; }
 
    friend std::ostream & operator << (std::ostream & str, const Hexahedron & t);

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};
 
static const uint8 Hexahedron20FacesID[6][8] = {
    {0,1,5,4,8,17,12,16},
    {1,2,6,5,9,18,13,17},
    {2,3,7,6,10,19,14,18},
    {3,0,4,7,11,16,15,19},
    {0,3,2,1,11,10,9,8},
    {4,5,6,7,12,13,14,15},
};
 
static const uint8 Hex20NodeSplit[20][2] = {
    {0,0},{1,1},{2,2},{3,3},{4,4},{5,5},{6,6},{7,7},
    {0,1},{1,2},{2,3},{3,0},
    {4,5},{5,6},{6,7},{7,4},
    {0,4},{1,5},{2,6},{3,7}
};


class DLLEXPORT Hexahedron20: public Hexahedron {
public:
    Hexahedron20(const std::vector < Node * > & nodes);
 
    virtual ~Hexahedron20();
 
    virtual uint rtti() const { return MESH_HEXAHEDRON20_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};


 
static const uint8 TriPrismFacesID[5][4] = {
    {1, 2, 5, 4},        // r
    {2, 0, 3, 5},        // r
    {0, 1, 4, 3},        // l
    {3, 4, 5, 255},      // l
    {0, 2, 1, 255},      // r
};
 
 
class DLLEXPORT TriPrism : public Cell {
public:
    TriPrism(const std::vector < Node * > & nodes);
 
    virtual ~TriPrism();
 
    virtual uint dim() const { return 3; }
 
    virtual uint rtti() const { return MESH_TRIPRISM_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    virtual uint neighborCellCount() const { return 5; }
 
    friend std::ostream & operator << (std::ostream & str, const Hexahedron & t);

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};
 
static const uint8 Prism15NodeSplit[15][2] = {
    {0,0},{1,1},{2,2},{3,3},{4,4},{5,5},
    {0,1},{1,2},{2,0},
    {3,4},{4,5},{5,3},
    {0,3},{1,4},{2,5}
};


 
class DLLEXPORT TriPrism15 : public TriPrism {
public:
    TriPrism15(const std::vector < Node * > & nodes);
 
    virtual ~TriPrism15();
 
    virtual uint rtti() const { return MESH_TRIPRISM15_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};


 
static const uint8 PyramidFacesID[5][4] = {
    {1, 2, 5, 255},    // l
    {2, 3, 5, 255},    // l
    {0, 5, 3, 255},    // l
    {0, 1, 5, 255},    // l
    {0, 3, 2, 1},      // r
};
 
class DLLEXPORT Pyramid : public Cell {
public:
    Pyramid(const std::vector < Node * > & nodes);
 
    virtual ~Pyramid();
 
    virtual uint dim() const { return 3; }
 
    virtual uint rtti() const { return MESH_PYRAMID_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    virtual uint neighborCellCount() const { return 5; }

    virtual std::vector < Node * > boundaryNodes(Index i) const;
 
protected:
};
 
//*! VTK, Flaherty, Gimli count: 1-2-3-4, 5(1-2), 6(2-3), 7(3-1), 8(1-4), 9(2-4), 10(3-4)* //
static const uint8 Pyramid13NodeSplit[13][2] = {
    {0,0},{1,1},{2,2},{3,3},{4,4},
    {0,1},{1,2},{2,3},{3,0},
    {0,4},{1,4},{2,4},{3,4}
};


 
class DLLEXPORT Pyramid13 : public Pyramid {
public:
    Pyramid13(const std::vector < Node * > & nodes);
 
    virtual ~Pyramid13();
 
    virtual uint rtti() const { return MESH_PYRAMID13_RTTI; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
protected:
};
 
DLLEXPORT std::ostream & operator << (std::ostream & str, const std::set < GIMLI::MeshEntity * > & ents);
 
DLLEXPORT std::ostream & operator << (std::ostream & str, const std::set < GIMLI::Boundary * > & bounds);
 
} // namespace GIMLI
 
#endif // MESHENTITIES__H