shape.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_Shape__H
#define _GIMLI_Shape__H
 
#include "gimli.h"
#include "polynomial.h"
#include "curvefitting.h"
 
#ifndef PYGIMLI_CAST // fails because of boost threads and clang problems
    #if USE_BOOST_THREAD
        #include <boost/thread.hpp>
        extern boost::mutex ShapeFunctionWriteCacheMutex__;
        //#error AVOID BOOST
    #else
        #include <mutex>
        extern std::mutex ShapeFunctionWriteCacheMutex__;
    #endif
#endif
 
namespace GIMLI{

DLLEXPORT double tetVolume(const RVector3 & p0, const RVector3 & p1,
                           const RVector3 & p2, const RVector3 & p3);

DLLEXPORT double triSize(const RVector3 & p0, const RVector3 & p1,
                         const RVector3 & p2);
 
 
DLLEXPORT std::vector < PolynomialFunction < double > >
createPolynomialShapeFunctions(const std::vector < RVector3 > & pnts,
                               uint dim, uint nCoeff,
                               bool pascale, bool serendipity,
                               const RVector & startVector);
 
 
template < class Ent > std::vector < PolynomialFunction < double > >
    createPolynomialShapeFunctions(const Ent & ent, uint nCoeff,
                                   bool pascale, bool serendipity,
                                   const RVector & startVector){
// __MS(ent)
    std::vector < RVector3 > pnts;
    for (Index i = 0; i < ent.nodeCount(); i ++){
        pnts.push_back(ent.rst(i));
    }
 
    return createPolynomialShapeFunctions(pnts, ent.dim(), nCoeff, pascale,
                                          serendipity, startVector);
}
 
// no default arg here .. pygimli@win64 linker bug
template < class Ent > std::vector < PolynomialFunction < double > >
    createPolynomialShapeFunctions(const Ent & ent, uint nCoeff,
                                   bool pascale, bool serendipity){
    RVector start(0);
    return createPolynomialShapeFunctions(ent, nCoeff, pascale, serendipity, start);
}
 
class DLLEXPORT ShapeFunctionCache : public Singleton< ShapeFunctionCache > {
public:
    friend class Singleton< ShapeFunctionCache >;
 
    template < class Ent > const std::vector < PolynomialFunction < double > > &
    shapeFunctions(const Ent & e) const {
 
        std::map < uint8, std::vector < PolynomialFunction < double > > >::const_iterator it = shapeFunctions_.find(e.rtti());
 
        if (it == shapeFunctions_.end()){
            this->createShapeFunctions_(e);
            it = shapeFunctions_.find(e.rtti());
        }
 
        return (*it).second;
    }
 
    template < class Ent > const std::vector < PolynomialFunction < double > > &
    deriveShapeFunctions(const Ent & e , uint dim) const {
        std::map < uint8, std::vector < std::vector < PolynomialFunction < double > > > >::const_iterator it = dShapeFunctions_.find(e.rtti());
 
        if (it == dShapeFunctions_.end()) {
            this->createShapeFunctions_(e);
            it = dShapeFunctions_.find(e.rtti());
        }
        return (*it).second[dim];
    }

    void clear() {
        shapeFunctions_.clear();
        dShapeFunctions_.clear();
    }
 
    std::vector< RMatrix3 > & RMatrix3Cache();
    std::vector< RMatrix > & RMatrixCache(uint rtti);
 
    RMatrix3 & cachedRMatrix3(uint i);
    RMatrix & cachedRMatrix(uint rtti, uint i);
private:

    template < class Ent > void createShapeFunctions_(const Ent & e) const {
        std::vector < PolynomialFunction < double > > N = e.createShapeFunctions();
 
        #if USE_BOOST_THREAD
        //#ifdef WIN32_LEAN_AND_MEAN
        //        __MS("pls check missing mutex")
            //boost::mutex::scoped_lock lock(ShapeFunctionWriteCacheMutex__);
        //#else
            #ifndef PYGIMLI_CAST // fails because of boost threads and clang problems
            boost::mutex::scoped_lock lock(ShapeFunctionWriteCacheMutex__);
            #endif
        //#endif
 
        #else
        // __MS("lock deactiated " << ShapeFunctionWriteCacheMutex__)
 
            std::unique_lock < std::mutex > lock(ShapeFunctionWriteCacheMutex__);
        #endif
 
 
        shapeFunctions_[e.rtti()] = N;
        dShapeFunctions_[e.rtti()] = std::vector < std::vector < PolynomialFunction < double > > >();
 
        dShapeFunctions_[e.rtti()].push_back(std::vector < PolynomialFunction < double > >());
        dShapeFunctions_[e.rtti()].push_back(std::vector < PolynomialFunction < double > >());
        dShapeFunctions_[e.rtti()].push_back(std::vector < PolynomialFunction < double > >());
 
        for (uint i = 0; i < N.size(); i ++){
            dShapeFunctions_[e.rtti()][0].push_back(N[i].derive(0));
            dShapeFunctions_[e.rtti()][1].push_back(N[i].derive(1));
            dShapeFunctions_[e.rtti()][2].push_back(N[i].derive(2));
        }
    }

