ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/visitors/ZconsVisitor.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/visitors/ZconsVisitor.cpp (file contents):
Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC

# Line 1 | Line 1
1 < #include "visitors/ZconsVisitor.hpp"
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   #include <cmath>
43 + #include "visitors/ZconsVisitor.hpp"
44 + #include "primitives/Molecule.hpp"
45 + #include "utils/StringUtils.hpp"
46 + namespace oopse {
47  
48 < ZConsVisitor::ZConsVisitor(SimInfo* info) : BaseVisitor(), zconsReader(NULL){
5 <  GenericData* data;
6 <  DoubleData* tolerance;  
7 <  ZConsParaData* zConsParaData;
8 <  StringData* filename;
9 <  DoubleData* sampleTime;
10 <  vector<ZConsParaItem>* parameters;
48 > ZConsVisitor::ZConsVisitor(SimInfo* info) : BaseVisitor(), info_(info), zconsReader_(NULL){
49  
50 <  this->info = info;
51 <  visitorName = "ZConsVisitor";
52 <  
15 <  //retrieve tolerance for z-constraint molecuels
16 <  data = info->getProperty(ZCONSTOL_ID);
50 >    visitorName = "ZConsVisitor";
51 >    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
52 >    Globals* simParam = info_->getSimParams();
53  
54 <  if (!data){
55 <    cerr << "Can not get zconstraint tolerance from SimInfo" << endl;
20 <    haveZcons = false;  
21 <    return;
22 <  }
23 <  else{
24 <    tolerance = dynamic_cast<DoubleData*>(data);
25 <
26 <    if (!tolerance){
27 <      cerr << "Can not get zconstraint tolerance  from SimInfo" << endl;
28 <      haveZcons = false;    
29 <      return;
54 >    if (simParam->haveZconstraintTime()){
55 >        zconsTime_ = simParam->getZconsTime();
56      }
57      else{
58 <      zconsTol = tolerance->getData();
58 >        sprintf(painCave.errMsg,
59 >        "ZConstraint error: If you use a ZConstraint,\n"
60 >        "\tyou must set zconsTime.\n");
61 >        painCave.isFatal = 1;
62 >        simError();
63      }
34  }
64  
65 <  //retrieve sample time of z-contraint
66 <  data = info->getProperty(ZCONSTIME_ID);
38 <
39 <  if (!data){
40 <    cerr << "Can not get zcons time  from SimInfo" << endl;
41 <    haveZcons = false;    
42 <    return;
43 <  }
44 <  else{
45 <    sampleTime = dynamic_cast<DoubleData*>(data);
46 <
47 <    if (!sampleTime){
48 <      cerr << "Can not get zcons time  from SimInfo" << endl;
49 <      haveZcons = false;    
50 <      return;
65 >    if (simParam->haveZconsTol()){
66 >        zconsTol_ = simParam->getZconsTol();
67      }
68      else{
69 <     zconsTime = sampleTime->getData();
69 >        zconsTol_ = 0.01;
70 >        sprintf(painCave.errMsg,
71 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
72 >            "\tOOPSE will use a default value of %f.\n"
73 >            "\tTo set the tolerance, use the zconsTol variable.\n",
74 >            zconsTol_);
75 >        painCave.isFatal = 0;
76 >        simError();      
77 >    }    
78 >        
79 >    int nZconstraints = simParam->getNzConstraints();
80 >    ZconStamp** stamp = simParam->getZconStamp();
81 >    for (int i = 0; i < nZconstraints; i++){
82 >        int zmolIndex = stamp[i]->getMolIndex();
83 >        zmolStates_.insert(std::make_pair(zmolIndex, zsMoving));
84      }
55  }
85  
57  //retrieve index of z-constraint molecules
58  data = info->getProperty(ZCONSPARADATA_ID);
59  if (!data){
60    cerr << "Can not get index of z-constraint molecules from SimInfo" << endl;
61    haveZcons = false;
62    return;
63  }
64  else{
65    zConsParaData = dynamic_cast<ZConsParaData*>(data);
86  
87 <    if (!zConsParaData){
88 <      cerr << "Can not get index of z-constraint molecules from SimInfo" << endl;
89 <      haveZcons = false;
90 <      return;
87 >    //fill zatomToZmol_ array
88 >    /** @todo only works for single version now*/
89 >    std::map<int, ZConsState>::iterator j;
90 >    for (j = zmolStates_.begin(); j != zmolStates_.end(); ++j) {
91 >        Molecule* mol = info_->getMoleculeByGlobalIndex(j->first);
92 >        assert(mol != NULL);
93 >        Molecule::AtomIterator ai;
94 >        Atom* at;
95 >        for (at = mol->beginAtom(ai); at != NULL; at = mol->nextAtom(ai)) {
96 >            zatomToZmol_.insert(std::make_pair(at->getGlobalIndex(), mol->getGlobalIndex()));
97 >        }
98      }
72    else{
73      vector<ZConsParaItem>::iterator i;
74      Molecule* mol;
99  
100 <      parameters = zConsParaData->getData();
101 <      for(i = parameters->begin(); i != parameters->end(); ++i){
100 >    zconsFilename_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
101 >    
102 >    zconsReader_ = new ZConsReader(info);
103  
104 <        mol = findZconsMol(i->zconsIndex);
105 <        if(mol != NULL)
81 <          zconsMol.push_back(mol);
82 <
83 <      }
84 <
85 <      if(zconsMol.size() < 1){
86 <        cerr << "number of zconstraint molecules is less than one" << endl;
87 <        haveZcons = false;
88 <        return;
89 <      }
90 <      
91 <    }//end if (!zConsParaData)
92 <
93 <  }//end  if (!data  
104 >    if (zconsReader_->hasNextFrame())
105 >        zconsReader_->readNextFrame();
106    
95  //retrieve output filename of z force
96  data = info->getProperty(ZCONSFILENAME_ID);
97  if (!data){
98    cerr << "Can not get filename of z-constraint output from SimInfo" << endl;
99    haveZcons = false;
100    return;
101  }
102  else{
103    filename = dynamic_cast<StringData*>(data);
104
105    if (!filename){
106      cerr << "Can not get filename of z-constraint output from SimInfo" << endl;
107      haveZcons = false;
108      return;
109    }
110    else{
111      zconsFilename = filename->getData();
112    }
113  }
114
115  zconsReader = new ZConsReader(info);
116
117  if (zconsReader->hasNextFrame())
118  zconsReader->readNextFrame();
119  
120  haveZcons = true;
107   }
108  
109   ZConsVisitor::~ZConsVisitor(){
110 <  if(!zconsReader)
111 <    delete zconsReader;
110 >  if(!zconsReader_)
111 >    delete zconsReader_;
112    
113   }
114  
115   void ZConsVisitor::visit(Atom* atom){
116 <  string prefix;
117 <  if(isZconstraint(atom->getIndex(), prefix))
118 <    internalVisit(atom, prefix);
116 >    std::string prefix;
117 >    if(isZconstraint(atom->getGlobalIndex(), prefix))
118 >        internalVisit(atom, prefix);
119   }
120  
121   void ZConsVisitor::visit(DirectionalAtom* datom){
122 <  string prefix;
137 <  
138 <  if(isZconstraint(datom->getIndex(), prefix))
139 <    internalVisit(datom, prefix);
140 < }
122 >    std::string prefix;
123  
124 < void ZConsVisitor::visit(RigidBody* rb){
125 <  string prefix;
144 <  vector<Atom*> atoms;
145 <  
146 <  atoms = rb->getAtoms();
147 <    
148 <  if(isZconstraint(atoms[0]->getIndex(), prefix))
149 <    internalVisit(rb, prefix);
124 >    if(isZconstraint(datom->getGlobalIndex(), prefix))
125 >        internalVisit(datom, prefix);
126   }
127  
128 < Molecule* ZConsVisitor::findZconsMol(int index){
129 <  Molecule* mols;
130 <  mols = info->molecules;
155 <  for(int i =0; i < info->n_mol; i++)
156 <    if(index == mols[i].getMyIndex())
157 <      return &mols[i];
128 > void ZConsVisitor::visit(RigidBody* rb){
129 >    std::string prefix;
130 >    std::vector<Atom*> atoms;
131  
132 <  return NULL;
132 >    atoms = rb->getAtoms();
133 >
134 >    if(isZconstraint(atoms[0]->getGlobalIndex(), prefix))
135 >        internalVisit(rb, prefix);
136   }
137  
138   void ZConsVisitor::update(){
139 <  Molecule* mol;
140 <  double com[3];
141 <  ZConsState state;
142 <  Atom** atoms;
143 <  int natoms;
139 >    Molecule* mol;
140 >    Vector3d com;
141 >    ZConsState state;
142 >    std::map<int, ZConsState>::iterator i;
143 >    for ( i = zmolStates_.begin(); i != zmolStates_.end(); ++i) {
144 >        i->second = zsMoving;
145 >    }
146 >    
147 >    readZconsFile(currSnapshot_->getTime());
148  
149 <  getZconsPos(info->currentTime);
150 <  
151 <  for(size_t i = 0; i < zconsMol.size(); i++){
152 <    zconsMol[i]->getCOM(com);
149 >    const std::vector<ZconsData>& fixedZmolData = zconsReader_->getFixedZMolData();
150 >    std::vector<ZconsData>::const_iterator j;
151 >    for (j = fixedZmolData.begin(); j != fixedZmolData.end(); ++j) {
152 >        std::map<int, ZConsState>::iterator k = zmolStates_.find(j->zmolIndex);
153 >        assert(k != zmolStates_.end());
154 >        k->second = zsFixed;
155 >    }
156 >    
157 > }
158  
159 <    if(fabs(com[2] - zconsPos[i]) < zconsTol)
160 <      state = zsFixed;
161 <    else
162 <      state = zsMoving;
163 <
164 <    //update the state of zconstraint atom
165 <    natoms = zconsMol[i]->getNAtoms();
166 <    atoms = zconsMol[i]->getMyAtoms();    
167 <    for(int j =0; j < natoms; j++)
168 <      zconsState[atoms[j]->getIndex()] = state;
184 <  }
185 <  
159 > void ZConsVisitor::readZconsFile(double time) {
160 >    double tempTime;
161 >    while(zconsReader_->hasNextFrame()){
162 >        tempTime = zconsReader_->getCurTime();
163 >        if(tempTime >= time) {
164 >            return;
165 >        }
166 >        
167 >        zconsReader_->readNextFrame();
168 >    }
169   }
170  
171 < void ZConsVisitor::getZconsPos(double time){
171 > void ZConsVisitor::internalVisit(StuntDouble* sd, const std::string& prefix){
172 >    GenericData* data;
173 >    AtomData* atomData;
174 >    AtomInfo* atomInfo;
175 >    std::vector<AtomInfo*>::iterator iter;
176  
177 <  vector<double> tempPos;
178 <  vector<double> prevPos;
179 <  double tempTime;
177 >    //if there is not atom data, just skip it
178 >    data = sd->getPropertyByName("ATOMDATA");
179 >    if(data != NULL){
180 >        atomData = dynamic_cast<AtomData*>(data);  
181 >        if(atomData == NULL)
182 >            return;
183 >        }
184 >    else
185 >        return;
186  
187 <  while(true){
188 <    tempTime = zconsReader->getCurTime();
196 <
197 <    zconsPos = zconsReader->getZConsPos();
198 <    
199 <    if(tempTime >= time)
200 <      return;
201 <    
202 <    zconsReader->readNextFrame();
203 <
204 <  }
205 <
187 >    for(atomInfo  = atomData->beginAtomInfo(iter); atomInfo; atomInfo = atomData->nextAtomInfo(iter))
188 >        (atomInfo->AtomType).insert(0, prefix);
189   }
190  
208 void ZConsVisitor::internalVisit(StuntDouble* sd, const string& prefix){
209  GenericData* data;
210  AtomData* atomData;
211  AtomInfo* atomInfo;
212  vector<AtomInfo*>::iterator iter;
191  
192 <  
193 <  //if there is not atom data, just skip it
194 <  data = sd->getProperty("ATOMDATA");
195 <  if(data != NULL){
196 <    atomData = dynamic_cast<AtomData*>(data);  
197 <    if(atomData == NULL)
198 <      return;
221 <  }
222 <  else
223 <    return;
192 > bool ZConsVisitor::isZconstraint(int atomIndex, std::string& prefix){
193 >    std::string prefixString[] = {"ZF", "ZM"};
194 >    std::map<int, int>::iterator i = zatomToZmol_.find(atomIndex);
195 >    if (i ==  zatomToZmol_.end() ){
196 >        prefix = "";
197 >        return false;
198 >    } else {
199  
200 <  for(atomInfo  = atomData->beginAtomInfo(iter); atomInfo; atomInfo = atomData->nextAtomInfo(iter))
201 <    (atomInfo->AtomType).insert(0, prefix);
200 >        std::map<int, ZConsState>::iterator j = zmolStates_.find(i->second);
201 >        assert(j !=zmolStates_.end());
202 >        prefix = prefixString[j->second];
203 >        return true;
204 >    }
205   }
206  
207 + const std::string ZConsVisitor::toString(){
208 +    char buffer[65535];
209 +    std::string result;
210  
211 < bool ZConsVisitor::isZconstraint(int index, string& prefix){
212 <  map<int, ZConsState>::iterator i;
232 <  string prefixString[] = {"ZF", "ZM"};
233 <  
234 <  i = zconsState.find(index);
235 <  if(i !=zconsState.end()){
236 <    prefix = prefixString[(*i).second];
237 <    return true;
238 <  }
239 <  else{
240 <    prefix = "";
241 <    return false;
242 <  }
243 < }
211 >    sprintf(buffer ,"------------------------------------------------------------------\n");
212 >    result += buffer;
213  
214 < const string ZConsVisitor::toString(){
215 <  char buffer[65535];
247 <  string result;
248 <  
249 <  sprintf(buffer ,"------------------------------------------------------------------\n");
250 <  result += buffer;
214 >    sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
215 >    result += buffer;
216  
217 <  sprintf(buffer ,"Visitor name: %s\n", visitorName.c_str());
218 <  result += buffer;
254 <  
255 <  sprintf(buffer , "number of zconstraint molecule: %d\n", zconsMol.size());
256 <  result += buffer;
217 >    sprintf(buffer , "number of zconstraint molecule: %d\n", zmolStates_.size());
218 >    result += buffer;
219  
220 <  sprintf(buffer , "zconstraint tolerance = %lf\n", zconsTol);
221 <  result += buffer;
220 >    sprintf(buffer , "zconstraint tolerance = %lf\n", zconsTol_);
221 >    result += buffer;
222  
223 <  sprintf(buffer , "zconstraint sample time = %lf\n", zconsTime);
224 <  result += buffer;
223 >    sprintf(buffer , "zconstraint sample time = %lf\n", zconsTime_);
224 >    result += buffer;
225  
226 <  sprintf(buffer , "zconstraint output filename = %s\n", zconsFilename.c_str());
265 <  result += buffer;
266 <  
267 <  for(size_t i = 0; i < zconsMol.size(); ++i){
268 <    sprintf(buffer ,"zconstraint molecule[%d] = %d\n", i, zconsMol[i]->getMyIndex());
226 >    sprintf(buffer , "zconstraint output filename = %s\n", zconsFilename_.c_str());
227      result += buffer;
228  
229 <  }
230 <
229 >    std::map<int, ZConsState>::iterator i;
230 >    int j = 0;
231 >    for ( i = zmolStates_.begin(); i != zmolStates_.end(); ++i) {
232 >        sprintf(buffer ,"zconstraint molecule[%d] = %d\n", j++, i->first);
233 >        result += buffer;
234 >    }
235    
236 <  sprintf(buffer ,"------------------------------------------------------------------\n");
237 <  result += buffer;
236 >    sprintf(buffer ,"------------------------------------------------------------------\n");
237 >    result += buffer;
238  
239 <  return result;
239 >    return result;
240   }
241 +
242 +
243 + }//namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines