--- trunk/src/applications/hydrodynamics/Hydro.cpp	2006/02/22 20:35:16	891
+++ trunk/src/applications/hydrodynamics/Hydro.cpp	2006/08/30 21:13:57	1027
@@ -47,19 +47,25 @@
 #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
 #include "applications/hydrodynamics/HydrodynamicsModelCreator.hpp"
 #include "applications/hydrodynamics/HydrodynamicsModelFactory.hpp"
+#include "applications/hydrodynamics/AnalyticalModel.hpp"
 #include "applications/hydrodynamics/BeadModel.hpp"
 #include "applications/hydrodynamics/RoughShell.hpp"
+#include "applications/hydrodynamics/ShapeBuilder.hpp"
 #include "brains/Register.hpp"
 #include "brains/SimCreator.hpp"
 #include "brains/SimInfo.hpp"
-
+#include "utils/StringUtils.hpp"
+#include "utils/simError.h"
+#include "utils/MemoryUtils.hpp"
 using namespace oopse;
 
-/** Register different hydrodynamics models */
+struct SDShape{
+    StuntDouble* sd;
+    Shape* shape;
+};
 void registerHydrodynamicsModels();
+void writeHydroProps(std::ostream& os);
 
-bool calcHydrodynamicsProp(const std::string& modelType, Molecule* mol, const DynamicProperty& param, const std::string& prefix,double viscosity);
-
 int main(int argc, char* argv[]){
   //register force fields
   registerForceFields();    
@@ -83,76 +89,112 @@ int main(int argc, char* argv[]){
     exit(1);
   }
   
-  mdFileName = dumpFileName;
-  mdFileName = mdFileName.substr(0, mdFileName.rfind(".")) + ".md";
-
   if (args_info.output_given){
     prefix = args_info.output_arg;
   } else {
-    prefix = "hydro";    
+    prefix = getPrefix(dumpFileName);    
   }
-
-  DynamicProperty param;
-  if (args_info.sigma_given) {
-    param.insert(DynamicProperty::value_type("Sigma", args_info.sigma_arg));
-  }
-  
-  
+  std::string outputFilename = prefix + ".diff";
+    
   //parse md file and set up the system
   SimCreator creator;
-  SimInfo* info = creator.createSim(mdFileName, true);
+  SimInfo* info = creator.createSim(dumpFileName, true);
     
   SimInfo::MoleculeIterator mi;
   Molecule* mol;
-  Molecule::RigidBodyIterator ri;
-  RigidBody* rb;
-    //update atoms of rigidbody
-  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
-      
-      //change the positions of atoms which belong to the rigidbodies
-      for (rb = mol->beginRigidBody(ri); rb != NULL; rb = mol->nextRigidBody(ri)) {
-        rb->updateAtoms();
-      }
+  Molecule::IntegrableObjectIterator  ii;
+  StuntDouble* integrableObject;
+  Mat3x3d identMat;
+  identMat(0,0) = 1.0;
+  identMat(1,1) = 1.0;
+  identMat(2,2) = 1.0;
+
+  Globals* simParams = info->getSimParams();
+  RealType temperature;
+  RealType viscosity;
+
+  if (simParams->haveViscosity()) {
+    viscosity = simParams->getViscosity();
+  } else {
+    sprintf(painCave.errMsg, "viscosity must be set\n");
+    painCave.isFatal = 1;
+    simError();  
   }
-    
-  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
-    calcHydrodynamicsProp(args_info.model_arg, mol, param, prefix, args_info.viscosity_arg);
+
+  if (simParams->haveTargetTemp()) {
+    temperature = simParams->getTargetTemp();
+  } else {
+    sprintf(painCave.errMsg, "target temperature must be set\n");
+    painCave.isFatal = 1;
+    simError();  
   }
+ 
+  std::map<std::string, SDShape> uniqueStuntDoubles;
   
-  delete info;
-   
-}
+  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
+      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
+	   integrableObject = mol->nextIntegrableObject(ii)) {
+          if (uniqueStuntDoubles.find(integrableObject->getType()) ==  uniqueStuntDoubles.end()) {
 
-void registerHydrodynamicsModels() {
-    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
-    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
+            integrableObject->setPos(V3Zero);
+            integrableObject->setA(identMat);
+            if (integrableObject->isRigidBody()) {
+                RigidBody* rb = static_cast<RigidBody*>(integrableObject);
+                rb->updateAtoms();
+            }
 
-}
+            SDShape tmp;
+            tmp.shape = ShapeBuilder::createShape(integrableObject);
+            tmp.sd = integrableObject;    
+            uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));
 
-bool calcHydrodynamicsProp(const std::string& modelType, Molecule* mol, const DynamicProperty& param, const std::string& prefix,double viscosity) {
-    HydrodynamicsModel* hydroModel = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(modelType, mol, param);
-    bool ret = false;
-    if (hydroModel == NULL) {
-        std::cout << "Integrator Factory can not create " << modelType <<std::endl;
-    }    
+          }
+	}
+  }
 
-    if (hydroModel->calcHydrodyanmicsProps(viscosity)) {
-        ret = true;
-        std::stringstream outputDiffTensor;
-        outputDiffTensor << prefix << "_" << mol->getType() << ".diff";
-        std::ofstream ofs;
-        ofs.open(outputDiffTensor.str().c_str());
-        hydroModel->writeDiffCenterAndDiffTensor(ofs);
-        ofs.close();
+
+  
+  std::ofstream outputDiff(outputFilename.c_str());
+  std::map<std::string, SDShape>::iterator si;
+  for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
+      HydrodynamicsModel* model;
+      Shape* shape = si->second.shape;
+      StuntDouble* sd = si->second.sd;;
+      if (args_info.model_given) {  
+        model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
+      } else if (shape->hasAnalyticalSolution()) {
+        model = new AnalyticalModel(sd, info);
+      } else {
+        model = new BeadModel(sd, info);
+      }
+
+        model->init();
         
+        std::ofstream ofs;
         std::stringstream outputBeads;
-        outputBeads << prefix << "_" << mol->getType() << ".xyz";
-        ofs.open(outputBeads.str().c_str());
-        hydroModel->writeBeads(ofs);
+        outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
+        ofs.open(outputBeads.str().c_str());        
+        model->writeBeads(ofs);
         ofs.close();
-    }
 
-    delete hydroModel;
+        //if beads option is turned on, skip the calculation
+        if (!args_info.beads_flag) {
+          model->calcHydroProps(shape, viscosity, temperature);
+          model->writeHydroProps(outputDiff);
+        }
+        
+        delete model;
+  }
 
-    return ret;    
+
+  //MemoryUtils::deletePointers(shapes);
+  delete info;
+   
 }
+
+void registerHydrodynamicsModels() {
+    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
+    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
+    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
+
+}