---

simSimulator.h

//////////////////////////////////////////////////////////////////////
// Copyright (c) 2002-2003 David Pritchard <drpritch@alumni.uwaterloo.ca>
// 
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License
// as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later
// version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
// 
// You should have received a copy of the GNU Lesser General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

#ifndef freecloth_sim_simSimulator_h
#define freecloth_sim_simSimulator_h

#ifndef freecloth_sim_package_h
#include <freecloth/simulator/package.h>
#endif

#ifndef freecloth_sim_simMatrix_h
#include <freecloth/simulator/simMatrix.h>
#endif

#ifndef freecloth_sim_simVector_h
#include <freecloth/simulator/simVector.h>
#endif

#ifndef freecloth_base_algorithm
#include <freecloth/base/algorithm>
#endif

#ifndef freecloth_base_list
#include <freecloth/base/list>
#endif

#ifndef freecloth_base_baTime_h
#include <freecloth/base/baTime.h>
#endif

#ifndef freecloth_resmgt_rcShdPtr_h
#include <freecloth/resmgt/rcShdPtr.h>
#endif

#ifndef freecloth_resmgt_rcBase_h
#include <freecloth/resmgt/rcBase.h>
#endif

#ifndef freecloth_geom_geMesh_h
#include <freecloth/geom/geMesh.h>
#endif

#ifndef freecloth_geom_geMeshWingedEdge_h
#include <freecloth/geom/geMeshWingedEdge.h>
#endif

FREECLOTH_NAMESPACE_START

////////////////////////////////////////////////////////////////////////////////
// FORWARD DECLARATIONS

class GeMatrix3;

////////////////////////////////////////////////////////////////////////////////
/*!
 * \class SimSimulator freecloth/simulator/simSimulator.h
 * \brief An implementation of Baraff & Witkin's cloth simulation algorithm.
 * 
 * Using the simulation is fairly straightforward. The client supplies an
 * initial mesh, which must have disc topology. The client also specifies the
 * density of the cloth, which indirectly defines the mass of each triangle of
 * the mesh. The client can also specify constraints on the cloth motion,
 * constraining both the position and the velocity of individual nodes of the
 * mesh.
 *
 * Time starts at zero, and is advanced by a fixed timestep by calls to step().
 * The client can retrieve the mesh after each timestep. The simulation can
 * be rewound to the beginning by calling rewind().
 *
 * There are many parameters for the simulation algorithm. See [BarWit98] or
 * the description in SimSimulator::Params for details.
 *
 * References:
 * - [BarWit98] D. Baraff and A. Witkin. Large Steps in Cloth Simulation.
 *    SIGGRAPH Conference Proceedings, 1998, 43-54.
 *    http://www-2.cs.cmu.edu/~baraff/papers/sig98.pdf
 * - [Pri02] D. Pritchard. "Implementing Baraff & Witkin's Cloth Simulation".
 *    Unpublished, 2002.
 *    http://freecloth.enigmati.ca/docs/report-chaps/
 * - [Mac00] D. Macri. Real-Time Cloth. Game Developers Conference, 2000.
 *    http://www.gdconf.com/archives/2000/macri.doc
 */

// FIXME: the mesh should be allowed to have cylindrical or spherical
// topology, not just disc topology.
//
// FIXME: should really implement this as a Bridge pattern, to hide the
// implementation from clients better.
class SimSimulator : public RCBase
{
public:
    // ----- classes -----
    ////////////////////////////////////////////////////////////////////////////
    /*!
     * \class Params freecloth/simulator/simSimulator.h
     * \brief Parameters for cloth simulator, as described in [BarWit98].
     *
     * The stretch, shear, and bend parameters control the strength of the
     * associated internal forces. For example, a higher _k_stretch value
     * makes the internal stretch forces stronger, and hence makes the cloth
     * harder to stretch. The bend force can be controlled independently in
     * u and v. Generally, stretch forces should be very, very strong,
     * and bend forces should be very weak.
     *
     * The damp parameter controls the strength of damping forces. A higher
     * damping value will make the cloth come to rest quicker, while a lower
     * value will allow it to oscillate back and forth before coming to rest.
     *
     * The b_u and b_v values allow the client to force stretch/compression
     * of the cloth. They're not very useful at present; they should really be
     * spatially varying to be useful.
     *
     * The g parameter is the constant for acceleration due to gravity,
     * by default 9.81 m/s^2.
     */
    class Params
    {
    public:
        // ----- member functions -----

