poromechanics_face_load_control_module_process.hpp

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Ignasi de Pouplana
//
 
 
#if !defined(KRATOS_POROMECHANICS_FACE_LOAD_CONTROL_MODULE_PROCESS )
#define  KRATOS_POROMECHANICS_FACE_LOAD_CONTROL_MODULE_PROCESS
 
#include "includes/table.h"
#include "includes/kratos_flags.h"
#include "includes/kratos_parameters.h"
#include "processes/process.h"
 
#include "poromechanics_application_variables.h"
 
namespace Kratos
{
 
 
class PoromechanicsFaceLoadControlModuleProcess : public Process
{
 
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(PoromechanicsFaceLoadControlModuleProcess);
 
    typedef Table<double,double> TableType;
 
 
    PoromechanicsFaceLoadControlModuleProcess(
        ModelPart& rModelPart,
        Parameters rParameters
        ) : Process(Flags()) , mrModelPart(rModelPart)
    {
        KRATOS_TRY
 
        Parameters default_parameters( R"(
            {
                "model_part_name":"MODEL_PART_NAME",
                "imposed_direction" : 0,
                "face_load_table_id" : 0,
                "initial_velocity" : 0.0,
                "limit_velocity" : 1.0,
                "velocity_factor" : 1.0,
                "initial_stiffness" : 1.0e7,
                "force_increment_tolerance": 100.0,
                "update_stiffness": true,
                "force_averaging_time": 1.0e-5
            }  )" );
 
        rParameters.ValidateAndAssignDefaults(default_parameters);
 
        double ImposedDirection = rParameters["imposed_direction"].GetInt();
 
        if(ImposedDirection == 2)
        {
            mDisplacementName = "DISPLACEMENT" + std::string("_Z");
            mReactionName = "REACTION" + std::string("_Z");
            mAverageReactionName = "AVERAGE_REACTION" + std::string("_Z");
            mTargetReactionName = "TARGET_REACTION" + std::string("_Z");
            mLoadingVelocityName = "LOADING_VELOCITY" + std::string("_Z");
        }
        else if(ImposedDirection == 1)
        {
            mDisplacementName = "DISPLACEMENT" + std::string("_Y");
            mReactionName = "REACTION" + std::string("_Y");
            mAverageReactionName = "AVERAGE_REACTION" + std::string("_Y");
            mTargetReactionName = "TARGET_REACTION" + std::string("_Y");
            mLoadingVelocityName = "LOADING_VELOCITY" + std::string("_Y");
        }
        else
        {
            mDisplacementName = "DISPLACEMENT" + std::string("_X");
            mReactionName = "REACTION" + std::string("_X");
            mAverageReactionName = "AVERAGE_REACTION" + std::string("_X");
            mTargetReactionName = "TARGET_REACTION" + std::string("_X");
            mLoadingVelocityName = "LOADING_VELOCITY" + std::string("_X");
        }
 
        mFaceLoadTableId = rParameters["face_load_table_id"].GetInt();
        mVelocity = rParameters["initial_velocity"].GetDouble();
        mLimitVelocity = rParameters["limit_velocity"].GetDouble();
        mVelocityFactor = rParameters["velocity_factor"].GetDouble();
        mForceIncrementTolerance = rParameters["force_increment_tolerance"].GetDouble();
        mUpdateStiffness = rParameters["update_stiffness"].GetBool();
        mReactionForceOld = 0.0;
        mStiffness = rParameters["initial_stiffness"].GetDouble();
        mForceAveragingTime = rParameters["force_averaging_time"].GetDouble(); // If set equal to the time step, no average is performed
        mVectorOfLastForces.resize(0);
 
        // Initialize Variables
        const int NNodes = static_cast<int>(mrModelPart.Nodes().size());
        ModelPart::NodesContainerType::iterator it_begin = mrModelPart.NodesBegin();
        array_1d<double,3> zero_vector = ZeroVector(3);
        #pragma omp parallel for
        for(int i = 0; i<NNodes; i++) {
            ModelPart::NodesContainerType::iterator it = it_begin + i;
            it->SetValue(AVERAGE_REACTION,zero_vector);
            it->SetValue(TARGET_REACTION,zero_vector);
            it->SetValue(LOADING_VELOCITY,zero_vector);
        }
 
        KRATOS_CATCH("");
    }
 
 
    ~PoromechanicsFaceLoadControlModuleProcess() override {}
 
 
    void Execute() override
    {
    }
 
