System.h

Go to the documentation of this file.

#ifndef SimTK_SimTKCOMMON_SYSTEM_H_
#define SimTK_SimTKCOMMON_SYSTEM_H_

/* -------------------------------------------------------------------------- *
 *                       Simbody(tm): SimTKcommon                             *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK biosimulation toolkit originating from           *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2006-13 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors: Peter Eastman                                                *
 *                                                                            *
 * 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.                                             *
 * -------------------------------------------------------------------------- */

#include "SimTKcommon/basics.h"
#include "SimTKcommon/Simmatrix.h"
#include "SimTKcommon/internal/State.h"
#include "SimTKcommon/internal/Subsystem.h"
#include "SimTKcommon/internal/SubsystemGuts.h"

#include <cassert>

namespace SimTK {

class DecorativeGeometry;
class DefaultSystemSubsystem;
class ScheduledEventHandler;
class ScheduledEventReporter;
class TriggeredEventHandler;
class TriggeredEventReporter;
class RealizeOptions;
class RealizeResults;
class ProjectOptions;
class ProjectResults;

//==============================================================================
//                                 SYSTEM
//==============================================================================
class SimTK_SimTKCOMMON_EXPORT System {
public:
class Guts; // local; name is System::Guts


//------------------------------------------------------------------------------
System& setUpDirection(const CoordinateDirection& up);

System& setUseUniformBackground(bool useUniformBackground);

System& setDefaultTimeScale(Real tc);

System& setDefaultLengthScale(Real lc);

void setHasTimeAdvancedEvents(bool); // default=false

CoordinateDirection getUpDirection() const;
bool getUseUniformBackground() const;
Real getDefaultTimeScale() const;
Real getDefaultLengthScale() const;
bool hasTimeAdvancedEvents() const;
//------------------------------------------------------------------------------
inline void addEventHandler(ScheduledEventHandler* handler);
inline void addEventHandler(TriggeredEventHandler* handler);
inline void addEventReporter(ScheduledEventReporter* handler) const;
inline void addEventReporter(TriggeredEventReporter* handler) const;
//------------------------------------------------------------------------------
const State& realizeTopology() const;

const State& getDefaultState() const;
State&       updDefaultState();

void realizeModel(State& state) const;

void realize(const State& state, Stage stage = Stage::HighestRuntime) const;
//------------------------------------------------------------------------------
void project(State& state, Real accuracy=-1) const;

void projectQ(State& state, Real accuracy=-1) const;

void projectU(State& state, Real accuracy=-1) const;

void projectQ(State& state, Vector& qErrEst, 
             const ProjectOptions& options, ProjectResults& results) const;

void projectU(State& state, Vector& uErrEst, 
             const ProjectOptions& options, ProjectResults& results) const;


void prescribe(State& state) const {
    realize(state, Stage::Time);
    prescribeQ(state);
    realize(state, Stage::Position);
    prescribeU(state);
}

bool prescribeQ(State& state) const;

bool prescribeU(State& state) const;

void getFreeQIndex(const State& state, Array_<SystemQIndex>& freeQs) const;

void getFreeUIndex(const State& state, Array_<SystemUIndex>& freeUs) const;
//------------------------------------------------------------------------------
void handleEvents(State&                        state, 
                  Event::Cause                  cause, 
                  const Array_<EventId>&        eventIds,
                  const HandleEventsOptions&    options,
                  HandleEventsResults&          results) const;
    
void reportEvents(const State&                  state, 
                  Event::Cause                  cause, 
                  const Array_<EventId>&        eventIds) const;

void calcEventTriggerInfo(const State&              state,
                          Array_<EventTriggerInfo>& triggerInfo) const;

void calcTimeOfNextScheduledEvent(const State&      state, 
                                  Real&             tNextEvent, 
                                  Array_<EventId>&  eventIds, 
                                  bool              includeCurrentTime) const;

void calcTimeOfNextScheduledReport(const State&     state, 
                                   Real&            tNextEvent, 
                                   Array_<EventId>& eventIds, 
                                   bool             includeCurrentTime) const;
//------------------------------------------------------------------------------
void relax(State& state, Stage stage, Real accuracy=-1) const;
//------------------------------------------------------------------------------
void multiplyByN(const State& state, const Vector& u, 
                 Vector& dq) const;
void multiplyByNTranspose(const State& state, const Vector& fq, 
                          Vector& fu) const;
void multiplyByNPInv(const State& state, const Vector& dq, 
                     Vector& u) const;
void multiplyByNPInvTranspose(const State& state, const Vector& fu, 
                              Vector& fq) const;
//------------------------------------------------------------------------------
void resetAllCountersToZero();

    // Realization

int getNumRealizationsOfThisStage(Stage) const;

int getNumRealizeCalls() const;

    // Prescribed motion

int getNumPrescribeQCalls() const;
int getNumPrescribeUCalls() const;

    // Projection

int getNumProjectQCalls() const;
int getNumFailedProjectQCalls() const;
int getNumQProjections() const;
int getNumQErrorEstimateProjections() const;

int getNumProjectUCalls() const;
int getNumFailedProjectUCalls() const;
int getNumUProjections() const;
int getNumUErrorEstimateProjections() const;

