explicit_builder.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Ruben Zorrilla
//
//
 
#if !defined(KRATOS_EXPLICIT_BUILDER)
#define  KRATOS_EXPLICIT_BUILDER
 
// System includes
#include <set>
#include <unordered_set>
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/model_part.h"
#include "utilities/parallel_utilities.h"
#include "utilities/constraint_utilities.h"
#include "includes/kratos_parameters.h"
#include "factories/factory.h"
#include "utilities/atomic_utilities.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
 
template<class TSparseSpace, class TDenseSpace >
class ExplicitBuilder
{
public:
 
    typedef std::size_t SizeType;
 
    typedef std::size_t IndexType;
 
    typedef typename TSparseSpace::DataType TDataType;
 
    typedef typename TSparseSpace::MatrixType TSystemMatrixType;
 
    typedef typename TSparseSpace::VectorType TSystemVectorType;
 
    typedef typename TSparseSpace::MatrixPointerType TSystemMatrixPointerType;
 
    typedef typename TSparseSpace::VectorPointerType TSystemVectorPointerType;
 
    typedef typename TDenseSpace::MatrixType LocalSystemMatrixType;
 
    typedef typename TDenseSpace::VectorType LocalSystemVectorType;
 
    typedef ModelPart::DofType DofType;
 
    typedef ModelPart::DofsArrayType DofsArrayType;
 
    typedef ModelPart::DofsVectorType DofsVectorType;
 
    typedef typename std::unordered_set<DofType::Pointer, DofPointerHasher> DofSetType;
 
    typedef ModelPart::NodesContainerType NodesArrayType;
    typedef ModelPart::ElementsContainerType ElementsArrayType;
    typedef ModelPart::ConditionsContainerType ConditionsArrayType;
 
    typedef PointerVectorSet<Element, IndexedObject> ElementsContainerType;
 
    typedef ExplicitBuilder<TSparseSpace, TDenseSpace> ClassType;
 
    KRATOS_CLASS_POINTER_DEFINITION(ExplicitBuilder);
 
 
    explicit ExplicitBuilder(Parameters ThisParameters)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    explicit ExplicitBuilder() = default;
 
    virtual ~ExplicitBuilder() = default;
 
 
    virtual typename ClassType::Pointer Create(Parameters ThisParameters) const
    {
        return Kratos::make_shared<ClassType>(ThisParameters);
    }
 
 
 
 
    bool GetCalculateReactionsFlag() const
    {
        return mCalculateReactionsFlag;
    }
 
    void SetCalculateReactionsFlag(bool CalculateReactionsFlag)
    {
        mCalculateReactionsFlag = CalculateReactionsFlag;
    }
 
    bool GetDofSetIsInitializedFlag() const
    {
        return mDofSetIsInitialized;
    }
 
    void SetDofSetIsInitializedFlag(bool DofSetIsInitialized)
    {
        mDofSetIsInitialized = DofSetIsInitialized;
    }
 
    bool GetResetDofSetFlag() const
    {
        return mResetDofSetFlag;
    }
 
    void SetResetDofSetFlag(bool ResetDofSetFlag)
    {
        mResetDofSetFlag = ResetDofSetFlag;
    }
 
    bool GetResetLumpedMassVectorFlag() const
    {
        return mResetLumpedMassVectorFlag;
    }
 
    void SetResetLumpedMassVectorFlag(bool ResetLumpedMassVectorFlag)
    {
        mResetLumpedMassVectorFlag = ResetLumpedMassVectorFlag;
    }
 
    unsigned int GetEquationSystemSize() const
    {
        return mEquationSystemSize;
    }
 
    DofsArrayType& GetDofSet()
    {
        return mDofSet;
    }
 
    const DofsArrayType& GetDofSet() const
    {
        return mDofSet;
    }
 
    TSystemVectorPointerType& pGetLumpedMassMatrixVector()
    {
        return mpLumpedMassVector;
    }
 
