GEO_PrimHexahedron.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME:        GEO_PrimHexahedron.h (GEO Library, C++)
 *
 * COMMENTS:    This is the hexahedron primitive type.
 */

#pragma once

#ifndef __GEO_PrimHexahedron__
#define __GEO_PrimHexahedron__

#include "GEO_API.h"
#include "GEO_Primitive.h"
#include "GEO_Detail.h"
#include "GEO_VolumeElementBase.h"
#include <GA/GA_Defines.h>
#include <UT/UT_Array.h>
#include <SYS/SYS_Inline.h>

class GA_Detail;
class GEO_ConvertParms;

// So it doesn't look like a magic number
#define PT_PER_HEX      8
#define FACE_PER_HEX    6
#define EDGE_PER_HEX    12

class GEO_API GEO_PrimHexahedron : public GEO_VolumeElementBase
{
public:
    /// NOTE: This constructor should only be called via GU_PrimitiveFactory.
    SYS_FORCE_INLINE
    GEO_PrimHexahedron(GA_Detail *d, GA_Offset offset = GA_INVALID_OFFSET)
        : GEO_VolumeElementBase(d, offset)
    {}

    /// This constructor is for making a representation of a hexahedron
    /// on the stack, so that you can call GEO_PrimHexahedron functions on
    /// the hexahedron without needing the detail to allocate one.
    SYS_FORCE_INLINE
    GEO_PrimHexahedron(GA_Detail *gdp, GA_Offset offset, const GA_OffsetListRef &vertex_list)
        : GEO_VolumeElementBase(gdp, offset)
    {
        UT_ASSERT_P(vertex_list.size() == PT_PER_HEX);
        myVertexList = vertex_list;
    }

    const GA_PrimitiveDefinition &getTypeDef() const override;

    /// Report approximate memory usage.
    int64       getMemoryUsage() const override;

    /// Count memory usage using a UT_MemoryCounter in order to count
    /// shared memory correctly.
    /// NOTE: This should always include sizeof(*this).
    void        countMemory(UT_MemoryCounter &counter) const override;

    void        reverse() override;
    UT_Vector3D  computeNormalD() const override;
    UT_Vector3   computeNormal() const override;

    // Compute the location of the breakpoint. Return 0 if OK, else -1.
    fpreal      calcVolume(const UT_Vector3 &refpt) const override;
    fpreal      calcArea() const override;
    fpreal      calcPerimeter() const override;

    /// @{
    /// Save/Load vertex list to a JSON stream
    bool                saveVertexArray(UT_JSONWriter &w,
                                const GA_SaveMap &map) const;
    bool                loadVertexArray(UT_JSONParser &p,
                                const GA_LoadMap &map);
    /// @}

    bool        isDegenerate() const override;

    // Take the whole set of points into consideration when applying the
    // point removal operation to this primitive. The method returns 0 if
    // successful, -1 if it failed because it would have become degenerate,
    // and -2 if it failed because it would have had to remove the primitive
    // altogether.
    int          detachPoints (GA_PointGroup &grp) override;

    /// Before a point is deleted, all primitives using the point will be
    /// notified.  The method should return "false" if it's impossible to
    /// delete the point.  Otherwise, the vertices should be removed.
    GA_DereferenceStatus dereferencePoint(
        GA_Offset point, bool dry_run=false) override;
    GA_DereferenceStatus dereferencePoints(
        const GA_RangeMemberQuery &pt_q, bool dry_run=false) override;

    //
    // Hexahedron layout
    //
    //   Z+  Y+
    //     \ | 
    //      \|
    //  X+---0
    //
    //  Point numbers: bit 0: X, bit 1: Y, bit 2: Z
    //  Face: Axis X:0, Y:2, Z:4, direction: 0/1
    //
    //  7---6    //  0 - right
    //  |\  |\   //  1 - left
    //  | 3---2  //  2 - bottom
    //  5-|-4 |  //  3 - top
    //   \|  \|  //  4 - front
    //    1---0  //  5 - back
    //
    // Each face starts with the least vertex #.

    SYS_FORCE_INLINE GA_Offset fastVertexOffset(GA_Size index) const
    { return getFastVertexOffset(index); }
    static GA_Size              rawVertexCount(const GA_Detail *gdp = nullptr, GA_Offset = GA_INVALID_OFFSET) { return PT_PER_HEX; }
    SYS_FORCE_INLINE GA_Size getFastVertexCount() const
    { return PT_PER_HEX; }
    SYS_FORCE_INLINE GA_Offset getFastVertexOffset(GA_Size index) const
    {
        UT_ASSERT_P(index >= 0 && index < PT_PER_HEX);
        return getVertexOffset(index);
    }

