vtf-logo

fsi/sfc-amroc/TubeCJBurnFrac/src/ShellManagerSpecific.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00004 // 
00005 //                                   Fehmi Cirak
00006 //                        California Institute of Technology
00007 //                           (C) 2004 All Rights Reserved
00008 //
00009 // <LicenseText>
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/fragmentation/ShellManagerFragmented.h"
00024 #include "shells/parallel/ShellManagerParallel.h"
00025 #include "shells/utilities/MakeUniqueName.h"
00026 #include "shells/utilities/PropertiesParser.h"
00027 #include "shells/driverCC/SVertexFunctors.h"
00028 
00029 #include "BasicSensorSpecific.h"          
00030 #include "ShellEnforceBCFunctor.h"
00031 
00032 
00033 typedef shells::ShellManagerParallel<shells::ShellManagerFragmented>  SMPF; 
00034 
00035 
00036 class ShellManagerSpecific : public SMPF {  
00037 private:
00038     typedef specific::BasicSensor<shells::SVertexDisplacement> BasicSensorSp;
00039     
00040     void applyTorsionToDisplacements();
00041   
00042 public:
00043     ShellManagerSpecific(const std::string& controlFileName, MPI_Comm solidComm,
00044                          int numFluidNodes, int firstFluidNode) : 
00045         SMPF(controlFileName, solidComm), SubIterations(1),
00046         _torsionAngle(0.0), _tubeLength(std::numeric_limits<double>::max()) {
00047         
00048 #ifdef DEBUG_PRINT
00049         int myRank;
00050         std::ostringstream obuf;
00051         MPI_Comm_rank(solidComm, &myRank);
00052         obuf << "S" << myRank << ".log" << static_cast<char>(0);
00053         oflog.open(obuf.str().c_str());     
00054         _olog = new std::ostream(oflog.rdbuf());
00055 #endif
00056 
00057         computeMassPrepareAdvance();
00058         
00059 #ifdef DEBUG_PRINT_ELC
00060         (*_olog << "*** LagrangianComm: " << numFluidNodes << "  " << firstFluidNode ).flush();
00061 #endif  
00062         // instantiate an elc object for data exchange
00063         _elcLag = new elc::LagrangianComm<DIM, double>(MPI_COMM_WORLD, solidComm, numFluidNodes,
00064                                                        firstFluidNode, elc::GlobalIdentifiers);
00065 #ifdef DEBUG_PRINT_ELC
00066         (*_olog << " created.\n").flush();
00067 #endif  
00068 
00069         // read data from input file
00070         utilities::PropertiesParser *prop = new utilities::PropertiesParser;
00071         prop->registerPropertiesVar("SubIterations", SubIterations);
00072         prop->registerPropertiesVar("torsionAngle", _torsionAngle);
00073         prop->registerPropertiesVar("tubeLength", _tubeLength);
00074 
00075         std::ifstream inputFile("./shellInput.dat");
00076         if (!inputFile.is_open()) assert(false);
00077         prop->readValues(inputFile);
00078         
00079         delete prop;
00080         inputFile.close();
00081 
00082         applyTorsionToDisplacements();
00083         
00084 #ifdef DEBUG_PRINT
00085         (*_olog << "*** ShellManagerSpecific created.\n").flush();
00086 #endif
00087 
00088         if (SubIterations<=0) 
00089           SubIterations=1;
00090     }
00091 
00092 
00093     ~ShellManagerSpecific() {
00094 //      if (_sensor!=NULL) delete _sensor;
00095         if (_elcLag!=NULL) delete _elcLag;
00096         if (_olog!=NULL) delete _olog;
00097         if (_gnuplotFile!=NULL) {
00098             _gnuplotFile->close();
00099             delete _gnuplotFile;
00100         }
00101     }
00102 
00103    
00104     void sendBoundaryReceivePressure() {
00105         double *elCoordinates = NULL;
00106         double *elVelocities = NULL;
00107         int *elGlobalNodeIDs = NULL;
00108         int elNumNodes;
00109         int *elConnectivity = NULL;
00110         int elNumElements;
00111         
00112         SMPF::decode(&elCoordinates, &elVelocities, &elGlobalNodeIDs,
00113                      &elNumNodes, &elConnectivity, &elNumElements);
00114         
00115         // allocate space for pressure and initialize
00116         _elPressures.resize(elNumElements);       
00117         std::fill(_elPressures.begin(), _elPressures.end(), 0.0);
00118         
00119 #ifdef DEBUG_PRINT_ELC
00120         ( *_olog << "*** ShellManagerSpecific::sendBoundaryReceivePressure" ).flush();
00121 #endif
00122         
00123         _elcLag->sendMesh(elNumNodes, reinterpret_cast<void*>(elGlobalNodeIDs), 
00124                           reinterpret_cast<void*>(elCoordinates), 
00125                           reinterpret_cast<void*>(elVelocities), 
00126                           elNumElements, reinterpret_cast<void*>(elConnectivity));
00127         _elcLag->waitForMesh();
00128         _elcLag->receivePressure(elNumElements, reinterpret_cast<void*>(&(_elPressures[0])));
00129         _elcLag->waitForPressure();
00130 
00131         // Necessary if solid domain larger than fluid domain
00132         for (unsigned int n=0; n<_elPressures.size(); n++) {
00133             if (_elPressures[n] == std::numeric_limits<double>::max()) {
00134                 _elPressures[n] = 0.0;
00135             }
00136         }
00137 
00138         encodePressure(&(_elPressures[0]), _elPressures.size(), element);
00139         
00140 #ifdef DEBUG_PRINT_ELC
00141         ( *_olog << " done.\n" ).flush();
00142 #endif
00143     }
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         // boundary conditions functor
00154         shells::ShellEnforceBCFunctor enforceBC(0.0, _tubeLength, _torsionAngle); 
00155 
00156         mShell()->iterateOverVertices(enforceBC);
00157         predictAndEnforceBC();
00158 
00159         sendBoundaryReceivePressure();
00160         internalExternalForces();       
00161 
00162         correct();      
00163         incrementCurrentTimeAndStep();
00164 
00165         for (int n=1; n<SubIterations; n++) {
00166             
00167           mShell()->iterateOverVertices(enforceBC);
00168           predictAndEnforceBC();         
00169 
00170           internalExternalForces();     
00171 
00172           correct();    
00173           incrementCurrentTimeAndStep();
00174         }
00175 
00176         dt = stableTimeStep()*SubIterations;
00177         t = getCurrentTime();
00178         
00179         // print sensors
00180         const int stepNum  = getCurrentStepNum();
00181         std::for_each(_sensors.begin(), _sensors.end(), 
00182                       std::bind2nd(std::mem_fun(&BasicSensorSp::printData), t));
00183     }
00184 
00185 
00186     void addSensors(double xyz[3], bool restart=false) {
00187         // open a file and instantiate a sensor
00188         const int myRank = communicatorRank();
00189         std::string nameGnuplot = utilities::makeUniqueName(_sensors.size(), myRank, "./dispSensor-", "dat");
00190 
00191         std::ofstream *gnuplotFile;
00192         if (restart) {
00193             gnuplotFile = new std::ofstream(nameGnuplot.c_str(), std::ofstream::app);
00194         } else {
00195             gnuplotFile = new std::ofstream(nameGnuplot.c_str());
00196         }
00197         BasicSensorSp *s = new BasicSensorSp(*gnuplotFile, ShellManagerFragmented::mShell(), xyz);
00198         
00199         _sensors.push_back(s);
00200     }
00201 
00202 
00203     void instrumentWithSensors(bool restart=false) {
00204         const double elementLength = 0.002528;
00205 //      const double eps = elementLength;
00206         const double eps = 1.0e-8;
00207         const double radius = 0.020195;
00208         const double begin = 0.1969;
00209         const double deltaZ = 2.*elementLength;
00210 
00211         double xyz0L[3] = { eps, radius, begin-deltaZ};
00212         double xyz0R[3] = {-eps, radius, begin-deltaZ};
00213         addSensors(xyz0L, restart);
00214         addSensors(xyz0R, restart);
00215 
00216         double xyz1L[3] = { eps, radius, begin-2.*deltaZ};
00217         double xyz1R[3] = {-eps, radius, begin-2.*deltaZ};
00218         addSensors(xyz1L, restart);
00219         addSensors(xyz1R, restart);
00220 
00221         double xyz2L[3] = { eps, radius, begin-3.*deltaZ};
00222         double xyz2R[3] = {-eps, radius, begin-3.*deltaZ};
00223         addSensors(xyz2L, restart);
00224         addSensors(xyz2R, restart);
00225 
00226         const double begin2 = 0.2601;
00227         double xyz3L[3] = { eps, radius, begin2+deltaZ};
00228         double xyz3R[3] = {-eps, radius, begin2+deltaZ};
00229         addSensors(xyz3L, restart);
00230         addSensors(xyz3R, restart);
00231 
00232         double xyz4L[3] = { eps, radius, begin2+2.*deltaZ};
00233         double xyz4R[3] = {-eps, radius, begin2+2.*deltaZ};
00234         addSensors(xyz4L, restart);
00235         addSensors(xyz4R, restart);
00236         
00237         double xyz5L[3] = { eps, radius, begin2+3.*deltaZ};
00238         double xyz5R[3] = {-eps, radius, begin2+3.*deltaZ};
00239         addSensors(xyz5L, restart);
00240         addSensors(xyz5R, restart);     
00241     }
00242 
00243 
00244          
00245     void initialize(double& t, double& dt){
00246         sendBoundaryReceivePressure();
00247         dt = stableTimeStep()*SubIterations;
00248         t = getCurrentTime();
00249         
00250         instrumentWithSensors(false);   
00251         
00252         return;
00253     }
00254 
00255  
00256     void output() {
00257         printDataParallel(true);
00258 
00259 #ifdef DEBUG_PRINT
00260         ( *_olog << " Printing interface mesh and pressure.\n" ).flush();
00261         printIFaceMeshPressureParallel();
00262 #endif
00263     }
00264 
00265 
00266     int nSteps() {
00267         return getCurrentStepNum()/SubIterations;
00268     }
00269 
00270 
00271     void checkpointing() {
00272         checkPointingParallel();
00273     }
00274 
00275     void restartSp(double& t, double& dt) {
00276         restartParallel();
00277         t = getCurrentTime();
00278         dt = stableTimeStep()*SubIterations;
00279 
00280         shells::ShellEnforceBCFunctor enforceBC(0.0, _tubeLength, _torsionAngle);       
00281         mShell()->iterateOverVertices(enforceBC);
00282 
00283         sendBoundaryReceivePressure();
00284 
00285         instrumentWithSensors(true);
00286 
00287         return;
00288     }
00289 
00290 
00291 // copy and assignment constructors - not implemented
00292 private:
00293     ShellManagerSpecific(const ShellManagerSpecific &);
00294     const ShellManagerSpecific & operator=(const ShellManagerSpecific &);   
00295 
00296 private:
00297     elc::LagrangianComm<DIM, double>     *_elcLag;
00298 
00299     std::ostream                         *_olog;
00300     std::ofstream                         oflog; 
00301 
00302     std::ofstream                        *_gnuplotFile;
00303         
00304     std::vector<double>                  _elPressures;
00305 
00306     std::vector<BasicSensorSp *>         _sensors;
00307 
00308     int                                   SubIterations;
00309 
00310     double                               _torsionAngle;
00311     double                               _tubeLength; 
00312 };
00313 
00314 
00315 // implementation
00316 void ShellManagerSpecific::applyTorsionToDisplacements() 
00317 {
00318     typedef std::vector<double> DCont;
00319     typedef std::back_insert_iterator<DCont> DContIns;
00320     
00321     // collect vertex coordinates
00322     DCont coor;
00323     shells::SVertexCollector<DContIns, shells::SVertexCoordinate> getCoor;
00324     mShell()->iterateOverVertices(std::bind2nd(getCoor, 
00325                                                std::back_inserter(coor)));              
00326     
00327     // modify displacements
00328     // we assume that the smallest z coordinate is zero
00329     // and the largest z coordinate is _tubeLength
00330     DCont disp;
00331     const double pi = 3.14159265358979323846; 
00332     const unsigned numNodes = static_cast<unsigned>(coor.size()/3);
00333     for (unsigned i=0; i<numNodes; ++i) {
00334         const double xcoor = coor[i*3];
00335         const double ycoor = coor[i*3+1];
00336         const double zcoor = coor[i*3+2];
00337         
00338         const double rad = std::sqrt(xcoor*xcoor+ycoor*ycoor);
00339         
00340         double phi = std::atan2(ycoor, xcoor);
00341         if (ycoor<0) phi =  2.0*pi+phi;
00342 
00343         double angle = _torsionAngle/_tubeLength*zcoor;
00344         phi += angle;
00345         const double xcoorN = rad*std::cos(phi);
00346         const double ycoorN = rad*std::sin(phi);
00347         
00348         disp.push_back(xcoorN-xcoor);
00349         disp.push_back(ycoorN-ycoor);
00350         disp.push_back(0.0);
00351     }
00352     
00353     assert(coor.size()==disp.size());
00354     
00355     // distribute vertex displacements
00356     shells::SVertexDistributor<DCont::iterator, shells::SVertexDisplacement> 
00357         setDisp(disp.begin(), disp.end());    
00358     mShell()->iterateOverVertices(setDisp);
00359     
00360     return;
00361 }
00362 
00363 #endif

Generated on Fri Aug 24 13:02:34 2007 for Virtual Test Facility Coupled Applications by  doxygen 1.4.7