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, "Error closing %s\n", filename_.c_str());
            painCave.isFatal = 1;            
            simError();
        }

        MemoryUtils::deleteVectorOfPointer(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,
                    "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,
                        "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,
                            "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) {
    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, "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,
                        "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, "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, "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, "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,
                                "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,
                                "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, "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,
                    "Not enough Tokens.\n");
            painCave.isFatal = 1;
            simError();
    }

    std::string name = tokenizer.nextToken();

    if (name != integrableObject->getType()) {
        
    }

    pos[0] = tokenizer.nextTokenAsDouble();
    pos[1] = tokenizer.nextTokenAsDouble();
    pos[2] = tokenizer.nextTokenAsDouble();
    integrableObject->setPos(pos);
    
    vel[0] = tokenizer.nextTokenAsDouble();
    vel[1] = tokenizer.nextTokenAsDouble();
    vel[2] = tokenizer.nextTokenAsDouble();
    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,
                    "initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
            painCave.isFatal = 1;
            simError();
            
        } 

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