    TSystemVectorType& GetLumpedMassMatrixVector()
    {
        KRATOS_ERROR_IF_NOT(mpLumpedMassVector) << "Lumped mass matrix vector is not initialized!" << std::endl;
        return (*mpLumpedMassVector);
    }
 
    virtual void BuildRHS(ModelPart& rModelPart)
    {
        KRATOS_TRY
 
        // Build the Right Hand Side without Dirichlet conditions
        // We skip setting the Dirichlet nodes residual to zero for the sake of efficiency
        // Note that this is not required as the Dirichlet conditions are set in the strategy
        BuildRHSNoDirichlet(rModelPart);
 
        KRATOS_CATCH("")
    }
 
    virtual void BuildRHSNoDirichlet(ModelPart& rModelPart)
    {
        KRATOS_TRY
 
        // Initialize the reaction (residual)
        InitializeDofSetReactions();
 
        // Gets the array of elements, conditions and constraints from the modeler
        const auto &r_elements_array = rModelPart.Elements();
        const auto &r_conditions_array = rModelPart.Conditions();
        const int n_elems = static_cast<int>(r_elements_array.size());
        const int n_conds = static_cast<int>(r_conditions_array.size());
 
        const auto& r_process_info = rModelPart.GetProcessInfo();
 
#pragma omp parallel firstprivate(n_elems, n_conds)
        {
#pragma omp for schedule(guided, 512) nowait
            // Assemble all elements
            for (int i_elem = 0; i_elem < n_elems; ++i_elem) {
                auto it_elem = r_elements_array.begin() + i_elem;
                // If the element is active
                if (it_elem->IsActive()) {
                    // Calculate elemental explicit residual contribution
                    // The explicit builder and solver assumes that the residual contribution is assembled in the REACTION variables
                    it_elem->AddExplicitContribution(r_process_info);
                }
            }
 
            // Assemble all conditions
#pragma omp for schedule(guided, 512)
            for (int i_cond = 0; i_cond < n_conds; ++i_cond) {
                auto it_cond = r_conditions_array.begin() + i_cond;
                // If the condition is active
                if (it_cond->IsActive()) {
                    // Calculate condition explicit residual contribution
                    // The explicit builder and solver assumes that the residual contribution is assembled in the REACTION variables
                    it_cond->AddExplicitContribution(r_process_info);
                }
            }
        }
 
        KRATOS_CATCH("")
    }
 
    virtual void ApplyConstraints(ModelPart& rModelPart)
    {
        // First we reset the slave dofs
        ConstraintUtilities::ResetSlaveDofs(rModelPart);
 
        // Now we apply the constraints
        ConstraintUtilities::ApplyConstraints(rModelPart);
    }
 
    virtual void Initialize(ModelPart& rModelPart)
    {
        if (!mDofSetIsInitialized) {
            // Initialize the DOF set and the equation ids
            this->SetUpDofSet(rModelPart);
            this->SetUpDofSetEquationIds();
            // Set up the lumped mass vector
            this->SetUpLumpedMassVector(rModelPart);
        } else if (!mLumpedMassVectorIsInitialized) {
            KRATOS_WARNING("ExplicitBuilder") << "Calling Initialize() with already initialized DOF set. Initializing lumped mass vector." << std::endl;;
            // Only set up the lumped mass vector
            this->SetUpLumpedMassVector(rModelPart);
        } else {
            KRATOS_WARNING("ExplicitBuilder") << "Calling Initialize() with already initialized DOF set and lumped mass vector." << std::endl;;
        }
    }
 
    virtual void InitializeSolutionStep(ModelPart& rModelPart)
    {
        // Check the operations that are required to be done
        if (mResetDofSetFlag) {
            // If required (e.g. topology changes) reset the DOF set
            // Note that we also set lumped mass vector in this case
            this->SetUpDofSet(rModelPart);
            this->SetUpDofSetEquationIds();
            this->SetUpLumpedMassVector(rModelPart);
        } else if (mResetLumpedMassVectorFlag) {
            // Only reset the lumped mass vector
            this->SetUpLumpedMassVector(rModelPart);
        }
 
        // Initialize the reactions (residual)
        this->InitializeDofSetReactions();
    }
 
