/
world.h
191 lines (154 loc) · 4.68 KB
/
world.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#ifndef WORLD_H
#define WORLD_H
#include "eigenIncludes.h"
// include elastic rod class
#include "elasticRod.h"
// include force classes
#include "elasticStretchingForce.h"
#include "elasticBendingForce.h"
#include "elasticTwistingForce.h"
#include "externalGravityForce.h"
#include "inertialForce.h"
// include external force
#include "RegularizedStokeslet.h" // RSS-based viscous force
#include "dampingForce.h" // MLRFT-based viscous force
#include "ResistiveForce.h" // RFT-based viscous force
#include "externalContactForce.h"
// include time stepper
#include "timeStepper.h"
// include input file and option
#include "setInput.h"
// include cppflow
#include "cppflow/cppflow.h"
// clamp function is copied from BASim code (Columbia Univ.)
/** Clamps scalar to the range [min,max]. */
template <typename T>
T clamp(T &scalar, T &min, T &max) // template <typename T> inline T clamp( const T& scalar, const T& min, const T& max)
{
if (scalar < min)
return min;
if (scalar > max)
return max;
return scalar;
}
class world
{
public:
world();
world(setInput &m_inputData);
~world();
void setRodStepper();
void updateTimeStep();
int simulationRunning();
int numPoints();
double getScaledCoordinate(int j);
double getCurrentTime();
double getTotalTime();
std::vector<double> Ct, Cp, Cz;
bool isRender();
// file output
void OpenFile(ofstream &outfile);
void CloseFile(ofstream &outfile);
void CoutData(ofstream &outfile);
private:
// Physical parameters
double RodLength;
double helixradius, helixpitchRatio, helixpitch;
double rodRadiusRatio, rodRadius;
double deltaLengthInput;
int numVertices;
double youngM;
double Poisson;
double shearM;
double deltaTime;
double totalTime;
double density;
Vector3d gVector;
double viscosity;
double epsilon;
double distf;
double nTurn;
double flagellaLength, contourRatio;
string modelname;
// Damping force parameters
double eta_per;
double eta_par;
double headSize, omega;
double C_translation, C_rotation; // numerical prefactor to multiply the force against head rotation and translation
bool translation; // option for MLRFT-translation model load.
double tol, stol;
int maxIter; // maximum number of iterations
int maxIterContact;
double characteristicForce;
double forceTol;
// Geometry
MatrixXd vertices;
double currentTime;
// Rod
elasticRod *rod;
// set up the time stepper
timeStepper *stepper;
double *totalForce;
// declare the forces
elasticStretchingForce *m_stretchForce;
elasticBendingForce *m_bendingForce;
elasticTwistingForce *m_twistingForce;
bool useRSS, useRFT; // if true, RSS will be used. RFT will be used, if both false, MLRFT used.
dampingForce *m_dampingForce; // Resistive force theory based model
resistiveForce *m_resistiveForce;
RegularizedStokeslet *m_RegularizedStokeslet; // Regularized Stokeslet Segment (RSS) based model
bool includeContact; // if true, contact will be handled. Otherwise, it will be ignored.
externalContactForce *m_externalContactForce; // contact force
inertialForce *m_inertialForce; // inertial force
externalGravityForce *m_gravityForce; // gravity
int Nstep;
int timeStep;
int iter;
void rodGeometry();
void rodBoundaryCondition();
// Variables about angular velocity
double deltaTwist; // actuation
string inputName; // name of file containing omega
string inputName2; // name of file containing coefficient based RSS
std::vector<double> timeSeries, omegaSeries;
int numOmegaPoints, currentOmegaIndex;
void ReadOmegaData();
void ReadCoefficientData(elasticRod &m_rod);
bool render; // should the OpenGL rendering be included?
bool saveData; // should data be written to a file?
void updateEachRod();
double axisLength, axisLengthInput;
double ClosestPtSegmentSegment(const Vector3d &p1, const Vector3d &q1, const Vector3d &p2, const Vector3d &q2, double &s, double &t, Vector3d &c1, Vector3d &c2);
double smallDist2;
void prepareForContact();
int ne, nv;
int contactNum;
double minDistace;
VectorXd reactionForce;
void computeReactionForce();
// input to NN
std::vector<float> inputNN;
int numVert;
VectorXd s1_NN, s2_NN, Vnorm_NN;
VectorXd VelocityT_NN, VelocityP_NN;
VectorXd ForceT_RSS, ForceN_RSS, ForceZ_RSS;
Vector3d vP_NN, t_NN, u_N, u_NN, p1_NN, p2_NN, p_NN, Vel_NN;
// denormalizing the output coefs
VectorXd Ct_e, Cp_e, Cz_e, Ct_norm, Cp_norm, Cz_norm;
vector<float> output_coef;
Vector3d xCurrent1;
Vector3d xCurrent2;
Vector3d rcurrent;
Vector3d Torque;
double Torquenorm, Forcehead, normTorque, Forcenorm;
Vector3d f_0_rod_0;
Vector3d f_1_rod_0;
Vector3d f_2_rod_0;
Vector3d f_rod_0;
Vector3d f_0_rod_1;
Vector3d f_1_rod_1;
Vector3d f_2_rod_1;
Vector3d f_rod_1;
Vector3d nextpos;
};
#endif