src/io/DumpReader.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#define _LARGEFILE_SOURCE64
#define _FILE_OFFSET_BITS 64

#include <sys/types.h>
#include <sys/stat.h>

#include <iostream>
#include <math.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/simError.h"
#include "utils/MemoryUtils.hpp"
#include "utils/StringTokenizer.hpp"

#ifdef IS_MPI

#include <mpi.h>
#define TAKE_THIS_TAG_CHAR 0
#define TAKE_THIS_TAG_INT 1

#endif // is_mpi


namespace oopse {
  
  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
    : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
    
#ifdef IS_MPI
    
      if (worldRank == 0) {
#endif
      
      inFile_ = new std::ifstream(filename_.c_str());
      
        if (inFile_->fail()) {
          sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
          painCave.isFatal = 1;
          simError();
        }
      
#ifdef IS_MPI
      
      }
    
      strcpy(checkPointMsg, "Dump file opened for reading successfully.");
      MPIcheckPoint();
    
#endif
    
      return;
    }
  
  DumpReader::~DumpReader() {
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif
      
      delete inFile_;
      
#ifdef IS_MPI
      
    }
    
    strcpy(checkPointMsg, "Dump file closed successfully.");
    MPIcheckPoint();
    
#endif
    
    return;
  }
  
  int DumpReader::getNFrames(void) {
    
    if (!isScanned_)
      scanFile();
    
    return nframes_;
  }
  
  void DumpReader::scanFile(void) {
    int i, j;
    int lineNum = 0;
    char readBuffer[maxBufferSize];
    std::streampos  currPos;
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      
      inFile_->seekg (0, std::ios::beg);
      

      currPos = inFile_->tellg();
      inFile_->getline(readBuffer, sizeof(readBuffer));
      lineNum++;
      
      if (inFile_->eof()) {
        sprintf(painCave.errMsg,
                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                filename_.c_str(),
                lineNum);
        painCave.isFatal = 1;
        simError();
      }
      
      while (!inFile_->eof()) {
        framePos_.push_back(currPos);
        
        i = atoi(readBuffer);
        
        inFile_->getline(readBuffer, sizeof(readBuffer));
        lineNum++;
        
        if (inFile_->eof()) {
          sprintf(painCave.errMsg,
                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                  filename_.c_str(),
                  lineNum);
          painCave.isFatal = 1;
          simError();
        }
        
        for(j = 0; j < i; j++) {
          inFile_->getline(readBuffer, sizeof(readBuffer));
          lineNum++;
          
          if (inFile_->eof()) {
            sprintf(painCave.errMsg,
                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
                    " with atom %d\n", filename_.c_str(),
                    lineNum,
                    j);
            
            painCave.isFatal = 1;
            simError();
          }
        }
        
        currPos = inFile_->tellg();
        inFile_->getline(readBuffer, sizeof(readBuffer));
        lineNum++;
      }

      inFile_->seekg (0, std::ios::beg);
      
      nframes_ = framePos_.size();
#ifdef IS_MPI
    }
    
    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
    
    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
    MPIcheckPoint();
    