    ShapeFunctionCache(){}
    virtual ~ShapeFunctionCache(){}
    ShapeFunctionCache(const ShapeFunctionCache &){};
    void operator = (const ShapeFunctionCache &){};
 
protected:

    mutable std::map < uint8, std::vector < PolynomialFunction < double > > > shapeFunctions_;

    mutable std::map < uint8, std::vector< std::vector < PolynomialFunction < double > > > > dShapeFunctions_;
 
    mutable std::vector< RMatrix3 > _rMatrix3Cache;
    mutable std::map< uint, std::vector< RMatrix > > _rMatrixCache;
 
};
 
static const double NodeCoordinates[1][3] = {
    {0.0, 0.0, 0.0}
};


class DLLEXPORT Shape {
public:
    Shape(MeshEntity * ent);

    virtual ~Shape();

    virtual int rtti() const = 0;
 
    virtual int dim() const = 0;

    virtual std::string name() const { return "Shape"; }

    Index nodeCount() const { return this->nodeCount_; }

    Node & node(Index i) const;
 
    // /*! Return a reference to the i-th \ref Node of this shape. */
    // Node & node(Index i);

    void setNodesPtr(const std::vector< Node * > & n) {
        this->nodeVector_ = &n;
    }

    const std::vector< Node * > & nodes() const { return * this->nodeVector_; }

    // std::vector < Node * > & nodes() { return nodeVector_ ; }
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const ;

    void createJacobian(RMatrix3 & J) const;
 
    RMatrix3 createJacobian() const;

    const RMatrix3 & invJacobian() const;

    virtual RVector N(const RVector3 & L) const;

    virtual void N(const RVector3 & L, RVector & ret) const;

    virtual void dNdrst(const RVector3 & rst, RMatrix & MdNdrst) const;
 
    virtual RMatrix dNdrst(const RVector3 & L) const;

    virtual void rst2xyz(const RVector3 & rst, RVector3 & xyz) const;

    virtual RVector3 xyz(const RVector3 & rst) const;

    virtual void xyz2rst(const RVector3 & xyz, RVector3 & rst) const;

    virtual RVector3 rst(const RVector3 & xyz) const;

    virtual RVector3 rst(Index i) const;

    inline double drstdxyz(uint rstI, uint xyzJ) const {
//         return invJacobian()[rstI][xyzJ];}
        return invJacobian()[rstI * 3 + xyzJ];}

    virtual bool isInside(const RVector3 & xyz, bool verbose=false) const;

    virtual bool isInside(const RVector3 & xyz, RVector & sf,
                          bool verbose=false) const;

    virtual bool touch(const RVector3 & pos, double tol=1e-6, bool verbose=false) const;
 
    /*!* Return true if the ray intersects the shape.
     * On boundary means inside too. The intersection position is stored in pos.
     * */
    virtual bool intersectRay(const RVector3 & start, const RVector3 & dir,
                              RVector3 & pos){
        __MS(rtti() << " " << name())
        THROW_TO_IMPL
        return false;
    };

    double domainSize() const;

    RVector3 center() const;

    virtual RVector3 norm() const;

    double h() const;

    virtual Plane plane() const;

    void changed();
 
    double jacobianDeterminant() const { return det(this->createJacobian()); }
 
    inline void resizeNodeSize_(Index n) { this->nodeCount_ = n; }
 
 
protected:

    virtual double domainSize_() const { return 0.0; }
 
    Index nodeCount_;
 
    mutable double domSize_;
    mutable bool hasDomSize_;
    mutable double _h;
 
    mutable RMatrix3 invJacobian_;
 