    void ExecuteInitialize() override
    {
        KRATOS_TRY;
 
        typedef Variable<double> component_type;
        const component_type& disp_component = KratosComponents< component_type >::Get(mDisplacementName);
        const component_type& loading_vel_component = KratosComponents< component_type >::Get(mLoadingVelocityName);
 
        const int NNodes = static_cast<int>(mrModelPart.Nodes().size());
        ModelPart::NodesContainerType::iterator it_begin = mrModelPart.NodesBegin();
 
        #pragma omp parallel for
        for(int i = 0; i<NNodes; i++) {
            ModelPart::NodesContainerType::iterator it = it_begin + i;
            it->Fix(disp_component);
            it->FastGetSolutionStepValue(disp_component) = 0.0;
            it->GetValue(loading_vel_component) = mVelocity;
        }
 
        KRATOS_CATCH("");
    }
 
    void ExecuteInitializeSolutionStep() override
    {
        KRATOS_TRY;
 
        const double CurrentTime = mrModelPart.GetProcessInfo()[TIME];
 
        const int NNodes = static_cast<int>(mrModelPart.Nodes().size());
        ModelPart::NodesContainerType::iterator it_begin = mrModelPart.NodesBegin();
 
        double ReactionForce = CalculateReactionForce();
        ReactionForce = UpdateVectorOfHistoricalForcesAndComputeNewAverage(ReactionForce);
 
        // Update K if required
        const double DeltaTime = mrModelPart.GetProcessInfo()[DELTA_TIME];
        if(mUpdateStiffness == true) {
            mStiffness = EstimateStiffness(ReactionForce,DeltaTime);
        }
        mReactionForceOld = ReactionForce;
 
        // Update velocity
        const double NextTargetReaction = this->GetTargetReaction(CurrentTime+DeltaTime);
        const double df_target = NextTargetReaction - ReactionForce;
        double delta_velocity = df_target/(mStiffness * DeltaTime) - mVelocity;
 
        if(std::abs(df_target) < mForceIncrementTolerance) {
            delta_velocity = -mVelocity;
        }
 
        mVelocity += mVelocityFactor * delta_velocity;
 
        if(std::abs(mVelocity) > std::abs(mLimitVelocity)) {
            if(mVelocity >= 0.0) {
                mVelocity = std::abs(mLimitVelocity);
            } else {
                mVelocity = - std::abs(mLimitVelocity);
            }
        }
 
        typedef Variable<double> component_type;
        const component_type& disp_component = KratosComponents< component_type >::Get(mDisplacementName);
 
        // Update Imposed displacement
        #pragma omp parallel for
        for(int i = 0; i<NNodes; i++) {
            ModelPart::NodesContainerType::iterator it = it_begin + i;
            it->FastGetSolutionStepValue(disp_component) += mVelocity * DeltaTime;
        }
 
        KRATOS_CATCH("");
    }
 
 
    void ExecuteFinalizeSolutionStep() override
    {
        const double CurrentTime = mrModelPart.GetProcessInfo()[TIME];
 
        const int NNodes = static_cast<int>(mrModelPart.Nodes().size());
        ModelPart::NodesContainerType::iterator it_begin = mrModelPart.NodesBegin();
 
        double ReactionForce = CalculateReactionForce();
 
        typedef Variable<double> component_type;
        const component_type& target_reaction_component = KratosComponents< component_type >::Get(mTargetReactionName);
        const component_type& average_reaction_component = KratosComponents< component_type >::Get(mAverageReactionName);
        const component_type& loading_vel_component = KratosComponents< component_type >::Get(mLoadingVelocityName);
 
        // Update Imposed displacement
        #pragma omp parallel for
        for(int i = 0; i<NNodes; i++) {
            ModelPart::NodesContainerType::iterator it = it_begin + i;
            // Save calculated velocity and reaction for print
            it->GetValue(target_reaction_component) = this->GetTargetReaction(CurrentTime);
            it->GetValue(average_reaction_component) = ReactionForce;
            it->GetValue(loading_vel_component) = mVelocity;
        }
    }
 
