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(); 
       
      RealType 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) { 
    RealType currTime; 
    Mat3x3d hmat; 
    RealType chi; 
    RealType 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:	963
Committed:	Wed May 17 21:51:42 2006 UTC (18 years, 11 months ago) by tim
Original Path:	*trunk/src/io/DumpReader.cpp*
File size:	19405 byte(s)
Log Message:	Adding single precision capabilities to c++ side
#	User	Rev	Content
1	chrisfen	721	/*
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	tim	963	RealType qlen = q.length();
542	chrisfen	721	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	tim	963	RealType currTime;
569	chrisfen	721	Mat3x3d hmat;
570	tim	963	RealType chi;
571			RealType integralOfChiDt;
572	chrisfen	721	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