Homec4science
Diffusion lammps (lammps-icms) Browse

lammps/lib/poems/workspace.cpp
237359ca79aelammps-icms

Last Change
Raw File

workspace.cpp
View Options

/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: workspace.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "workspace.h"
#include "system.h"
#include "solver.h"
#include "SystemProcessor.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
using namespace std;
void Workspace::allocateNewSystem() {
currentIndex++;
if(currentIndex < maxAlloc)
{
system[currentIndex].system = new System;
}
else
{
maxAlloc = (maxAlloc + 1) * 2;
SysData * tempPtrSys = new SysData[maxAlloc];
for(int i = 0; i < currentIndex; i++)
{
tempPtrSys[i] = system[i];
}
delete [] system;
system = tempPtrSys;
system[currentIndex].system = new System;
}
}
Workspace::Workspace(){
currentIndex = -1;
maxAlloc = 0;
system = NULL;
}
Workspace::~Workspace(){
for(int i = 0; i <= currentIndex; i++){
delete system[i].system;
}
delete [] system;
}
bool Workspace::LoadFile(char* filename){
bool ans;
ifstream file;
file.open(filename, ifstream::in);
if( !file.is_open() ){
cerr << "File '" << filename << "' not found." << endl;
return false;
}
allocateNewSystem();
ans = system[currentIndex].system->ReadIn(file);
file.close();
return ans;
}
void Workspace::SetLammpsValues(double dtv, double dthalf, double tempcon){
Thalf = dthalf;
Tfull = dtv;
ConFac = tempcon;
}
bool Workspace::MakeSystem(int& nbody, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space, int &njoint, int **&jointbody, double **&xjoint, int& nfree, int*freelist, double dthalf, double dtv, double tempcon, double KE){
SetLammpsValues(dtv, dthalf, tempcon);
if(njoint){
SystemProcessor sys;
sys.processArray(jointbody, njoint);
List<POEMSChain> * results = sys.getSystemData();
int numsyschains = results->GetNumElements();
int headvalue = 0;
List<POEMSChain> * newresults = results;
ListElement<POEMSChain> * tempNode = results->GetHeadElement();
int stop = 1;
int counter = 1;
for(int n = 1; n<=numsyschains; n++){
while(stop){
if ( (*(tempNode->value->listOfNodes.GetHeadElement()->value) == (headvalue+1) ) || (*(tempNode->value->listOfNodes.GetTailElement()->value) == (headvalue+1) ) ) {
newresults->Append(tempNode->value);
headvalue = headvalue + tempNode->value->listOfNodes.GetNumElements();
tempNode = results->GetHeadElement();
stop = 0;
counter ++;
}
else{
tempNode = tempNode->next;
}
}
stop=1;
}
ListElement<POEMSChain> * TNode = newresults->GetHeadElement();
ListElement<POEMSChain> * TTNode = newresults->GetHeadElement();
for(int kk=1; kk<=numsyschains; kk++){
if(kk!=numsyschains)
TTNode = TNode->next;
newresults->Remove(TNode);
if(kk!=numsyschains)
TNode = TTNode;
}
ListElement<POEMSChain> * NodeValue = newresults->GetHeadElement();
int count = 0;
int * array;
int ** arrayFromChain;
int numElementsInSystem;
int ttk = 0;
while(NodeValue != NULL) {
array = new int[NodeValue->value->listOfNodes.GetNumElements()];
arrayFromChain = NodeValue->value->listOfNodes.CreateArray();
numElementsInSystem = NodeValue->value->listOfNodes.GetNumElements();
for(int counter = 0; counter < numElementsInSystem; counter++){
array[counter] = *arrayFromChain[counter];
}
SetKE(1,KE);
allocateNewSystem();
system[currentIndex].system->Create_System_LAMMPS(nbody,masstotal,inertia,xcm,xjoint,vcm,omega,ex_space,ey_space,ez_space, numElementsInSystem, array, count);
system[currentIndex].solver = ONSOLVER;
ttk = ttk + 1;
count = ttk;
delete [] array;
delete [] arrayFromChain;
NodeValue= NodeValue->next;
}
}
if(nfree){
MakeDegenerateSystem(nfree,freelist,masstotal,inertia,xcm,vcm,omega,ex_space,ey_space,ez_space);
}
return true;
}
bool Workspace::MakeDegenerateSystem(int& nfree, int*freelist, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space){
allocateNewSystem();
system[currentIndex].system->Create_DegenerateSystem(nfree,freelist,masstotal,inertia,xcm,vcm,omega,ex_space,ey_space,ez_space);
system[currentIndex].solver = ONSOLVER;
return true;
}
bool Workspace::SaveFile(char* filename, int index){
if(index < 0){
index = currentIndex;
}
ofstream file;
file.open(filename, ofstream::out);
if( !file.is_open() ){
cerr << "File '" << filename << "' could not be opened." << endl;
return false;
}
if(index >= 0 && index <= currentIndex){
system[index].system->WriteOut(file);
}
else {
cerr << "Error, requested system index " << index << ", minimum index 0 and maximum index " << currentIndex << endl;
}
file.close();
return true;
}
System* Workspace::GetSystem(int index){
if(index <= currentIndex){
if(index >= 0){
return system[index].system;
}
else{
return system[currentIndex].system;
}
}
else{
return NULL;
}
}
void Workspace::AddSolver(Solver* s, int index){
if(currentIndex >= index){
if(index >= 0){
system[index].solver = (int)s->GetSolverType();
}
else{
system[currentIndex].solver = (int)s->GetSolverType();
}
}
else{
cout << "Error adding solver to index " << index << endl;
}
}
int Workspace::getNumberOfSystems(){
return currentIndex + 1;
}
void Workspace::SetKE(int temp, double SysKE){
KE_val = SysKE;
FirstTime =temp;
}
void Workspace::LobattoOne(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space){
int numsys = currentIndex;
int numbodies;
double time = 0.0;
int * mappings;
double SysKE=0.0;
for (int i = 0; i <= numsys; i++){
mappings = system[i].system->GetMappings();
numbodies = system[i].system->GetNumBodies() - 1;
Matrix FF(6,numbodies);
FF.Zeros();
for (int j=1; j<=numbodies; j++){
FF(1,j) = torque[mappings[j - 1]-1][0]*ConFac;
FF(2,j) = torque[mappings[j - 1]-1][1]*ConFac;
FF(3,j) = torque[mappings[j - 1]-1][2]*ConFac;
FF(4,j) = fcm[mappings[j - 1]-1][0]*ConFac;
FF(5,j) = fcm[mappings[j - 1]-1][1]*ConFac;
FF(6,j) = fcm[mappings[j - 1]-1][2]*ConFac;
}
//------------------------------------//
// Get a solver and solve that system.
Solver * theSolver = Solver::GetSolver((SolverType)system[i].solver);
theSolver->SetSystem(system[i].system);
theSolver->Solve(time, FF);
theSolver->Solve(time, FF);
ColMatrix tempx = *((*theSolver).GetState());
ColMatrix tempv = *((*theSolver).GetStateDerivative());
ColMatrix tempa = *((*theSolver).GetStateDerivativeDerivative());
for(int numint =0 ; numint<3; numint++){
theSolver->Solve(time, FF);
tempa = *((*theSolver).GetStateDerivativeDerivative());
*((*theSolver).GetStateDerivative())= tempv + Thalf*tempa;
}
ColMatrix TempV= *((*theSolver).GetStateDerivative());
*((*theSolver).GetState())= tempx + Tfull*TempV;
int numjoints = system[i].system->joints.GetNumElements();
for(int k = 0; k < numjoints; k++){
system[i].system->joints(k)->ForwardKinematics();
}
for(int k = 0; k < numbodies; k++){
Vect3 temp1 =system[i].system->bodies(k+1)->r;
Vect3 temp2 =system[i].system->bodies(k+1)->v;
Vect3 temp3 =system[i].system->bodies(k+1)->omega;
Mat3x3 temp4 =system[i].system->bodies(k+1)->n_C_k;
for(int m = 0; m < 3; m++){
xcm[mappings[k]-1][m] = temp1(m+1);
vcm[mappings[k]-1][m] = temp2(m+1);
omega[mappings[k]-1][m] = temp3(m+1);
ex_space[mappings[k]-1][m] = temp4(m+1,1);
ey_space[mappings[k]-1][m] = temp4(m+1,2);
ez_space[mappings[k]-1][m] = temp4(m+1,3);
}
SysKE = SysKE + system[i].system->bodies(k+1)->KE;
}
delete theSolver;
}
}
void Workspace::LobattoTwo(double **&vcm,double **&omega,double **&torque, double **&fcm){
int numsys = currentIndex;
int numbodies;
double time = 0.0;
int * mappings;
double SysKE =0.0;
for (int i = 0; i <= numsys; i++){
mappings = system[i].system->GetMappings();
numbodies = system[i].system->GetNumBodies() - 1;
Matrix FF(6,numbodies);
for (int j=1; j<=numbodies; j++){
FF(1,j) = torque[mappings[j - 1]-1][0]*ConFac;
FF(2,j) = torque[mappings[j - 1]-1][1]*ConFac;
FF(3,j) = torque[mappings[j - 1]-1][2]*ConFac;
FF(4,j) = fcm[mappings[j - 1]-1][0]*ConFac;
FF(5,j) = fcm[mappings[j - 1]-1][1]*ConFac;
FF(6,j) = fcm[mappings[j - 1]-1][2]*ConFac;
}
//------------------------------------//
// Get a solver and solve that system.
Solver * theSolver = Solver::GetSolver((SolverType)system[i].solver);
theSolver->SetSystem(system[i].system);
theSolver->Solve(time, FF);
ColMatrix tempv = *((*theSolver).GetStateDerivative());
ColMatrix tempa = *((*theSolver).GetStateDerivativeDerivative());
*((*theSolver).GetStateDerivative()) = tempv + Thalf*tempa;
int numjoints = system[i].system->joints.GetNumElements();
for(int k = 0; k < numjoints; k++){
system[i].system->joints(k)->ForwardKinematics();
}
for(int k = 0; k < numbodies; k++){
Vect3 temp1 =system[i].system->bodies(k+1)->r;
Vect3 temp2 =system[i].system->bodies(k+1)->v;
Vect3 temp3 =system[i].system->bodies(k+1)->omega;
SysKE = SysKE + system[i].system->bodies(k+1)->KE;
for(int m = 0; m < 3; m++){
vcm[mappings[k]-1][m] = temp2(m+1);
omega[mappings[k]-1][m] = temp3(m+1);
}
}
delete theSolver;
}
}
void Workspace::RKStep(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space){
double a[6];
double b[6][6];
double c[6];
//double e[6];
a[1] = 0.2;
a[2] = 0.3;
a[3] = 0.6;
a[4] = 1.0;
a[5] = 0.875;
b[1][0] = 0.2;
b[2][0] = 0.075; b[2][1] = 0.225;
b[3][0] = 0.3; b[3][1] = -0.9; b[3][2] = 1.2;
b[4][0] = -11.0/54.0; b[4][1] = 2.5; b[4][2] = -70.0/27.0; b[4][3] = 35.0/27.0;
b[5][0] = 1631.0/55296.0; b[5][1] = 175.0/512.0; b[5][2] = 575.0/13824.0; b[5][3] = 44275.0/110592.0; b[5][4] = 253.0/4096.0;
c[0] = 37.0/378.0;
c[1] = 0.0;
c[2] = 250.0/621.0;
c[3] = 125.0/594.0;
c[4] = 0.0;
c[5] = 512.0/1771.0;
int numsys = currentIndex;
int numbodies;
double time = 0.0;
int * mappings;
double SysKE =0.0;
for (int i = 0; i <= numsys; i++){
mappings = system[i].system->GetMappings();
numbodies = system[i].system->GetNumBodies() - 1;
Matrix FF(6,numbodies);
for (int j=1; j<=numbodies; j++){
FF(1,j) = ConFac*torque[mappings[j - 1]-1][0];
FF(2,j) = ConFac*torque[mappings[j - 1]-1][1];
FF(3,j) = ConFac*torque[mappings[j - 1]-1][2];
FF(4,j) = ConFac*fcm[mappings[j - 1]-1][0];
FF(5,j) = ConFac*fcm[mappings[j - 1]-1][1];
FF(6,j) = ConFac*fcm[mappings[j - 1]-1][2];
}
//------------------------------------//
// Get a solver and solve that system.
Solver * theSolver = Solver::GetSolver((SolverType)system[i].solver);
theSolver->SetSystem(system[i].system);
theSolver->Solve(time, FF);
ColMatrix initial_x;
ColMatrix initial_xdot;
ColMatMap* x;
ColMatMap* xdot;
ColMatMap* xdotdot;
x = theSolver->GetState();
xdot = theSolver->GetStateDerivative();
xdotdot=theSolver->GetStateDerivativeDerivative();
initial_x = *x;
initial_xdot = *xdot;
ColMatrix f[6];
ColMatrix ff[6];
ff[0] = initial_xdot;
f[0] = *xdotdot;
for(int ii=1;ii<6;ii++){
time = a[ii] * Tfull;
(*x) = initial_x;
(*xdot) = initial_xdot;
for(int j=0;j<ii;j++){
(*x) = (*x) + (b[ii][j]*Tfull)*ff[j];
(*xdot) = (*xdot) + (b[ii][j]*Tfull)*f[j];
}
theSolver->Solve(time,FF);
f[ii] = (*xdotdot);
ff[ii] = (*xdot);
}
(*x) = initial_x + (c[0]*Tfull)*ff[0] + (c[2]*Tfull)*ff[2] + (c[3]*Tfull)*ff[3] + (c[5]*Tfull)*ff[5];
(*xdot) = initial_xdot + (c[0]*Tfull)*f[0] + (c[2]*Tfull)*f[2] + (c[3]*Tfull)*f[3] + (c[5]*Tfull)*f[5];
int numjoints = system[i].system->joints.GetNumElements();
for(int k = 0; k < numjoints; k++){
system[i].system->joints(k)->ForwardKinematics();
}
for(int k = 0; k < numbodies; k++){
Vect3 temp1 =system[i].system->bodies(k+1)->r;
Vect3 temp2 =system[i].system->bodies(k+1)->v;
Vect3 temp3 =system[i].system->bodies(k+1)->omega;
Mat3x3 temp4 =system[i].system->bodies(k+1)->n_C_k;
SysKE = SysKE + system[i].system->bodies(k+1)->KE;
for(int m = 0; m < 3; m++){
xcm[mappings[k]-1][m] = temp1(m+1);
vcm[mappings[k]-1][m] = temp2(m+1);
omega[mappings[k]-1][m] = temp3(m+1);
ex_space[mappings[k]-1][m] = temp4(m+1,1);
ey_space[mappings[k]-1][m] = temp4(m+1,2);
ez_space[mappings[k]-1][m] = temp4(m+1,3);
}
}
delete theSolver;
}
}
void Workspace::WriteFile(char * filename){
int numsys = currentIndex;
int numbodies;
for (int i = 0; i <= numsys; i++){
numbodies = system[i].system->GetNumBodies() - 1;
ofstream outfile;
outfile.open(filename,ofstream::out | ios::app);
outfile << numbodies<<endl;
outfile << "Atoms "<<endl;
for(int k = 0; k < numbodies; k++){
Vect3 temp1 =system[i].system->bodies(k+1)->r;
outfile<<1<<"\t"<<temp1(1)<<"\t"<<temp1(2)<<"\t"<<temp1(3)<<endl;
}
outfile.close();
}
}
c4science · Help