    /// Return the topology of the hexahedron
    static GA_Size      rawFaceCount(const GA_Detail *gdp = nullptr, GA_Offset offset = GA_INVALID_OFFSET) { return FACE_PER_HEX; }
    GA_Size     getFaceCount() const override { return rawFaceCount(); }

    static int rawFaceIndices(const GA_Detail *gdp, GA_Offset offset, GA_Size faceno, UT_Array<int> &vtxlist)
    {
        vtxlist.clear();

        const int *f = fastFaceIndices(faceno);
        vtxlist.append(f[0]);
        vtxlist.append(f[1]);
        vtxlist.append(f[2]);
        vtxlist.append(f[3]);

        return vtxlist.entries();
    }
    int getFaceIndices(GA_Size faceno, UT_Array<int> &vtxlist) const override
    { return rawFaceIndices(nullptr, GA_INVALID_OFFSET, faceno, vtxlist); }

    static int rawFaceIndexCount(const GA_Detail *gdp, GA_Offset offset, GA_Size faceno)
    { return 4; }
    int getFaceIndexCount(GA_Size faceno) const override
    { return rawFaceIndexCount(nullptr, GA_INVALID_OFFSET, faceno); }

    SYS_FORCE_INLINE
    static const int *fastFaceIndices(GA_Size faceno)
    {
        static const int hex_faceidx[FACE_PER_HEX][4] =
        {
            { 0, 2, 6, 4 },
            { 1, 5, 7, 3 },
            { 0, 4, 5, 1 },
            { 2, 3, 7, 6 },
            { 0, 1, 3, 2 },
            { 4, 6, 7, 5 },
        };

        return hex_faceidx[faceno];
    }

    /// Return the barycentric coordinate of the face.  Convertes from
    /// the natural parameterization of the face:
    /// s = vertex0 to vertex1
    /// t = vertex0 to vertex3
    /// (s, t) -> (u, v, w)
    /// (0, 0) -> (u0, v0, w0)
    /// (1, 1) -> (u1, v1, w1)
    static UT_Vector3 remapFaceCoords(GA_Size faceno, float s, float t);

    /// Get the vertex offsets of the specified quad.
    void getQuadVertexOffsets(
        GA_Size quad, GA_Offset &v0, GA_Offset &v1, GA_Offset &v2, GA_Offset &v3) const
    {
        const int *indices = fastFaceIndices(quad);
        v0 = fastVertexOffset(indices[0]);
        v1 = fastVertexOffset(indices[1]);
        v2 = fastVertexOffset(indices[2]);
        v3 = fastVertexOffset(indices[3]);
    }

    static GA_Size rawEdgeCount(const GA_Detail *gdp = nullptr, GA_Offset = GA_INVALID_OFFSET)
    { return EDGE_PER_HEX; }
    GA_Size     getEdgeCount() const override
    { return rawEdgeCount(); }
    static void rawEdgeIndices(const GA_Detail *gdp, GA_Offset offset, GA_Size edgeno, int &e0, int &e1)
    {
        static const int                edgeidx[EDGE_PER_HEX][2] =
        {
            { 0, 1 },
            { 0, 2 },
            { 0, 4 },
            { 1, 3 },
            { 1, 5 },
            { 2, 3 },
            { 2, 6 },
            { 3, 7 },
            { 4, 5 },
            { 4, 6 },
            { 5, 7 },
            { 6, 7 },
        };

        e0 = edgeidx[edgeno][0];
        e1 = edgeidx[edgeno][1];
    }
    void        getEdgeIndices(GA_Size edgeno, int &e0, int &e1) const override
    { return rawEdgeIndices(nullptr, GA_INVALID_OFFSET, edgeno, e0, e1); }

    GA_Offset findSharedHex(GA_Size faceno) const;
    SYS_FORCE_INLINE
    bool isFaceShared(GA_Size faceno) const
    { return GAisValid(findSharedHex(faceno)); }

    GA_Offset findSharedFace(GA_Size faceno) const override
    { return findSharedHex(faceno); }

    // These methods find the index of the given vertex (vtx or the vertex
    // that contains a pt). Return -1 if not found. 
    GA_Size     findVertex(GA_Offset vtx) const
                        {
                            for (GA_Size i = 0; i < PT_PER_HEX; i++)
                                if (fastVertexOffset(i) == vtx) return i;
                            return -1;
                        }

    GA_Size     find(GA_Offset pt) const
                        { 
                            for (GA_Size i = 0; i < PT_PER_HEX; i++)
                                if (vertexPoint(i) == pt) return i;
                            return -1;
                        }

