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

Go to the documentation of this file.

// -*- C++ -*-
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 
//                                   Fehmi Cirak
//                        California Institute of Technology
//                           (C) 2004 All Rights Reserved
//
// <LicenseText>
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#ifndef SHELLMANAGERSPECIFIC_H
#define SHELLMANAGERSPECIFIC_H
#include <vector>
#include <fstream>
#include <string>
#include <limits>

#include "mpi.h"
#include "elc/LagrangianComm.h"

#include "shells/fragmentation/ShellManagerFragmented.h"
#include "shells/parallel/ShellManagerParallel.h"
#include "shells/utilities/MakeUniqueName.h"
#include "shells/utilities/PropertiesParser.h"
#include "shells/driverCC/SVertexFunctors.h"

#include "BasicSensorSpecific.h"          
#include "ShellEnforceBCFunctor.h"


typedef shells::ShellManagerParallel<shells::ShellManagerFragmented>  SMPF; 


class ShellManagerSpecific : public SMPF {  
private:
    typedef specific::BasicSensor<shells::SVertexDisplacement> BasicSensorSp;
    
    void applyTorsionToDisplacements();
  
public:
    ShellManagerSpecific(const std::string& controlFileName, MPI_Comm solidComm,
                         int numFluidNodes, int firstFluidNode) : 
        SMPF(controlFileName, solidComm), SubIterations(1),
        _torsionAngle(0.0), _tubeLength(std::numeric_limits<double>::max()) {
        
#ifdef DEBUG_PRINT
        int myRank;
        std::ostringstream obuf;
        MPI_Comm_rank(solidComm, &myRank);
        obuf << "S" << myRank << ".log" << static_cast<char>(0);
        oflog.open(obuf.str().c_str());     
        _olog = new std::ostream(oflog.rdbuf());
#endif

        computeMassPrepareAdvance();
        
#ifdef DEBUG_PRINT_ELC
        (*_olog << "*** LagrangianComm: " << numFluidNodes << "  " << firstFluidNode ).flush();
#endif  
        // instantiate an elc object for data exchange
        _elcLag = new elc::LagrangianComm<DIM, double>(MPI_COMM_WORLD, solidComm, numFluidNodes,
                                                       firstFluidNode, elc::GlobalIdentifiers);
#ifdef DEBUG_PRINT_ELC
        (*_olog << " created.\n").flush();
#endif  

        // read data from input file
        utilities::PropertiesParser *prop = new utilities::PropertiesParser;
        prop->registerPropertiesVar("SubIterations", SubIterations);
        prop->registerPropertiesVar("torsionAngle", _torsionAngle);
        prop->registerPropertiesVar("tubeLength", _tubeLength);

        std::ifstream inputFile("./shellInput.dat");
        if (!inputFile.is_open()) assert(false);
        prop->readValues(inputFile);
        
        delete prop;
        inputFile.close();

        applyTorsionToDisplacements();
        
#ifdef DEBUG_PRINT
        (*_olog << "*** ShellManagerSpecific created.\n").flush();
#endif

        if (SubIterations<=0) 
          SubIterations=1;
    }