#endif // is_mpi
    
    isScanned_ = true;
  }
  
  void DumpReader::readFrame(int whichFrame) {
    if (!isScanned_)
      scanFile();
        
    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
    
    if (storageLayout & DataStorage::dslPosition) {
      needPos_ = true;
    } else {
      needPos_ = false;
    }
    
    if (storageLayout & DataStorage::dslVelocity) {
      needVel_ = true;
    } else {
      needVel_ = false;
    }
    
    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
      needQuaternion_ = true;
    } else {
      needQuaternion_ = false;
    }
    
    if (storageLayout & DataStorage::dslAngularMomentum) {
      needAngMom_ = true;
    } else {
      needAngMom_ = false;    
    }
    
    readSet(whichFrame);
  }
  
  void DumpReader::readSet(int whichFrame) {
    int i;
    int nTotObjs;                  // the number of atoms
    char read_buffer[maxBufferSize];  //the line buffer for reading
    char * eof_test;               // ptr to see when we reach the end of the file
    
    Molecule* mol;
    StuntDouble* integrableObject;
    SimInfo::MoleculeIterator mi;
    Molecule::IntegrableObjectIterator ii;
    
#ifndef IS_MPI
    inFile_->clear(); 
    inFile_->seekg(framePos_[whichFrame]);
        
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
      sprintf(painCave.errMsg,
              "DumpReader error: error reading 1st line of \"%s\"\n",
              filename_.c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    nTotObjs = atoi(read_buffer);
    
    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
      sprintf(painCave.errMsg,
              "DumpReader error. %s nIntegrable, %d, "
              "does not match the meta-data file's nIntegrable, %d.\n",
              filename_.c_str(),
              nTotObjs,
              info_->getNGlobalIntegrableObjects());
      
      painCave.isFatal = 1;
      simError();
    }
    
    //read the box mat from the comment line
    
    
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
              filename_.c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
    
    //parse dump lines
    
    i = 0;
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
      
      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
           integrableObject = mol->nextIntegrableObject(ii)) {           
        
        
        if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
          sprintf(painCave.errMsg,
                  "DumpReader Error: error in reading file %s\n"
                  "natoms  = %d; index = %d\n"
                  "error reading the line from the file.\n",
                  filename_.c_str(),
                  nTotObjs,
                  i);
          
          painCave.isFatal = 1;
          simError();
        }
        
        parseDumpLine(read_buffer, integrableObject);
        i++;
      }
    }
    
    // MPI Section of code..........
    
#else //IS_MPI
    
    // first thing first, suspend fatalities.
    int masterNode = 0;
    int nCurObj;
    painCave.isEventLoop = 1;
    
    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
    int haveError;
    
    MPI_Status istatus;
    int nitems;
    
    nTotObjs = info_->getNGlobalIntegrableObjects();
    haveError = 0;
    
    if (worldRank == masterNode) {
      inFile_->clear();            
      inFile_->seekg(framePos_[whichFrame]);
      
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
                filename_.c_str());
        painCave.isFatal = 1;
        simError();
      }
      
      nitems = atoi(read_buffer);
      
      // Check to see that the number of integrable objects in the
      // intial configuration file is the same as derived from the
      // meta-data file.
      
      if (nTotObjs != nitems) {
        sprintf(painCave.errMsg,
                "DumpReader Error. %s nIntegrable, %d, "
                "does not match the meta-data file's nIntegrable, %d.\n",
                filename_.c_str(),
                nTotObjs,
                info_->getNGlobalIntegrableObjects());
        
        painCave.isFatal = 1;
        simError();
      }
      
      //read the boxMat from the comment line
      
      
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
                filename_.c_str());
        painCave.isFatal = 1;
        simError();
      }
      
      //Every single processor will parse the comment line by itself
      //By using this way, we might lose some efficiency, but if we want to add
      //more parameters into comment line, we only need to modify function
      //parseCommentLine
      
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
      
      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
        int which_node = info_->getMolToProc(i);
        
        if (which_node == masterNode) {
          //molecules belong to master node
          
          mol = info_->getMoleculeByGlobalIndex(i);
          
          if (mol == NULL) {
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
            painCave.isFatal = 1;
            simError();
          }
          
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
               integrableObject = mol->nextIntegrableObject(ii)){
            
            
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
              sprintf(painCave.errMsg,
                      "DumpReader Error: error in reading file %s\n"
                      "natoms  = %d; index = %d\n"
                      "error reading the line from the file.\n",
                      filename_.c_str(),
                      nTotObjs,
                      i);
              
              painCave.isFatal = 1;
              simError();
            }
            
            parseDumpLine(read_buffer, integrableObject);
          }
        } else {
          //molecule belongs to slave nodes
          
          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
                   MPI_COMM_WORLD, &istatus);
          
          for(int j = 0; j < nCurObj; j++) {
            
            
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
              sprintf(painCave.errMsg,
                      "DumpReader Error: error in reading file %s\n"
                      "natoms  = %d; index = %d\n"
                      "error reading the line from the file.\n",
                      filename_.c_str(),
                      nTotObjs,
                      i);
              
              painCave.isFatal = 1;
              simError();
            }
            
            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
          }
        }
      }
    } else {
      //actions taken at slave nodes
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
      
      /**@todo*/
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
      
      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
        int which_node = info_->getMolToProc(i);
        
        if (which_node == worldRank) {
          //molecule with global index i belongs to this processor
          
          mol = info_->getMoleculeByGlobalIndex(i);
          if (mol == NULL) {
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
            painCave.isFatal = 1;
            simError();
          }
          
          nCurObj = mol->getNIntegrableObjects();
          
          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
                   MPI_COMM_WORLD);
          
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
               integrableObject = mol->nextIntegrableObject(ii)){
            
            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
            
            parseDumpLine(read_buffer, integrableObject);
          }
          
        }
        
      }
      
    }
    
