00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef SHELLMANAGERSPECIFIC_H
00014 #define SHELLMANAGERSPECIFIC_H
00015 #include <vector>
00016 #include <fstream>
00017 #include <string>
00018 #include <limits>
00019
00020 #include "mpi.h"
00021 #include "elc/LagrangianComm.h"
00022
00023 #include "shells/driverCC/ShellManagerBasic.h"
00024 #include "shells/parallel/ShellManagerParallel.h"
00025 #include "shells/utilities/BasicSensor.h"
00026 #include "shells/utilities/MakeUniqueName.h"
00027 #include "shells/utilities/PropertiesParser.h"
00028 #include "shells/driverCC/SVertexFunctors.h"
00029 #include "shells/driverCC/PrescribeVarFunctor.h"
00030 #include <shells/driverCC/StableTimeStep.h>
00031
00032 #include "ShellEnforceBCFunctor.h"
00033
00034 typedef shells::ShellManagerParallel<shells::ShellManagerBasic> SMPF;
00035
00036
00037 class ShellManagerSpecific : public SMPF {
00038 private:
00039 typedef utilities::BasicSensor<shells::SVertexDisplacement> BasicSensorSp;
00040
00041 public:
00042 ShellManagerSpecific(const std::string& controlFileName, MPI_Comm solidComm,
00043 int numFluidNodes, int firstFluidNode) :
00044 SMPF(controlFileName, solidComm), _timeStepReduction(1.), _subIterations(1),
00045 _sensoroutput(0) {
00046
00047 #ifdef DEBUG_PRINT
00048 int myRank;
00049 std::ostringstream obuf;
00050 MPI_Comm_rank(solidComm, &myRank);
00051 obuf << "S" << myRank << ".log" << static_cast<char>(0);
00052 oflog.open(obuf.str().c_str());
00053 _olog = new std::ostream(oflog.rdbuf());
00054 #endif
00055
00056 computeMassPrepareAdvance();
00057
00058 #ifdef DEBUG_PRINT_ELC
00059 (*_olog << "*** LagrangianComm: " << numFluidNodes << " " << firstFluidNode ).flush();
00060 #endif
00061
00062 _elcLag = new elc::LagrangianComm<DIM, double>(MPI_COMM_WORLD, solidComm, numFluidNodes,
00063 firstFluidNode, elc::GlobalIdentifiers);
00064 #ifdef DEBUG_PRINT_ELC
00065 (*_olog << " created.\n").flush();
00066 #endif
00067
00068
00069 utilities::PropertiesParser *prop = new utilities::PropertiesParser;
00070 prop->registerPropertiesVar("timeStepReduction", _timeStepReduction);
00071 prop->registerPropertiesVar("SubIterations", _subIterations);
00072 prop->registerPropertiesVar("SensorOutput", _sensoroutput);
00073 std::ifstream inputFile("./shellInput.dat");
00074 if (!inputFile.is_open()) assert(false);
00075 prop->readValues(inputFile);
00076 delete prop;
00077 inputFile.close();
00078
00079 #ifdef DEBUG_PRINT
00080 (*_olog << "*** ShellManagerSpecific created.\n").flush();
00081 #endif
00082
00083 if (_subIterations<=0)
00084 _subIterations=1;
00085 }
00086
00087
00088 ~ShellManagerSpecific() {
00089 if (_elcLag!=NULL) delete _elcLag;
00090 if (_olog!=NULL) delete _olog;
00091 if (_gnuplotFile!=NULL) {
00092 _gnuplotFile->close();
00093 delete _gnuplotFile;
00094 }
00095 }
00096
00097
00098 void sendBoundaryReceivePressure() {
00099 double *elCoordinates = NULL;
00100 double *elVelocities = NULL;
00101 int *elGlobalNodeIDs = NULL;
00102 int elNumNodes;
00103 int *elConnectivity = NULL;
00104 int elNumElements;
00105
00106 SMPF::decode(&elCoordinates, &elVelocities, &elGlobalNodeIDs,
00107 &elNumNodes, &elConnectivity, &elNumElements);
00108
00109
00110 _elPressures.resize(elNumElements);
00111 std::fill(_elPressures.begin(), _elPressures.end(), 0.0);
00112
00113 #ifdef DEBUG_PRINT_ELC
00114 ( *_olog << "*** ShellManagerSpecific::sendBoundaryReceivePressure" ).flush();
00115 #endif
00116
00117 _elcLag->sendMesh(elNumNodes, reinterpret_cast<void*>(elGlobalNodeIDs),
00118 reinterpret_cast<void*>(elCoordinates),
00119 reinterpret_cast<void*>(elVelocities),
00120 elNumElements, reinterpret_cast<void*>(elConnectivity));
00121 _elcLag->waitForMesh();
00122 _elcLag->receivePressure(elNumElements, reinterpret_cast<void*>(&(_elPressures[0])));
00123 _elcLag->waitForPressure();
00124
00125
00126 for (unsigned int n=0; n<_elPressures.size(); n++) {
00127 if (_elPressures[n] == std::numeric_limits<double>::max()) {
00128 _elPressures[n] = 0.0;
00129 }
00130 }
00131
00132 encodePressure(&(_elPressures[0]), _elPressures.size(), element);
00133
00134 #ifdef DEBUG_PRINT_ELC
00135 ( *_olog << " done.\n" ).flush();
00136 #endif
00137 }
00138
00139
00140 virtual double stableTimeStep() {
00141 return _timeStepReduction*shells::computeStableTimeStep(SMPF::mShell(),
00142 SMPF::Thickness(),
00143 shells::MShell::active);
00144 }
00145
00146 virtual void advanceSp(double& t, double& dt){
00147 dt /= _subIterations;
00148 setTimeStep(dt);
00149
00150 #ifdef DEBUG_PRINT
00151 (*_olog << "*** advancing in time.\n").flush();
00152 #endif
00153
00154
00155 shells::ShellEnforceBCFunctor enforceBC;
00156 mShell()->iterateOverVertices(enforceBC);
00157
00158 predictAndEnforceBC();
00159 sendBoundaryReceivePressure();
00160 internalExternalForces();
00161 correct();
00162 incrementCurrentTimeAndStep();
00163
00164 for (int n=1; n<_subIterations; n++) {
00165 mShell()->iterateOverVertices(enforceBC);
00166 predictAndEnforceBC();
00167 internalExternalForces();
00168 correct();
00169 incrementCurrentTimeAndStep();
00170 }
00171
00172 dt = stableTimeStep()*_subIterations;
00173 t = getCurrentTime();
00174
00175
00176 if (_sensoroutput)
00177 std::for_each(_sensors.begin(), _sensors.end(),
00178 std::bind2nd(std::mem_fun(&BasicSensorSp::printData), t));
00179 }
00180
00181
00182 void addSensors(double xyz[3], bool restart=false) {
00183
00184 const int myRank = communicatorRank();
00185 std::string nameGnuplot = utilities::makeUniqueName(_sensors.size(), myRank, "./dispSensor-", "dat");
00186
00187 std::ofstream *gnuplotFile;
00188 if (restart) {
00189 gnuplotFile = new std::ofstream(nameGnuplot.c_str(), std::ofstream::app);
00190 } else {
00191 gnuplotFile = new std::ofstream(nameGnuplot.c_str());
00192 }
00193 BasicSensorSp *s = new BasicSensorSp(*gnuplotFile, ShellManagerBasic::mShell(), xyz);
00194
00195 _sensors.push_back(s);
00196 }
00197
00198
00199 void instrumentWithSensors(bool restart=false) {
00200 double xyz[3] = { 0., 0., 0.};
00201 addSensors(xyz, restart);
00202
00203 double r = 0.032;
00204 for (int n=-20; n<=20; n++) {
00205 xyz[1] = r/20.*n;
00206 addSensors(xyz, restart);
00207 }
00208
00209 xyz[1] = 0.;
00210 for (int n=-20; n<=20; n++) {
00211 xyz[2] = r/20.*n;
00212 addSensors(xyz, restart);
00213 }
00214
00215 double pi4 = std::atan(1.);
00216 for (int n=-20; n<=20; n++) {
00217 xyz[1] = r/20.*n*std::cos(pi4);
00218 xyz[2] = r/20.*n*std::sin(pi4);
00219 addSensors(xyz, restart);
00220 }
00221
00222 for (int n=-20; n<=20; n++) {
00223 xyz[1] = r/20.*n*std::cos(-pi4);
00224 xyz[2] = r/20.*n*std::sin(-pi4);
00225 addSensors(xyz, restart);
00226 }
00227 }
00228
00229
00230 void initialize(double& t, double& dt){
00231 sendBoundaryReceivePressure();
00232 dt = stableTimeStep()*_subIterations;
00233 t = getCurrentTime();
00234 if (_sensoroutput)
00235 instrumentWithSensors(false);
00236 }
00237
00238
00239 void output() {
00240 printDataParallel(true);
00241
00242 #ifdef DEBUG_PRINT
00243 ( *_olog << " Printing interface mesh and pressure.\n" ).flush();
00244 printIFaceMeshPressureParallel();
00245 #endif
00246 }
00247
00248
00249 int nSteps() {
00250 return getCurrentStepNum()/_subIterations;
00251 }
00252
00253
00254 void checkpointing() {
00255 checkPointingParallel();
00256 }
00257
00258 void restartSp(double& t, double& dt) {
00259 restartParallel();
00260 t = getCurrentTime();
00261 dt = stableTimeStep()*_subIterations;
00262
00263 shells::ShellEnforceBCFunctor enforceBC;
00264 mShell()->iterateOverVertices(enforceBC);
00265
00266 sendBoundaryReceivePressure();
00267
00268 if (_sensoroutput)
00269 instrumentWithSensors(true);
00270 }
00271
00272
00273
00274 private:
00275 ShellManagerSpecific(const ShellManagerSpecific &);
00276 const ShellManagerSpecific & operator=(const ShellManagerSpecific &);
00277
00278 private:
00279 elc::LagrangianComm<DIM, double> *_elcLag;
00280
00281 std::ostream *_olog;
00282 std::ofstream oflog;
00283
00284 std::ofstream *_gnuplotFile;
00285
00286 std::vector<double> _elPressures;
00287
00288 double _timeStepReduction;
00289
00290 int _subIterations;
00291
00292 int _sensoroutput;
00293
00294 std::vector<BasicSensorSp *> _sensors;
00295 };
00296
00297 #endif