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

Go to the documentation of this file.

// -*- C++ -*-
//
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//  
//                           Ralf Deiterding, Fehmi Cirak 
//                        California Institute of Technology
//                           (C) 2006 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/BasicSensor.h"
#include "shells/utilities/MakeUniqueName.h"
#include "shells/utilities/PropertiesParser.h"
#include "shells/driverCC/SVertexFunctors.h"
#include "shells/driverCC/PrescribeVarFunctor.h"
#include <shells/driverCC/StableTimeStep.h>

#include "ShellEnforceBCFunctor.h"

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


class ShellManagerSpecific : public SMPF {  
private:
    typedef utilities::BasicSensor<shells::SVertexDisplacement> BasicSensorSp;
    
public:
    ShellManagerSpecific(const std::string& controlFileName, MPI_Comm solidComm,
                         int numFluidNodes, int firstFluidNode) : 
        SMPF(controlFileName, solidComm), _timeStepReduction(1.), _subIterations(1), 
        _sensoroutput(0) {
        
#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 _subIterations from input file
        utilities::PropertiesParser *prop = new utilities::PropertiesParser;
        prop->registerPropertiesVar("timeStepReduction", _timeStepReduction);
        prop->registerPropertiesVar("SubIterations", _subIterations);
        prop->registerPropertiesVar("SensorOutput", _sensoroutput);
        std::ifstream inputFile("./shellInput.dat");
        if (!inputFile.is_open()) assert(false);
        prop->readValues(inputFile);    
        delete prop;
        inputFile.close();
        
#ifdef DEBUG_PRINT
        (*_olog << "*** ShellManagerSpecific created.\n").flush();
#endif

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


    ~ShellManagerSpecific() {
        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 double stableTimeStep() {
      return _timeStepReduction*shells::computeStableTimeStep(SMPF::mShell(), 
                                                              SMPF::Thickness(), 
                                                              shells::MShell::active);
    }

    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; 
        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
        if (_sensoroutput) 
          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) {
        double xyz[3] = { 0., 0., 0.};
        addSensors(xyz, restart);

        double r = 0.032;
        for (int n=-20; n<=20; n++) {
          xyz[1] = r/20.*n;
          addSensors(xyz, restart);
        }

        xyz[1] = 0.;
        for (int n=-20; n<=20; n++) {
          xyz[2] = r/20.*n;
          addSensors(xyz, restart);
        }

        double pi4 = std::atan(1.);
        for (int n=-20; n<=20; n++) {
          xyz[1] = r/20.*n*std::cos(pi4);
          xyz[2] = r/20.*n*std::sin(pi4);
          addSensors(xyz, restart);
        }

        for (int n=-20; n<=20; n++) {
          xyz[1] = r/20.*n*std::cos(-pi4);
          xyz[2] = r/20.*n*std::sin(-pi4);
          addSensors(xyz, restart);
        }
    }


    void initialize(double& t, double& dt){
        sendBoundaryReceivePressure();
        dt = stableTimeStep()*_subIterations;
        t = getCurrentTime();
        if (_sensoroutput) 
          instrumentWithSensors(false); 
    }

 
    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;        
        mShell()->iterateOverVertices(enforceBC);

        sendBoundaryReceivePressure();

        if (_sensoroutput) 
          instrumentWithSensors(true);
    }

    
// 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;

    double                               _timeStepReduction;

    int                                  _subIterations;

    int                                  _sensoroutput;

    std::vector<BasicSensorSp *>         _sensors;
};

#endif

---