1 /* -------------------------------------------------------------------------- *
2 * OpenSim: mainFatigue.cpp *
3 * -------------------------------------------------------------------------- *
4 * The OpenSim API is a toolkit for musculoskeletal modeling and simulation. *
5 * See http://opensim.stanford.edu and the NOTICE file for more information. *
6 * OpenSim is developed at Stanford University and supported by the US *
7 * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA *
8 * through the Warrior Web program. *
9 * *
10 * Copyright (c) 2005-2017 Stanford University and the Authors *
11 * Author(s): Peter Loan, Ajay Seth, Ayman Habib *
12 * *
13 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
14 * not use this file except in compliance with the License. You may obtain a *
15 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
16 * *
17 * Unless required by applicable law or agreed to in writing, software *
18 * distributed under the License is distributed on an "AS IS" BASIS, *
19 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied *
20 * See the License for the specific language governing permissions and *
21 * limitations under the License. *
22 * -------------------------------------------------------------------------- */
23
24 /*
25 * Below is an example of an OpenSim main() routine.
26 * This program is a forward simulation of a tug-of-war between two muscles
27 * pulling on a block. One of the muscles fatigues and the other does not.
28 */
29
30 //=============================================================================
31 //=============================================================================
32 #include <OpenSim/OpenSim.h>
33 #include "FatigableMuscle.h"
34 #include <OpenSim/Common/IO.h>
35 #include "OpenSim/Common/STOFileAdapter.h"
36
37 using namespace OpenSim;
38 using namespace SimTK;
39
40 //_____________________________________________________________________________
41 /**
42 * Run a simulation of a sliding block being pulled by two muscle
43 */
main()44 int main()
45 {
46 std::clock_t startTime = std::clock();
47
48 try {
49 ///////////////////////////////////////////////
50 // DEFINE THE SIMULATION START AND END TIMES //
51 ///////////////////////////////////////////////
52 // Define the initial and final simulation times
53 double initialTime = 0.0;
54 double finalTime = 10.0;
55
56 ///////////////////////////////////////////
57 // DEFINE BODIES AND JOINTS OF THE MODEL //
58 ///////////////////////////////////////////
59 // Create an OpenSim model and set its name
60 Model osimModel;
61 osimModel.setName("tugOfWar");
62
63 // GROUND FRAME
64
65 // Get a reference to the model's ground body
66 Ground& ground = osimModel.updGround();
67
68 // Add display geometry to the ground to visualize in the GUI
69 ground.attachGeometry(new Mesh("ground.vtp"));
70 ground.attachGeometry(new Mesh("anchor1.vtp"));
71 ground.attachGeometry(new Mesh("anchor2.vtp"));
72
73 // BLOCK BODY
74
75 // Specify properties of a 20 kg, 10cm length block body
76 double blockMass = 20.0, blockSideLength = 0.1;
77 Vec3 blockMassCenter(0);
78 Inertia blockInertia = blockMass*Inertia::brick(blockSideLength,
79 blockSideLength, blockSideLength);
80
81 // Create a new block body with the specified properties
82 OpenSim::Body *block = new OpenSim::Body("block", blockMass,
83 blockMassCenter, blockInertia);
84
85 // Add display geometry to the block to visualize in the GUI
86 block->attachGeometry(new Mesh("block.vtp"));
87
88 // FREE JOINT
89
90 // Create a new free joint with 6 degrees-of-freedom (coordinates)
91 // between the block and ground bodies
92 double halfLength = blockSideLength/2.0;
93 Vec3 locationInParent(0, halfLength, 0), orientationInParent(0);
94 Vec3 locationInBody(0, halfLength, 0), orientationInBody(0);
95 FreeJoint *blockToGround = new FreeJoint("blockToGround", ground,
96 locationInParent, orientationInParent,
97 *block, locationInBody, orientationInBody);
98
99 // Set the angle and position ranges for the free (6-degree-of-freedom)
100 // joint between the block and ground frames.
101 double angleRange[2] = {-SimTK::Pi/2, SimTK::Pi/2};
102 double positionRange[2] = {-1, 1};
103 blockToGround->updCoordinate(FreeJoint::Coord::Rotation1X).setRange(angleRange);
104 blockToGround->updCoordinate(FreeJoint::Coord::Rotation2Y).setRange(angleRange);
105 blockToGround->updCoordinate(FreeJoint::Coord::Rotation3Z).setRange(angleRange);
106 blockToGround->updCoordinate(FreeJoint::Coord::TranslationX).setRange(positionRange);
107 blockToGround->updCoordinate(FreeJoint::Coord::TranslationY).setRange(positionRange);
108 blockToGround->updCoordinate(FreeJoint::Coord::TranslationZ).setRange(positionRange);
109
110 // Add the block body to the model
111 osimModel.addBody(block);
112 osimModel.addJoint(blockToGround);
113
114 ///////////////////////////////////////
115 // DEFINE FORCES ACTING ON THE MODEL //
116 ///////////////////////////////////////
117 // MUSCLE FORCES
118
119 // Create two new muscles
120 double maxIsometricForce = 1000.0, optimalFiberLength = 0.2,
121 tendonSlackLength = 0.1, pennationAngle = 0.0,
122 fatigueFactor = 0.30, recoveryFactor = 0.20;
123
124 // fatigable muscle (Millard2012EquilibriumMuscle with fatigue)
125 FatigableMuscle* fatigable = new FatigableMuscle("fatigable",
126 maxIsometricForce, optimalFiberLength, tendonSlackLength,
127 pennationAngle, fatigueFactor, recoveryFactor);
128
129 // original muscle model (muscle without fatigue)
130 Millard2012EquilibriumMuscle* original =
131 new Millard2012EquilibriumMuscle("original",
132 maxIsometricForce, optimalFiberLength, tendonSlackLength,
133 pennationAngle);
134
135 // Define the path of the muscles
136 fatigable->addNewPathPoint("fatigable-point1", ground,
137 Vec3(0.0, halfLength, -0.35));
138 fatigable->addNewPathPoint("fatigable-point2", *block,
139 Vec3(0.0, halfLength, -halfLength));
140
141 original->addNewPathPoint("original-point1", ground,
142 Vec3(0.0, halfLength, 0.35));
143 original->addNewPathPoint("original-point2", *block,
144 Vec3(0.0, halfLength, halfLength));
145
146 // Define the default states for the two muscles
147 // Activation
148 fatigable->setDefaultActivation(0.01);
149 original->setDefaultActivation(0.01);
150 // Fiber length
151 fatigable->setDefaultFiberLength(optimalFiberLength);
152 original->setDefaultFiberLength(optimalFiberLength);
153
154 // Add the two muscles (as forces) to the model
155 osimModel.addForce(fatigable);
156 osimModel.addForce(original);
157
158 ///////////////////////////////////
159 // DEFINE CONTROLS FOR THE MODEL //
160 ///////////////////////////////////
161 // Create a prescribed controller that simply supplies controls as
162 // a function of time.
163 // For muscles, controls are normalized stoor-neuron excitations
164 PrescribedController *muscleController = new PrescribedController();
165 muscleController->setActuators(osimModel.updActuators());
166
167 // Set the prescribed muscle controller to use the same muscle control function for each muscle
168 muscleController->prescribeControlForActuator("fatigable", new Constant(1.0));
169 muscleController->prescribeControlForActuator("original", new Constant(1.0));
170
171 // Add the muscle controller to the model
172 osimModel.addController(muscleController);
173
174 // Add a Muscle analysis
175 MuscleAnalysis* muscAnalysis = new MuscleAnalysis(&osimModel);
176 Array<std::string> coords(blockToGround->getCoordinate(FreeJoint::Coord::TranslationZ).getName(),1);
177 muscAnalysis->setCoordinates(coords);
178 muscAnalysis->setComputeMoments(false);
179 osimModel.addAnalysis(muscAnalysis);
180
181 // Turn on the visualizer to view the simulation run live.
182 osimModel.setUseVisualizer(false);
183
184 //////////////////////////
185 // PERFORM A SIMULATION //
186 //////////////////////////
187
188 // Initialize the system and get the state
189 SimTK::State& si = osimModel.initSystem();
190
191 // Init coords to 0 and lock the rotational degrees of freedom so the block doesn't twist
192 CoordinateSet& coordinates = osimModel.updCoordinateSet();
193 coordinates[0].setValue(si, 0);
194 coordinates[1].setValue(si, 0);
195 coordinates[2].setValue(si, 0);
196 coordinates[3].setValue(si, 0);
197 coordinates[4].setValue(si, 0);
198 coordinates[5].setValue(si, 0);
199 coordinates[0].setLocked(si, true);
200 coordinates[1].setLocked(si, true);
201 coordinates[2].setLocked(si, true);
202 // Last coordinate (index 5) is the Z translation of the block
203 coordinates[4].setLocked(si, true);
204
205 // Compute initial conditions for muscles
206 osimModel.equilibrateMuscles(si);
207
208 // Create the force reporter
209 ForceReporter* reporter = new ForceReporter(&osimModel);
210 osimModel.updAnalysisSet().adoptAndAppend(reporter);
211 // Create the manager
212 Manager manager(osimModel);
213 manager.setIntegratorAccuracy(1.0e-6);
214
215 // Print out details of the model
216 osimModel.printDetailedInfo(si, std::cout);
217
218 // Integrate from initial time to final time
219 si.setTime(initialTime);
220 manager.initialize(si);
221 std::cout<<"\nIntegrating from "<<initialTime<<" to "<<finalTime<<std::endl;
222 manager.integrate(finalTime);
223
224 //////////////////////////////
225 // SAVE THE RESULTS TO FILE //
226 //////////////////////////////
227
228 // Save the simulation results
229 // Save the states
230 auto statesTable = manager.getStatesTable();
231 STOFileAdapter_<double>::write(statesTable,
232 "tugOfWar_fatigue_states.sto");
233
234 auto forcesTable = reporter->getForcesTable();
235 STOFileAdapter_<double>::write(forcesTable,
236 "tugOfWar_fatigue_forces.sto");
237
238 // Save the muscle analysis results
239 IO::makeDir("MuscleAnalysisResults");
240 muscAnalysis->printResults("fatigue", "MuscleAnalysisResults");
241
242 // To print (serialize) the latest connections of the model, it is
243 // necessary to finalizeConnections() first.
244 osimModel.finalizeConnections();
245 // Save the OpenSim model to a file
246 osimModel.print("tugOfWar_fatigue_model.osim");
247 }
248 catch (const std::exception& ex)
249 {
250 std::cout << ex.what() << std::endl;
251 return 1;
252 }
253 catch (...)
254 {
255 std::cout << "UNRECOGNIZED EXCEPTION" << std::endl;
256 return 1;
257 }
258
259 std::cout << "main() routine time = " << 1.e3*(std::clock()-startTime)/CLOCKS_PER_SEC << "ms\n";
260
261 std::cout << "OpenSim example completed successfully.\n";
262 return 0;
263 }
264