    // const std::vector< Node * > & nodes()
    const std::vector < Node * > * nodeVector_;
};
 
DLLEXPORT std::ostream & operator << (std::ostream & str, const Shape & c);
 
class DLLEXPORT NodeShape : public Shape{
public:
    NodeShape(MeshEntity * ent) : Shape(ent) { resizeNodeSize_(1); }
 
    virtual ~NodeShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_NODE_RTTI; }
 
    virtual int dim() const { return 0; }
 
    virtual std::string name() const { return "NodeShape"; }
 
    virtual RVector3 rst(Index i) const;
 
    virtual double domainSize_() const { return 1.0; }
 
    virtual RVector3 norm() const;
 
protected:
    //virtual double jacobianDeterminant_() const { return 0.0; }
};
 
static const double EdgeCoordinates[2][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0}
};
 
class DLLEXPORT EdgeShape : public Shape {
public:
    EdgeShape(MeshEntity * ent) : Shape(ent) { resizeNodeSize_(2); }
 
    virtual ~EdgeShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_EDGE_RTTI; }
 
    virtual int dim() const { return 1; }
 
    virtual std::string name() const { return "EdgeShape"; }

    virtual RVector3 rst(Index i) const;
 
    /*!* Return true if the ray intersects the shape.
     * On boundary means inside too. The intersection position is stored in pos.
     * */
    virtual bool intersectRay(const RVector3 & start, const RVector3 & dir,
                              RVector3 & pos);

    virtual bool touch(const RVector3 & pos, double tol=1e-6, bool verbose=false) const;
 
    double length() const;
 
    virtual RVector3 norm() const;
 
protected:
 
    virtual double domainSize_() const { return length(); }
};
 
static const double TriCoordinates[3][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {0.0, 1.0, 0.0}
};

class DLLEXPORT TriangleShape : public Shape {
public:
 
    TriangleShape(MeshEntity * ent) : Shape(ent) { resizeNodeSize_(3); }
 
    virtual ~TriangleShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_TRIANGLE_RTTI; }
 
    virtual int dim() const { return 2; }
 
    virtual std::string name() const { return "TriangleShape"; }

    virtual RVector3 rst(Index i) const;

    virtual void xyz2rst(const RVector3 & pos, RVector3 & rst) const;
 
    double area() const;
 
    virtual RVector3 norm() const;
 
    /*!* Return true if the ray intersects the shape.
     * On boundary means inside too. The intersection position is stored in pos.
     * */
    virtual bool intersectRay(const RVector3 & start, const RVector3 & dir,
                              RVector3 & pos);
 
protected:

    virtual double domainSize_() const { return area(); }
};
 
 
static const double QuadCoordinates[4][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {1.0, 1.0, 0.0},
    {0.0, 1.0, 0.0}
};


class DLLEXPORT QuadrangleShape : public Shape {
public:
 
    QuadrangleShape(MeshEntity * ent) : Shape(ent) { resizeNodeSize_(4); }
 
    virtual ~QuadrangleShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_QUADRANGLE_RTTI; }
 
    virtual int dim() const { return 2; }
 
    virtual std::string name() const { return "QuadrangleShape"; }

    virtual RVector3 rst(Index i) const;
 
    double area() const;
 
//     virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
//     /*! See Shape::N. */
//     virtual void N(const RVector3 & L, RVector & n) const;
 
//     /*! See Shape::dNdrst. */
//     virtual RMatrix dNdrst(const RVector3 & rst) const;
//
    virtual RVector3 norm() const;
 
    virtual bool intersectRay(const RVector3 & start, const RVector3 & dir,
                              RVector3 & pos);
 
protected:
    virtual double domainSize_() const { return area(); }
};
 
 
class DLLEXPORT PolygonShape : public Shape {
public:
    PolygonShape(MeshEntity * ent);
 
    virtual ~PolygonShape();
 
    virtual int rtti() const { return MESH_SHAPE_POLYGON_FACE_RTTI; }
 
    virtual int dim() const { return 3; }
 
    virtual std::string name() const { return "PolygonFaceShape"; }
 
    virtual RVector3 norm() const;
 
    virtual RVector3 rst(Index i) const;
 
    virtual bool isInside(const RVector3 & xyz, bool verbose) const;
 
protected:
};
 
 
static const double TetCoordinates[4][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {0.0, 1.0, 0.0},
    {0.0, 0.0, 1.0}
};