    void        setVertexPoint(int i, GA_Offset pt)
                        {
                            if (i < PT_PER_HEX)
                                wireVertex(fastVertexOffset(i), pt);
                        }

    const GA_PrimitiveJSON      *getJSON() const override;

    class Face
    {
    public:
        Face(const GEO_PrimHexahedron &tet,int face)
            : myHex(tet)
            , myI(face)
            , myF(GEO_PrimHexahedron::fastFaceIndices(myI))
        {}
        inline GA_Size getFastVertexCount() const
        { return 4; }
        inline GA_Size getVertexCount() const
        { return 4; }
        inline GA_Offset getFastVertexOffset(GA_Size i) const
        { return myHex.getFastVertexOffset(myF[i]); }
        inline GA_Offset getVertexOffset(GA_Size i) const
        { return getFastVertexOffset(i); }
        inline GA_Offset getPointOffset(GA_Size i) const
        { return myHex.getPointOffset(myF[i]); }
        inline UT_Vector3 getPos3(GA_Size i) const
        { return myHex.getPos3(myF[i]); }
        inline bool isClosed() const
        { return true; }
    private:
        const GEO_PrimHexahedron &myHex;
        int myI;
        const int *const myF;
    };

    /// NOTE: Do not use this under normal circumstances!
    ///       This is for if a tet has been added using
    ///       GA_Detail::appendPrimitive[Block] and vertices have
    ///       to be added one at a time, e.g. from VEX.
    ///       It will return GA_INVALID_OFFSET and won't add a vertex
    ///       if all PT_PER_HEX have already been assigned.
    GA_Offset unsafeAppendVertex();

    void normal(NormalComp &output) const override;
    void normal(NormalCompD &output) const override;

    // Conversion Methods
    GEO_Primitive *convert(GEO_ConvertParms &parms,
                           GA_PointGroup *usedpts = nullptr) override;
    GEO_Primitive *convertNew(GEO_ConvertParms &parms) override;

    /// Builds a single hexahedron primitive, optionally creating new points
    /// and wiring them to its vertices.
    static GEO_PrimHexahedron *build(GEO_Detail *gdp, bool appendPts = true);

    static void rawComputeInteriorPointWeights(
        const GA_Detail *gdp, GA_Offset offset,
        UT_Array<GA_Offset> &vtxlist, UT_Array<float> &weightlist,
        fpreal u, fpreal v, fpreal w);
    void computeInteriorPointWeights(
        UT_Array<GA_Offset> &vtxlist, UT_Array<float> &weightlist,
        fpreal u, fpreal v, fpreal w) const override final
    { rawComputeInteriorPointWeights(getParent(), getMapOffset(),
            vtxlist, weightlist, u, v, w); }

    /// Trilinear Hex
    ///
    /// Cube deformed with tri-linear deformation    
    class TLHex
    {
    public:
        TLHex(const UT_Vector3 pos[8])
        {
            myCoefficients[0] = pos[0];
            myCoefficients[1] = pos[1] - pos[0];
            myCoefficients[2] = pos[2] - pos[0];
            myCoefficients[3] = pos[3] - pos[1] - pos[2] + pos[0];
            myCoefficients[4] = pos[4] - pos[0];
            myCoefficients[5] = pos[5] - pos[4] - pos[1] + pos[0];
            myCoefficients[6] = pos[6] - pos[4] - pos[2] + pos[0];
            myCoefficients[7] = pos[7] - pos[6] - pos[5] + pos[4] - pos[3] + pos[2] + pos[1] - pos[0];

            UT_BoundingBox bbox;
            bbox.initBounds(pos[0]);
            for(int i = 1; i < 8; ++i)
                bbox.enlargeBounds(pos[i]);
            myTol2 = 1e-6 * bbox.size().length2();
        }

        // convert position relative to a unit lattice cell to worldspace
        UT_Vector3 unitToHex(const UT_Vector3 &u) const
        {
            return  myCoefficients[0] + unitToHexTranslated(u);
        }