#endif
    
  }
  
  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
    
    Vector3d pos;  // position place holders
    Vector3d vel;  // velocity placeholders
    Quat4d q;    // the quaternions
    Vector3d ji;   // angular velocity placeholders;
    StringTokenizer tokenizer(line);
    int nTokens;
    
    nTokens = tokenizer.countTokens();
    
    if (nTokens < 14) {
      sprintf(painCave.errMsg,
              "DumpReader Error: Not enough Tokens.\n%s\n", line);
      painCave.isFatal = 1;
      simError();
    }
    
    std::string name = tokenizer.nextToken();
    
    if (name != integrableObject->getType()) {
      
      sprintf(painCave.errMsg,
              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
      painCave.isFatal = 1;
      simError();        
    }
    
    pos[0] = tokenizer.nextTokenAsDouble();
    pos[1] = tokenizer.nextTokenAsDouble();
    pos[2] = tokenizer.nextTokenAsDouble();
    if (needPos_) {
      integrableObject->setPos(pos);
    }
    
    vel[0] = tokenizer.nextTokenAsDouble();
    vel[1] = tokenizer.nextTokenAsDouble();
    vel[2] = tokenizer.nextTokenAsDouble();
    if (needVel_) {
      integrableObject->setVel(vel);
    }
    
    if (integrableObject->isDirectional()) {
      
      q[0] = tokenizer.nextTokenAsDouble();
      q[1] = tokenizer.nextTokenAsDouble();
      q[2] = tokenizer.nextTokenAsDouble();
      q[3] = tokenizer.nextTokenAsDouble();
      
      double qlen = q.length();
      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
        
        sprintf(painCave.errMsg,
                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0\n");
        painCave.isFatal = 1;
        simError();
        
      } 
      
      q.normalize();
      if (needQuaternion_) {           
        integrableObject->setQ(q);
      }
      
      ji[0] = tokenizer.nextTokenAsDouble();
      ji[1] = tokenizer.nextTokenAsDouble();
      ji[2] = tokenizer.nextTokenAsDouble();
      if (needAngMom_) {
        integrableObject->setJ(ji);
      }
    }
    
  }
  
  
  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
    double currTime;
    Mat3x3d hmat;
    double chi;
    double integralOfChiDt;
    Mat3x3d eta;
    
    StringTokenizer tokenizer(line);
    int nTokens;
    
    nTokens = tokenizer.countTokens();
    
    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
    if (nTokens < 10) {
      sprintf(painCave.errMsg,
              "DumpReader Error: Not enough tokens in comment line: %s", line);
      painCave.isFatal = 1;
      simError();   
    }
    
    //read current time
    currTime = tokenizer.nextTokenAsDouble();
    s->setTime(currTime);
    
    //read h-matrix
    hmat(0, 0) = tokenizer.nextTokenAsDouble();
    hmat(0, 1) = tokenizer.nextTokenAsDouble();
    hmat(0, 2) = tokenizer.nextTokenAsDouble();
    hmat(1, 0) = tokenizer.nextTokenAsDouble();
    hmat(1, 1) = tokenizer.nextTokenAsDouble();
    hmat(1, 2) = tokenizer.nextTokenAsDouble();
    hmat(2, 0) = tokenizer.nextTokenAsDouble();
    hmat(2, 1) = tokenizer.nextTokenAsDouble();
    hmat(2, 2) = tokenizer.nextTokenAsDouble();
    s->setHmat(hmat);
    
    //read chi and integralOfChidt, they should apprear in pair
    if (tokenizer.countTokens() >= 2) {
      chi = tokenizer.nextTokenAsDouble();
      integralOfChiDt = tokenizer.nextTokenAsDouble();            
      
      s->setChi(chi);
      s->setIntegralOfChiDt(integralOfChiDt);
    }
    
    //read eta (eta is 3x3 matrix)
    if (tokenizer.countTokens() >= 9) {
      eta(0, 0) = tokenizer.nextTokenAsDouble();
      eta(0, 1) = tokenizer.nextTokenAsDouble();
      eta(0, 2) = tokenizer.nextTokenAsDouble();
      eta(1, 0) = tokenizer.nextTokenAsDouble();
      eta(1, 1) = tokenizer.nextTokenAsDouble();
      eta(1, 2) = tokenizer.nextTokenAsDouble();
      eta(2, 0) = tokenizer.nextTokenAsDouble();
      eta(2, 1) = tokenizer.nextTokenAsDouble();
      eta(2, 2) = tokenizer.nextTokenAsDouble();      
      
      s->setEta(eta);
    }
    
    
  }
  
}//end namespace oopse
Revision:	2334
Committed:	Wed Sep 28 16:32:30 2005 UTC (18 years, 9 months ago) by gezelter
File size:	19397 byte(s)
Log Message:	tracking down quaternion bug
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42	#define _LARGEFILE_SOURCE64
43	#define _FILE_OFFSET_BITS 64
44
45	#include <sys/types.h>
46	#include <sys/stat.h>
47
48	#include <iostream>
49	#include <math.h>
50
51	#include <stdio.h>
52	#include <stdlib.h>
53	#include <string.h>
54
55	#include "io/DumpReader.hpp"
56	#include "primitives/Molecule.hpp"
57	#include "utils/simError.h"
58	#include "utils/MemoryUtils.hpp"
59	#include "utils/StringTokenizer.hpp"
60
61	#ifdef IS_MPI
62
63	#include <mpi.h>
64	#define TAKE_THIS_TAG_CHAR 0
65	#define TAKE_THIS_TAG_INT 1
66
67	#endif // is_mpi
68
69
70	namespace oopse {
71
72	DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73	: info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74
75	#ifdef IS_MPI
76
77	if (worldRank == 0) {
78	#endif
79
80	inFile_ = new std::ifstream(filename_.c_str());
81
82	if (inFile_->fail()) {
83	sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84	painCave.isFatal = 1;
85	simError();
86	}
87
88	#ifdef IS_MPI
89
90	}
91
92	strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93	MPIcheckPoint();
94
95	#endif
96
97	return;
98	}
99
100	DumpReader::~DumpReader() {
101
102	#ifdef IS_MPI
103
104	if (worldRank == 0) {
105	#endif
106
107	delete inFile_;
108
109	#ifdef IS_MPI
110
111	}
112
113	strcpy(checkPointMsg, "Dump file closed successfully.");
114	MPIcheckPoint();
115
116	#endif
117
118	return;
119	}
120
121	int DumpReader::getNFrames(void) {
122
123	if (!isScanned_)
124	scanFile();
125
126	return nframes_;
127	}
128
129	void DumpReader::scanFile(void) {
130	int i, j;
131	int lineNum = 0;
132	char readBuffer[maxBufferSize];
133	std::streampos currPos;
134
135	#ifdef IS_MPI
136
137	if (worldRank == 0) {
138	#endif // is_mpi
139
140	inFile_->seekg (0, std::ios::beg);
141
142
143	currPos = inFile_->tellg();
144	inFile_->getline(readBuffer, sizeof(readBuffer));
145	lineNum++;
146
147	if (inFile_->eof()) {
148	sprintf(painCave.errMsg,
149	"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
150	filename_.c_str(),
151	lineNum);
152	painCave.isFatal = 1;
153	simError();
154	}
155
156	while (!inFile_->eof()) {
157	framePos_.push_back(currPos);
158
159	i = atoi(readBuffer);
160
161	inFile_->getline(readBuffer, sizeof(readBuffer));
162	lineNum++;
163
164	if (inFile_->eof()) {
165	sprintf(painCave.errMsg,
166	"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
167	filename_.c_str(),
168	lineNum);
169	painCave.isFatal = 1;
170	simError();
171	}
172
173	for(j = 0; j < i; j++) {
174	inFile_->getline(readBuffer, sizeof(readBuffer));
175	lineNum++;
176
177	if (inFile_->eof()) {
178	sprintf(painCave.errMsg,
179	"DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
180	" with atom %d\n", filename_.c_str(),
181	lineNum,
182	j);
183
184	painCave.isFatal = 1;
185	simError();
186	}
187	}
188
189	currPos = inFile_->tellg();
190	inFile_->getline(readBuffer, sizeof(readBuffer));
191	lineNum++;
192	}
193
194	inFile_->seekg (0, std::ios::beg);
195
196	nframes_ = framePos_.size();
197	#ifdef IS_MPI
198	}
199
200	MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
201
202	strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
203	MPIcheckPoint();
204
205	#endif // is_mpi
206
207	isScanned_ = true;
208	}
209
210	void DumpReader::readFrame(int whichFrame) {
211	if (!isScanned_)
212	scanFile();
213
214	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
215
216	if (storageLayout & DataStorage::dslPosition) {
217	needPos_ = true;
218	} else {
219	needPos_ = false;
220	}
221
222	if (storageLayout & DataStorage::dslVelocity) {
223	needVel_ = true;
224	} else {
225	needVel_ = false;
226	}
227
228	if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
229	needQuaternion_ = true;
230	} else {
231	needQuaternion_ = false;
232	}
233
234	if (storageLayout & DataStorage::dslAngularMomentum) {
235	needAngMom_ = true;
236	} else {
237	needAngMom_ = false;
238	}
239
240	readSet(whichFrame);
241	}
242
243	void DumpReader::readSet(int whichFrame) {
244	int i;
245	int nTotObjs; // the number of atoms
246	char read_buffer[maxBufferSize]; //the line buffer for reading
247	char * eof_test; // ptr to see when we reach the end of the file
248
249	Molecule* mol;
250	StuntDouble* integrableObject;
251	SimInfo::MoleculeIterator mi;
252	Molecule::IntegrableObjectIterator ii;
253
254	#ifndef IS_MPI
255	inFile_->clear();
256	inFile_->seekg(framePos_[whichFrame]);
257
258	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
259	sprintf(painCave.errMsg,
260	"DumpReader error: error reading 1st line of \"%s\"\n",
261	filename_.c_str());
262	painCave.isFatal = 1;
263	simError();
264	}
265
266	nTotObjs = atoi(read_buffer);
267
268	if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
269	sprintf(painCave.errMsg,
270	"DumpReader error. %s nIntegrable, %d, "
271	"does not match the meta-data file's nIntegrable, %d.\n",
272	filename_.c_str(),
273	nTotObjs,
274	info_->getNGlobalIntegrableObjects());
275
276	painCave.isFatal = 1;
277	simError();
278	}
279
280	//read the box mat from the comment line
281
282
283	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
284	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
285	filename_.c_str());
286	painCave.isFatal = 1;
287	simError();
288	}
289
290	parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
291
292	//parse dump lines
293
294	i = 0;
295	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
296
297	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
298	integrableObject = mol->nextIntegrableObject(ii)) {
299
300
301
302	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
303	sprintf(painCave.errMsg,
304	"DumpReader Error: error in reading file %s\n"
305	"natoms = %d; index = %d\n"
306	"error reading the line from the file.\n",
307	filename_.c_str(),
308	nTotObjs,
309	i);
310
311	painCave.isFatal = 1;
312	simError();
313	}
314
315	parseDumpLine(read_buffer, integrableObject);
316	i++;
317	}
318	}
319
320	// MPI Section of code..........
321
322	#else //IS_MPI
323
324	// first thing first, suspend fatalities.
325	int masterNode = 0;
326	int nCurObj;
327	painCave.isEventLoop = 1;
328
329	int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
330	int haveError;
331
332	MPI_Status istatus;
333	int nitems;
334
335	nTotObjs = info_->getNGlobalIntegrableObjects();
336	haveError = 0;
337
338	if (worldRank == masterNode) {
339	inFile_->clear();
340	inFile_->seekg(framePos_[whichFrame]);
341
342	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
343	sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
344	filename_.c_str());
345	painCave.isFatal = 1;
346	simError();
347	}
348
349	nitems = atoi(read_buffer);
350
351	// Check to see that the number of integrable objects in the
352	// intial configuration file is the same as derived from the
353	// meta-data file.
354
355	if (nTotObjs != nitems) {
356	sprintf(painCave.errMsg,
357	"DumpReader Error. %s nIntegrable, %d, "
358	"does not match the meta-data file's nIntegrable, %d.\n",
359	filename_.c_str(),
360	nTotObjs,
361	info_->getNGlobalIntegrableObjects());
362
363	painCave.isFatal = 1;
364	simError();
365	}
366
367	//read the boxMat from the comment line
368
369
370
371	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
372	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
373	filename_.c_str());
374	painCave.isFatal = 1;
375	simError();
376	}
377
378	//Every single processor will parse the comment line by itself
379	//By using this way, we might lose some efficiency, but if we want to add
380	//more parameters into comment line, we only need to modify function
381	//parseCommentLine
382
383	MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
384	parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
385
386	for(i = 0; i < info_->getNGlobalMolecules(); i++) {
387	int which_node = info_->getMolToProc(i);
388
389	if (which_node == masterNode) {
390	//molecules belong to master node
391
392	mol = info_->getMoleculeByGlobalIndex(i);
393
394	if (mol == NULL) {
395	sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
396	painCave.isFatal = 1;
397	simError();
398	}
399
400	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
401	integrableObject = mol->nextIntegrableObject(ii)){
402
403
404
405	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
406	sprintf(painCave.errMsg,
407	"DumpReader Error: error in reading file %s\n"
408	"natoms = %d; index = %d\n"
409	"error reading the line from the file.\n",
410	filename_.c_str(),
411	nTotObjs,
412	i);
413
414	painCave.isFatal = 1;
415	simError();
416	}
417
418	parseDumpLine(read_buffer, integrableObject);
419	}
420	} else {
421	//molecule belongs to slave nodes
422
423	MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
424	MPI_COMM_WORLD, &istatus);
425
426	for(int j = 0; j < nCurObj; j++) {
427
428
429	if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
430	sprintf(painCave.errMsg,
431	"DumpReader Error: error in reading file %s\n"
432	"natoms = %d; index = %d\n"
433	"error reading the line from the file.\n",
434	filename_.c_str(),
435	nTotObjs,
436	i);
437
438	painCave.isFatal = 1;
439	simError();
440	}
441
442	MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
443	TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
444	}
445	}
446	}
447	} else {
448	//actions taken at slave nodes
449	MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
450
451	/*@todo/
452	parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
453
454	for(i = 0; i < info_->getNGlobalMolecules(); i++) {
455	int which_node = info_->getMolToProc(i);
456
457	if (which_node == worldRank) {
458	//molecule with global index i belongs to this processor
459
460	mol = info_->getMoleculeByGlobalIndex(i);
461	if (mol == NULL) {
462	sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
463	painCave.isFatal = 1;
464	simError();
465	}
466
467	nCurObj = mol->getNIntegrableObjects();
468
469	MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
470	MPI_COMM_WORLD);
471
472	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
473	integrableObject = mol->nextIntegrableObject(ii)){
474
475	MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
476	TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
477
478	parseDumpLine(read_buffer, integrableObject);
479	}
480
481	}
482
483	}
484
485	}
486
487	#endif
488
489	}
490
491	void DumpReader::parseDumpLine(char line, StuntDouble integrableObject) {
492
493	Vector3d pos; // position place holders
494	Vector3d vel; // velocity placeholders
495	Quat4d q; // the quaternions
496	Vector3d ji; // angular velocity placeholders;
497	StringTokenizer tokenizer(line);
498	int nTokens;
499
500	nTokens = tokenizer.countTokens();
501
502	if (nTokens < 14) {
503	sprintf(painCave.errMsg,
504	"DumpReader Error: Not enough Tokens.\n%s\n", line);
505	painCave.isFatal = 1;
506	simError();
507	}
508
509	std::string name = tokenizer.nextToken();
510
511	if (name != integrableObject->getType()) {
512
513	sprintf(painCave.errMsg,
514	"DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
515	name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
516	painCave.isFatal = 1;
517	simError();
518	}
519
520	pos[0] = tokenizer.nextTokenAsDouble();
521	pos[1] = tokenizer.nextTokenAsDouble();
522	pos[2] = tokenizer.nextTokenAsDouble();
523	if (needPos_) {
524	integrableObject->setPos(pos);
525	}
526
527	vel[0] = tokenizer.nextTokenAsDouble();
528	vel[1] = tokenizer.nextTokenAsDouble();
529	vel[2] = tokenizer.nextTokenAsDouble();
530	if (needVel_) {
531	integrableObject->setVel(vel);
532	}
533
534	if (integrableObject->isDirectional()) {
535
536	q[0] = tokenizer.nextTokenAsDouble();
537	q[1] = tokenizer.nextTokenAsDouble();
538	q[2] = tokenizer.nextTokenAsDouble();
539	q[3] = tokenizer.nextTokenAsDouble();
540
541	double qlen = q.length();
542	if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
543
544	sprintf(painCave.errMsg,
545	"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0\n");
546	painCave.isFatal = 1;
547	simError();
548
549	}
550
551	q.normalize();
552	if (needQuaternion_) {
553	integrableObject->setQ(q);
554	}
555
556	ji[0] = tokenizer.nextTokenAsDouble();
557	ji[1] = tokenizer.nextTokenAsDouble();
558	ji[2] = tokenizer.nextTokenAsDouble();
559	if (needAngMom_) {
560	integrableObject->setJ(ji);
561	}
562	}
563
564	}
565
566
567	void DumpReader::parseCommentLine(char* line, Snapshot* s) {
568	double currTime;
569	Mat3x3d hmat;
570	double chi;
571	double integralOfChiDt;
572	Mat3x3d eta;
573
574	StringTokenizer tokenizer(line);
575	int nTokens;
576
577	nTokens = tokenizer.countTokens();
578
579	//comment line should at least contain 10 tokens: current time(1 token) and h-matrix(9 tokens)
580	if (nTokens < 10) {
581	sprintf(painCave.errMsg,
582	"DumpReader Error: Not enough tokens in comment line: %s", line);
583	painCave.isFatal = 1;
584	simError();
585	}
586
587	//read current time
588	currTime = tokenizer.nextTokenAsDouble();
589	s->setTime(currTime);
590
591	//read h-matrix
592	hmat(0, 0) = tokenizer.nextTokenAsDouble();
593	hmat(0, 1) = tokenizer.nextTokenAsDouble();
594	hmat(0, 2) = tokenizer.nextTokenAsDouble();
595	hmat(1, 0) = tokenizer.nextTokenAsDouble();
596	hmat(1, 1) = tokenizer.nextTokenAsDouble();
597	hmat(1, 2) = tokenizer.nextTokenAsDouble();
598	hmat(2, 0) = tokenizer.nextTokenAsDouble();
599	hmat(2, 1) = tokenizer.nextTokenAsDouble();
600	hmat(2, 2) = tokenizer.nextTokenAsDouble();
601	s->setHmat(hmat);
602
603	//read chi and integralOfChidt, they should apprear in pair
604	if (tokenizer.countTokens() >= 2) {
605	chi = tokenizer.nextTokenAsDouble();
606	integralOfChiDt = tokenizer.nextTokenAsDouble();
607
608	s->setChi(chi);
609	s->setIntegralOfChiDt(integralOfChiDt);
610	}
611
612	//read eta (eta is 3x3 matrix)
613	if (tokenizer.countTokens() >= 9) {
614	eta(0, 0) = tokenizer.nextTokenAsDouble();
615	eta(0, 1) = tokenizer.nextTokenAsDouble();
616	eta(0, 2) = tokenizer.nextTokenAsDouble();
617	eta(1, 0) = tokenizer.nextTokenAsDouble();
618	eta(1, 1) = tokenizer.nextTokenAsDouble();
619	eta(1, 2) = tokenizer.nextTokenAsDouble();
620	eta(2, 0) = tokenizer.nextTokenAsDouble();
621	eta(2, 1) = tokenizer.nextTokenAsDouble();
622	eta(2, 2) = tokenizer.nextTokenAsDouble();
623
624	s->setEta(eta);
625	}
626
627
628	}
629
630	}//end namespace oopse