simple_mortar_mapper_process.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#pragma once
 
// System includes
#include <unordered_map>
 
// External includes
 
// Project includes
#include "processes/process.h"
#include "includes/kratos_parameters.h"
#include "includes/model_part.h"
#include "includes/mortar_classes.h"
#include "spaces/ublas_space.h"
#include "linear_solvers/linear_solver.h"
#include "utilities/atomic_utilities.h"
#include "utilities/exact_mortar_segmentation_utility.h"
 
/* Tree structures */
#include "spatial_containers/point_object.h"
#include "spatial_containers/spatial_containers.h" // kd-tree
 
namespace Kratos
{
 
 
    using SizeType = std::size_t;
 
 
 
 
template< const SizeType TDim, const SizeType TNumNodes, class TVarType, const SizeType TNumNodesMaster = TNumNodes>
class KRATOS_API(KRATOS_CORE) SimpleMortarMapperProcess
        : public Process
{
public:
 
    // DEFINITION OF FLAGS TO CONTROL THE BEHAVIOUR
    KRATOS_DEFINE_LOCAL_FLAG(AVERAGE_NORMAL);                      
    KRATOS_DEFINE_LOCAL_FLAG(DISCONTINOUS_INTERFACE);              
    KRATOS_DEFINE_LOCAL_FLAG(ORIGIN_IS_HISTORICAL);                
    KRATOS_DEFINE_LOCAL_FLAG(DESTINATION_IS_HISTORICAL);           
    KRATOS_DEFINE_LOCAL_FLAG(ORIGIN_SKIN_IS_CONDITION_BASED);      
    KRATOS_DEFINE_LOCAL_FLAG(DESTINATION_SKIN_IS_CONDITION_BASED); 
    KRATOS_DEFINE_LOCAL_FLAG(CONSIDER_TESELLATION);                
 
    KRATOS_CLASS_POINTER_DEFINITION(SimpleMortarMapperProcess);
 
    using PointType = Point;
 
    using GeometryType = Geometry<Node>;
 
    using GeometryPointType = Geometry<PointType>;
 
    using IntegrationMethod = GeometryData::IntegrationMethod;
 
    using LineType = Line2D2<PointType>;
 
    using TriangleType = Triangle3D3<PointType>;
 
    using DecompositionType = typename std::conditional<TDim == 2, LineType, TriangleType>::type;
 
    using SparseSpaceType = UblasSpace<double, CompressedMatrix, Vector>;
 
    using LocalSpaceType = UblasSpace<double, Matrix, Vector>;
 
    using MatrixType = typename SparseSpaceType::MatrixType;
 
    using VectorType = typename SparseSpaceType::VectorType;
 
    using LinearSolverType = LinearSolver<SparseSpaceType, LocalSpaceType>;
 
    using IndexType = std::size_t;
 
    using IntMap = std::unordered_map<IndexType, IndexType>;
 
    using BoundedMatrixType = BoundedMatrix<double, TNumNodes, TNumNodes>;
 
    using PointMapperType = PointObject<GeometricalObject>;
    using PointTypePointer = typename PointMapperType::Pointer;
    using PointVector = std::vector<PointTypePointer>;
 
    // KDtree definitions
    using BucketType = Bucket<3ul, PointMapperType, PointVector>;
    using KDTreeType = Tree<KDTreePartition<BucketType>>;
 
    using MortarKinematicVariablesType = MortarKinematicVariables<TNumNodes, TNumNodesMaster>;
    using MortarOperatorType = MortarOperator<TNumNodes, TNumNodesMaster>;
    using DualLagrangeMultiplierOperatorsType = DualLagrangeMultiplierOperators<TNumNodes, TNumNodesMaster>;
    using ExactMortarIntegrationUtilityType = ExactMortarIntegrationUtility<TDim, TNumNodes, false, TNumNodesMaster>;
 
    struct TLS {
        MortarKinematicVariablesType this_kinematic_variables;    // Create and initialize condition variables:
        MortarOperatorType this_mortar_operators;                 // Create the mortar operators
        ExactMortarIntegrationUtilityType integration_utility;    // We call the exact integration utility
    };
 
 
    SimpleMortarMapperProcess(
        ModelPart& rOriginModelPart,
        ModelPart& rDestinationModelPart,
        Parameters ThisParameters = Parameters(R"({})" ),
        LinearSolverType::Pointer pThisLinearSolver = nullptr
        );
 
    SimpleMortarMapperProcess(
        ModelPart& rOriginModelPart,
        ModelPart& rDestinationModelPart,
        TVarType& rThisVariable,
        Parameters ThisParameters = Parameters(R"({})" ),
        LinearSolverType::Pointer pThisLinearSolver = nullptr
        );
 
    SimpleMortarMapperProcess(
        ModelPart& rOriginModelPart,
        ModelPart& rDestinationModelPart,
        TVarType& rOriginVariable,
        TVarType& rDestinationVariable,
        Parameters ThisParameters = Parameters(R"({})" ),
        LinearSolverType::Pointer pThisLinearSolver = nullptr
        );
 
    ~SimpleMortarMapperProcess() override = default;
 
 
 
 
 
 
    void operator()()
    {
        Execute();
    }
 
 
    void Execute() override;
 
    void ExecuteInitializeSolutionStep() override;
 
    void Map(
        TVarType& rOriginVariable,
        TVarType& rDestinationVariable,
        const Flags Flag = Flags()
        )
    {
        // Reassign the variables
        mpOriginVariable = &rOriginVariable;
        mpDestinationVariable = &rDestinationVariable;
 
        // Execute the process
        Execute();
    }
 
    const Parameters GetDefaultParameters() const override;
 
 
 
 
    std::string Info() const override
    {
        return "SimpleMortarMapperProcess";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "SimpleMortarMapperProcess";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
    }
 
 
private:
 
    ModelPart& mOriginModelPart;                  
    ModelPart& mDestinationModelPart;             
    const TVarType* mpOriginVariable;             
    const TVarType* mpDestinationVariable;        
 
    double mMappingCoefficient = 1.0;             
    Flags mOptions;                               
 
    unsigned int mEchoLevel = 0;                  
    Parameters mThisParameters;                   
 
    LinearSolverType::Pointer mpThisLinearSolver; 
 
 
 
    void CheckAndPerformSearch();
 
    void ResetNodalArea();
 
    double GetReferenceArea();
 
    void AssemblyMortarOperators(
        const std::vector<array_1d<PointType,TDim>>& rGeometricalObjectsPointSlave,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const array_1d<double, 3>& rMasterNormal,
        MortarKinematicVariablesType& rThisKinematicVariables,
        MortarOperatorType& rThisMortarOperators,
        const IntegrationMethod& rThisIntegrationMethod,
        const BoundedMatrixType Ae = IdentityMatrix(TNumNodes)
        );
 
    static inline BoundedMatrixType CalculateAe(
        GeometryType& rSlaveGeometry,
        MortarKinematicVariablesType& rThisKinematicVariables,
        std::vector<array_1d<PointType,TDim>>& rConditionsPointsSlave,
        const IntegrationMethod& rThisIntegrationMethod
        );
 
    static inline BoundedMatrixType InvertDiagonalMatrix(const BoundedMatrixType& rInputMatrix);
 
    static inline void InvertDiagonalMatrix(
        const BoundedMatrixType& rInputMatrix,
        BoundedMatrixType& rInvertedMatrix
        );
 
    void LumpMatrix(BoundedMatrixType& rInputMatrix);
 
    void GetSystemSize(SizeType& rSizeSystem);
 
    void CreateSlaveConectivityDatabase(
        SizeType& rSizeSystem,
        IntMap& rConectivityDatabase,
        IntMap& rInverseConectivityDatabase
        );
 
    IntegrationMethod GetIntegrationMethod();
 
    bool CheckWholeVector(std::vector<bool>& rVectorToCheck);
 
    void ComputeResidualMatrix(
        Matrix& rResidualMatrix,
        const GeometryType& rSlaveGeometry,
        const GeometryType& rMasterGeometry,
        const MortarOperatorType& rThisMortarOperators
        );
 
    void AssembleRHSAndLHS(
        MatrixType& rA,
        std::vector<VectorType>& rb,
        const SizeType VariableSize,
        const Matrix& rResidualMatrix,
        const GeometryType& rSlaveGeometry,
        IntMap& rInverseConectivityDatabase,
        const MortarOperatorType& rThisMortarOperators
        );
 
    void AssembleRHS(
        std::vector<VectorType>& rb,
        const SizeType VariableSize,
        const Matrix& rResidualMatrix,
        const GeometryType& rSlaveGeometry,
        IntMap& rInverseConectivityDatabase
        );
 
    void ExecuteExplicitMapping();
 
    void ExecuteImplicitMapping();
 
    template<class TClassType, bool TImplicit = false>
    void PerformMortarOperations(
        MatrixType& rA,
        std::vector<VectorType>& rb,
        IntMap& rInverseConectivityDatabase,
        typename TClassType::Pointer pIndexesPairs,
        GeometricalObject& rGeometricalObject,
        ExactMortarIntegrationUtilityType& rIntegrationUtility,
        MortarKinematicVariablesType& rThisKineticVariables,
        MortarOperatorType& rThisMortarOperators,
        const IndexType Iteration
        )
    {
        // The root model part
        ModelPart& r_root_model_part = mOriginModelPart.GetRootModelPart();
 
        // Getting the auxiliar variable
        const TVarType& r_aux_variable = KratosComponents<TVarType>::Get(MortarUtilities::GetAuxiliarVariable<TVarType>());
 
        // Indexes of the pair to be removed
        std::vector<IndexType> indexes_to_remove, geometrical_objects_to_erase;
 
        // Getting discontinous factor
        const double discontinous_interface_factor = mOptions.Is(DISCONTINOUS_INTERFACE) ? mThisParameters["discontinous_interface_factor"].GetDouble() : 1.0;
 
        // Declare auxiliar coordinates
        GeometryType::CoordinatesArrayType aux_coords;
 
        // Geometrical values
        auto& r_slave_geometry = rGeometricalObject.GetGeometry();
        r_slave_geometry.PointLocalCoordinates(aux_coords, r_slave_geometry.Center());
        const array_1d<double, 3> slave_normal = r_slave_geometry.UnitNormal(aux_coords);
 
        // The model part as const to avoid race conditions
        const auto& r_const_origin_model_part = mOriginModelPart;
        const auto& r_const_destination_model_part = mDestinationModelPart;
 
        for (auto it_pair = pIndexesPairs->begin(); it_pair != pIndexesPairs->end(); ++it_pair ) {
            const IndexType master_id = pIndexesPairs->GetId(it_pair); // MASTER
 
            const auto& r_master_geometry = mOptions.Is(ORIGIN_SKIN_IS_CONDITION_BASED) ? r_const_origin_model_part.GetCondition(master_id).GetGeometry() : r_const_origin_model_part.GetElement(master_id).GetGeometry();
            r_master_geometry.PointLocalCoordinates(aux_coords, r_master_geometry.Center());
            const array_1d<double, 3> master_normal = r_master_geometry.UnitNormal(aux_coords);
 
            const IntegrationMethod& r_integration_method = GetIntegrationMethod();
 
            // Reading integration points
            std::vector<array_1d<PointType,TDim>> geometrical_objects_points_slave; // These are the segmentation points, with this points it is possible to create the lines or triangles used on the mapping
            const bool is_inside = rIntegrationUtility.GetExactIntegration(r_slave_geometry, slave_normal, r_master_geometry, master_normal, geometrical_objects_points_slave);
 
            if (is_inside) {
                // Initialize general variables for the current master element
                rThisKineticVariables.Initialize();
 
                // Initialize the mortar operators
                rThisMortarOperators.Initialize();
 
                const BoundedMatrixType Ae = CalculateAe(r_slave_geometry, rThisKineticVariables, geometrical_objects_points_slave, r_integration_method);
 
                AssemblyMortarOperators( geometrical_objects_points_slave, r_slave_geometry, r_master_geometry,master_normal, rThisKineticVariables, rThisMortarOperators, r_integration_method, Ae);
 
                /* We compute the residual */
                const IndexType size_to_compute = MortarUtilities::SizeToCompute<TDim, TVarType>();
                Matrix residual_matrix(TNumNodes, size_to_compute);
                ComputeResidualMatrix(residual_matrix, r_slave_geometry, r_master_geometry, rThisMortarOperators);
 
                if (!TImplicit) {
                    MortarUtilities::AddValue<TVarType, MortarUtilitiesSettings::SaveAsNonHistoricalVariable>(r_slave_geometry, r_aux_variable, residual_matrix);
                }
 
                // We check if DOperator is diagonal
                if (mEchoLevel > 1) {
                    BoundedMatrixType aux_copy_D = rThisMortarOperators.DOperator;
                    LumpMatrix(aux_copy_D);
                    const BoundedMatrixType aux_diff = aux_copy_D - rThisMortarOperators.DOperator;
                    const double norm_diff = norm_frobenius(aux_diff);
                    if (norm_diff > 1.0e-4)
                        KRATOS_WARNING("D OPERATOR") << " THE MORTAR OPERATOR D IS NOT DIAGONAL" << std::endl;
                    if (mEchoLevel == 3) {
                        KRATOS_WATCH(norm_diff);
                        KRATOS_WATCH(rThisMortarOperators.DOperator);
                    }
                }
 
                if (Iteration == 0) { // Just assembled the first iteration
                    if constexpr (TImplicit) {
                        /* We compute the residual and assemble */
                        const SizeType variable_size = MortarUtilities::SizeToCompute<TDim, TVarType>();
                        AssembleRHSAndLHS(rA, rb, variable_size, residual_matrix, r_slave_geometry, rInverseConectivityDatabase, rThisMortarOperators);
                    } else {
                        for (IndexType i_node = 0; i_node < TNumNodes; ++i_node) {
                            double& r_nodal_area = r_slave_geometry[i_node].GetValue(NODAL_AREA);
                            AtomicAdd(r_nodal_area, rThisMortarOperators.DOperator(i_node, i_node));
                        }
                        // In case of discontinous interface we add contribution to near nodes
                        if (mOptions.Is(DISCONTINOUS_INTERFACE)) {
                            const double element_length = r_slave_geometry.Length();
 
                            // Iterating over nodes
                            for (IndexType i_node = 0; i_node < TNumNodes; ++i_node) {
                                const double nodal_area_contribution = rThisMortarOperators.DOperator(i_node, i_node);
 
                                // The original node coordinates
                                const auto& r_slave_node = r_slave_geometry[i_node];
 
                                // Iterating over other paired geometrical objects
                                const auto& r_index_masp_master = mOptions.Is(ORIGIN_SKIN_IS_CONDITION_BASED) ? r_const_origin_model_part.GetCondition(master_id).GetValue(INDEX_SET) : r_const_origin_model_part.GetElement(master_id).GetValue(INDEX_SET);
                                for (auto it_master_pair = r_index_masp_master->begin(); it_master_pair != r_index_masp_master->end(); ++it_master_pair ) {
 
                                    const IndexType auxiliar_slave_id = r_index_masp_master->GetId(it_master_pair);
                                    if (rGeometricalObject.Id() != auxiliar_slave_id) {
                                        GeometryType& r_auxiliar_slave_geometry =  const_cast<GeometryType&>(mOptions.Is(DESTINATION_SKIN_IS_CONDITION_BASED) ? r_const_destination_model_part.GetCondition(auxiliar_slave_id).GetGeometry() : r_const_destination_model_part.GetElement(auxiliar_slave_id).GetGeometry());
 
                                        for (IndexType j_node = 0; j_node < TNumNodes; ++j_node) {
                                            // The auxiliary node distance
                                            auto& r_auxiliary_slave_node = r_auxiliar_slave_geometry[j_node];
                                            const double distance = r_auxiliary_slave_node.Distance(r_slave_node);
                                            const double contribution_coeff = 1.0/std::pow((1.0 + distance/(discontinous_interface_factor * element_length)), 2);
 
                                            double& r_nodal_area = r_auxiliary_slave_node.GetValue(NODAL_AREA);
                                            AtomicAdd(r_nodal_area, contribution_coeff * nodal_area_contribution);
                                        }
                                    }
                                }
                            }
                        }
                    }
                } else if constexpr (TImplicit) {
                    const SizeType variable_size = MortarUtilities::SizeToCompute<TDim, TVarType>();
                    AssembleRHS(rb, variable_size, residual_matrix, r_slave_geometry, rInverseConectivityDatabase);
                }
            } else { // NOTE: The geometrical object considered maybe is to tight
                indexes_to_remove.push_back(master_id);
                const IndexType other_id = pIndexesPairs->GetOtherId(it_pair);
                if (std::is_same<TClassType, IndexMap>::value && other_id != 0) {
                    geometrical_objects_to_erase.push_back(other_id);
                }
            }
        }
 
        // Clear indexes
        for (IndexType i_to_remove = 0; i_to_remove < indexes_to_remove.size(); ++i_to_remove) {
            if (mOptions.Is(ORIGIN_SKIN_IS_CONDITION_BASED)) {
                for (auto& id : geometrical_objects_to_erase ) {
                    auto& r_cond = r_root_model_part.GetCondition(id);
                    r_cond.Set(TO_ERASE, true);
                }
            } else {
                for (auto& id : geometrical_objects_to_erase ) {
                    auto& r_elem = r_root_model_part.GetElement(id);
                    r_elem.Set(TO_ERASE, true);
                }
            }
            pIndexesPairs->RemoveId(indexes_to_remove[i_to_remove]);
        }
    }
 
    template<class TClassType>
    void ClearIndexes(
        typename TClassType::Pointer pIndexesPairs,
        GeometricalObject& rGeometricalObject,
        ExactMortarIntegrationUtilityType& rIntegrationUtility
        )
    {
        // The root model part
        ModelPart& r_root_model_part = mOriginModelPart.GetRootModelPart();
 
        // Indexes of the pair to be removed
        std::vector<IndexType> indexes_to_remove, geometrical_objects_to_erase;
 
        // Declare auxiliar coordinates
        GeometryType::CoordinatesArrayType aux_coords;
 
        // Geometrical values
        auto& r_slave_geometry = rGeometricalObject.GetGeometry();
        r_slave_geometry.PointLocalCoordinates(aux_coords, r_slave_geometry.Center());
        const array_1d<double, 3> slave_normal = r_slave_geometry.UnitNormal(aux_coords);
 
        // The model part as const to avoid race conditions
        const auto& r_const_origin_model_part = mOriginModelPart;
 
        for (auto it_pair = pIndexesPairs->begin(); it_pair != pIndexesPairs->end(); ++it_pair ) {
            const IndexType master_id = pIndexesPairs->GetId(it_pair); // MASTER
 
            const auto& r_master_geometry = mOptions.Is(ORIGIN_SKIN_IS_CONDITION_BASED) ? r_const_origin_model_part.GetCondition(master_id).GetGeometry() : r_const_origin_model_part.GetElement(master_id).GetGeometry();
            r_master_geometry.PointLocalCoordinates(aux_coords, r_master_geometry.Center());
            const array_1d<double, 3> master_normal = r_master_geometry.UnitNormal(aux_coords);
 
            // Reading integration points
            std::vector<array_1d<PointType,TDim>> geometrical_objects_points_slave; // These are the segmentation points, with this points it is possible to create the lines or triangles used on the mapping
            const bool is_inside = rIntegrationUtility.GetExactIntegration(r_slave_geometry, slave_normal, r_master_geometry, master_normal, geometrical_objects_points_slave);
 
            if (!is_inside) {
                indexes_to_remove.push_back(master_id);
                const IndexType other_id = pIndexesPairs->GetOtherId(it_pair);
                if (std::is_same<TClassType, IndexMap>::value && other_id != 0) {
                    geometrical_objects_to_erase.push_back(other_id);
                }
            }
        }
 
        // Clear indexes
        for (IndexType i_to_remove = 0; i_to_remove < indexes_to_remove.size(); ++i_to_remove) {
            if (mOptions.Is(ORIGIN_SKIN_IS_CONDITION_BASED)) {
                for (auto& id : geometrical_objects_to_erase ) {
                    auto& r_cond = r_root_model_part.GetCondition(id);
                    r_cond.Set(TO_ERASE, true);
                }
            } else {
                for (auto& id : geometrical_objects_to_erase ) {
                    auto& r_elem = r_root_model_part.GetElement(id);
                    r_elem.Set(TO_ERASE, true);
                }
            }
            pIndexesPairs->RemoveId(indexes_to_remove[i_to_remove]);
        }
    }
 
    template<class TEntity>
    void FillDatabase(
        TEntity& rGeometricalObject,
        KDTreeType& rTreePoints,
        const SizeType AllocationSize,
        const double SearchFactor
        )
    {
        // Initialize values
        PointVector points_found(AllocationSize);
 
        GeometryType& r_geometry = rGeometricalObject.GetGeometry();
        const Point center = r_geometry.Center();
 
        double radius = 0.0;
        for(IndexType i_node = 0; i_node < r_geometry.PointsNumber(); ++i_node)  {
            const array_1d<double, 3> aux_vector = center.Coordinates() - r_geometry[i_node].Coordinates();
            const double aux_value = inner_prod(aux_vector, aux_vector);
            if(aux_value > radius) radius = aux_value;
        }
 
        const double search_radius = SearchFactor * std::sqrt(radius);
 
        const SizeType number_points_found = rTreePoints.SearchInRadius(center, search_radius, points_found.begin(), AllocationSize);
 
        if (number_points_found > 0) {
            // In case of missing is created
            if (!rGeometricalObject.Has(INDEX_SET))
                rGeometricalObject.SetValue(INDEX_SET, Kratos::make_shared<IndexSet>());
 
            // Accessing to the index set
            IndexSet::Pointer indexes_set = rGeometricalObject.GetValue(INDEX_SET);
 
            for (IndexType i_point = 0; i_point < number_points_found; ++i_point ) {
                auto p_geometrical_object_master = points_found[i_point]->pGetObject();
                indexes_set->AddId(p_geometrical_object_master->Id());
            }
        }
    }
 
    void CreateInverseDatabase();
 
    void UpdateInterface();
 
 
 
 
 
    SimpleMortarMapperProcess& operator=(SimpleMortarMapperProcess const& rOther) = delete;
 
    //SimpleMortarMapperProcess(SimpleMortarMapperProcess const& rOther);
 
};// class SimpleMortarMapperProcess
 
 
 
 
 
}  // namespace Kratos.