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 891 by tim, Wed Feb 22 20:35:16 2006 UTC vs.
Revision 907 by tim, Fri Mar 17 23:51:54 2006 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 calcHydrodynamicsProp(HydrodynamicsModel* model, Shape* shape,double viscosity,  double temperature, std::ostream& os, const std::string& prefix);
68  
61 – bool calcHydrodynamicsProp(const std::string& modelType, Molecule* mol, const DynamicProperty& param, const std::string& prefix,double viscosity);
62 –
69   int main(int argc, char* argv[]){
70    //register force fields
71    registerForceFields();    
# Line 89 | Line 95 | int main(int argc, char* argv[]){
95    if (args_info.output_given){
96      prefix = args_info.output_arg;
97    } else {
98 <    prefix = "hydro";    
98 >    prefix = getPrefix(mdFileName);    
99    }
100 <
101 <  DynamicProperty param;
96 <  if (args_info.sigma_given) {
97 <    param.insert(DynamicProperty::value_type("Sigma", args_info.sigma_arg));
98 <  }
99 <  
100 <  
100 >  std::string outputFilename = prefix + ".diff";
101 >    
102    //parse md file and set up the system
103    SimCreator creator;
104    SimInfo* info = creator.createSim(mdFileName, true);
105      
106    SimInfo::MoleculeIterator mi;
107    Molecule* mol;
108 <  Molecule::RigidBodyIterator ri;
109 <  RigidBody* rb;
110 <    //update atoms of rigidbody
111 <  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
112 <      
113 <      //change the positions of atoms which belong to the rigidbodies
114 <      for (rb = mol->beginRigidBody(ri); rb != NULL; rb = mol->nextRigidBody(ri)) {
115 <        rb->updateAtoms();
116 <      }
108 >  Molecule::IntegrableObjectIterator  ii;
109 >  StuntDouble* integrableObject;
110 >  Mat3x3d identMat;
111 >  identMat(0,0) = 1.0;
112 >  identMat(1,1) = 1.0;
113 >  identMat(2,2) = 1.0;
114 >
115 >  Globals* simParams = info->getSimParams();
116 >  double temperature;
117 >  double viscosity;
118 >
119 >  if (simParams->haveViscosity()) {
120 >    viscosity = simParams->getViscosity();
121 >  } else {
122 >    sprintf(painCave.errMsg, "viscosity must be set\n");
123 >    painCave.isFatal = 1;
124 >    simError();  
125    }
126 <    
126 >
127 >  if (simParams->haveTargetTemp()) {
128 >    temperature = simParams->getTargetTemp();
129 >  } else {
130 >    sprintf(painCave.errMsg, "target temperaturemust be set\n");
131 >    painCave.isFatal = 1;
132 >    simError();  
133 >  }
134 >
135 >  std::map<std::string, SDShape> uniqueStuntDoubles;
136 >  
137    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
138 <    calcHydrodynamicsProp(args_info.model_arg, mol, param, prefix, args_info.viscosity_arg);
138 >      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
139 >           integrableObject = mol->nextIntegrableObject(ii)) {
140 >          if (uniqueStuntDoubles.find(integrableObject->getType()) ==  uniqueStuntDoubles.end()) {
141 >
142 >            SDShape tmp;
143 >            tmp.shape = ShapeBuilder::createShape(integrableObject);
144 >            tmp.sd = integrableObject;    
145 >            uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));
146 >            integrableObject->setPos(V3Zero);
147 >            integrableObject->setA(identMat);
148 >            if (integrableObject->isRigidBody()) {
149 >                RigidBody* rb = static_cast<RigidBody*>(integrableObject);
150 >                rb->updateAtoms();
151 >            }
152 >          }
153 >        }
154    }
155 +
156 +
157    
158 +  std::ofstream outputDiff(outputFilename.c_str());
159 +  std::map<std::string, SDShape>::iterator si;
160 +  for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
161 +      HydrodynamicsModel* model;
162 +      Shape* shape = si->second.shape;
163 +      StuntDouble* sd = si->second.sd;;
164 +      if (args_info.model_given) {  
165 +        model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
166 +      } else if (shape->hasAnalyticalSolution()) {
167 +        model = new AnalyticalModel(sd, info);
168 +      } else {
169 +        model = new BeadModel(sd, info);
170 +      }
171 +        calcHydrodynamicsProp(model, shape, viscosity, temperature, outputDiff, prefix);          
172 +        delete model;
173 +  }
174 +
175 +
176 +  //MemoryUtils::deletePointers(shapes);
177    delete info;
178    
179   }
# Line 126 | Line 181 | void registerHydrodynamicsModels() {
181   void registerHydrodynamicsModels() {
182      HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
183      HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
184 +    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
185  
186   }
187 + void calcHydrodynamicsProp(HydrodynamicsModel* model, Shape* shape,double viscosity,  double temperature, std::ostream& os, const std::string& prefix) {
188  
189 < bool calcHydrodynamicsProp(const std::string& modelType, Molecule* mol, const DynamicProperty& param, const std::string& prefix,double viscosity) {
190 <    HydrodynamicsModel* hydroModel = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(modelType, mol, param);
134 <    bool ret = false;
135 <    if (hydroModel == NULL) {
136 <        std::cout << "Integrator Factory can not create " << modelType <<std::endl;
137 <    }    
138 <
139 <    if (hydroModel->calcHydrodyanmicsProps(viscosity)) {
140 <        ret = true;
141 <        std::stringstream outputDiffTensor;
142 <        outputDiffTensor << prefix << "_" << mol->getType() << ".diff";
189 >       shape->calcHydroProps(model, viscosity, temperature);
190 >        model->writeHydroProps(os);
191          std::ofstream ofs;
144 –        ofs.open(outputDiffTensor.str().c_str());
145 –        hydroModel->writeDiffCenterAndDiffTensor(ofs);
146 –        ofs.close();
147 –        
192          std::stringstream outputBeads;
193 <        outputBeads << prefix << "_" << mol->getType() << ".xyz";
194 <        ofs.open(outputBeads.str().c_str());
195 <        hydroModel->writeBeads(ofs);
193 >        outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
194 >        ofs.open(outputBeads.str().c_str());        
195 >        model->writeBeads(ofs);
196          ofs.close();
153 –    }
154 –
155 –    delete hydroModel;
156 –
157 –    return ret;    
197   }

Diff Legend

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