Creating Your Own Analysis Part Two

Build a Body Position, Velocity, and Acceleration Analysis

The Analysis from the previous section outputs a body's position. We will now extend it to also output the body's velocity and accelerations. Below, we will show you snippets of the existing code that outputs positions. You should search for that code snippet and use it as a model to add the code necessary to output velocities and acceleration. Note that if you want to use the same set of column labels as are defined here for positions, you need to output the same number of columns as that for positions (6 per body).

The following steps outline the procedure:

Declare additional output storage and internal working arrays. The number of outputs has changed. Before we had one storage file with position data:

/** Storage for recording body positions. */
Storage _storePos; 
/** Internal work array to hold the computed positions. */
Array<double> _bodypos;

Add additional storage for velocities and accelerations and provide the needed working arrays for velocities and accelerations in the .h file.

Update the description of the Analysis in constructDescription() in the .cpp file.
Setup the storage for the velocity and acceleration results, following the example for positions:

setupStorage()
{
  // Positions
  _storePos.reset(0);
  _storePos.setName("Positions");
  _storePos.setDescription(getDescription());
  _storePos.setColumnLabels(getColumnLabels());
  …
}

Correctly size the working arrays :

setModel(Model& aModel)
{ …
  int numBodies = _bodyIndices.getSize();
  _bodypos.setSize(6*numBodies);
  …
}

An Analysis' record()method is the heart of the analysis. It collects and, if necessary, computes the data to output the results of an analysis. In this case, the analysis requires adding a calculation (call to the SimbodyEngine) to get the model accelerations.

/*
* Compute and record the results.
*
* This method, for the purpose of example, records the position and
* orientation of each body in the model. You can customize it
* to perform your analysis.
*
* @param aState Current state of the system.
*/ 

record(const SimTK::State& aState)
{
  … 
  // After setting the state of the model and applying forces
  // Compute the derivative of the multibody system (speeds and
  // accelerations) 

  // POSITION
	const BodySet& bodySet = _model->getBodySet();
	for(int i=0;i<_bodyIndices.getSize();i++) {
		const Body& body = bodySet.get(_bodyIndices[i]);
		SimTK::Vec3 com;
		body.getMassCenter(com);
		// GET POSITIONS AND EULER ANGLES
		_model->getSimbodyEngine().getPosition(s,body,com,vec);
		_model->getSimbodyEngine().getDirectionCosines(s,body,dirCos);
		_model->getSimbodyEngine().convertDirectionCosinesToAngles(dirCos,
			&angVec[0],&angVec[1],&angVec[2]);
		// CONVERT TO DEGREES?
		if(getInDegrees()) {
			angVec *= SimTK_RADIAN_TO_DEGREE;
		}			
		// FILL KINEMATICS ARRAY
		int I=6*i;
		memcpy(&_bodypos[I],&vec[0],3*sizeof(double));
		memcpy(&_bodypos[I+3],&angVec[0],3*sizeof(double));
	}
	_storePos.append(s.getTime(),_bodypos.getSize(),&_bodypos[0]);

// VELOCITY 
…

Repeat this process for velocities and accelerations. Check the Doxygen documents for calls to the SimbodyEngine to get velocities and accelerations.

Note: You need to first realize the velocity and acceleration state before the calculations. For example:

Realize state

_model->getMultibodySystem().realize(s, SimTK::Stage::Acceleration);

In begin() reset the storage objects at the specified time.

// RESET STORAGE 
storePos.reset(s.getTime());

An analysis is finalized by printing the results out to file:

/*
* Print results.
*
* The file names are constructed as
* aDir + "/" + aBaseName + "_" + ComponentName + aExtension
*
* @param aDir Directory in which the results reside.
* @param aBaseName Base file name.
* @param aDT Desired time interval between adjacent storage vectors.
*     Linear interpolation is used to print the data out at the
*     desired interval.
* @param aExtension File extension.
*
* @return 0 on success, -1 on error.
*/

printResults(const string &aBaseName,const string &aDir,double aDT,const string &aExtension)
{
  // POSITIONS
  _storePos.scaleTime(_model->getTimeNormConstant());
  Storage::printResult(&_storePos,aBaseName+"_"+getName()+"_pos",aDir,aDT,aExtension);

  // VELOCITIES
  …

Compile and debug in Visual Studio.
Build install when satisfied (this will overwrite the previous plugin in the OpenSim install directory if you do not change the name)
Repeat the process from the previous section to run and test.

Try plotting pelvis_y with respect of time. You should see a plot like the one below:

Next: Adding New Functionality

Previous: Creating Your Own Analysis Part One

Home: Scripting and Development | Developer's Guide | Creating Your Own Analysis