residual_based_bdf_displacement_scheme.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
 
#if !defined(KRATOS_RESIDUAL_BASED_BDF_DISPLACEMENT_SCHEME )
#define  KRATOS_RESIDUAL_BASED_BDF_DISPLACEMENT_SCHEME
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "solving_strategies/schemes/residual_based_bdf_scheme.h"
#include "includes/variables.h"
#include "includes/checks.h"
 
namespace Kratos
{
 
template<class TSparseSpace,  class TDenseSpace>
class ResidualBasedBDFDisplacementScheme
    : public ResidualBasedBDFScheme<TSparseSpace, TDenseSpace>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedBDFDisplacementScheme );
 
    typedef Scheme<TSparseSpace,TDenseSpace>                                  BaseType;
    typedef ResidualBasedImplicitTimeScheme<TSparseSpace,TDenseSpace> ImplicitBaseType;
    typedef ResidualBasedBDFScheme<TSparseSpace,TDenseSpace>               BDFBaseType;
    typedef ResidualBasedBDFDisplacementScheme<TSparseSpace, TDenseSpace>    ClassType;
 
    typedef typename BDFBaseType::TDataType                                  TDataType;
    typedef typename BDFBaseType::TSystemMatrixType                  TSystemMatrixType;
    typedef typename BDFBaseType::TSystemVectorType                  TSystemVectorType;
    typedef typename BDFBaseType::LocalSystemVectorType          LocalSystemVectorType;
    typedef typename BDFBaseType::LocalSystemMatrixType          LocalSystemMatrixType;
 
    typedef typename BDFBaseType::DofsArrayType                          DofsArrayType;
    typedef typename Element::DofsVectorType                            DofsVectorType;
 
    typedef ModelPart::NodesContainerType                               NodesArrayType;
    typedef ModelPart::ElementsContainerType                         ElementsArrayType;
    typedef ModelPart::ConditionsContainerType                     ConditionsArrayType;
 
    typedef double              ComponentType;
 
 
    explicit ResidualBasedBDFDisplacementScheme(Parameters ThisParameters)
        : BDFBaseType()
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    explicit ResidualBasedBDFDisplacementScheme(const std::size_t Order = 2)
        :BDFBaseType(Order)
    {
    }
 
    explicit ResidualBasedBDFDisplacementScheme(ResidualBasedBDFDisplacementScheme& rOther)
        :BDFBaseType(rOther)
    {
    }
 
    typename BaseType::Pointer Clone() override
    {
        return Kratos::make_shared<ResidualBasedBDFDisplacementScheme>(*this);
    }
 
    ~ResidualBasedBDFDisplacementScheme
    () override {}
 
 
 
    typename BaseType::Pointer Create(Parameters ThisParameters) const override
    {
        return Kratos::make_shared<ClassType>(ThisParameters);
    }
 
    void InitializeSolutionStep(
        ModelPart& rModelPart,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        ) override
    {
        KRATOS_TRY;
 
        // The current process info
        const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
 
        BDFBaseType::InitializeSolutionStep(rModelPart, rA, rDx, rb);
 
        // Updating time derivatives (nodally for efficiency)
        const auto it_node_begin = rModelPart.Nodes().begin();
 
        // Getting dimension
        KRATOS_WARNING_IF("ResidualBasedBDFDisplacementScheme", !r_current_process_info.Has(DOMAIN_SIZE)) << "DOMAIN_SIZE not defined. Please define DOMAIN_SIZE. 3D case will be assumed" << std::endl;
        const std::size_t dimension = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3;
 
        // Getting position
        const int velpos = it_node_begin->HasDofFor(VELOCITY_X) ? static_cast<int>(it_node_begin->GetDofPosition(VELOCITY_X)) : -1;
        const int accelpos = it_node_begin->HasDofFor(ACCELERATION_X) ? static_cast<int>(it_node_begin->GetDofPosition(ACCELERATION_X)) : -1;
 
        std::array<bool, 3> fixed = {false, false, false};
        const std::array<const Variable<ComponentType>*, 3> disp_components = {&DISPLACEMENT_X, &DISPLACEMENT_Y, &DISPLACEMENT_Z};
        const std::array<const Variable<ComponentType>*, 3> vel_components = {&VELOCITY_X, &VELOCITY_Y, &VELOCITY_Z};
        const std::array<const Variable<ComponentType>*, 3> accel_components = {&ACCELERATION_X, &ACCELERATION_Y, &ACCELERATION_Z};
 
        block_for_each(rModelPart.Nodes(), fixed, [&](Node& rNode, std::array<bool,3>& rFixedTLS){
            for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                rFixedTLS[i_dim] = false;
            }
 
            if (accelpos > -1) {
                for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                    if (rNode.GetDof(*accel_components[i_dim], accelpos + i_dim).IsFixed()) {
                        rNode.Fix(*disp_components[i_dim]);
                        rFixedTLS[i_dim] = true;
                    }
                }
            }
            if (velpos > -1) {
                for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                    if (rNode.GetDof(*vel_components[i_dim], velpos + i_dim).IsFixed() && !rFixedTLS[i_dim]) {
                        rNode.Fix(*disp_components[i_dim]);
                    }
                }
            }
        });
 
        KRATOS_CATCH("ResidualBasedBDFDisplacementScheme.InitializeSolutionStep");
    }
 
    void Predict(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        TSystemMatrixType& A,
        TSystemVectorType& Dx,
        TSystemVectorType& b
        ) override
    {
        KRATOS_TRY;
 
        // The current process info
        const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
        const double delta_time = r_current_process_info[DELTA_TIME];
 
        // Updating time derivatives (nodally for efficiency)
        const int num_nodes = static_cast<int>( rModelPart.Nodes().size() );
        const auto it_node_begin = rModelPart.Nodes().begin();
 
        // Getting position
        KRATOS_ERROR_IF_NOT(it_node_begin->HasDofFor(DISPLACEMENT_X)) << "ResidualBasedBDFDisplacementScheme:: DISPLACEMENT is not added" << std::endl;
        const int disppos = it_node_begin->GetDofPosition(DISPLACEMENT_X);
        const int velpos = it_node_begin->HasDofFor(VELOCITY_X) ? static_cast<int>(it_node_begin->GetDofPosition(VELOCITY_X)) : -1;
        const int accelpos = it_node_begin->HasDofFor(ACCELERATION_X) ? static_cast<int>(it_node_begin->GetDofPosition(ACCELERATION_X)) : -1;
 
        // Getting dimension
        KRATOS_WARNING_IF("ResidualBasedBDFDisplacementScheme", !r_current_process_info.Has(DOMAIN_SIZE)) << "DOMAIN_SIZE not defined. Please define DOMAIN_SIZE. 3D case will be assumed" << std::endl;
        const std::size_t dimension = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3;
 
        // Auxiliar variables
        std::array<bool, 3> predicted = {false, false, false};
        const std::array<const Variable<ComponentType>*, 3> disp_components = {&DISPLACEMENT_X, &DISPLACEMENT_Y, &DISPLACEMENT_Z};
        const std::array<const Variable<ComponentType>*, 3> vel_components = {&VELOCITY_X, &VELOCITY_Y, &VELOCITY_Z};
        const std::array<const Variable<ComponentType>*, 3> accel_components = {&ACCELERATION_X, &ACCELERATION_Y, &ACCELERATION_Z};
 
        IndexPartition<int>(num_nodes).for_each(predicted, [&](int Index, std::array<bool, 3>& rPredictedTLS){
            auto it_node = it_node_begin + Index;
 
            for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                rPredictedTLS[i_dim] = false;
            }
 
            const array_1d<double, 3>& dot2un1 = it_node->FastGetSolutionStepValue(ACCELERATION, 1);
            const array_1d<double, 3>& dotun1 = it_node->FastGetSolutionStepValue(VELOCITY,     1);
            const array_1d<double, 3>& un1 = it_node->FastGetSolutionStepValue(DISPLACEMENT, 1);
            const array_1d<double, 3>& dot2un0 = it_node->FastGetSolutionStepValue(ACCELERATION);
            array_1d<double, 3>& dotun0 = it_node->FastGetSolutionStepValue(VELOCITY);
            array_1d<double, 3>& un0 = it_node->FastGetSolutionStepValue(DISPLACEMENT);
 
            if (accelpos > -1) {
                for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                    if (it_node->GetDof(*accel_components[i_dim], accelpos + i_dim).IsFixed()) {
                        dotun0[i_dim] = dot2un0[i_dim];
                        for (std::size_t i_order = 1; i_order < this->mOrder + 1; ++i_order)
                            dotun0[i_dim] -= this->mBDF[i_order] * it_node->FastGetSolutionStepValue(*vel_components[i_dim], i_order);
                        dotun0[i_dim] /= this->mBDF[i_dim];
 
                        un0[i_dim] = dotun0[i_dim];
                        for (std::size_t i_order = 1; i_order < this->mOrder + 1; ++i_order)
                            un0[i_dim] -= this->mBDF[i_order] * it_node->FastGetSolutionStepValue(*disp_components[i_dim], i_order);
                        un0[i_dim] /= this->mBDF[i_dim];
                        rPredictedTLS[i_dim] = true;
                    }
                }
            }
            if (velpos > -1) {
                for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                    if (it_node->GetDof(*vel_components[i_dim], velpos + i_dim).IsFixed() && !rPredictedTLS[i_dim]) {
                        un0[i_dim] = dotun0[i_dim];
                        for (std::size_t i_order = 1; i_order < this->mOrder + 1; ++i_order)
                            un0[i_dim] -= this->mBDF[i_order] * it_node->FastGetSolutionStepValue(*disp_components[i_dim], i_order);
                        un0[i_dim] /= this->mBDF[i_dim];
                        rPredictedTLS[i_dim] = true;
                    }
                }
            }
            for (std::size_t i_dim = 0; i_dim < dimension; ++i_dim) {
                if (!it_node->GetDof(*disp_components[i_dim], disppos + i_dim).IsFixed() && !rPredictedTLS[i_dim]) {
                    un0[i_dim] = un1[i_dim] + delta_time * dotun1[i_dim] + 0.5 * std::pow(delta_time, 2) * dot2un1[i_dim];
                }
            }
 
            // Updating time derivatives
            UpdateFirstDerivative(it_node);
            UpdateSecondDerivative(it_node);
        });
 
        KRATOS_CATCH( "" );
    }
 
    int Check(const ModelPart& rModelPart) const override
    {
        KRATOS_TRY;
 
        const int err = BDFBaseType::Check(rModelPart);
        if(err!=0) return err;
 
        // Check that variables are correctly allocated
        for(auto& rnode : rModelPart.Nodes()) {
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(DISPLACEMENT,rnode)
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(VELOCITY,rnode)
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(ACCELERATION,rnode)
 
            KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_X, rnode)
            KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_Y, rnode)
            KRATOS_CHECK_DOF_IN_NODE(DISPLACEMENT_Z, rnode)
        }
 
        KRATOS_CATCH( "" );
 
        return 0;
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"               : "bdf_displacement_scheme",
            "integration_order"  : 2
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BDFBaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "bdf_displacement_scheme";
    }
 
 
 
 
    std::string Info() const override
    {
        return "ResidualBasedBDFDisplacementScheme";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
 
 
 
    inline void UpdateFirstDerivative(NodesArrayType::iterator itNode) override
    {
        array_1d<double, 3>& dotun0 = itNode->FastGetSolutionStepValue(VELOCITY);
        noalias(dotun0) = this->mBDF[0] * itNode->FastGetSolutionStepValue(DISPLACEMENT);
        for (std::size_t i_order = 1; i_order < this->mOrder + 1; ++i_order)
            noalias(dotun0) += this->mBDF[i_order] * itNode->FastGetSolutionStepValue(DISPLACEMENT, i_order);
    }
 
    inline void UpdateSecondDerivative(NodesArrayType::iterator itNode) override
    {
        array_1d<double, 3>& dot2un0 = itNode->FastGetSolutionStepValue(ACCELERATION);
        noalias(dot2un0) = this->mBDF[0] * itNode->FastGetSolutionStepValue(VELOCITY);
        for (std::size_t i_order = 1; i_order < this->mOrder + 1; ++i_order)
            noalias(dot2un0) += this->mBDF[i_order] * itNode->FastGetSolutionStepValue(VELOCITY, i_order);
    }
 
 
 
 
private:
 
 
 
 
 
 
}; /* Class ResidualBasedBDFDisplacementScheme */
}  /* namespace Kratos.*/
 
#endif /* KRATOS_RESIDUAL_BASED_BDF_DISPLACEMENT_SCHEME defined */