    ~ShellManagerSpecific() {
//      if (_sensor!=NULL) delete _sensor;
        if (_elcLag!=NULL) delete _elcLag;
        if (_olog!=NULL) delete _olog;
        if (_gnuplotFile!=NULL) {
            _gnuplotFile->close();
            delete _gnuplotFile;
        }
    }

   
    void sendBoundaryReceivePressure() {
        double *elCoordinates = NULL;
        double *elVelocities = NULL;
        int *elGlobalNodeIDs = NULL;
        int elNumNodes;
        int *elConnectivity = NULL;
        int elNumElements;
        
        SMPF::decode(&elCoordinates, &elVelocities, &elGlobalNodeIDs,
                     &elNumNodes, &elConnectivity, &elNumElements);
        
        // allocate space for pressure and initialize
        _elPressures.resize(elNumElements);       
        std::fill(_elPressures.begin(), _elPressures.end(), 0.0);
        
#ifdef DEBUG_PRINT_ELC
        ( *_olog << "*** ShellManagerSpecific::sendBoundaryReceivePressure" ).flush();
#endif
        
        _elcLag->sendMesh(elNumNodes, reinterpret_cast<void*>(elGlobalNodeIDs), 
                          reinterpret_cast<void*>(elCoordinates), 
                          reinterpret_cast<void*>(elVelocities), 
                          elNumElements, reinterpret_cast<void*>(elConnectivity));
        _elcLag->waitForMesh();
        _elcLag->receivePressure(elNumElements, reinterpret_cast<void*>(&(_elPressures[0])));
        _elcLag->waitForPressure();

        // Necessary if solid domain larger than fluid domain
        for (unsigned int n=0; n<_elPressures.size(); n++) {
            if (_elPressures[n] == std::numeric_limits<double>::max()) {
                _elPressures[n] = 0.0;
            }
        }

        encodePressure(&(_elPressures[0]), _elPressures.size(), element);
        
#ifdef DEBUG_PRINT_ELC
        ( *_olog << " done.\n" ).flush();
#endif
    }


    virtual void advanceSp(double& t, double& dt){
        dt /= SubIterations;
        setTimeStep(dt);

#ifdef DEBUG_PRINT
        (*_olog << "*** advancing in time.\n").flush();
#endif
        // boundary conditions functor
        shells::ShellEnforceBCFunctor enforceBC(0.0, _tubeLength, _torsionAngle); 

        mShell()->iterateOverVertices(enforceBC);
        predictAndEnforceBC();

        sendBoundaryReceivePressure();
        internalExternalForces();       

        correct();      
        incrementCurrentTimeAndStep();

        for (int n=1; n<SubIterations; n++) {
            
          mShell()->iterateOverVertices(enforceBC);
          predictAndEnforceBC();         

          internalExternalForces();     

          correct();    
          incrementCurrentTimeAndStep();
        }

        dt = stableTimeStep()*SubIterations;
        t = getCurrentTime();
        
        // print sensors
        const int stepNum  = getCurrentStepNum();
        std::for_each(_sensors.begin(), _sensors.end(), 
                      std::bind2nd(std::mem_fun(&BasicSensorSp::printData), t));
    }


    void addSensors(double xyz[3], bool restart=false) {
        // open a file and instantiate a sensor
        const int myRank = communicatorRank();
        std::string nameGnuplot = utilities::makeUniqueName(_sensors.size(), myRank, "./dispSensor-", "dat");

        std::ofstream *gnuplotFile;
        if (restart) {
            gnuplotFile = new std::ofstream(nameGnuplot.c_str(), std::ofstream::app);
        } else {
            gnuplotFile = new std::ofstream(nameGnuplot.c_str());
        }
        BasicSensorSp *s = new BasicSensorSp(*gnuplotFile, ShellManagerFragmented::mShell(), xyz);
        
        _sensors.push_back(s);
    }


