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_ = fopen(filename_.c_str(), "r");

        if (inFile_ == NULL) {
            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

        int error;
        error = fclose(inFile_);

        if (error) {
            sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
            painCave.isFatal = 1;            
            simError();
        }

        MemoryUtils::deletePointers(framePos_);

#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];
  fpos_t * currPos;

#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi

        rewind(inFile_);

        currPos = new fpos_t;
        fgetpos(inFile_, currPos);
        fgets(readBuffer, sizeof(readBuffer), inFile_);
        lineNum++;

        if (feof(inFile_)) {
            sprintf(painCave.errMsg,
                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                    filename_.c_str(),
                    lineNum);
            painCave.isFatal = 1;
            simError();
        }

        while (!feof(inFile_)) {
            framePos_.push_back(currPos);

            i = atoi(readBuffer);

            fgets(readBuffer, sizeof(readBuffer), inFile_);
            lineNum++;

            if (feof(inFile_)) {
                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++) {
                fgets(readBuffer, sizeof(readBuffer), inFile_);
                lineNum++;

                if (feof(inFile_)) {
                    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 = new fpos_t;
            fgetpos(inFile_, currPos);
            fgets(readBuffer, sizeof(readBuffer), inFile_);
            lineNum++;
        }

        delete currPos;
        rewind(inFile_);
        
        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) {
    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

    fsetpos(inFile_, framePos_[whichFrame]);
    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

    if (eof_test == NULL) {
        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

    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

    if (eof_test == NULL) {
        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

    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {

        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {           

            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

            if (eof_test == NULL) {
                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);
            
            }
    }

    // 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) {
        fsetpos(inFile_, framePos_[whichFrame]);

        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

        if (eof_test == NULL) {
            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

        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

        if (eof_test == NULL) {
            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)){
                        
                    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

                    if (eof_test == NULL) {
                        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++) {
                    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);

                    if (eof_test == NULL) {
                        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:	2082
Committed:	Mon Mar 7 22:39:33 2005 UTC (19 years, 4 months ago) by tim
File size:	19673 byte(s)
Log Message:	Fixing a bug in BitSet.cpp
#	User	Rev	Content
1	gezelter	1930	/*
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	gezelter	1490	#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	gezelter	1930	#include "io/DumpReader.hpp"
56			#include "primitives/Molecule.hpp"
57	tim	1492	#include "utils/simError.h"
58	gezelter	1930	#include "utils/MemoryUtils.hpp"
59			#include "utils/StringTokenizer.hpp"
60	gezelter	1490
61			#ifdef IS_MPI
62	gezelter	1930
63	gezelter	1490	#include <mpi.h>
64			#define TAKE_THIS_TAG_CHAR 0
65			#define TAKE_THIS_TAG_INT 1
66	gezelter	1930
67	gezelter	1490	#endif // is_mpi
68
69
70	gezelter	1930	namespace oopse {
71	gezelter	1490
72	gezelter	1930	DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73			: info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74	gezelter	1490
75			#ifdef IS_MPI
76
77	gezelter	1930	if (worldRank == 0) {
78	gezelter	1490	#endif
79
80	gezelter	1930	inFile_ = fopen(filename_.c_str(), "r");
81
82			if (inFile_ == NULL) {
83	tim	1957	sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84	gezelter	1930	painCave.isFatal = 1;
85			simError();
86			}
87
88	gezelter	1490	#ifdef IS_MPI
89
90	gezelter	1930	}
91	gezelter	1490
92	gezelter	1930	strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93			MPIcheckPoint();
94
95	gezelter	1490	#endif
96
97	gezelter	1930	return;
98	gezelter	1490	}
99
100	gezelter	1930	DumpReader::~DumpReader() {
101	gezelter	1490
102	gezelter	1930	#ifdef IS_MPI
103	gezelter	1490
104	gezelter	1930	if (worldRank == 0) {
105			#endif
106	gezelter	1490
107	gezelter	1930	int error;
108			error = fclose(inFile_);
109	gezelter	1490
110	gezelter	1930	if (error) {
111	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
112	gezelter	1930	painCave.isFatal = 1;
113			simError();
114			}
115
116	tim	2082	MemoryUtils::deletePointers(framePos_);
117	gezelter	1930
118	gezelter	1490	#ifdef IS_MPI
119	gezelter	1930
120	gezelter	1490	}
121
122	gezelter	1930	strcpy(checkPointMsg, "Dump file closed successfully.");
123			MPIcheckPoint();
124	gezelter	1490
125	gezelter	1930	#endif
126	gezelter	1490
127	gezelter	1930	return;
128	gezelter	1490	}
129
130	gezelter	1930	int DumpReader::getNFrames(void) {
131	gezelter	1490
132	gezelter	1930	if (!isScanned_)
133			scanFile();
134	gezelter	1490
135	gezelter	1930	return nframes_;
136	gezelter	1490	}
137
138	gezelter	1930	void DumpReader::scanFile(void) {
139			int i, j;
140			int lineNum = 0;
141			char readBuffer[maxBufferSize];
142			fpos_t * currPos;
143	gezelter	1490
144			#ifdef IS_MPI
145
146	gezelter	1930	if (worldRank == 0) {
147			#endif // is_mpi
148	gezelter	1490
149	gezelter	1930	rewind(inFile_);
150	gezelter	1490
151	gezelter	1930	currPos = new fpos_t;
152			fgetpos(inFile_, currPos);
153			fgets(readBuffer, sizeof(readBuffer), inFile_);
154			lineNum++;
155	gezelter	1490
156	gezelter	1930	if (feof(inFile_)) {
157			sprintf(painCave.errMsg,
158	tim	1958	"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
159	gezelter	1930	filename_.c_str(),
160			lineNum);
161			painCave.isFatal = 1;
162			simError();
163			}
164	gezelter	1490
165	gezelter	1930	while (!feof(inFile_)) {
166			framePos_.push_back(currPos);
167	gezelter	1490
168	gezelter	1930	i = atoi(readBuffer);
169	gezelter	1490
170	gezelter	1930	fgets(readBuffer, sizeof(readBuffer), inFile_);
171			lineNum++;
172	gezelter	1490
173	gezelter	1930	if (feof(inFile_)) {
174			sprintf(painCave.errMsg,
175	tim	1958	"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
176	gezelter	1930	filename_.c_str(),
177			lineNum);
178			painCave.isFatal = 1;
179			simError();
180			}
181	gezelter	1490
182	gezelter	1930	for(j = 0; j < i; j++) {
183			fgets(readBuffer, sizeof(readBuffer), inFile_);
184			lineNum++;
185	gezelter	1490
186	gezelter	1930	if (feof(inFile_)) {
187			sprintf(painCave.errMsg,
188	tim	1958	"DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
189	gezelter	1930	" with atom %d\n", filename_.c_str(),
190			lineNum,
191			j);
192	gezelter	1490
193	gezelter	1930	painCave.isFatal = 1;
194			simError();
195			}
196			}
197	gezelter	1490
198	gezelter	1930	currPos = new fpos_t;
199			fgetpos(inFile_, currPos);
200			fgets(readBuffer, sizeof(readBuffer), inFile_);
201			lineNum++;
202			}
203	gezelter	1490
204	gezelter	1930	delete currPos;
205			rewind(inFile_);
206
207			nframes_ = framePos_.size();
208			#ifdef IS_MPI
209	gezelter	1490	}
210
211	gezelter	1930	MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
212	gezelter	1490
213	gezelter	1930	strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
214			MPIcheckPoint();
215	gezelter	1490
216	gezelter	1930	#endif // is_mpi
217	gezelter	1490
218	gezelter	1930	isScanned_ = true;
219			}
220	gezelter	1490
221	gezelter	1930	void DumpReader::readFrame(int whichFrame) {
222	tim	2002	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
223
224			if (storageLayout & DataStorage::dslPosition) {
225			needPos_ = true;
226			} else {
227			needPos_ = false;
228			}
229
230			if (storageLayout & DataStorage::dslVelocity) {
231			needVel_ = true;
232			} else {
233			needVel_ = false;
234			}
235
236			if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
237			needQuaternion_ = true;
238			} else {
239			needQuaternion_ = false;
240			}
241
242			if (storageLayout & DataStorage::dslAngularMomentum) {
243			needAngMom_ = true;
244			} else {
245			needAngMom_ = false;
246			}
247
248	gezelter	1930	readSet(whichFrame);
249			}
250	gezelter	1490
251	gezelter	1930	void DumpReader::readSet(int whichFrame) {
252			int i;
253			int nTotObjs; // the number of atoms
254			char read_buffer[maxBufferSize]; //the line buffer for reading
255			char * eof_test; // ptr to see when we reach the end of the file
256	gezelter	1490
257	gezelter	1930	Molecule* mol;
258			StuntDouble* integrableObject;
259			SimInfo::MoleculeIterator mi;
260			Molecule::IntegrableObjectIterator ii;
261	gezelter	1490
262	gezelter	1930	#ifndef IS_MPI
263	gezelter	1490
264	gezelter	1930	fsetpos(inFile_, framePos_[whichFrame]);
265			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
266	gezelter	1490
267	gezelter	1930	if (eof_test == NULL) {
268			sprintf(painCave.errMsg,
269			"DumpReader error: error reading 1st line of \"%s\"\n",
270			filename_.c_str());
271			painCave.isFatal = 1;
272			simError();
273	gezelter	1490	}
274
275	gezelter	1930	nTotObjs = atoi(read_buffer);
276	gezelter	1490
277	gezelter	1930	if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
278			sprintf(painCave.errMsg,
279			"DumpReader error. %s nIntegrable, %d, "
280			"does not match the meta-data file's nIntegrable, %d.\n",
281			filename_.c_str(),
282			nTotObjs,
283	tim	1957	info_->getNGlobalIntegrableObjects());
284	gezelter	1490
285	gezelter	1930	painCave.isFatal = 1;
286			simError();
287	gezelter	1490	}
288
289	gezelter	1930	//read the box mat from the comment line
290	gezelter	1490
291	gezelter	1930	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
292	gezelter	1490
293	gezelter	1930	if (eof_test == NULL) {
294	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
295	gezelter	1930	filename_.c_str());
296			painCave.isFatal = 1;
297	gezelter	1490	simError();
298	gezelter	1930	}
299	gezelter	1490
300	gezelter	1930	parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
301	gezelter	1490
302	gezelter	1930	//parse dump lines
303	gezelter	1490
304	gezelter	1930	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
305	gezelter	1490
306	gezelter	1930	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
307			integrableObject = mol->nextIntegrableObject(ii)) {
308	gezelter	1490
309	gezelter	1930	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
310	gezelter	1490
311	gezelter	1930	if (eof_test == NULL) {
312			sprintf(painCave.errMsg,
313	tim	1958	"DumpReader Error: error in reading file %s\n"
314	gezelter	1930	"natoms = %d; index = %d\n"
315			"error reading the line from the file.\n",
316			filename_.c_str(),
317			nTotObjs,
318			i);
319	gezelter	1490
320	gezelter	1930	painCave.isFatal = 1;
321			simError();
322			}
323
324			parseDumpLine(read_buffer, integrableObject);
325
326			}
327	gezelter	1490	}
328
329	gezelter	1930	// MPI Section of code..........
330	gezelter	1490
331	gezelter	1930	#else //IS_MPI
332	gezelter	1490
333	gezelter	1930	// first thing first, suspend fatalities.
334			int masterNode = 0;
335			int nCurObj;
336			painCave.isEventLoop = 1;
337	gezelter	1490
338	gezelter	1930	int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
339			int haveError;
340	gezelter	1490
341	gezelter	1930	MPI_Status istatus;
342			int nitems;
343	gezelter	1490
344	gezelter	1930	nTotObjs = info_->getNGlobalIntegrableObjects();
345			haveError = 0;
346	gezelter	1490
347	gezelter	1930	if (worldRank == masterNode) {
348			fsetpos(inFile_, framePos_[whichFrame]);
349	gezelter	1490
350	gezelter	1930	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
351	gezelter	1490
352	gezelter	1930	if (eof_test == NULL) {
353	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
354	gezelter	1930	filename_.c_str());
355			painCave.isFatal = 1;
356			simError();
357	gezelter	1490	}
358
359	gezelter	1930	nitems = atoi(read_buffer);
360	gezelter	1490
361	gezelter	1930	// Check to see that the number of integrable objects in the
362			// intial configuration file is the same as derived from the
363			// meta-data file.
364	gezelter	1490
365	gezelter	1930	if (nTotObjs != nitems) {
366			sprintf(painCave.errMsg,
367			"DumpReader Error. %s nIntegrable, %d, "
368			"does not match the meta-data file's nIntegrable, %d.\n",
369			filename_.c_str(),
370			nTotObjs,
371			info_->getNGlobalIntegrableObjects());
372	gezelter	1490
373	gezelter	1930	painCave.isFatal = 1;
374			simError();
375			}
376	gezelter	1490
377	gezelter	1930	//read the boxMat from the comment line
378	gezelter	1490
379	gezelter	1930	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
380	gezelter	1490
381	gezelter	1930	if (eof_test == NULL) {
382	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
383	gezelter	1930	filename_.c_str());
384			painCave.isFatal = 1;
385			simError();
386			}
387	gezelter	1490
388	gezelter	1930	//Every single processor will parse the comment line by itself
389			//By using this way, we might lose some efficiency, but if we want to add
390			//more parameters into comment line, we only need to modify function
391			//parseCommentLine
392	gezelter	1490
393	gezelter	1930	MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
394			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
395	gezelter	1490
396	gezelter	1930	for(i = 0; i < info_->getNGlobalMolecules(); i++) {
397			int which_node = info_->getMolToProc(i);
398	gezelter	1490
399	gezelter	1930	if (which_node == masterNode) {
400			//molecules belong to master node
401	gezelter	1490
402	gezelter	1930	mol = info_->getMoleculeByGlobalIndex(i);
403	gezelter	1490
404	gezelter	1930	if (mol == NULL) {
405	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
406	gezelter	1930	painCave.isFatal = 1;
407			simError();
408			}
409	gezelter	1490
410	gezelter	1930	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
411			integrableObject = mol->nextIntegrableObject(ii)){
412
413			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
414	gezelter	1490
415	gezelter	1930	if (eof_test == NULL) {
416			sprintf(painCave.errMsg,
417	tim	1958	"DumpReader Error: error in reading file %s\n"
418	gezelter	1930	"natoms = %d; index = %d\n"
419			"error reading the line from the file.\n",
420			filename_.c_str(),
421			nTotObjs,
422			i);
423	gezelter	1490
424	gezelter	1930	painCave.isFatal = 1;
425			simError();
426			}
427	gezelter	1490
428	gezelter	1930	parseDumpLine(read_buffer, integrableObject);
429			}
430			} else {
431			//molecule belongs to slave nodes
432	gezelter	1490
433	gezelter	1930	MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
434			MPI_COMM_WORLD, &istatus);
435	gezelter	1490
436	gezelter	1930	for(int j = 0; j < nCurObj; j++) {
437			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
438	gezelter	1490
439	gezelter	1930	if (eof_test == NULL) {
440			sprintf(painCave.errMsg,
441	tim	1958	"DumpReader Error: error in reading file %s\n"
442	gezelter	1930	"natoms = %d; index = %d\n"
443			"error reading the line from the file.\n",
444			filename_.c_str(),
445			nTotObjs,
446			i);
447	gezelter	1490
448	gezelter	1930	painCave.isFatal = 1;
449			simError();
450			}
451
452			MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
453			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
454			}
455			}
456			}
457			} else {
458			//actions taken at slave nodes
459			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
460	gezelter	1490
461	gezelter	1930	/*@todo/
462			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
463	gezelter	1490
464	gezelter	1930	for(i = 0; i < info_->getNGlobalMolecules(); i++) {
465			int which_node = info_->getMolToProc(i);
466	gezelter	1490
467	gezelter	1930	if (which_node == worldRank) {
468			//molecule with global index i belongs to this processor
469
470			mol = info_->getMoleculeByGlobalIndex(i);
471			if (mol == NULL) {
472	tim	1958	sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
473	gezelter	1930	painCave.isFatal = 1;
474			simError();
475			}
476
477			nCurObj = mol->getNIntegrableObjects();
478	gezelter	1490
479	gezelter	1930	MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
480			MPI_COMM_WORLD);
481	gezelter	1490
482	gezelter	1930	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
483			integrableObject = mol->nextIntegrableObject(ii)){
484
485			MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
486			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
487	gezelter	1490
488	gezelter	1930	parseDumpLine(read_buffer, integrableObject);
489			}
490
491			}
492
493			}
494
495	gezelter	1490	}
496
497	gezelter	1930	#endif
498	gezelter	1490
499	gezelter	1930	}
500	gezelter	1490
501	gezelter	1930	void DumpReader::parseDumpLine(char line, StuntDouble integrableObject) {
502	gezelter	1490
503	gezelter	1930	Vector3d pos; // position place holders
504			Vector3d vel; // velocity placeholders
505			Quat4d q; // the quaternions
506			Vector3d ji; // angular velocity placeholders;
507			StringTokenizer tokenizer(line);
508			int nTokens;
509
510			nTokens = tokenizer.countTokens();
511	gezelter	1490
512	gezelter	1930	if (nTokens < 14) {
513			sprintf(painCave.errMsg,
514	tim	2065	"DumpReader Error: Not enough Tokens.\n%s\n", line);
515	gezelter	1930	painCave.isFatal = 1;
516			simError();
517	gezelter	1490	}
518
519	gezelter	1930	std::string name = tokenizer.nextToken();
520
521			if (name != integrableObject->getType()) {
522	tim	1958
523			sprintf(painCave.errMsg,
524			"DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
525			name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
526			painCave.isFatal = 1;
527			simError();
528	gezelter	1490	}
529
530	gezelter	1930	pos[0] = tokenizer.nextTokenAsDouble();
531			pos[1] = tokenizer.nextTokenAsDouble();
532			pos[2] = tokenizer.nextTokenAsDouble();
533	tim	2002	if (needPos_) {
534			integrableObject->setPos(pos);
535			}
536	gezelter	1930
537			vel[0] = tokenizer.nextTokenAsDouble();
538			vel[1] = tokenizer.nextTokenAsDouble();
539			vel[2] = tokenizer.nextTokenAsDouble();
540	tim	2002	if (needVel_) {
541			integrableObject->setVel(vel);
542			}
543
544	gezelter	1930	if (integrableObject->isDirectional()) {
545
546			q[0] = tokenizer.nextTokenAsDouble();
547			q[1] = tokenizer.nextTokenAsDouble();
548			q[2] = tokenizer.nextTokenAsDouble();
549			q[3] = tokenizer.nextTokenAsDouble();
550	gezelter	1490
551	gezelter	1930	double qlen = q.length();
552			if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
553
554			sprintf(painCave.errMsg,
555	tim	1958	"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
556	gezelter	1930	painCave.isFatal = 1;
557			simError();
558
559			}
560	gezelter	1490
561	gezelter	1930	q.normalize();
562	tim	2002	if (needQuaternion_) {
563			integrableObject->setQ(q);
564			}
565
566	gezelter	1930	ji[0] = tokenizer.nextTokenAsDouble();
567			ji[1] = tokenizer.nextTokenAsDouble();
568			ji[2] = tokenizer.nextTokenAsDouble();
569	tim	2002	if (needAngMom_) {
570			integrableObject->setJ(ji);
571			}
572	gezelter	1490	}
573
574			}
575
576
577	gezelter	1930	void DumpReader::parseCommentLine(char* line, Snapshot* s) {
578			double currTime;
579			Mat3x3d hmat;
580			double chi;
581			double integralOfChiDt;
582			Mat3x3d eta;
583	gezelter	1490
584	gezelter	1930	StringTokenizer tokenizer(line);
585			int nTokens;
586	gezelter	1490
587	gezelter	1930	nTokens = tokenizer.countTokens();
588	gezelter	1490
589	gezelter	1930	//comment line should at least contain 10 tokens: current time(1 token) and h-matrix(9 tokens)
590			if (nTokens < 10) {
591			sprintf(painCave.errMsg,
592	tim	1958	"DumpReader Error: Not enough tokens in comment line: %s", line);
593	gezelter	1930	painCave.isFatal = 1;
594			simError();
595	gezelter	1490	}
596
597	gezelter	1930	//read current time
598			currTime = tokenizer.nextTokenAsDouble();
599			s->setTime(currTime);
600
601			//read h-matrix
602			hmat(0, 0) = tokenizer.nextTokenAsDouble();
603			hmat(0, 1) = tokenizer.nextTokenAsDouble();
604			hmat(0, 2) = tokenizer.nextTokenAsDouble();
605			hmat(1, 0) = tokenizer.nextTokenAsDouble();
606			hmat(1, 1) = tokenizer.nextTokenAsDouble();
607			hmat(1, 2) = tokenizer.nextTokenAsDouble();
608			hmat(2, 0) = tokenizer.nextTokenAsDouble();
609			hmat(2, 1) = tokenizer.nextTokenAsDouble();
610			hmat(2, 2) = tokenizer.nextTokenAsDouble();
611			s->setHmat(hmat);
612
613	tim	2002	//read chi and integralOfChidt, they should apprear in pair
614	gezelter	1930	if (tokenizer.countTokens() >= 2) {
615			chi = tokenizer.nextTokenAsDouble();
616			integralOfChiDt = tokenizer.nextTokenAsDouble();
617	gezelter	1490
618	gezelter	1930	s->setChi(chi);
619			s->setIntegralOfChiDt(integralOfChiDt);
620	gezelter	1490	}
621
622	gezelter	1930	//read eta (eta is 3x3 matrix)
623			if (tokenizer.countTokens() >= 9) {
624			eta(0, 0) = tokenizer.nextTokenAsDouble();
625			eta(0, 1) = tokenizer.nextTokenAsDouble();
626			eta(0, 2) = tokenizer.nextTokenAsDouble();
627			eta(1, 0) = tokenizer.nextTokenAsDouble();
628			eta(1, 1) = tokenizer.nextTokenAsDouble();
629			eta(1, 2) = tokenizer.nextTokenAsDouble();
630			eta(2, 0) = tokenizer.nextTokenAsDouble();
631			eta(2, 1) = tokenizer.nextTokenAsDouble();
632			eta(2, 2) = tokenizer.nextTokenAsDouble();
633
634			s->setEta(eta);
635	gezelter	1490	}
636	gezelter	1930
637	gezelter	1490
638			}
639
640	gezelter	1930	}//end namespace oopse