        // MSVC bug requires this to be inline
        // Maya values: .5, .35, .01 (unknown scaling...)
        inline Params()
          : _k_stretch( 5000.f ),
            _k_shear( 500.f ),
            _k_bend_u( .00001f ),
            _k_bend_v( .00001f ),
            _k_damp( .2f ),
            _k_drag( .1f ),
            _b_u( 1.f ),
            _b_v( 1.f ),
            _g( 9.81f )
        {
        }

        // ----- data members -----
        // TODO: Maya has separate u/v control for this
        Float           _k_stretch;
        Float           _k_shear;
        Float           _k_bend_u, _k_bend_v;
        Float           _k_damp;
        //! Drag constant.
        //! Technically, this is 0.5 * C_D * rho, where rho is density of
        //! the medium and C_D is the drag coefficient. We combine this into
        //! a single constant for simplicity.
        Float           _k_drag;
        //! Stretch constant
        // FIXME: should be spatially varying, via texture.
        Float           _b_u,_b_v;
        //! Gravitational acceleration
        Float           _g;
    };

    // ----- types and enumerations -----

    //! The individual classes of forces.
    enum ForceType {
        F_STRETCH,
        F_SHEAR,
        F_BEND,
        F_GRAVITY,
        F_DRAG,

        NB_FORCES
    };

    // ----- member functions -----

    explicit SimSimulator( const GeMesh& initialMesh );

    const GeMesh& getMesh() const;
    const RCShdPtr<GeMesh>& getMeshPtr() const;

    //! Rewind simulation time to zero, and return the mesh to its initial
    //! state.
    void rewind();

    //! Advance the simulation by one timestep.
    //! step() is equivalent to calling
    //! - preSubSteps()
    //! - repeated subStep() calls
    //! - postSubSteps()
    void step();

    //@{
    //! Finer grained simulation control - allows the client more control over
    //! the computation. See step() for details.
    bool subStepsDone() const;
    void preSubSteps();
    void subStep();
    void postSubSteps();
    //@}

    //! Retrieve the current simulation time.
    BaTime::Instant getTime() const;

    //@{
    //! Mutator
    void setTimestep( Float );
    void setParams( const Params& );
    void setDensity( Float rho );
    void setPCGTolerance( Float );
    //@}

    //@{
    //! Accessor
    Float getTimestep() const;
    const Params& getParams() const;
    Float getDensity() const;
    Float getPCGTolerance() const;
    //@}

    //! Remove all constraints on all vertices.
    void removeAllConstraints();
    //! Constrain movement of the vertex to the plane defined by the given
    //! normal.
    void setPosConstraintPlane( GeMesh::VertexId, const GeVector& );
    //! Constrain movement of the vertex to the given line.
    void setPosConstraintLine( GeMesh::VertexId, const GeVector& );
    //! Completely constrain the vertex, allowing no movement.
    void setPosConstraintFull( GeMesh::VertexId );

    //! Force the vertex to move with the specified velocity. This overrides
    //! any position constraints.
    void setVelConstraint( GeMesh::VertexId, const GeVector& );

    //@{
    //! Only exposed for debugging purposes.
    GeVector getVelocity( GeMesh::VertexId vid ) const;
    GeVector getForce( GeMesh::VertexId vid ) const;
    Float getEnergy( ForceType ) const;
    Float getTriEnergy( ForceType, GeMesh::FaceId ) const;
    // The final two will only work in debug builds.
    GeVector getForce( ForceType type, GeMesh::VertexId vid ) const;
    GeVector getDampingForce( ForceType type, GeMesh::VertexId vid ) const;

    void setMesh( const GeMesh& );
    //@}

private:

    // ----- types and enumerations -----

    // FIXME: implement custom matrix classes for these some time, if it
    // helps the performance much.
    typedef SimMatrix Matrix;
    typedef SimMatrix SymMatrix;
    typedef SimMatrix TridiagMatrix;