    std::string Info() const override
    {
        return "PoromechanicsFaceLoadControlModuleProcess";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "PoromechanicsFaceLoadControlModuleProcess";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
    }
 
 
protected:
 
 
    ModelPart& mrModelPart;
    // unsigned int mImposedDirection;
    std::string mDisplacementName;
    std::string mReactionName;
    std::string mAverageReactionName;
    std::string mTargetReactionName;
    std::string mLoadingVelocityName;
    unsigned int mFaceLoadTableId;
    double mVelocity;
    double mLimitVelocity;
    double mVelocityFactor;
    double mReactionForceOld;
    double mForceIncrementTolerance;
    double mStiffness;
    bool mUpdateStiffness;
    std::vector<double> mVectorOfLastForces;
    double mForceAveragingTime;
 
 
    virtual double GetTargetReaction(const double& time) {
 
        // Calculate face_area
        const int NCons = static_cast<int>(mrModelPart.Conditions().size());
        ModelPart::ConditionsContainerType::iterator con_begin = mrModelPart.ConditionsBegin();
        double FaceArea = 0.0;
        #pragma omp parallel for reduction(+:FaceArea)
        for(int i = 0; i < NCons; i++)
        {
            ModelPart::ConditionsContainerType::iterator itCond = con_begin + i;
            FaceArea += itCond->GetGeometry().Area();
        }
 
        // Get Target Face Load from table
        TableType::Pointer pFaceLoadTable = mrModelPart.pGetTable(mFaceLoadTableId);
        const double face_load = pFaceLoadTable->GetValue(time);
 
        return face_load*FaceArea;
    }
 
    virtual double UpdateVectorOfHistoricalForcesAndComputeNewAverage(const double& last_reaction) {
        KRATOS_TRY;
        int length_of_vector = mVectorOfLastForces.size();
        if (length_of_vector == 0) { //only the first time
            int number_of_steps_for_force_averaging = (int) (mForceAveragingTime / mrModelPart.GetProcessInfo()[DELTA_TIME]);
            if(number_of_steps_for_force_averaging < 1) number_of_steps_for_force_averaging = 1;
            mVectorOfLastForces.resize(number_of_steps_for_force_averaging);
            KRATOS_INFO("Poro control module") << " 'number_of_steps_for_force_averaging' is "<< number_of_steps_for_force_averaging << std::endl;
        }
 
        length_of_vector = mVectorOfLastForces.size();
 
        if(length_of_vector > 1) {
            for(int i=1; i<length_of_vector; i++) {
                mVectorOfLastForces[i-1] = mVectorOfLastForces[i];
            }
        }
        mVectorOfLastForces[length_of_vector-1] = last_reaction;
 
        double average = 0.0;
        for(int i=0; i<length_of_vector; i++) {
            average += mVectorOfLastForces[i];
        }
        average /= (double) length_of_vector;
        return average;
 
        KRATOS_CATCH("");
    }
 
    virtual double CalculateReactionForce() {
 
        typedef Variable<double> component_type;
        const component_type& reaction_component = KratosComponents< component_type >::Get(mReactionName);
 
        const int NNodes = static_cast<int>(mrModelPart.Nodes().size());
        ModelPart::NodesContainerType::iterator it_begin = mrModelPart.NodesBegin();
 
        // Calculate Average ReactionForce
        double FaceReaction = 0.0;
        #pragma omp parallel for reduction(+:FaceReaction)
        for(int i = 0; i<NNodes; i++){
            ModelPart::NodesContainerType::iterator it = it_begin + i;
            FaceReaction += it->FastGetSolutionStepValue(reaction_component);
        }
 
        return FaceReaction/NNodes;
    }
 
    virtual double EstimateStiffness(const double& rReactionForce, const double& rDeltaTime) {
        double K_estimated = mStiffness;
        if(std::abs(mVelocity) > 1.0e-12 && std::abs(rReactionForce-mReactionForceOld) > mForceIncrementTolerance) {
            K_estimated = std::abs((rReactionForce-mReactionForceOld)/(mVelocity * rDeltaTime));
        }
        return K_estimated;
    }
 
 
private:
 
    PoromechanicsFaceLoadControlModuleProcess& operator=(PoromechanicsFaceLoadControlModuleProcess const& rOther);
 
    //PoromechanicsFaceLoadControlModuleProcess(PoromechanicsFaceLoadControlModuleProcess const& rOther);
 
 
 
}; // Class PoromechanicsFaceLoadControlModuleProcess
 
inline std::istream& operator >> (std::istream& rIStream,
                                  PoromechanicsFaceLoadControlModuleProcess& rThis);
 
inline std::ostream& operator << (std::ostream& rOStream,
                                  const PoromechanicsFaceLoadControlModuleProcess& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
} // namespace Kratos.
 
#endif /* KRATOS_POROMECHANICS_FACE_LOAD_CONTROL_MODULE_PROCESS defined */