GU_Cosserat.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: GU_Cosserat.h ( CE Library, C++)
 *
 * COMMENTS: Functions to perform stable cosserat rods simulation.
 */

#pragma once

#include "GU_API.h"
#include <CE/CE_Array.h>
#include <UT/UT_Quaternion.h>
#include <UT/UT_Matrix3.h>
#include <variant>

class CE_Context;

namespace GU_Cosserat
{

class CachedItems
{
public:
    GU_API void loadKernels(
            CE_Context &context,
            int n,
            bool recompile);

    template <typename T>
    class Buffer
    {
    public:
        void initIfInternal(exint size)
        {
            if (std::holds_alternative<CE_Array<T>>(myBuffer))
                return std::get<CE_Array<T>>(myBuffer).init(size);
        };

        cl::Buffer buffer()
        {
            if (std::holds_alternative<CE_Array<T>>(myBuffer))
                return std::get<CE_Array<T>>(myBuffer).buffer();
            else
                return std::get<cl::Buffer>(myBuffer);
        }

        void setExternal(cl::Buffer buffer)
        {
            myBuffer = buffer;
        }

        void setInternal()
        {
            if (!std::holds_alternative<CE_Array<T>>(myBuffer))
                myBuffer = CE_Array<T>();
        }
        bool isExternal()
        {
            return std::holds_alternative<cl::Buffer>(myBuffer);
        }
    private:
        std::variant<CE_Array<T>, cl::Buffer> myBuffer;
    };

    class Buffers
    {
    public:
        GU_API void init(int n);

        Buffer<float> myPos;
        Buffer<float> myOrient;
        Buffer<float> myRestOrient;
        Buffer<float> myFlags;
        Buffer<float> myVel;
        Buffer<float> myDx;
        Buffer<float> myRefPos;
        Buffer<float> myRestLength;
        Buffer<float> myInvMass;
        Buffer<int> myClumpPtsIndex;
        Buffer<int> myClumpPts;
        Buffer<float> myClumpDists;

        // Ping-pong buffers: alternate read/write sides each iteration instead
        // of copying.  myPosSwapped/myOrientSwapped track which side is current.
        // GUcosseratFinalize normalises back to myPos/myOrient when needed.
        Buffer<float> myPIn;
        Buffer<float> myDxIn;
        Buffer<float> myOrientIn;

        bool myPosSwapped    = false;
        bool myOrientSwapped = false;

    } myBuffers;

    struct Kernels
    {
        cl::KernelFunctor myInitBuffer;
        cl::KernelFunctor myInitRestLength;
        cl::KernelFunctor myInitInvMass;
        cl::KernelFunctor myCopyBuffer;
        cl::Kernel myIterate;
        cl::KernelFunctor myUpdateOrient;
    } myKernels;

public:
    GU_API Buffers &getBuffers() { return myBuffers; }

    GU_API Kernels &getKernels() { return myKernels; }
};

enum class PointFlags
{
    HASPREV = 1,
    HASNEXT = 2,
    NEXTORIENT = 4,
    FIXORIENT = 8
};

struct SimSettings
{
    float myTimeStep = 1.0f;
    float myDamping = 1.0f;
    float myAccelerationRatio = 1.0f;
    float myKStretch = 1.0f;
    float myKBend = 1.0f;
    float myDensity = 1.0f;
    float myRootBias = 0.5f;
    float myRefPosMult = 1.0f;
    float myClumpStiffness = 1.0f;
    float myClumpDamping = 1.0f;
    bool myLockRoots = false;
};

int computeInitFlags(int curveptindex, int ncurvepts)
{
    int flags = 0;

    if (curveptindex > 0)
        flags |= (int)GU_Cosserat::PointFlags::HASPREV;

    if (curveptindex < ncurvepts-1)
    {
        flags |= (int)GU_Cosserat::PointFlags::HASNEXT;

        if (curveptindex < ncurvepts-2)
            flags |= (int)GU_Cosserat::PointFlags::NEXTORIENT;
    }

    return flags;
}

UT_Quaternion computeMinimalTwistOrient(
        const UT_Vector3 &lastpos,
        const UT_Vector3 &pos,
        const UT_Quaternion &lastorient)
{
    UT_Vector3F seg = pos - lastpos;
    if (seg.length2() < 1e-12f)
        return lastorient;
    UT_Matrix3F temp;
    UT_Vector3F xaxis = {0,0,1};
    temp.dihedral(xaxis, seg);
    UT_Quaternion qq;
    qq.updateFromRotationMatrix(temp);

    // minimize twist
    UT_QuaternionF invqq(qq);
    invqq.invert();
    UT_QuaternionF t = invqq * lastorient;
    t.x() = 0.0f;
    t.y() = 0.0f;
    t.normalize();

    return qq * t;
}
} // namespace GU_Cosserat

// by default this uses internal CE_Arrays for every buffer.
// call setExternal beforehand to use external cl::Buffer
GU_API void GUcosseratInit(
        CE_Context &context,
        int n,
        GU_Cosserat::SimSettings &sim_settings,
        GU_Cosserat::CachedItems &cached_items,
        bool recompile);

// by default this uses internal CE_Arrays for every buffer.
// call setExternal beforehand to use external cl::Buffer
GU_API void GUcosseratInitRefPos(
        CE_Context &context,
        int n,
        GU_Cosserat::CachedItems &cached_items);

GU_API void GUcosseratIter(
        CE_Context &context,
        int n,
        GU_Cosserat::SimSettings &sim_settings,
        GU_Cosserat::CachedItems &cached_items);

GU_API void GUcosseratUpdateOrient(
        CE_Context &context,
        int n,
        GU_Cosserat::SimSettings &sim_settings,
        GU_Cosserat::CachedItems &cached_items);

// Symmetrize a CSR clump-constraint graph so every directed edge (a→b) has a
// matching reverse edge (b→a).  The initial build from a nearest-neighbour
// search is asymmetric (maxconstraints=3 by default), causing one-sided forces
// in the Jacobi solver.  Call once after building the initial arrays.
GU_API void GUcosseratSymmetrizeClumpConstraints(
        UT_Array<int>   &clumppts,
        UT_Array<int>   &clumpptsindex,
        UT_Array<float> &clumpdists,
        exint            npts);

// Must be called once after the iteration loop.  For even iteration counts
// this is a no-op; for odd counts it copies the ping-pong buffer back to
// myPos/myOrient so callers always read results from the primary buffers.
GU_API void GUcosseratFinalize(
        CE_Context &context,
        int n,
        GU_Cosserat::CachedItems &cached_items);