gsPreCICE.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsConfig.h>
#include <precice/SolverInterface.hpp>

namespace gismo {

template<class T>
class gsPreCICE
{
public:

    gsPreCICE()
    :
    m_interface("participantName", "configurationFileName", 0, 1)
    {
        // #ifdef GISMO_WITH_MPI --> rank and size
        // m_interface = precice::SolverInterface( "participantName", "configurationFileName", 0, 1 ); // last arguments rank, size
    }

    gsPreCICE(std::string participantName, std::string configurationFileName)
    :
    m_interface(participantName, configurationFileName, 0, 1) // last arguments rank, size
    {
    }

    bool isCouplingOngoing() const { return m_interface.isCouplingOngoing(); }
    bool isReadDataAvailable() const { return m_interface.isReadDataAvailable(); }
    bool isWriteDataRequired(real_t computedTimestepLength) const { return m_interface.isWriteDataRequired(computedTimestepLength); }
    bool isTimeWindowComplete() const { return m_interface.isTimeWindowComplete(); }
    bool isActionRequired(const std::string &action) const { return m_interface.isActionRequired(action); }
    void markActionFulfilled(const std::string &action) { m_interface.markActionFulfilled(action); }

    // TODO: These functions are precice constants and are preferably called outside of the class
    const std::string actionWriteInitialData() { return precice::constants::actionWriteInitialData(); }
    const std::string actionWriteIterationCheckpoint() { return precice::constants::actionWriteIterationCheckpoint(); }
    const std::string actionReadIterationCheckpoint() { return precice::constants::actionReadIterationCheckpoint(); }


    T initialize()
    {
        m_precicedt = m_interface.initialize();
        return m_precicedt;
    }


    void addMesh(const std::string & meshName, const gsMatrix<T> & points)
    {
        m_meshNames.push_back(meshName);
        m_meshIDs.push_back(m_interface.getMeshID(m_meshNames.back()));

        gsMatrix<T> pointsTranspose = points;
        pointsTranspose.blockTransposeInPlace(1);
        m_sizes.push_back(points.cols());
        m_positions.push_back(pointsTranspose.data());

        index_t * vertexIDs;
        m_interface.setMeshVertices(m_meshIDs.back(),m_sizes.back(),m_positions.back(),vertexIDs);
        m_vertexIDs.push_back(vertexIDs);
    }

    T advance(T dt)
    {
        m_precicedt = m_interface.advance(dt);
        return m_precicedt;
    }

    void finalize()
    {
        m_interface.finalize();
    }

    void readBlockScalarData(const index_t & dataID, const index_t & size, const index_t * IDs, gsMatrix<T> & values) const
    {
        T * values_ptr = new T[size];
        m_interface.readBlockScalarData(dataID,size,IDs,values_ptr);

        values = Eigen::Map<typename gsMatrix<T>::Base>(values_ptr,1,size);
    }

    void readBlockScalarData(const index_t & meshID, const index_t & dataID, const gsMatrix<T> & coords, gsMatrix<T> & values) const
    {
        index_t * IDs = new index_t[coords.cols()];
        this->getMeshVertexIDsFromPositions(meshID,coords,IDs);

        this->readBlockScalarData(dataID,coords.cols(),IDs,values);
    }

    // void readBlockVectorData(const index_t & dataID, const index_t & size, const index_t * IDs, gsMatrix<T> & values) const
    // {
    //     T * values_ptr = new T[size];
    //     m_interface.readBlockVectorData(dataID,size,IDs,values_ptr);
    //     index_t rows =
    //     values = Eigen::Map<typename gsMatrix<T>::Base>(values_ptr,1,size);

    //     // CAST VALUES TO gsMATRIX
    // }

    // void readBlockVectorData(const index_t & meshID, const index_t & dataID, const gsMatrix<T> & coords, gsMatrix<T> & values) const
    // {
    //     index_t * IDs = new index_t[coords.cols()];
    //     this->getMeshVertexIDsFromPositions(meshID,coords,IDs);

    //     this->readBlockVectorData(dataID,coords.cols(),IDs,values);
    // }

    void writeBlockScalarData(const index_t & dataID, const index_t & size, const index_t * IDs, const gsMatrix<T> & values)
    {
        m_interface.writeBlockScalarData(dataID,size,IDs,values.data());
    }

    void writeBlockScalarData(const index_t & meshID, const index_t & dataID, const gsMatrix<T> & coords, const gsMatrix<T> & values)
    {
        index_t * IDs = new index_t[coords.cols()];
        this->getMeshVertexIDsFromPositions(meshID,coords,IDs);

        this->writeBlockScalarData(dataID,coords.cols(),IDs,values);
    }

    // void writeBlockVectorData(const index_t & dataID, const index_t & size, const index_t * IDs, const gsMatrix<T> & values) const
    // {
    //     m_interface.writeBlockVectorData(dataID,size,IDs,values.data());
    // }

    // void writeBlockVectorData(const index_t & meshID, const index_t & dataID, const gsMatrix<T> & coords, const gsMatrix<T> & values) const
    // {
    //     index_t * IDs = new index_t[coords.cols()];
    //     this->getMeshVertexIDsFromPositions(meshID,coords,IDs);

    //     this->writeBlockVectorData(dataID,coords.cols(),IDs,values);
    // }

    void getMeshVertexIDsFromPositions(const index_t & meshID, const gsMatrix<T> & coords, gsMatrix<index_t> & IDs) const
    {
        index_t * ID_ptr;
        this->getMeshVertexIDsFromPositions(meshID,coords,ID_ptr);
        for (index_t k=0; k!=coords.cols(); k++)
            gsDebugVar(ID_ptr[k]);

        IDs = Eigen::Map<typename gsMatrix<index_t>::Base>(ID_ptr,1,coords.cols());
    }

    void getMeshVertexIDsFromPositions(const index_t & meshID, const gsMatrix<T> & coords, index_t * ID_ptr) const
    {
        gsMatrix<T> coordsTranspose = coords;
        coordsTranspose.blockTransposeInPlace(1);
        T * positions = coordsTranspose.data();
        m_interface.getMeshVertexIDsFromPositions(meshID,coords.cols(),positions,ID_ptr);
    }

    index_t getMeshID(const std::string &dataName) const
    {
        return m_interface.getMeshID(dataName);
    }

    index_t getDataID(const std::string &dataName, int meshID) const
    {
        return m_interface.getDataID(dataName,meshID);
    }

    virtual std::ostream &print(std::ostream &os) const
    {
        os << precice::getVersionInformation()<<"\n\n";
        os << "Interface has the following meshes:\n";
        for (index_t k=0; k!=m_meshNames.size(); k++)
            gsInfo<<m_meshNames[k]<<"\t (ID="<<m_meshIDs[k]<<")\n";
        return os;
    }

private:
    std::vector<std::string> m_meshNames;
    std::vector<index_t> m_meshIDs;
    std::vector<index_t> m_sizes;
    std::vector<T *> m_positions;
    std::vector<index_t *> m_vertexIDs;
    precice::SolverInterface m_interface;
    T m_precicedt;
};

} //namespace gismo