1 |
|
#include <iostream> |
2 |
|
#include <stdlib.h> |
3 |
|
#include <math.h> |
4 |
< |
|
4 |
> |
#include "Rattle.hpp" |
5 |
> |
#include "Roll.hpp" |
6 |
|
#ifdef IS_MPI |
7 |
|
#include "mpiSimulation.hpp" |
8 |
|
#include <unistd.h> |
33 |
|
|
34 |
|
nAtoms = info->n_atoms; |
35 |
|
integrableObjects = info->integrableObjects; |
36 |
< |
|
36 |
> |
|
37 |
> |
consFramework = new RattleFramework(info); |
38 |
> |
|
39 |
> |
if(consFramework == NULL){ |
40 |
> |
sprintf(painCave.errMsg, |
41 |
> |
"Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); |
42 |
> |
painCave.isFatal = 1; |
43 |
> |
simError(); |
44 |
> |
} |
45 |
> |
|
46 |
> |
/* |
47 |
|
// check for constraints |
48 |
|
|
49 |
|
constrainedA = NULL; |
56 |
|
nConstrained = 0; |
57 |
|
|
58 |
|
checkConstraints(); |
59 |
< |
|
59 |
> |
*/ |
60 |
|
} |
61 |
|
|
62 |
|
template<typename T> Integrator<T>::~Integrator(){ |
63 |
+ |
if (consFramework != NULL) |
64 |
+ |
delete consFramework; |
65 |
+ |
/* |
66 |
|
if (nConstrained){ |
67 |
|
delete[] constrainedA; |
68 |
|
delete[] constrainedB; |
71 |
|
delete[] moved; |
72 |
|
delete[] oldPos; |
73 |
|
} |
74 |
+ |
*/ |
75 |
|
} |
76 |
|
|
77 |
+ |
/* |
78 |
|
template<typename T> void Integrator<T>::checkConstraints(void){ |
79 |
|
isConstrained = 0; |
80 |
|
|
109 |
|
if (constrained){ |
110 |
|
dummy_plug = theArray[j]->get_constraint(); |
111 |
|
temp_con[nConstrained].set_a(dummy_plug->get_a()); |
112 |
< |
temp_con[nConstrained].set_b(dummy_plug->get_b()); |
112 |
> |
temp_con[nConstrained].set_b(Dummy_plug->get_b()); |
113 |
|
temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr()); |
114 |
|
|
115 |
|
nConstrained++; |
167 |
|
|
168 |
|
delete[] temp_con; |
169 |
|
} |
170 |
+ |
*/ |
171 |
|
|
155 |
– |
|
172 |
|
template<typename T> void Integrator<T>::integrate(void){ |
173 |
|
|
174 |
|
double runTime = info->run_time; |
208 |
|
// initialize the forces before the first step |
209 |
|
|
210 |
|
calcForce(1, 1); |
211 |
< |
|
212 |
< |
if (nConstrained){ |
213 |
< |
preMove(); |
214 |
< |
constrainA(); |
215 |
< |
calcForce(1, 1); |
216 |
< |
constrainB(); |
201 |
< |
} |
211 |
> |
|
212 |
> |
//execute constraint algorithm to make sure at the very beginning the system is constrained |
213 |
> |
consFramework->doPreConstraint(); |
214 |
> |
consFramework->doConstrainA(); |
215 |
> |
calcForce(1, 1); |
216 |
> |
consFramework->doConstrainB(); |
217 |
|
|
218 |
|
if (info->setTemp){ |
219 |
|
thermalize(); |
270 |
|
|
271 |
|
if (info->getTime() >= currStatus){ |
272 |
|
statOut->writeStat(info->getTime()); |
258 |
– |
if (info->useSolidThermInt || info->useLiquidThermInt) |
259 |
– |
statOut->writeRaw(info->getTime()); |
273 |
|
calcPot = 0; |
274 |
|
calcStress = 0; |
275 |
|
currStatus += statusTime; |
309 |
|
startProfile(pro3); |
310 |
|
#endif //profile |
311 |
|
|
312 |
< |
preMove(); |
312 |
> |
//save old state (position, velocity etc) |
313 |
> |
consFramework->doPreConstraint(); |
314 |
|
|
315 |
|
#ifdef PROFILE |
316 |
|
endProfile(pro3); |
406 |
|
} |
407 |
|
} |
408 |
|
|
409 |
< |
if (nConstrained){ |
396 |
< |
constrainA(); |
397 |
< |
} |
409 |
> |
consFramework->doConstrainA(); |
410 |
|
} |
411 |
|
|
412 |
|
|
447 |
|
} |
448 |
|
} |
449 |
|
|
450 |
< |
if (nConstrained){ |
439 |
< |
constrainB(); |
440 |
< |
} |
450 |
> |
consFramework->doConstrainB(); |
451 |
|
} |
452 |
|
|
453 |
+ |
/* |
454 |
|
template<typename T> void Integrator<T>::preMove(void){ |
455 |
|
int i, j; |
456 |
|
double pos[3]; |
709 |
|
simError(); |
710 |
|
} |
711 |
|
} |
712 |
< |
|
712 |
> |
*/ |
713 |
|
template<typename T> void Integrator<T>::rotationPropagation |
714 |
|
( StuntDouble* sd, double ji[3] ){ |
715 |
|
|