ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/Hydro.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/Hydro.cpp (file contents):
Revision 892 by tim, Thu Feb 23 23:16:43 2006 UTC vs.
Revision 1177 by xsun, Tue Aug 14 17:40:33 2007 UTC

# Line 47 | Line 47
47   #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
48   #include "applications/hydrodynamics/HydrodynamicsModelCreator.hpp"
49   #include "applications/hydrodynamics/HydrodynamicsModelFactory.hpp"
50 + #include "applications/hydrodynamics/AnalyticalModel.hpp"
51   #include "applications/hydrodynamics/BeadModel.hpp"
52   #include "applications/hydrodynamics/RoughShell.hpp"
53 + #include "applications/hydrodynamics/ShapeBuilder.hpp"
54   #include "brains/Register.hpp"
55   #include "brains/SimCreator.hpp"
56   #include "brains/SimInfo.hpp"
57 <
57 > #include "utils/StringUtils.hpp"
58 > #include "utils/simError.h"
59 > #include "utils/MemoryUtils.hpp"
60   using namespace oopse;
61  
62 < /** Register different hydrodynamics models */
62 > struct SDShape{
63 >    StuntDouble* sd;
64 >    Shape* shape;
65 > };
66   void registerHydrodynamicsModels();
67 + void writeHydroProps(std::ostream& os);
68  
61 – bool calcHydrodynamicsProp(const std::string& modelType, StuntDouble* sd, const DynamicProperty& param, std::ostream& os,  const std::string& prefix);
62 –
69   int main(int argc, char* argv[]){
70    //register force fields
71    registerForceFields();    
# Line 83 | Line 89 | int main(int argc, char* argv[]){
89      exit(1);
90    }
91    
86 –  mdFileName = dumpFileName;
87 –  mdFileName = mdFileName.substr(0, mdFileName.rfind(".")) + ".md";
88 –
92    if (args_info.output_given){
93      prefix = args_info.output_arg;
94    } else {
95 <    prefix = "hydro";    
95 >    prefix = getPrefix(dumpFileName);    
96    }
97    std::string outputFilename = prefix + ".diff";
98 <
96 <  DynamicProperty param;
97 <  param.insert(DynamicProperty::value_type("Viscosity", args_info.viscosity_arg));
98 <  param.insert(DynamicProperty::value_type("Temperature", args_info.temperature_arg));
99 <  
100 <  if (args_info.sigma_given) {
101 <    param.insert(DynamicProperty::value_type("Sigma", args_info.sigma_arg));
102 <  }
103 <  
104 <  
98 >    
99    //parse md file and set up the system
100    SimCreator creator;
101 <  SimInfo* info = creator.createSim(mdFileName, true);
101 >  SimInfo* info = creator.createSim(dumpFileName, true);
102      
103    SimInfo::MoleculeIterator mi;
104    Molecule* mol;
# Line 114 | Line 108 | int main(int argc, char* argv[]){
108    identMat(0,0) = 1.0;
109    identMat(1,1) = 1.0;
110    identMat(2,2) = 1.0;
111 <  
111 >
112 >  Globals* simParams = info->getSimParams();
113 >  RealType temperature;
114 >  RealType viscosity;
115 >
116 >  if (simParams->haveViscosity()) {
117 >    viscosity = simParams->getViscosity();
118 >  } else {
119 >    sprintf(painCave.errMsg, "viscosity must be set\n");
120 >    painCave.isFatal = 1;
121 >    simError();  
122 >  }
123 >
124 >  if (simParams->haveTargetTemp()) {
125 >    temperature = simParams->getTargetTemp();
126 >  } else {
127 >    sprintf(painCave.errMsg, "target temperature must be set\n");
128 >    painCave.isFatal = 1;
129 >    simError();  
130 >  }
131  
132 <  std::map<std::string, StuntDouble*> uniqueStuntDoubles;
132 >  std::map<std::string, SDShape> uniqueStuntDoubles;
133    
134    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
135        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
136             integrableObject = mol->nextIntegrableObject(ii)) {
137            if (uniqueStuntDoubles.find(integrableObject->getType()) ==  uniqueStuntDoubles.end()) {
138 <            uniqueStuntDoubles.insert(std::map<std::string, StuntDouble*>::value_type(integrableObject->getType(), integrableObject));
138 >
139              integrableObject->setPos(V3Zero);
140              integrableObject->setA(identMat);
141              if (integrableObject->isRigidBody()) {
142                  RigidBody* rb = static_cast<RigidBody*>(integrableObject);
143                  rb->updateAtoms();
144              }
145 +
146 +            SDShape tmp;
147 +            tmp.shape = ShapeBuilder::createShape(integrableObject);
148 +            tmp.sd = integrableObject;    
149 +            uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));
150 +
151            }
152          }
153    }
154  
155 <  std::map<std::string, StuntDouble*>::iterator iter;
137 <  std::ofstream outputDiff(outputFilename.c_str());
138 <  for (iter = uniqueStuntDoubles.begin(); iter != uniqueStuntDoubles.end(); ++iter) {
139 <    calcHydrodynamicsProp(args_info.model_arg, iter->second, param, outputDiff, prefix);
140 <  }
155 >
156    
157 +  std::map<std::string, SDShape>::iterator si;
158 +  for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
159 +      HydrodynamicsModel* model;
160 +      Shape* shape = si->second.shape;
161 +      StuntDouble* sd = si->second.sd;;
162 +      if (args_info.model_given) {  
163 +        model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
164 +      } else if (shape->hasAnalyticalSolution()) {
165 +        model = new AnalyticalModel(sd, info);
166 +      } else {
167 +        model = new BeadModel(sd, info);
168 +      }
169 +
170 +        model->init();
171 +        
172 +        std::ofstream ofs;
173 +        std::stringstream outputBeads;
174 +        outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
175 +        ofs.open(outputBeads.str().c_str());        
176 +        model->writeBeads(ofs);
177 +        ofs.close();
178 +
179 +        //if beads option is turned on, skip the calculation
180 +        if (!args_info.beads_flag) {
181 +          model->calcHydroProps(shape, viscosity, temperature);
182 +          std::ofstream outputDiff;
183 +          outputDiff.open(outputFilename.c_str());
184 +          model->writeHydroProps(outputDiff);
185 +          outputDiff.close();
186 +        }
187 +        
188 +        delete model;
189 +  }
190 +
191 +
192 +  //MemoryUtils::deletePointers(shapes);
193    delete info;
194    
195   }
# Line 146 | Line 197 | void registerHydrodynamicsModels() {
197   void registerHydrodynamicsModels() {
198      HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
199      HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
200 +    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
201  
202   }
151 –
152 – bool calcHydrodynamicsProp(const std::string& modelType, StuntDouble* sd, const DynamicProperty& param, std::ostream& os, const std::string& prefix) {
153 –    HydrodynamicsModel* hydroModel = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(modelType, sd, param);
154 –    bool ret = false;
155 –    if (hydroModel == NULL) {
156 –        std::cout << "Integrator Factory can not create " << modelType <<std::endl;
157 –    }    
158 –
159 –    if (hydroModel->calcHydrodyanmicsProps()) {
160 –        ret = true;        
161 –        hydroModel->writeDiffCenterAndDiffTensor(os);
162 –
163 –        std::ofstream ofs;
164 –        std::stringstream outputBeads;
165 –        outputBeads << prefix << "_" << sd->getType() << ".xyz";
166 –        ofs.open(outputBeads.str().c_str());
167 –        hydroModel->writeBeads(ofs);
168 –        ofs.close();
169 –    }
170 –
171 –    delete hydroModel;
172 –
173 –    return ret;    
174 – }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines