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), needCOMprops_(false) { 
    
#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 lineNo = 0; 
    std::streampos prevPos;
    std::streampos  currPos; 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif // is_mpi 
      
      currPos = inFile_->tellg();
      prevPos = currPos;
      bool foundOpenSnapshotTag = false;
      bool foundClosedSnapshotTag = false;
      while(inFile_->getline(buffer, bufferSize)) {
        ++lineNo;
        
        std::string line = buffer;
        currPos = inFile_->tellg(); 
        if (line.find("<Snapshot>")!= std::string::npos) {
          if (foundOpenSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError();           
          }
          foundOpenSnapshotTag = true;
          foundClosedSnapshotTag = false;
          framePos_.push_back(prevPos);
          
        } else if (line.find("</Snapshot>") != std::string::npos){
          if (!foundOpenSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError(); 
          }
          
          if (foundClosedSnapshotTag) {
            sprintf(painCave.errMsg, 
                    "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo, 
                    filename_.c_str()); 
            painCave.isFatal = 1; 
            simError(); 
          }
          foundClosedSnapshotTag = true;
          foundOpenSnapshotTag = false;
        }
        prevPos = currPos;
      }
      
      // only found <Snapshot> for the last frame means the file is corrupted, we should discard
      // it and give a warning message
      if (foundOpenSnapshotTag) {
        sprintf(painCave.errMsg, 
                "DumpReader: last frame in %s is invalid\n", filename_.c_str()); 
        painCave.isFatal = 0; 
        simError();       
        framePos_.pop_back();
      }
      
      nframes_ = framePos_.size(); 
      
      if (nframes_ == 0) {
        sprintf(painCave.errMsg, 
                "DumpReader: %s does not contain a valid frame\n", filename_.c_str()); 
        painCave.isFatal = 1; 
        simError();      
      }
#ifdef IS_MPI 
    } 
     
    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD); 
    
#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); 

    if (needCOMprops_) {
      Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
      Vector3d com;
      Vector3d comvel;
      Vector3d comw;
      info_->getComAll(com, comvel);
      comw = info_->getAngularMomentum();
      s->setCOMprops(com, comvel, comw);      
    }

  } 
   
  void DumpReader::readSet(int whichFrame) {     
    std::string line;

#ifndef IS_MPI 

    inFile_->clear();  
    inFile_->seekg(framePos_[whichFrame]); 

    std::istream& inputStream = *inFile_;     

#else 
    int masterNode = 0;
    std::stringstream sstream;
    if (worldRank == masterNode) {
      std::string sendBuffer;

      inFile_->clear();  
      inFile_->seekg(framePos_[whichFrame]); 
      
      while (inFile_->getline(buffer, bufferSize)) {

        line = buffer;
        sendBuffer += line;
        sendBuffer += '\n';
        if (line.find("</Snapshot>") != std::string::npos) {
          break;
        }        
      }

      int sendBufferSize = sendBuffer.size();
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);     
      MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);     
      
      sstream.str(sendBuffer);
    } else {
      int sendBufferSize;
      MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);     
      char * recvBuffer = new char[sendBufferSize+1];
      MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);     
      sstream.str(recvBuffer);
    }      

    std::istream& inputStream = sstream;  