    void instrumentWithSensors(bool restart=false) {
        const double elementLength = 0.002528;
//      const double eps = elementLength;
        const double eps = 1.0e-8;
        const double radius = 0.020195;
        const double begin = 0.1969;
        const double deltaZ = 2.*elementLength;

        double xyz0L[3] = { eps, radius, begin-deltaZ};
        double xyz0R[3] = {-eps, radius, begin-deltaZ};
        addSensors(xyz0L, restart);
        addSensors(xyz0R, restart);

        double xyz1L[3] = { eps, radius, begin-2.*deltaZ};
        double xyz1R[3] = {-eps, radius, begin-2.*deltaZ};
        addSensors(xyz1L, restart);
        addSensors(xyz1R, restart);

        double xyz2L[3] = { eps, radius, begin-3.*deltaZ};
        double xyz2R[3] = {-eps, radius, begin-3.*deltaZ};
        addSensors(xyz2L, restart);
        addSensors(xyz2R, restart);

        const double begin2 = 0.2601;
        double xyz3L[3] = { eps, radius, begin2+deltaZ};
        double xyz3R[3] = {-eps, radius, begin2+deltaZ};
        addSensors(xyz3L, restart);
        addSensors(xyz3R, restart);

        double xyz4L[3] = { eps, radius, begin2+2.*deltaZ};
        double xyz4R[3] = {-eps, radius, begin2+2.*deltaZ};
        addSensors(xyz4L, restart);
        addSensors(xyz4R, restart);
        
        double xyz5L[3] = { eps, radius, begin2+3.*deltaZ};
        double xyz5R[3] = {-eps, radius, begin2+3.*deltaZ};
        addSensors(xyz5L, restart);
        addSensors(xyz5R, restart);     
    }


         
    void initialize(double& t, double& dt){
        sendBoundaryReceivePressure();
        dt = stableTimeStep()*SubIterations;
        t = getCurrentTime();
        
        instrumentWithSensors(false);   
        
        return;
    }

 
    void output() {
        printDataParallel(true);

#ifdef DEBUG_PRINT
        ( *_olog << " Printing interface mesh and pressure.\n" ).flush();
        printIFaceMeshPressureParallel();
#endif
    }


    int nSteps() {
        return getCurrentStepNum()/SubIterations;
    }


    void checkpointing() {
        checkPointingParallel();
    }

    void restartSp(double& t, double& dt) {
        restartParallel();
        t = getCurrentTime();
        dt = stableTimeStep()*SubIterations;

        shells::ShellEnforceBCFunctor enforceBC(0.0, _tubeLength, _torsionAngle);       
        mShell()->iterateOverVertices(enforceBC);

        sendBoundaryReceivePressure();

        instrumentWithSensors(true);

        return;
    }


// copy and assignment constructors - not implemented
private:
    ShellManagerSpecific(const ShellManagerSpecific &);
    const ShellManagerSpecific & operator=(const ShellManagerSpecific &);   

private:
    elc::LagrangianComm<DIM, double>     *_elcLag;

    std::ostream                         *_olog;
    std::ofstream                         oflog; 

    std::ofstream                        *_gnuplotFile;
        
    std::vector<double>                  _elPressures;

    std::vector<BasicSensorSp *>         _sensors;

    int                                   SubIterations;

    double                               _torsionAngle;
    double                               _tubeLength; 
};


// implementation
void ShellManagerSpecific::applyTorsionToDisplacements() 
{
    typedef std::vector<double> DCont;
    typedef std::back_insert_iterator<DCont> DContIns;
    
    // collect vertex coordinates
    DCont coor;
    shells::SVertexCollector<DContIns, shells::SVertexCoordinate> getCoor;
    mShell()->iterateOverVertices(std::bind2nd(getCoor, 
                                               std::back_inserter(coor)));              
    
    // modify displacements
    // we assume that the smallest z coordinate is zero
    // and the largest z coordinate is _tubeLength
    DCont disp;
    const double pi = 3.14159265358979323846; 
    const unsigned numNodes = static_cast<unsigned>(coor.size()/3);
    for (unsigned i=0; i<numNodes; ++i) {
        const double xcoor = coor[i*3];
        const double ycoor = coor[i*3+1];
        const double zcoor = coor[i*3+2];
        
        const double rad = std::sqrt(xcoor*xcoor+ycoor*ycoor);
        
        double phi = std::atan2(ycoor, xcoor);
        if (ycoor<0) phi =  2.0*pi+phi;

        double angle = _torsionAngle/_tubeLength*zcoor;
        phi += angle;
        const double xcoorN = rad*std::cos(phi);
        const double ycoorN = rad*std::sin(phi);
        
        disp.push_back(xcoorN-xcoor);
        disp.push_back(ycoorN-ycoor);
        disp.push_back(0.0);
    }
    
    assert(coor.size()==disp.size());
    
    // distribute vertex displacements
    shells::SVertexDistributor<DCont::iterator, shells::SVertexDisplacement> 
        setDisp(disp.begin(), disp.end());    
    mShell()->iterateOverVertices(setDisp);
    
    return;
}

#endif

---