    // Event handling and reporting

int getNumHandlerCallsThatChangedStage(Stage) const;

int getNumHandleEventCalls() const;

int getNumReportEventCalls() const;
//------------------------------------------------------------------------------
System() : guts(0) { }
System(const System&);
System& operator=(const System&);
~System();

const String& getName()    const;
const String& getVersion() const;

SubsystemIndex adoptSubsystem(Subsystem& child);

int getNumSubsystems() const;
const Subsystem& getSubsystem(SubsystemIndex)   const;
Subsystem& updSubsystem(SubsystemIndex);
const DefaultSystemSubsystem& getDefaultSubsystem() const;
DefaultSystemSubsystem& updDefaultSubsystem();

inline operator const Subsystem&() const; // implemented below
inline operator Subsystem&();

bool systemTopologyHasBeenRealized() const;

StageVersion getSystemTopologyCacheVersion() const;

void setSystemTopologyCacheVersion(StageVersion topoVersion) const;

void invalidateSystemTopologyCache() const;

void calcDecorativeGeometryAndAppend(const State&, Stage, 
                                     Array_<DecorativeGeometry>&) const;


bool isSameSystem(const System& otherSystem) const;


const Guts& getSystemGuts() const {assert(guts); return *guts;}
Guts&       updSystemGuts()       {assert(guts); return *guts;}

void adoptSystemGuts(System::Guts* g);

explicit System(System::Guts* g) : guts(g) { }
bool hasGuts() const {return guts!=0;}

bool isOwnerHandle() const;
bool isEmptyHandle() const;
private:
friend class Guts;
// This is the only data member in this class. Also, any class derived from
// System must have *NO* data members at all (data goes in the Guts class).
Guts*   guts;
};


class SimTK_SimTKCOMMON_EXPORT DefaultSystemSubsystem : public Subsystem {
public:
    explicit DefaultSystemSubsystem(System& sys);
    void addEventHandler(ScheduledEventHandler* handler);
    void addEventHandler(TriggeredEventHandler* handler);
    void addEventReporter(ScheduledEventReporter* handler) const;
    void addEventReporter(TriggeredEventReporter* handler) const;
    EventId createEventId(SubsystemIndex subsys, const State& state) const;
    void findSubsystemEventIds
       (SubsystemIndex subsys, const State& state, 
        const Array_<EventId>& allEvents, 
        Array_<EventId>& eventsForSubsystem) const;
  // don't let doxygen see this private class
    class Guts;
private:
    const Guts& getGuts() const;
    Guts& updGuts();
};

inline void System::addEventHandler(ScheduledEventHandler* handler)
{   updDefaultSubsystem().addEventHandler(handler); }
inline void System::addEventHandler(TriggeredEventHandler* handler)
{   updDefaultSubsystem().addEventHandler(handler); }
inline void System::addEventReporter(ScheduledEventReporter* handler) const
{   getDefaultSubsystem().addEventReporter(handler); }
inline void System::addEventReporter(TriggeredEventReporter* handler) const
{   getDefaultSubsystem().addEventReporter(handler); }

inline System::operator const Subsystem&() const {return getDefaultSubsystem();}
inline System::operator Subsystem&() {return updDefaultSubsystem();}


//==============================================================================
//                     PROJECT OPTIONS and PROJECT RESULTS
//==============================================================================
class ProjectOptions {
public:
    enum Option {
        None            = 0x0000,
        LocalOnly       = 0x0001,
        DontThrow       = 0x0002,
        UseInfinityNorm = 0x0004,
        ForceProjection = 0x0008,
        ForceFullNewton = 0x0010
    };

    ProjectOptions() {clear();}
    explicit ProjectOptions(Real accuracy) 
    {   clear(); setRequiredAccuracy(accuracy); }
    explicit ProjectOptions(Option opt)
    {   clear(); setOption(opt); }

    ProjectOptions& clear() 
    {   optionSet=0; setAccuracyDefaults(); return *this; }

    ProjectOptions& setRequiredAccuracy(Real accuracy) {
        requiredAccuracy = accuracy > 0 ? accuracy 
                                        : getDefaultRequiredAccuracy();
        return *this;
    }

    ProjectOptions& setOvershootFactor(Real overshoot) {
        assert(0 < overshoot && overshoot <= 1);
        desiredOvershoot = overshoot;
        return *this;
    }

    ProjectOptions& setProjectionLimit(Real limit) {
        assert(limit > 0);
        projectionLimit = limit;
        return *this;
    }

    ProjectOptions& clearOption(Option opt) 
    {   optionSet &= ~(unsigned)opt; return *this; }
    ProjectOptions& setOption  (Option opt) 
    {   optionSet |= (unsigned)opt; return *this; }

    Real getRequiredAccuracy()       const {return requiredAccuracy;}
    Real getOvershootFactor() const {return desiredOvershoot;}
    Real getProjectionLimit() const {return projectionLimit;}

