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, 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(); |
95 |
|
if (args_info.output_given){ |
96 |
|
prefix = args_info.output_arg; |
97 |
|
} else { |
98 |
< |
prefix = "hydro"; |
98 |
> |
prefix = getPrefix(mdFileName); |
99 |
|
} |
100 |
|
std::string outputFilename = prefix + ".diff"; |
101 |
< |
|
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 |
< |
|
101 |
> |
|
102 |
|
//parse md file and set up the system |
103 |
|
SimCreator creator; |
104 |
|
SimInfo* info = creator.createSim(mdFileName, true); |
111 |
|
identMat(0,0) = 1.0; |
112 |
|
identMat(1,1) = 1.0; |
113 |
|
identMat(2,2) = 1.0; |
114 |
< |
|
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 |
> |
|
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, StuntDouble*> uniqueStuntDoubles; |
135 |
> |
std::map<std::string, SDShape> uniqueStuntDoubles; |
136 |
|
|
137 |
|
for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
138 |
|
for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
139 |
|
integrableObject = mol->nextIntegrableObject(ii)) { |
140 |
|
if (uniqueStuntDoubles.find(integrableObject->getType()) == uniqueStuntDoubles.end()) { |
141 |
< |
uniqueStuntDoubles.insert(std::map<std::string, StuntDouble*>::value_type(integrableObject->getType(), integrableObject)); |
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()) { |
153 |
|
} |
154 |
|
} |
155 |
|
|
156 |
< |
std::map<std::string, StuntDouble*>::iterator iter; |
156 |
> |
|
157 |
> |
|
158 |
|
std::ofstream outputDiff(outputFilename.c_str()); |
159 |
< |
for (iter = uniqueStuntDoubles.begin(); iter != uniqueStuntDoubles.end(); ++iter) { |
160 |
< |
calcHydrodynamicsProp(args_info.model_arg, iter->second, param, outputDiff, prefix); |
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 |
< |
|
174 |
> |
|
175 |
> |
|
176 |
> |
//MemoryUtils::deletePointers(shapes); |
177 |
|
delete info; |
178 |
|
|
179 |
|
} |
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, StuntDouble* sd, const DynamicProperty& param, std::ostream& os, const std::string& prefix) { |
190 |
< |
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 |
< |
|
189 |
> |
shape->calcHydroProps(model, viscosity, temperature); |
190 |
> |
model->writeHydroProps(os); |
191 |
|
std::ofstream ofs; |
192 |
|
std::stringstream outputBeads; |
193 |
< |
outputBeads << prefix << "_" << sd->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(); |
169 |
– |
} |
170 |
– |
|
171 |
– |
delete hydroModel; |
172 |
– |
|
173 |
– |
return ret; |
197 |
|
} |