#endif

    inputStream.getline(buffer, bufferSize);

    line = buffer;
    if (line.find("<Snapshot>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: can not find <Snapshot>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    } 
    
    //read frameData
    readFrameProperties(inputStream);

    //read StuntDoubles
    readStuntDoubles(inputStream);     

    inputStream.getline(buffer, bufferSize);
    line = buffer;
    if (line.find("</Snapshot>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: can not find </Snapshot>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }        
   
  } 
   
  void DumpReader::parseDumpLine(const std::string& line) { 

       
    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 2) { 
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 

    int index = tokenizer.nextTokenAsInt();
 
    StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);

    if (integrableObject == NULL) {
      return;
    }
    std::string type = tokenizer.nextToken(); 
    int size = type.size();

    for(int i = 0; i < size; ++i) {
      switch(type[i]) {
        
        case 'p': {
            Vector3d pos;
            pos[0] = tokenizer.nextTokenAsDouble(); 
            pos[1] = tokenizer.nextTokenAsDouble(); 
            pos[2] = tokenizer.nextTokenAsDouble(); 
            if (needPos_) { 
              integrableObject->setPos(pos); 
            }             
            break;
        }
        case 'v' : {
            Vector3d vel;
            vel[0] = tokenizer.nextTokenAsDouble(); 
            vel[1] = tokenizer.nextTokenAsDouble(); 
            vel[2] = tokenizer.nextTokenAsDouble(); 
            if (needVel_) { 
              integrableObject->setVel(vel); 
            } 
            break;
        }

        case 'q' : {
           Quat4d q;
           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); 
             }               
           }            
           break;
        }  
        case 'j' : {
          Vector3d ji;
          if (integrableObject->isDirectional()) {
             ji[0] = tokenizer.nextTokenAsDouble(); 
             ji[1] = tokenizer.nextTokenAsDouble(); 
             ji[2] = tokenizer.nextTokenAsDouble(); 
             if (needAngMom_) { 
               integrableObject->setJ(ji); 
             } 
          }
          break;
        }  
        case 'f': {

          Vector3d force;
          force[0] = tokenizer.nextTokenAsDouble(); 
          force[1] = tokenizer.nextTokenAsDouble(); 
          force[2] = tokenizer.nextTokenAsDouble();           
          integrableObject->setFrc(force); 
          break;
        }
        case 't' : {

           Vector3d torque;
           torque[0] = tokenizer.nextTokenAsDouble(); 
           torque[1] = tokenizer.nextTokenAsDouble(); 
           torque[2] = tokenizer.nextTokenAsDouble();           
           integrableObject->setTrq(torque);          
           break;
        }
        default: {
               sprintf(painCave.errMsg, 
                       "DumpReader Error: %s is an unrecognized type\n", type.c_str()); 
               painCave.isFatal = 1; 
               simError(); 
          break;   
        }

      }
    }
     
  } 
   

  void  DumpReader::readStuntDoubles(std::istream& inputStream) {

    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);
    
    if (line.find("<StuntDoubles>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: Missing <StuntDoubles>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</StuntDoubles>") != std::string::npos) {
        break;
      }

      parseDumpLine(line);
    }
  
  }

  void DumpReader::readFrameProperties(std::istream& inputStream) {

    Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
    inputStream.getline(buffer, bufferSize);
    std::string line(buffer);

    if (line.find("<FrameData>") == std::string::npos) {
      sprintf(painCave.errMsg, 
              "DumpReader Error: Missing <FrameData>\n"); 
      painCave.isFatal = 1; 
      simError(); 
    }

    while(inputStream.getline(buffer, bufferSize)) {
      line = buffer;
      
      if(line.find("</FrameData>") != std::string::npos) {
        break;
      }
      
      StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
      if (!tokenizer.hasMoreTokens()) {
        sprintf(painCave.errMsg, 
                "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str()); 
        painCave.isFatal = 1; 
        simError();      
      }

      std::string propertyName = tokenizer.nextToken();
      if (propertyName == "Time") {
        RealType currTime = tokenizer.nextTokenAsDouble(); 
        s->setTime(currTime); 
      } else if (propertyName == "Hmat"){
        Mat3x3d hmat;
        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);      
      } else if (propertyName == "Thermostat") {
        RealType chi = tokenizer.nextTokenAsDouble();
        RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
        s->setChi(chi); 
        s->setIntegralOfChiDt(integralOfChiDt);        
     } else if (propertyName == "Barostat") {
        Mat3x3d eta;
        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); 
      } else {
        sprintf(painCave.errMsg, 
                "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str()); 
        painCave.isFatal = 0; 
        simError();         
      }
      
    }

  }

   
}//end namespace oopse 
Revision:	1104
Committed:	Fri Dec 29 21:43:18 2006 UTC (18 years, 4 months ago) by gezelter
Original Path:	*trunk/src/io/DumpReader.cpp*
File size:	16121 byte(s)
Log Message:	Adding the ability to compute Center of Mass properties on reading in a dump file. DumpReader must be told if it needs to do this, however.
#	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), needCOMprops_(false) {
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,
84	"DumpReader: Cannot open file: %s\n",
85	filename_.c_str());
86	painCave.isFatal = 1;
87	simError();
88	}
89
90	#ifdef IS_MPI
91
92	}
93
94	strcpy(checkPointMsg, "Dump file opened for reading successfully.");
95	MPIcheckPoint();
96
97	#endif
98
99	return;
100	}
101
102	DumpReader::~DumpReader() {
103
104	#ifdef IS_MPI
105
106	if (worldRank == 0) {
107	#endif
108
109	delete inFile_;
110
111	#ifdef IS_MPI
112
113	}
114
115	strcpy(checkPointMsg, "Dump file closed successfully.");
116	MPIcheckPoint();
117
118	#endif
119
120	return;
121	}
122
123	int DumpReader::getNFrames(void) {
124
125	if (!isScanned_)
126	scanFile();
127
128	return nframes_;
129	}
130
131	void DumpReader::scanFile(void) {
132	int lineNo = 0;
133	std::streampos prevPos;
134	std::streampos currPos;
135
136	#ifdef IS_MPI
137
138	if (worldRank == 0) {
139	#endif // is_mpi
140
141	currPos = inFile_->tellg();
142	prevPos = currPos;
143	bool foundOpenSnapshotTag = false;
144	bool foundClosedSnapshotTag = false;
145	while(inFile_->getline(buffer, bufferSize)) {
146	++lineNo;
147
148	std::string line = buffer;
149	currPos = inFile_->tellg();
150	if (line.find("<Snapshot>")!= std::string::npos) {
151	if (foundOpenSnapshotTag) {
152	sprintf(painCave.errMsg,
153	"DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
154	filename_.c_str());
155	painCave.isFatal = 1;
156	simError();
157	}
158	foundOpenSnapshotTag = true;
159	foundClosedSnapshotTag = false;
160	framePos_.push_back(prevPos);
161
162	} else if (line.find("</Snapshot>") != std::string::npos){
163	if (!foundOpenSnapshotTag) {
164	sprintf(painCave.errMsg,
165	"DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
166	filename_.c_str());
167	painCave.isFatal = 1;
168	simError();
169	}
170
171	if (foundClosedSnapshotTag) {
172	sprintf(painCave.errMsg,
173	"DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
174	filename_.c_str());
175	painCave.isFatal = 1;
176	simError();
177	}
178	foundClosedSnapshotTag = true;
179	foundOpenSnapshotTag = false;
180	}
181	prevPos = currPos;
182	}
183
184	// only found <Snapshot> for the last frame means the file is corrupted, we should discard
185	// it and give a warning message
186	if (foundOpenSnapshotTag) {
187	sprintf(painCave.errMsg,
188	"DumpReader: last frame in %s is invalid\n", filename_.c_str());
189	painCave.isFatal = 0;
190	simError();
191	framePos_.pop_back();
192	}
193
194	nframes_ = framePos_.size();
195
196	if (nframes_ == 0) {
197	sprintf(painCave.errMsg,
198	"DumpReader: %s does not contain a valid frame\n", filename_.c_str());
199	painCave.isFatal = 1;
200	simError();
201	}
202	#ifdef IS_MPI
203	}
204
205	MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
206
207	#endif // is_mpi
208
209	isScanned_ = true;
210	}
211
212	void DumpReader::readFrame(int whichFrame) {
213	if (!isScanned_)
214	scanFile();
215
216	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
217
218	if (storageLayout & DataStorage::dslPosition) {
219	needPos_ = true;
220	} else {
221	needPos_ = false;
222	}
223
224	if (storageLayout & DataStorage::dslVelocity) {
225	needVel_ = true;
226	} else {
227	needVel_ = false;
228	}
229
230	if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
231	needQuaternion_ = true;
232	} else {
233	needQuaternion_ = false;
234	}
235
236	if (storageLayout & DataStorage::dslAngularMomentum) {
237	needAngMom_ = true;
238	} else {
239	needAngMom_ = false;
240	}
241
242	readSet(whichFrame);
243
244	if (needCOMprops_) {
245	Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
246	Vector3d com;
247	Vector3d comvel;
248	Vector3d comw;
249	info_->getComAll(com, comvel);
250	comw = info_->getAngularMomentum();
251	s->setCOMprops(com, comvel, comw);
252	}
253
254	}
255
256	void DumpReader::readSet(int whichFrame) {
257	std::string line;
258
259	#ifndef IS_MPI
260
261	inFile_->clear();
262	inFile_->seekg(framePos_[whichFrame]);
263
264	std::istream& inputStream = *inFile_;
265
266	#else
267	int masterNode = 0;
268	std::stringstream sstream;
269	if (worldRank == masterNode) {
270	std::string sendBuffer;
271
272	inFile_->clear();
273	inFile_->seekg(framePos_[whichFrame]);
274
275	while (inFile_->getline(buffer, bufferSize)) {
276
277	line = buffer;
278	sendBuffer += line;
279	sendBuffer += '\n';
280	if (line.find("</Snapshot>") != std::string::npos) {
281	break;
282	}
283	}
284
285	int sendBufferSize = sendBuffer.size();
286	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
287	MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
288
289	sstream.str(sendBuffer);
290	} else {
291	int sendBufferSize;
292	MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
293	char * recvBuffer = new char[sendBufferSize+1];
294	MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
295	sstream.str(recvBuffer);
296	}
297
298	std::istream& inputStream = sstream;
299	#endif
300
301	inputStream.getline(buffer, bufferSize);
302
303	line = buffer;
304	if (line.find("<Snapshot>") == std::string::npos) {
305	sprintf(painCave.errMsg,
306	"DumpReader Error: can not find <Snapshot>\n");
307	painCave.isFatal = 1;
308	simError();
309	}
310
311	//read frameData
312	readFrameProperties(inputStream);
313
314	//read StuntDoubles
315	readStuntDoubles(inputStream);
316
317	inputStream.getline(buffer, bufferSize);
318	line = buffer;
319	if (line.find("</Snapshot>") == std::string::npos) {
320	sprintf(painCave.errMsg,
321	"DumpReader Error: can not find </Snapshot>\n");
322	painCave.isFatal = 1;
323	simError();
324	}
325
326	}
327
328	void DumpReader::parseDumpLine(const std::string& line) {
329
330
331	StringTokenizer tokenizer(line);
332	int nTokens;
333
334	nTokens = tokenizer.countTokens();
335
336	if (nTokens < 2) {
337	sprintf(painCave.errMsg,
338	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
339	painCave.isFatal = 1;
340	simError();
341	}
342
343	int index = tokenizer.nextTokenAsInt();
344
345	StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
346
347	if (integrableObject == NULL) {
348	return;
349	}
350	std::string type = tokenizer.nextToken();
351	int size = type.size();
352
353	for(int i = 0; i < size; ++i) {
354	switch(type[i]) {
355
356	case 'p': {
357	Vector3d pos;
358	pos[0] = tokenizer.nextTokenAsDouble();
359	pos[1] = tokenizer.nextTokenAsDouble();
360	pos[2] = tokenizer.nextTokenAsDouble();
361	if (needPos_) {
362	integrableObject->setPos(pos);
363	}
364	break;
365	}
366	case 'v' : {
367	Vector3d vel;
368	vel[0] = tokenizer.nextTokenAsDouble();
369	vel[1] = tokenizer.nextTokenAsDouble();
370	vel[2] = tokenizer.nextTokenAsDouble();
371	if (needVel_) {
372	integrableObject->setVel(vel);
373	}
374	break;
375	}
376
377	case 'q' : {
378	Quat4d q;
379	if (integrableObject->isDirectional()) {
380
381	q[0] = tokenizer.nextTokenAsDouble();
382	q[1] = tokenizer.nextTokenAsDouble();
383	q[2] = tokenizer.nextTokenAsDouble();
384	q[3] = tokenizer.nextTokenAsDouble();
385
386	RealType qlen = q.length();
387	if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
388
389	sprintf(painCave.errMsg,
390	"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
391	painCave.isFatal = 1;
392	simError();
393
394	}
395
396	q.normalize();
397	if (needQuaternion_) {
398	integrableObject->setQ(q);
399	}
400	}
401	break;
402	}
403	case 'j' : {
404	Vector3d ji;
405	if (integrableObject->isDirectional()) {
406	ji[0] = tokenizer.nextTokenAsDouble();
407	ji[1] = tokenizer.nextTokenAsDouble();
408	ji[2] = tokenizer.nextTokenAsDouble();
409	if (needAngMom_) {
410	integrableObject->setJ(ji);
411	}
412	}
413	break;
414	}
415	case 'f': {
416
417	Vector3d force;
418	force[0] = tokenizer.nextTokenAsDouble();
419	force[1] = tokenizer.nextTokenAsDouble();
420	force[2] = tokenizer.nextTokenAsDouble();
421	integrableObject->setFrc(force);
422	break;
423	}
424	case 't' : {
425
426	Vector3d torque;
427	torque[0] = tokenizer.nextTokenAsDouble();
428	torque[1] = tokenizer.nextTokenAsDouble();
429	torque[2] = tokenizer.nextTokenAsDouble();
430	integrableObject->setTrq(torque);
431	break;
432	}
433	default: {
434	sprintf(painCave.errMsg,
435	"DumpReader Error: %s is an unrecognized type\n", type.c_str());
436	painCave.isFatal = 1;
437	simError();
438	break;
439	}
440
441	}
442	}
443
444	}
445
446
447	void DumpReader::readStuntDoubles(std::istream& inputStream) {
448
449	inputStream.getline(buffer, bufferSize);
450	std::string line(buffer);
451
452	if (line.find("<StuntDoubles>") == std::string::npos) {
453	sprintf(painCave.errMsg,
454	"DumpReader Error: Missing <StuntDoubles>\n");
455	painCave.isFatal = 1;
456	simError();
457	}
458
459	while(inputStream.getline(buffer, bufferSize)) {
460	line = buffer;
461
462	if(line.find("</StuntDoubles>") != std::string::npos) {
463	break;
464	}
465
466	parseDumpLine(line);
467	}
468
469	}
470
471	void DumpReader::readFrameProperties(std::istream& inputStream) {
472
473	Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
474	inputStream.getline(buffer, bufferSize);
475	std::string line(buffer);
476
477	if (line.find("<FrameData>") == std::string::npos) {
478	sprintf(painCave.errMsg,
479	"DumpReader Error: Missing <FrameData>\n");
480	painCave.isFatal = 1;
481	simError();
482	}
483
484	while(inputStream.getline(buffer, bufferSize)) {
485	line = buffer;
486
487	if(line.find("</FrameData>") != std::string::npos) {
488	break;
489	}
490
491	StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
492	if (!tokenizer.hasMoreTokens()) {
493	sprintf(painCave.errMsg,
494	"DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
495	painCave.isFatal = 1;
496	simError();
497	}
498
499	std::string propertyName = tokenizer.nextToken();
500	if (propertyName == "Time") {
501	RealType currTime = tokenizer.nextTokenAsDouble();
502	s->setTime(currTime);
503	} else if (propertyName == "Hmat"){
504	Mat3x3d hmat;
505	hmat(0, 0) = tokenizer.nextTokenAsDouble();
506	hmat(0, 1) = tokenizer.nextTokenAsDouble();
507	hmat(0, 2) = tokenizer.nextTokenAsDouble();
508	hmat(1, 0) = tokenizer.nextTokenAsDouble();
509	hmat(1, 1) = tokenizer.nextTokenAsDouble();
510	hmat(1, 2) = tokenizer.nextTokenAsDouble();
511	hmat(2, 0) = tokenizer.nextTokenAsDouble();
512	hmat(2, 1) = tokenizer.nextTokenAsDouble();
513	hmat(2, 2) = tokenizer.nextTokenAsDouble();
514	s->setHmat(hmat);
515	} else if (propertyName == "Thermostat") {
516	RealType chi = tokenizer.nextTokenAsDouble();
517	RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
518	s->setChi(chi);
519	s->setIntegralOfChiDt(integralOfChiDt);
520	} else if (propertyName == "Barostat") {
521	Mat3x3d eta;
522	eta(0, 0) = tokenizer.nextTokenAsDouble();
523	eta(0, 1) = tokenizer.nextTokenAsDouble();
524	eta(0, 2) = tokenizer.nextTokenAsDouble();
525	eta(1, 0) = tokenizer.nextTokenAsDouble();
526	eta(1, 1) = tokenizer.nextTokenAsDouble();
527	eta(1, 2) = tokenizer.nextTokenAsDouble();
528	eta(2, 0) = tokenizer.nextTokenAsDouble();
529	eta(2, 1) = tokenizer.nextTokenAsDouble();
530	eta(2, 2) = tokenizer.nextTokenAsDouble();
531	s->setEta(eta);
532	} else {
533	sprintf(painCave.errMsg,
534	"DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
535	painCave.isFatal = 0;
536	simError();
537	}
538
539	}
540
541	}
542
543
544	}//end namespace oopse