    virtual void FinalizeSolutionStep(ModelPart& rModelPart)
    {
        // If required, calculate the reactions
        if (mCalculateReactionsFlag) {
            this->CalculateReactions();
        }
    }
 
    virtual void Clear()
    {
        this->mDofSet = DofsArrayType();
        this->mpLumpedMassVector.reset();
 
        KRATOS_INFO_IF("ExplicitBuilder", this->GetEchoLevel() > 0) << "Clear Function called" << std::endl;
    }
 
    virtual int Check(const ModelPart& rModelPart) const
    {
        KRATOS_TRY
 
        return 0;
 
        KRATOS_CATCH("");
    }
 
    virtual Parameters GetDefaultParameters() const
    {
        const Parameters default_parameters = Parameters(R"(
        {
            "name" : "explicit_builder"
        })");
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "explicit_builder";
    }
 
    void SetEchoLevel(int Level)
    {
        mEchoLevel = Level;
    }
 
    int GetEchoLevel() const
    {
        return mEchoLevel;
    }
 
 
 
 
 
 
    virtual std::string Info() const
    {
        return "ExplicitBuilder";
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const
    {
        rOStream << Info();
    }
 
    virtual void PrintData(std::ostream& rOStream) const
    {
        rOStream << Info();
    }
 
 
 
protected:
 
 
 
    DofsArrayType mDofSet; 
 
    TSystemVectorPointerType mpLumpedMassVector; // The lumped mass vector associated to the DOF set
 
    bool mResetDofSetFlag = false;  
 
    bool mResetLumpedMassVectorFlag = false;  
 
    bool mDofSetIsInitialized = false; 
 
    bool mLumpedMassVectorIsInitialized = false; 
 
    bool mCalculateReactionsFlag = false; 
 
    unsigned int mEquationSystemSize; 
 
    int mEchoLevel = 0;
 
 
 
 
    virtual void SetUpDofSet(const ModelPart& rModelPart)
    {
        KRATOS_TRY;
 
        KRATOS_INFO_IF("ExplicitBuilder", this->GetEchoLevel() > 1) << "Setting up the dofs" << std::endl;
 
        // Gets the array of elements, conditions and constraints from the modeler
        const auto &r_elements_array = rModelPart.Elements();
        const auto &r_conditions_array = rModelPart.Conditions();
        const auto &r_constraints_array = rModelPart.MasterSlaveConstraints();
        const int n_elems = static_cast<int>(r_elements_array.size());
        const int n_conds = static_cast<int>(r_conditions_array.size());
        const int n_constraints = static_cast<int>(r_constraints_array.size());
 
        // Global dof set
        DofSetType dof_global_set;
        dof_global_set.reserve(n_elems*20);
 
        // Auxiliary DOFs list
        DofsVectorType dof_list;
        DofsVectorType second_dof_list; // The second_dof_list is only used on constraints to include master/slave relations
 
#pragma omp parallel firstprivate(dof_list, second_dof_list)
        {
            const auto& r_process_info = rModelPart.GetProcessInfo();
 
            // We cleate the temporal set and we reserve some space on them
            DofSetType dofs_tmp_set;
            dofs_tmp_set.reserve(20000);
 
            // Get the DOFs list from each element and insert it in the temporary set
#pragma omp for schedule(guided, 512) nowait
            for (int i_elem = 0; i_elem < n_elems; ++i_elem) {
                const auto it_elem = r_elements_array.begin() + i_elem;
                it_elem->GetDofList(dof_list, r_process_info);
                dofs_tmp_set.insert(dof_list.begin(), dof_list.end());
            }
 
            // Get the DOFs list from each condition and insert it in the temporary set
#pragma omp for schedule(guided, 512) nowait
            for (int i_cond = 0; i_cond < n_conds; ++i_cond) {
                const auto it_cond = r_conditions_array.begin() + i_cond;
                it_cond->GetDofList(dof_list, r_process_info);
                dofs_tmp_set.insert(dof_list.begin(), dof_list.end());
            }
 
            // Get the DOFs list from each constraint and insert it in the temporary set
#pragma omp for  schedule(guided, 512) nowait
            for (int i_const = 0; i_const < n_constraints; ++i_const) {
                auto it_const = r_constraints_array.begin() + i_const;
                it_const->GetDofList(dof_list, second_dof_list, r_process_info);
                dofs_tmp_set.insert(dof_list.begin(), dof_list.end());
                dofs_tmp_set.insert(second_dof_list.begin(), second_dof_list.end());
            }
 
            // We merge all the sets in one thread
#pragma omp critical
            {
                dof_global_set.insert(dofs_tmp_set.begin(), dofs_tmp_set.end());
            }
        }
 
        KRATOS_INFO_IF("ExplicitBuilder", ( this->GetEchoLevel() > 2)) << "Initializing ordered array filling\n" << std::endl;
 
        // Ordering the global DOF set
        mDofSet = DofsArrayType();
        DofsArrayType temp_dof_set;
        temp_dof_set.reserve(dof_global_set.size());
        for (auto it_dof = dof_global_set.begin(); it_dof != dof_global_set.end(); ++it_dof) {
            temp_dof_set.push_back(*it_dof);
        }
        temp_dof_set.Sort();
        mDofSet = temp_dof_set;
        mEquationSystemSize = mDofSet.size();
 
        // DoFs set checks
        // Throws an exception if there are no Degrees Of Freedom involved in the analysis
        KRATOS_ERROR_IF(mDofSet.size() == 0) << "No degrees of freedom!" << std::endl;
 
        // Check if each DOF has a reaction. Note that the explicit residual is stored in these
        for (auto it_dof = mDofSet.begin(); it_dof != mDofSet.end(); ++it_dof) {
            KRATOS_ERROR_IF_NOT(it_dof->HasReaction()) << "Reaction variable not set for the following : " << std::endl
                << "Node : " << it_dof->Id() << std::endl
                << "Dof : " << (*it_dof) << std::endl << "Not possible to calculate reactions." << std::endl;
        }
 
        mDofSetIsInitialized = true;
 
        KRATOS_INFO_IF("ExplicitBuilder", ( this->GetEchoLevel() > 2)) << "Number of degrees of freedom:" << mDofSet.size() << std::endl;
        KRATOS_INFO_IF("ExplicitBuilder", ( this->GetEchoLevel() > 2 && rModelPart.GetCommunicator().MyPID() == 0)) << "Finished setting up the dofs" << std::endl;
        KRATOS_INFO_IF("ExplicitBuilder", ( this->GetEchoLevel() > 2)) << "End of setup dof set\n" << std::endl;
 
        KRATOS_CATCH("");
    }
 
    virtual void SetUpDofSetEquationIds()
    {
        // Firstly check that the DOF set is initialized
        KRATOS_ERROR_IF_NOT(mDofSetIsInitialized) << "Trying to set the equation ids. before initializing the DOF set. Please call the SetUpDofSet() before." << std::endl;
        KRATOS_ERROR_IF(mEquationSystemSize == 0) << "Trying to set the equation ids. in an empty DOF set (equation system size is 0)." << std::endl;
 
        // Loop the DOF set to assign the equation ids
        IndexPartition<int>(mEquationSystemSize).for_each(
            [&](int i_dof){
                auto it_dof = mDofSet.begin() + i_dof;
                it_dof->SetEquationId(i_dof);
            }
        );
    }
 
    virtual void SetUpLumpedMassVector(const ModelPart &rModelPart)
    {
        KRATOS_TRY;
 
        KRATOS_INFO_IF("ExplicitBuilder", this->GetEchoLevel() > 1) << "Setting up the lumped mass matrix vector" << std::endl;
 
        // Initialize the lumped mass matrix vector
        // Note that the lumped mass matrix vector size matches the dof set one
        mpLumpedMassVector = TSystemVectorPointerType(new TSystemVectorType(GetDofSet().size()));
        TDenseSpace::SetToZero(*mpLumpedMassVector);
 
        // Loop the elements to get the lumped mass vector
        LocalSystemVectorType elem_mass_vector;
        std::vector<std::size_t> elem_equation_id;
        const auto &r_elements_array = rModelPart.Elements();
        const auto &r_process_info = rModelPart.GetProcessInfo();
        const int n_elems = static_cast<int>(r_elements_array.size());
 
#pragma omp for private(elem_mass_vector) schedule(guided, 512) nowait
        for (int i_elem = 0; i_elem < n_elems; ++i_elem) {
            const auto it_elem = r_elements_array.begin() + i_elem;
 
            // Calculate the elemental lumped mass vector
            it_elem->CalculateLumpedMassVector(elem_mass_vector, r_process_info);
            it_elem->EquationIdVector(elem_equation_id, r_process_info);
 
            // Update value of lumped mass vector
            for (IndexType i = 0; i < elem_equation_id.size(); ++i) {
                AtomicAdd((*mpLumpedMassVector)[elem_equation_id[i]], elem_mass_vector(i));
            }
        }
 
        // Set the lumped mass vector flag as true
        mLumpedMassVectorIsInitialized = true;
 
        KRATOS_CATCH("");
    }
 
    virtual void InitializeDofSetReactions()
    {
        // Firstly check that the DOF set is initialized
        KRATOS_ERROR_IF_NOT(mDofSetIsInitialized) << "Trying to initialize the explicit residual but the DOFs set is not initialized yet." << std::endl;
        KRATOS_ERROR_IF(mEquationSystemSize == 0) << "Trying to set the equation ids. in an empty DOF set (equation system size is 0)." << std::endl;
 
        // Loop the reactions to initialize them to zero
        block_for_each(
            mDofSet,
            [](DofType& rDof){
                rDof.GetSolutionStepReactionValue() = 0.0;
            }
        );
    }
 
    virtual void CalculateReactions()
    {
        if (mCalculateReactionsFlag) {
            // Firstly check that the DOF set is initialized
            KRATOS_ERROR_IF_NOT(mDofSetIsInitialized) << "Trying to initialize the explicit residual but the DOFs set is not initialized yet." << std::endl;
            KRATOS_ERROR_IF(mEquationSystemSize == 0) << "Trying to set the equation ids. in an empty DOF set (equation system size is 0)." << std::endl;
 
            // Calculate the reactions as minus the current value
            // Note that we take advantage of the fact that Kratos always works with a residual based formulation
            // This means that the reactions are minus the residual. As we use the reaction as residual container
            // during the explicit resolution of the problem, the calculate reactions is as easy as switching the sign
            block_for_each(
                mDofSet,
                [](DofType& rDof){
                    auto& r_reaction_value = rDof.GetSolutionStepReactionValue();
                    r_reaction_value *= -1.0;
                }
            );
        }
    }
 
    virtual Parameters ValidateAndAssignParameters(
        Parameters ThisParameters,
        const Parameters DefaultParameters
        ) const
    {
        ThisParameters.ValidateAndAssignDefaults(DefaultParameters);
        return ThisParameters;
    }
 
    virtual void AssignSettings(const Parameters ThisParameters)
    {
    }
 
 
 
 
 
 
 
private:
 
    static std::vector<Internals::RegisteredPrototypeBase<ClassType>> msPrototypes;
 
 
 
 
 
 
 
 
 
 
 
 
}; /* Class ExplicitBuilder */
 
 
} /* namespace Kratos.*/
 
#endif /* KRATOS_EXPLICIT_BUILDER  defined */