    // ----- classes -----

    class CommonVars;
    class StretchVars;
    class ShearVars;
    class BendVars;
    
    ////////////////////////////////////////////////////////////////////////////
    /*!
     * \class ModPCGSolver freecloth/simulator/simSimulator.h
     *
     * Preconditioned modified conjugate gradient solver. Solves a linear
     * system Ax=b subject to constraints defined by S and z, as defined by
     * [BarWit98] in section 5.3.
     */
    class ModPCGSolver
    {
    public:
        // ----- member functions -----

        ModPCGSolver();
        void preStep();
        void step();
        bool done() const;
        const SimVector& result() const;

        //! Specify the tolerance.
        void setTolerance( Float );
        //! Access the tolerance
        Float getTolerance() const;

        // ----- data members -----
        
        std::vector<GeMatrix3> _S;
        SimVector       _z;

        // Final Ax=b system
        SymMatrix       _A;
        SimVector       _b;

    private:
        // ----- member functions -----

        void setupPreconditioner();
        void setupCG();
        SimVector filter( const SimVector& ) const;
        void filterInPlace( SimVector& ) const;

        // ----- data members -----
        
        SimMatrix       _P, _Pinv;
        bool            _done;
        UInt32          _nbSteps;
        Float           _tolerance;

        SimVector       _x, _r, _c;
        Float           _delta0, _deltaNew;
        
        // Temporaries within step()
        SimVector       _q, _s;
    };


    // ----- member functions -----

    //! Fill the mass matrix (_M) using the initial mesh's triangles areas
    //! and the density parameter, as per [BarWit98] section 2.2.
    void setupMass();

    //! Calculate values common to the stretch and shear conditions.
    void calcStretchShear(
        const GeMesh::FaceWrapper& face
    );
    //! Calculate the stretch condition and its derivatives.
    void calcStretch(
        const GeMesh::FaceWrapper& face,
        const CommonVars &cv,
        const GePoint x[3],
        const GePoint tp[3],
        const GeVector v0[ 3 ]
    );
    //! Calculate the shear condition and its derivatives.
    void calcShear(
        const GeMesh::FaceWrapper& face,
        const CommonVars &cv,
        const GePoint x[3],
        const GePoint tp[3],
        const GeVector v0[ 3 ]
    );
    //! Calculate the bend condition and its derivatives.
    void calcBend(
        const GeMeshWingedEdge::HalfEdgeWrapper& edge
    );
    
    // ----- data members -----
    
    // FIXME: initial mesh should be RCShdPtr<const GeMesh>
    RCShdPtr<GeMesh> _initialMesh;
    RCShdPtr<GeMeshWingedEdge> _initialMeshWingedEdge;
    RCShdPtr<GeMesh> _mesh;
    
    Params          _params;
    //! Density of cloth ( kg / m^2 )
    Float           _rho;
    //! Timestep
    Float           _h;
    //! Current time
    Float           _time;


    //! Variables from the [BarWit98], as defined in equations (6) and (14)
    TridiagMatrix   _M;
    //! Sum of diagonals of _M
    Float           _totalMass;
    SimVector       _f0;
    SymMatrix       _df_dx;
    SymMatrix       _df_dv;
    SimVector       _v0;
    SimVector       _z0;

    //@{
    //! For debugging
    SimVector       _f0i[ NB_FORCES ];
    SimVector       _d0i[ NB_FORCES ];
    SymMatrix       _df_dxi[ NB_FORCES ];
    // df_dv is always zero...
    // Damping forces from the [BarWit98], as defined in section 4.5
    SymMatrix       _dd_dvi[ NB_FORCES ];
    SymMatrix       _dd_dxi[ NB_FORCES ];
    //@}
    
    //@{
    //! Total energy of cloth
    Float           _fenergy[ NB_FORCES ];
    Float           _venergy;
    Float           _energy;
    std::vector<Float> _trienergy[ NB_FORCES ];
    //@}

    ModPCGSolver    _modPCG;
};

FREECLOTH_NAMESPACE_END

#endif

---