    bool isOptionSet(Option opt) const {return (optionSet&(unsigned)opt) != 0;}

    static Real getDefaultRequiredAccuracy() {return Real(1e-4);}
    static Real getDefaultOvershootFactor()  {return Real(0.1);} //i.e., 1e-5

    // Set operators: not, or, and, set difference
    ProjectOptions& operator|=(const ProjectOptions& opts) 
    {   optionSet |= opts.optionSet; return *this; }
    ProjectOptions& operator&=(const ProjectOptions& opts) 
    {   optionSet &= opts.optionSet; return *this; }
    ProjectOptions& operator-=(const ProjectOptions& opts) 
    {   optionSet &= ~opts.optionSet; return *this; }

    ProjectOptions& operator|=(Option opt) {setOption(opt); return *this;}
    ProjectOptions& operator-=(Option opt) {clearOption(opt); return *this;}

private:
    Real     requiredAccuracy;
    Real     desiredOvershoot; // try for accuracy*overshoot
    Real     projectionLimit;  // abort if initial norm is worse than this
    unsigned optionSet;

    void setAccuracyDefaults() {
        requiredAccuracy = getDefaultRequiredAccuracy();
        desiredOvershoot = getDefaultOvershootFactor(); 
        projectionLimit  = Infinity; // we'll try from however far away
    }
};

class ProjectResults {
public:
    ProjectResults() {clear();}

    enum Status {
        Invalid                 = -1,
        Succeeded               = 0,
        FailedToAchieveAccuracy = 1,
        FailedToConverge        = 2    
    };

    ProjectResults& clear() {
        m_exitStatus = Invalid;
        m_anyChangeMade = m_projectionLimitExceeded = false;
        m_numIterations = 0;
        m_worstError = -1;
        m_normOnEntrance = m_normOnExit = NaN;
        return *this;
    }
    bool    isValid()        const {return m_exitStatus != Invalid;}
    Status  getExitStatus()  const {return m_exitStatus;}

    bool getAnyChangeMade()  const {assert(isValid());return m_anyChangeMade;}
    int  getNumIterations()  const {assert(isValid());return m_numIterations;}
    Real getNormOnEntrance() const {assert(isValid());return m_normOnEntrance;}
    Real getNormOnExit()     const {assert(isValid());return m_normOnExit;}
    int  getWorstErrorOnEntrance()    const 
    {   assert(isValid());return m_worstError; }
    bool getProjectionLimitExceeded() const 
    {   assert(isValid());return m_projectionLimitExceeded; }

    ProjectResults& setExitStatus(Status status) 
    {   m_exitStatus=status; return *this; }
    ProjectResults& setAnyChangeMade(bool changeMade) 
    {   m_anyChangeMade=changeMade; return *this; }
    ProjectResults& setProjectionLimitExceeded(bool limitExceeded) 
    {   m_projectionLimitExceeded=limitExceeded; return *this; }
    ProjectResults& setNumIterations(int numIterations) 
    {   m_numIterations=numIterations; return *this; }
    ProjectResults& setNormOnEntrance(Real norm, int worstError) 
    {   m_normOnEntrance=norm; m_worstError=worstError; return *this; }
    ProjectResults& setNormOnExit(Real norm) 
    {   m_normOnExit=norm; return *this; }
private:
    Status  m_exitStatus;
    bool    m_anyChangeMade;
    bool    m_projectionLimitExceeded;
    int     m_numIterations;
    int     m_worstError;       // index of worst error on entrance
    Real    m_normOnEntrance;   // in selected rms or infinity norm
    Real    m_normOnExit;
};



//==============================================================================
//                    REALIZE OPTIONS and REALIZE RESULTS 
//==============================================================================
class RealizeOptions {
    unsigned int optionSet;
    explicit RealizeOptions(unsigned o) : optionSet(o) { }
public:

    enum Option {
        None      = 0x00,
        DontThrow = 0x01
    };


    RealizeOptions() : optionSet(0) { }

    // This is an implicit conversion
    RealizeOptions(Option opt) : optionSet((unsigned)opt) { }

    // Implicit conversion to bool when needed
    operator bool() const {return optionSet != 0;}
    bool isEmpty() const {return optionSet==0;}

    bool isOptionSet(Option opt) const {return (optionSet&(unsigned)opt) != 0;}
    void clear() {optionSet=0;}
    void clearOption(Option opt) {optionSet &= ~(unsigned)opt;}
    void setOption  (Option opt) {optionSet |= (unsigned)opt;}

    // Set operators: or, and
    RealizeOptions& operator|=(RealizeOptions opts) {optionSet |= opts.optionSet; return *this;}
    RealizeOptions& operator&=(RealizeOptions opts) {optionSet &= opts.optionSet; return *this;}

    RealizeOptions& operator|=(Option opt) {setOption(opt); return *this;}
    RealizeOptions& operator-=(Option opt) {clearOption(opt); return *this;}
};

class RealizeResults {
};




} // namespace SimTK

#endif // SimTK_SimTKCOMMON_SYSTEM_H_