class DLLEXPORT TetrahedronShape : public Shape {
public:
    TetrahedronShape(MeshEntity * ent) : Shape(ent) { resizeNodeSize_(4); }
 
    virtual ~TetrahedronShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_TETRAHEDRON_RTTI; }
 
    virtual int dim() const { return 3; }
 
    virtual std::string name() const { return "TetrahedronShape"; }

    virtual RVector3 rst(Index i) const;

    void xyz2rst(const RVector3 & pos, RVector3 & rst) const;
 
//     /*! See Shape::N. */
//     virtual void N(const RVector3 & L, RVector & n) const;
//
//     /*! See Shape::dNdrst. */
//     virtual RMatrix dNdrst(const RVector3 & rst) const;
 
    double volume() const;
 
protected:
 
    virtual double domainSize_() const { return volume(); }
};
 
static const int HexahedronSplit5TetID[5][4] = {
    {1, 4, 5, 6},
    {3, 6, 7, 4},
    {1, 0, 4, 3},
    {1, 2, 3, 6},
    {1, 4, 6, 3}
};
 
static const int HexahedronSplit6TetID[6][4] = {
    {0, 1, 2, 6},
    {0, 2, 3, 6},
    {0, 1, 6, 5},
    {0, 4, 5, 6},
    {0, 3, 7, 6},
    {0, 4, 6, 7}
};
 
static const double HexCoordinates[8][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {1.0, 1.0, 0.0},
    {0.0, 1.0, 0.0},
    {0.0, 0.0, 1.0},
    {1.0, 0.0, 1.0},
    {1.0, 1.0, 1.0},
    {0.0, 1.0, 1.0}
};


 
class DLLEXPORT HexahedronShape : public Shape {
public:
    HexahedronShape(MeshEntity * ent) : Shape(ent){ resizeNodeSize_(8); }
 
    virtual ~HexahedronShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_HEXAHEDRON_RTTI; }
 
    virtual int dim() const { return 3; }

    virtual RVector3 rst(Index i) const;
 
    virtual std::string name() const { return "HexahedronShape"; }
 
    double volume() const;

    virtual bool enforcePositiveDirection();
 
//     /*! See Shape::N. */
//     virtual void N(const RVector3 & L, RVector & n) const;
//
//     /*! See Shape::dNdrst. */
//     virtual RMatrix dNdrst(const RVector3 & rst) const;
 
protected:
 
    virtual double domainSize_() const { return volume(); }
};
 
static const uint8 TriPrimSplit3TetID[3][4] = {
    {0, 5, 1, 2},
    {3, 1, 5, 4},
    {3, 5, 1, 0}
};
 
static const double PrismCoordinates[6][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {0.0, 1.0, 0.0},
    {0.0, 0.0, 1.0},
    {1.0, 0.0, 1.0},
    {0.0, 1.0, 1.0}
};


 
class DLLEXPORT TriPrismShape : public Shape {
public:
    TriPrismShape(MeshEntity * ent) : Shape(ent){ resizeNodeSize_(6); }
 
    virtual ~TriPrismShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_TRIPRISM_RTTI; }
 
    virtual int dim() const { return 3; }
 
    virtual std::string name() const { return "TriagonalPrismShape"; }

    virtual RVector3 rst(Index i) const;
 
    virtual std::vector < PolynomialFunction < double > > createShapeFunctions() const;
 
    double volume() const;
 
protected:
 
    virtual double domainSize_() const { return volume(); }
};
 
static const double PyramidCoordinates[5][3] = {
    {0.0, 0.0, 0.0},
    {1.0, 0.0, 0.0},
    {1.0, 1.0, 0.0},
    {0.0, 1.0, 0.0},
    {0.0, 0.0, 0.1}
};


class DLLEXPORT PyramidShape : public Shape {
public:
    PyramidShape(MeshEntity * ent) : Shape(ent){ resizeNodeSize_(5); }
 
    virtual ~PyramidShape(){ }
 
    virtual int rtti() const { return MESH_SHAPE_PYRAMID_RTTI; }
 
    virtual int dim() const { return 3; }

    virtual RVector3 rst(Index i) const;
 
//     /*! See Shape::N. */
//     virtual void N(const RVector3 & L, RVector & n) const;
//
//     /*! See Shape::dNdrst. */
//     virtual RMatrix dNdrst(const RVector3 & rst) const;
 
};
 
} // namespace GIMLI
 
#endif // _GIMLI_SHAPE__H