80 |
|
|
81 |
|
const unsigned int BufferSize = 65535; |
82 |
|
char buffer[BufferSize]; |
83 |
+ |
Mat3x3d Ddtt; |
84 |
+ |
Mat3x3d Ddtr; |
85 |
+ |
Mat3x3d Ddrr; |
86 |
|
while (ifs.getline(buffer, BufferSize)) { |
87 |
|
StringTokenizer tokenizer(buffer); |
88 |
|
HydroProp currProp; |
89 |
< |
if (tokenizer.countTokens() >= 67) { |
89 |
> |
if (tokenizer.countTokens() >= 40) { |
90 |
|
std::string atomName = tokenizer.nextToken(); |
91 |
< |
currProp.cod[0] = tokenizer.nextTokenAsDouble(); |
92 |
< |
currProp.cod[1] = tokenizer.nextTokenAsDouble(); |
93 |
< |
currProp.cod[2] = tokenizer.nextTokenAsDouble(); |
91 |
> |
currProp.cor[0] = tokenizer.nextTokenAsDouble(); |
92 |
> |
currProp.cor[1] = tokenizer.nextTokenAsDouble(); |
93 |
> |
currProp.cor[2] = tokenizer.nextTokenAsDouble(); |
94 |
|
|
92 |
– |
currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble(); |
93 |
– |
currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble(); |
94 |
– |
currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble(); |
95 |
– |
currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble(); |
96 |
– |
currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble(); |
97 |
– |
currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble(); |
98 |
– |
currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble(); |
99 |
– |
currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble(); |
100 |
– |
currProp.Ddtt(2,2) = tokenizer.nextTokenAsDouble(); |
95 |
|
|
96 |
< |
currProp.Ddtr(0,0) = tokenizer.nextTokenAsDouble(); |
97 |
< |
currProp.Ddtr(0,1) = tokenizer.nextTokenAsDouble(); |
98 |
< |
currProp.Ddtr(0,2) = tokenizer.nextTokenAsDouble(); |
99 |
< |
currProp.Ddtr(1,0) = tokenizer.nextTokenAsDouble(); |
100 |
< |
currProp.Ddtr(1,1) = tokenizer.nextTokenAsDouble(); |
101 |
< |
currProp.Ddtr(1,2) = tokenizer.nextTokenAsDouble(); |
102 |
< |
currProp.Ddtr(2,0) = tokenizer.nextTokenAsDouble(); |
103 |
< |
currProp.Ddtr(2,1) = tokenizer.nextTokenAsDouble(); |
104 |
< |
currProp.Ddtr(2,2) = tokenizer.nextTokenAsDouble(); |
96 |
> |
Ddtt(0,0) = tokenizer.nextTokenAsDouble(); |
97 |
> |
Ddtt(0,1) = tokenizer.nextTokenAsDouble(); |
98 |
> |
Ddtt(0,2) = tokenizer.nextTokenAsDouble(); |
99 |
> |
Ddtt(1,0) = tokenizer.nextTokenAsDouble(); |
100 |
> |
Ddtt(1,1) = tokenizer.nextTokenAsDouble(); |
101 |
> |
Ddtt(1,2) = tokenizer.nextTokenAsDouble(); |
102 |
> |
Ddtt(2,0) = tokenizer.nextTokenAsDouble(); |
103 |
> |
Ddtt(2,1) = tokenizer.nextTokenAsDouble(); |
104 |
> |
Ddtt(2,2) = tokenizer.nextTokenAsDouble(); |
105 |
|
|
106 |
< |
currProp.Ddrr(0,0) = tokenizer.nextTokenAsDouble(); |
107 |
< |
currProp.Ddrr(0,1) = tokenizer.nextTokenAsDouble(); |
108 |
< |
currProp.Ddrr(0,2) = tokenizer.nextTokenAsDouble(); |
109 |
< |
currProp.Ddrr(1,0) = tokenizer.nextTokenAsDouble(); |
110 |
< |
currProp.Ddrr(1,1) = tokenizer.nextTokenAsDouble(); |
111 |
< |
currProp.Ddrr(1,2) = tokenizer.nextTokenAsDouble(); |
112 |
< |
currProp.Ddrr(2,0) = tokenizer.nextTokenAsDouble(); |
113 |
< |
currProp.Ddrr(2,1) = tokenizer.nextTokenAsDouble(); |
114 |
< |
currProp.Ddrr(2,2) = tokenizer.nextTokenAsDouble(); |
106 |
> |
Ddtr(0,0) = tokenizer.nextTokenAsDouble(); |
107 |
> |
Ddtr(0,1) = tokenizer.nextTokenAsDouble(); |
108 |
> |
Ddtr(0,2) = tokenizer.nextTokenAsDouble(); |
109 |
> |
Ddtr(1,0) = tokenizer.nextTokenAsDouble(); |
110 |
> |
Ddtr(1,1) = tokenizer.nextTokenAsDouble(); |
111 |
> |
Ddtr(1,2) = tokenizer.nextTokenAsDouble(); |
112 |
> |
Ddtr(2,0) = tokenizer.nextTokenAsDouble(); |
113 |
> |
Ddtr(2,1) = tokenizer.nextTokenAsDouble(); |
114 |
> |
Ddtr(2,2) = tokenizer.nextTokenAsDouble(); |
115 |
|
|
116 |
< |
currProp.Xidtt(0,0) = tokenizer.nextTokenAsDouble(); |
117 |
< |
currProp.Xidtt(0,1) = tokenizer.nextTokenAsDouble(); |
118 |
< |
currProp.Xidtt(0,2) = tokenizer.nextTokenAsDouble(); |
119 |
< |
currProp.Xidtt(1,0) = tokenizer.nextTokenAsDouble(); |
120 |
< |
currProp.Xidtt(1,1) = tokenizer.nextTokenAsDouble(); |
121 |
< |
currProp.Xidtt(1,2) = tokenizer.nextTokenAsDouble(); |
122 |
< |
currProp.Xidtt(2,0) = tokenizer.nextTokenAsDouble(); |
123 |
< |
currProp.Xidtt(2,1) = tokenizer.nextTokenAsDouble(); |
124 |
< |
currProp.Xidtt(2,2) = tokenizer.nextTokenAsDouble(); |
116 |
> |
Ddrr(0,0) = tokenizer.nextTokenAsDouble(); |
117 |
> |
Ddrr(0,1) = tokenizer.nextTokenAsDouble(); |
118 |
> |
Ddrr(0,2) = tokenizer.nextTokenAsDouble(); |
119 |
> |
Ddrr(1,0) = tokenizer.nextTokenAsDouble(); |
120 |
> |
Ddrr(1,1) = tokenizer.nextTokenAsDouble(); |
121 |
> |
Ddrr(1,2) = tokenizer.nextTokenAsDouble(); |
122 |
> |
Ddrr(2,0) = tokenizer.nextTokenAsDouble(); |
123 |
> |
Ddrr(2,1) = tokenizer.nextTokenAsDouble(); |
124 |
> |
Ddrr(2,2) = tokenizer.nextTokenAsDouble(); |
125 |
|
|
126 |
< |
currProp.Xidrt(0,0) = tokenizer.nextTokenAsDouble(); |
127 |
< |
currProp.Xidrt(0,1) = tokenizer.nextTokenAsDouble(); |
128 |
< |
currProp.Xidrt(0,2) = tokenizer.nextTokenAsDouble(); |
129 |
< |
currProp.Xidrt(1,0) = tokenizer.nextTokenAsDouble(); |
130 |
< |
currProp.Xidrt(1,1) = tokenizer.nextTokenAsDouble(); |
131 |
< |
currProp.Xidrt(1,2) = tokenizer.nextTokenAsDouble(); |
132 |
< |
currProp.Xidrt(2,0) = tokenizer.nextTokenAsDouble(); |
133 |
< |
currProp.Xidrt(2,1) = tokenizer.nextTokenAsDouble(); |
134 |
< |
currProp.Xidrt(2,2) = tokenizer.nextTokenAsDouble(); |
141 |
< |
|
142 |
< |
currProp.Xidtr(0,0) = tokenizer.nextTokenAsDouble(); |
143 |
< |
currProp.Xidtr(0,1) = tokenizer.nextTokenAsDouble(); |
144 |
< |
currProp.Xidtr(0,2) = tokenizer.nextTokenAsDouble(); |
145 |
< |
currProp.Xidtr(1,0) = tokenizer.nextTokenAsDouble(); |
146 |
< |
currProp.Xidtr(1,1) = tokenizer.nextTokenAsDouble(); |
147 |
< |
currProp.Xidtr(1,2) = tokenizer.nextTokenAsDouble(); |
148 |
< |
currProp.Xidtr(2,0) = tokenizer.nextTokenAsDouble(); |
149 |
< |
currProp.Xidtr(2,1) = tokenizer.nextTokenAsDouble(); |
150 |
< |
currProp.Xidtr(2,2) = tokenizer.nextTokenAsDouble(); |
126 |
> |
currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble(); |
127 |
> |
currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble(); |
128 |
> |
currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble(); |
129 |
> |
currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble(); |
130 |
> |
currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble(); |
131 |
> |
currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble(); |
132 |
> |
currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble(); |
133 |
> |
currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble(); |
134 |
> |
currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble(); |
135 |
|
|
136 |
< |
currProp.Xidrr(0,0) = tokenizer.nextTokenAsDouble(); |
137 |
< |
currProp.Xidrr(0,1) = tokenizer.nextTokenAsDouble(); |
138 |
< |
currProp.Xidrr(0,2) = tokenizer.nextTokenAsDouble(); |
139 |
< |
currProp.Xidrr(1,0) = tokenizer.nextTokenAsDouble(); |
140 |
< |
currProp.Xidrr(1,1) = tokenizer.nextTokenAsDouble(); |
141 |
< |
currProp.Xidrr(1,2) = tokenizer.nextTokenAsDouble(); |
142 |
< |
currProp.Xidrr(2,0) = tokenizer.nextTokenAsDouble(); |
143 |
< |
currProp.Xidrr(2,1) = tokenizer.nextTokenAsDouble(); |
144 |
< |
currProp.Xidrr(2,2) = tokenizer.nextTokenAsDouble(); |
136 |
> |
currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble(); |
137 |
> |
currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble(); |
138 |
> |
currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble(); |
139 |
> |
currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble(); |
140 |
> |
currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble(); |
141 |
> |
currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble(); |
142 |
> |
currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble(); |
143 |
> |
currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble(); |
144 |
> |
currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble(); |
145 |
> |
|
146 |
> |
currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble(); |
147 |
> |
currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble(); |
148 |
> |
currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble(); |
149 |
> |
currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble(); |
150 |
> |
currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble(); |
151 |
> |
currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble(); |
152 |
> |
currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble(); |
153 |
> |
currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble(); |
154 |
> |
currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble(); |
155 |
> |
|
156 |
> |
currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble(); |
157 |
> |
currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble(); |
158 |
> |
currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble(); |
159 |
> |
currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble(); |
160 |
> |
currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble(); |
161 |
> |
currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble(); |
162 |
> |
currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble(); |
163 |
> |
currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble(); |
164 |
> |
currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble(); |
165 |
> |
|
166 |
> |
SquareMatrix<double, 6> Xir; |
167 |
> |
Xir.setSubMatrix(0, 0, currProp.Xirtt); |
168 |
> |
Xir.setSubMatrix(0, 3, currProp.Xirrt); |
169 |
> |
Xir.setSubMatrix(3, 0, currProp.Xirtr); |
170 |
> |
Xir.setSubMatrix(3, 3, currProp.Xirrr); |
171 |
> |
CholeskyDecomposition(Xir, currProp.S); |
172 |
> |
|
173 |
|
props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp)); |
174 |
|
} |
175 |
|
} |
215 |
|
omega[2] = angMom[2] /I(2, 2); |
216 |
|
} |
217 |
|
|
218 |
< |
//apply friction force and torque at center of diffusion |
218 |
> |
//apply friction force and torque at center of resistance |
219 |
|
A = integrableObject->getA(); |
220 |
|
Atrans = A.transpose(); |
221 |
< |
Vector3d rcd = Atrans * hydroProps_[index].cod; |
222 |
< |
Vector3d vcd = vel + cross(omega, rcd); |
223 |
< |
vcd = A* vcd; |
224 |
< |
Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega); |
225 |
< |
frictionForce = Atrans*frictionForce; |
226 |
< |
integrableObject->addFrc(frictionForce); |
227 |
< |
Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega); |
228 |
< |
frictionTorque = Atrans*frictionTorque; |
229 |
< |
integrableObject->addTrq(frictionTorque+ cross(rcd, frictionForce)); |
230 |
< |
|
231 |
< |
//apply random force and torque at center of diffustion |
232 |
< |
Vector3d randomForce; |
233 |
< |
Vector3d randomTorque; |
234 |
< |
genRandomForceAndTorque(randomForce, randomTorque, index, variance_); |
235 |
< |
randomForce = Atrans*randomForce; |
236 |
< |
randomTorque = Atrans* randomTorque; |
237 |
< |
integrableObject->addFrc(randomForce); |
238 |
< |
integrableObject->addTrq(randomTorque + cross(rcd, randomForce )); |
239 |
< |
|
221 |
> |
Vector3d rcr = Atrans * hydroProps_[index].cor; |
222 |
> |
Vector3d vcdLab = vel + cross(omega, rcr); |
223 |
> |
Vector3d vcdBody = A* vcdLab; |
224 |
> |
Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega); |
225 |
> |
Vector3d frictionForceLab = Atrans*frictionForceBody; |
226 |
> |
integrableObject->addFrc(frictionForceLab); |
227 |
> |
Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega); |
228 |
> |
Vector3d frictionTorqueLab = Atrans*frictionTorqueBody; |
229 |
> |
integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab)); |
230 |
> |
|
231 |
> |
//apply random force and torque at center of resistance |
232 |
> |
Vector3d randomForceBody; |
233 |
> |
Vector3d randomTorqueBody; |
234 |
> |
genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_); |
235 |
> |
Vector3d randomForceLab = Atrans*randomForceBody; |
236 |
> |
Vector3d randomTorqueLab = Atrans* randomTorqueBody; |
237 |
> |
integrableObject->addFrc(randomForceLab); |
238 |
> |
integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab )); |
239 |
> |
|
240 |
|
} else { |
241 |
|
//spheric atom |
242 |
< |
Vector3d frictionForce = -(hydroProps_[index].Xidtt *vel); |
242 |
> |
Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel); |
243 |
|
Vector3d randomForce; |
244 |
|
Vector3d randomTorque; |
245 |
|
genRandomForceAndTorque(randomForce, randomTorque, index, variance_); |
246 |
|
|
235 |
– |
//randomForce /= OOPSEConstant::energyConvert; |
236 |
– |
//randomTorque /= OOPSEConstant::energyConvert; |
247 |
|
integrableObject->addFrc(frictionForce+randomForce); |
248 |
|
} |
249 |
|
|
259 |
|
} |
260 |
|
|
261 |
|
void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) { |
252 |
– |
/* |
253 |
– |
SquareMatrix<double, 6> Dd; |
254 |
– |
SquareMatrix<double, 6> S; |
255 |
– |
Vector<double, 6> Z; |
256 |
– |
Vector<double, 6> generalForce; |
257 |
– |
Dd.setSubMatrix(0, 0, hydroProps_[index].Ddtt); |
258 |
– |
Dd.setSubMatrix(0, 3, hydroProps_[index].Ddtr.transpose()); |
259 |
– |
Dd.setSubMatrix(3, 0, hydroProps_[index].Ddtr); |
260 |
– |
Dd.setSubMatrix(3, 3, hydroProps_[index].Ddrr); |
261 |
– |
CholeskyDecomposition(Dd, S); |
262 |
– |
*/ |
262 |
|
|
263 |
< |
SquareMatrix<double, 6> Xid; |
265 |
< |
SquareMatrix<double, 6> S; |
263 |
> |
|
264 |
|
Vector<double, 6> Z; |
265 |
|
Vector<double, 6> generalForce; |
268 |
– |
Xid.setSubMatrix(0, 0, hydroProps_[index].Xidtt); |
269 |
– |
Xid.setSubMatrix(0, 3, hydroProps_[index].Xidrt); |
270 |
– |
Xid.setSubMatrix(3, 0, hydroProps_[index].Xidtr); |
271 |
– |
Xid.setSubMatrix(3, 3, hydroProps_[index].Xidrr); |
272 |
– |
CholeskyDecomposition(Xid, S); |
266 |
|
|
274 |
– |
/* |
275 |
– |
Xid *= variance; |
276 |
– |
Z[0] = randNumGen_.randNorm(0, 1.0); |
277 |
– |
Z[1] = randNumGen_.randNorm(0, 1.0); |
278 |
– |
Z[2] = randNumGen_.randNorm(0, 1.0); |
279 |
– |
Z[3] = randNumGen_.randNorm(0, 1.0); |
280 |
– |
Z[4] = randNumGen_.randNorm(0, 1.0); |
281 |
– |
Z[5] = randNumGen_.randNorm(0, 1.0); |
282 |
– |
*/ |
267 |
|
|
268 |
|
Z[0] = randNumGen_.randNorm(0, variance); |
269 |
|
Z[1] = randNumGen_.randNorm(0, variance); |
273 |
|
Z[5] = randNumGen_.randNorm(0, variance); |
274 |
|
|
275 |
|
|
276 |
< |
generalForce = S*Z; |
276 |
> |
generalForce = hydroProps_[index].S*Z; |
277 |
|
|
278 |
|
force[0] = generalForce[0]; |
279 |
|
force[1] = generalForce[1]; |