        // convert a position in worldspace to one relative to a unit lattice cell
        bool hexToUnit(UT_Vector3 &u, const UT_Vector3 &p) const
        {
            UT_Vector3 pos = p - myCoefficients[0];

            constexpr static float GRID_LOWER_TOL = 1e-4;
            constexpr static float GRID_UPPER_TOL = GRID_LOWER_TOL + 1;

            constexpr static int maxIterations = 100;

            u.assign(0.5, 0.5, 0.5);
            for (int i = 0; i < maxIterations; ++i)
            {
                UT_Matrix3 J = getJacobian(u);
                auto det = J.determinant();
                if (det == 0)
                    return false;

                auto divDet = 1/det;

                auto J00 = J(0,0);
                auto J01 = J(0,1);
                auto J02 = J(0,2);
                auto J10 = J(1,0);
                auto J11 = J(1,1);
                auto J12 = J(1,2);

                UT_Vector3 val = pos - unitToHexTranslated(u);
                auto cx = val[0];
                auto cy = val[1];
                auto cz = val[2];

                UT_Vector3 delta;
                J(0,0) = cx; J(0,1) = cy; J(0,2) = cz;
                delta.x() = J.determinant() * divDet;

                J(0,0) = J00; J(0,1) = J01; J(0,2) = J02;
                J(1,0) = cx; J(1,1) = cy; J(1,2) = cz;
                delta.y() = J.determinant() * divDet;

                J(1,0) = J10; J(1,1) = J11; J(1,2) = J12;
                J(2,0) = cx; J(2,1) = cy; J(2,2) = cz;
                delta.z() = J.determinant() * divDet;

                u += delta;
                if (delta.length2() <= myTol2)
                {
                    for (int j = 0; j < 3; ++j)
                    {
                        auto &v = u(j);
                        if (!SYSisFinite(v) || v < GRID_LOWER_TOL || v >= GRID_UPPER_TOL)
                            return false;
                    }
                    return true;
                }
            }

            return false;
        }

        // given a position relative to a unit lattice cell, return the Jacobian
        // of the transformation applied to produce the deformed position
        UT_Matrix3 getJacobian(const UT_Vector3 &u) const
        {
            auto dfd0 = myCoefficients[1] +
                        myCoefficients[3] * u[1] +
                        myCoefficients[5] * u[2] +
                        myCoefficients[7] * u[1] * u[2];

            auto dfd1 = myCoefficients[2] +
                        myCoefficients[3] * u[0] +
                        myCoefficients[6] * u[2] +
                        myCoefficients[7] * u[0] * u[2];

            auto dfd2 = myCoefficients[4] +
                        myCoefficients[5] * u[0] +
                        myCoefficients[6] * u[1] +
                        myCoefficients[7] * u[0] * u[1];

            return UT_Matrix3(dfd0, dfd1, dfd2);
        }

    private:
        // convert position relative to a unit lattice cell to worldspace
        // translated to put corner 0 at the origin
        UT_Vector3 unitToHexTranslated(const UT_Vector3 &u) const
        {
            return  myCoefficients[1] * u[0] +
                    myCoefficients[2] * u[1] +
                    myCoefficients[3] * u[0] * u[1] +
                    myCoefficients[4] * u[2] +
                    myCoefficients[5] * u[0] * u[2] +
                    myCoefficients[6] * u[1] * u[2] +
                    myCoefficients[7] * u[0] * u[1] * u[2];
        }

        UT_Vector3 myCoefficients[8];
        float myTol2;
    };

protected:
    void allocateVertices();

    template <typename T>
    UT_Vector3T<T> internalComputeNormal() const;
    template <typename T>
    fpreal      internalCalcVolume(const UT_Vector3T<T> &refpt) const;
    template <typename T>
    fpreal      internalCalcArea() const;
    template <typename T>
    fpreal      internalCalcPerimeter() const;

    GA_DECLARE_INTRINSICS(override)

    // Check the validity of the data. Meant to be called especially at loading
    // time. The method returns 1 if OK and 0 if trouble.
    virtual bool validate() const;

    /// @warning vertexPoint() doesn't check the bounds.  Use with caution.
    GA_Offset vertexPoint(GA_Size i) const
    { return getDetail().vertexPoint(fastVertexOffset(i)); }

    // Evaluates the position at domain point (u,v,w) in the interior of the
    // geometry, rather than on the perimeter of the poly.
    // (Note that u, v, w don't have to be converted to unit domain.  We assume
    // that the barycentric coords or bilinear interpolants are already 
    // correctly scaled).
    // Edge 0-1 is the U, 0-2 the V, and 0-4 the W.
    bool evaluateInteriorPointRefMap(
        GA_Offset result_vtx, GA_AttributeRefMap &map,
        fpreal u, fpreal v, fpreal w = 0) const override;
    int  evaluateInteriorPointV4(
        UT_Vector4 &pos, fpreal u, fpreal v, fpreal w = 0) const override;

    static GA_PrimitiveDefinition *theDefinition;
    friend class GU_PrimitiveFactory;
private:
};

using GU_PrimHexahedron = GEO_PrimHexahedron;

#undef PT_PER_HEX
#undef FACE_PER_HEX
#undef EDGE_PER_HEX

#endif