| 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 | 
  | 
 |