normal_calculation_utils.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Pooyan Dadvand
//                   Riccardo Rossi
//                   Vicente Mataix Ferrandiz
//                   Ruben Zorrilla
//
 
#if !defined(KRATOS_NORMAL_CALCULATION_UTILS )
#define  KRATOS_NORMAL_CALCULATION_UTILS
 
// System includes
 
// External includes
 
// Project includes
#include "includes/model_part.h"
#include "utilities/math_utils.h"
#include "utilities/variable_utils.h"
#include "includes/deprecated_variables.h"
 
namespace Kratos
{
 
 
 
 
 
 
class KRATOS_API(KRATOS_CORE) NormalCalculationUtils
{
public:
 
    typedef std::size_t IndexType;
 
    typedef std::size_t SizeType;
 
    // Node definitions
    typedef Node NodeType;
 
    typedef Geometry<NodeType> GeometryType;
 
    typedef ModelPart::ConditionType ConditionType;
 
    typedef ModelPart::ConditionsContainerType ConditionsArrayType;
 
    using NormalVariableType = Variable<array_1d<double,3>>;
 
 
 
 
    template<class TContainerType>
    void CalculateNormalsInContainer(
        ModelPart& rModelPart,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    template<class TContainerType, bool TIsHistorical = true>
    void CalculateNormals(
        ModelPart& rModelPart,
        const bool EnforceGenericGeometryAlgorithm = false,
        const bool ConsiderUnitNormal = false,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    template<class TContainerType, bool TIsHistorical = true>
    void CalculateUnitNormals(
        ModelPart& rModelPart,
        const bool EnforceGenericGeometryAlgorithm = false,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateOnSimplex(
        ConditionsArrayType& rConditions,
        const std::size_t Dimension,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateOnSimplexNonHistorical(
        ConditionsArrayType& rConditions,
        const std::size_t Dimension,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateNormalShapeDerivativesOnSimplex(
        ConditionsArrayType& rConditions,
        const std::size_t Dimension
    );
 
    void CalculateOnSimplex(
        ModelPart& rModelPart,
        const std::size_t Dimension,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateOnSimplexNonHistorical(
        ModelPart& rModelPart,
        const std::size_t Dimension,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateOnSimplex(
        ModelPart& rModelPart,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void CalculateOnSimplexNonHistorical(
        ModelPart& rModelPart,
        const NormalVariableType& rNormalVariable = NORMAL
        );
 
    void SwapNormals(ModelPart& rModelPart);
 
    template<class TValueType>
    void CalculateOnSimplex(
        ModelPart& rModelPart,
        const std::size_t Dimension,
        const Variable<TValueType>& rVariable,
        const TValueType Zero,
        const NormalVariableType& rNormalVariable = NORMAL
        )
    {
        KRATOS_TRY;
 
        // Reset normals
        //TODO: This can be parallel
        const array_1d<double,3> ZeroNormal(3,0.0);
        for(ModelPart::NodesContainerType::iterator it = rModelPart.NodesBegin(); it !=rModelPart.NodesEnd(); it++) {
            noalias(it->FastGetSolutionStepValue(rNormalVariable)) = ZeroNormal;
        }
 
        // Calculate new condition normals, using only conditions with rVariable == rValue
        array_1d<double,3> An(3,0.0);
 
        if ( Dimension == 2 ) {
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal2D(*itCond,An,rNormalVariable);
            }
        } else if ( Dimension == 3 ) {
            array_1d<double,3> v1(3,0.0);
            array_1d<double,3> v2(3,0.0);
 
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal3D(*itCond,An,v1,v2,rNormalVariable);
            }
        }
 
        // Transfer normals to nodes
        for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
            Condition::GeometryType& rGeom = itCond->GetGeometry();
            const double Coef = 1.0 / rGeom.PointsNumber();
            const auto& r_normal = itCond->GetValue(rNormalVariable);
            for ( Condition::GeometryType::iterator itNode = rGeom.begin(); itNode != rGeom.end(); ++itNode)
                noalias(itNode->FastGetSolutionStepValue(rNormalVariable)) += r_normal * Coef;
        }
 
        // For MPI: correct values on partition boundaries
        rModelPart.GetCommunicator().AssembleCurrentData(rNormalVariable);
 
        KRATOS_CATCH("");
    }
 
 
    template<class TValueType>
    void CalculateOnSimplex(
        ModelPart& rModelPart,
        const std::size_t Dimension,
        const Variable<TValueType>& rVariable
        )
    {
        CalculateOnSimplex(rModelPart,Dimension,rVariable,TValueType(),NORMAL);
    }
 
    template<class TValueType>
    void CalculateOnSimplex(
        ModelPart& rModelPart,
        const std::size_t Dimension,
        const Variable<TValueType>& rVariable,
        const TValueType Zero,
        const double rAlpha,
        const NormalVariableType& rNormalVariable = NORMAL
        )
    {
        KRATOS_TRY;
 
        // Reset normals
        //TODO: This can be parallel
        const array_1d<double,3> ZeroNormal(3,0.0);
        for(ModelPart::NodesContainerType::iterator it =  rModelPart.NodesBegin(); it !=rModelPart.NodesEnd(); it++) {
            noalias(it->FastGetSolutionStepValue(NORMAL)) = ZeroNormal;
            it->FastGetSolutionStepValue(NODAL_PAUX) = 0.0;
        }
 
        // Calculate new condition normals, using only conditions with rVariable == rValue
        array_1d<double,3> An(3,0.0);
 
        if ( Dimension == 2 ) {
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal2D(*itCond,An,rNormalVariable);
            }
        } else if ( Dimension == 3 ) {
            array_1d<double,3> v1(3,0.0);
            array_1d<double,3> v2(3,0.0);
 
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal3D(*itCond,An,v1,v2,rNormalVariable);
            }
        }
 
        // Loop over nodes to set normals
        for(ModelPart::NodesContainerType::iterator it =  rModelPart.NodesBegin(); it !=rModelPart.NodesEnd(); it++) {
            std::vector< array_1d<double,3> > N_Mat;
            N_Mat.reserve(10);
            double nodal_area = 0.0;
 
            GlobalPointersVector<Condition >& ng_cond = it->GetValue(NEIGHBOUR_CONDITIONS);
 
            if(ng_cond.size() != 0) {
                for(GlobalPointersVector<Condition >::iterator ic = ng_cond.begin(); ic!=ng_cond.end(); ic++) {
                    Condition::GeometryType& pGeom = ic->GetGeometry();
                    const auto& rNormal = ic->GetValue(rNormalVariable);
                    const double Coef = 1.0 / pGeom.PointsNumber();
                    double norm_normal = norm_2( rNormal );
 
                    if(norm_normal != 0.0) {
                        nodal_area += Coef * norm_normal;
 
                        if(N_Mat.size() == 0.0) {
                            N_Mat.push_back( rNormal * Coef );
                        } else {
                            int added = 0;
                            for(unsigned int ii=0; ii<N_Mat.size();++ii) {
                                const array_1d<double,3>& temp_normal = N_Mat[ii];
                                double norm_temp = norm_2( temp_normal );
 
                                double cos_alpha=temp_normal[0]*rNormal[0] + temp_normal[1]*rNormal[1] +temp_normal[2]*rNormal[2];
                                cos_alpha /= (norm_temp*norm_normal);
 
                                if( cos_alpha > std::cos(0.017453293*rAlpha)) {
                                    N_Mat[ii] += rNormal * Coef;
                                    added = 1;
                                }
                            }
 
                            if(!added)
                                N_Mat.push_back( rNormal*Coef );
                        }
                    }
                }
            }
 
            // Compute NORMAL and mark
            array_1d<double,3> sum_Normal(3,0.0);
 
            for(unsigned int ii=0; ii<N_Mat.size(); ++ii) {
                sum_Normal += N_Mat[ii];
            }
 
            noalias( it->FastGetSolutionStepValue(rNormalVariable) ) = sum_Normal;
            it->FastGetSolutionStepValue(NODAL_PAUX) = nodal_area;
            //assign IS_SLIP = 0 for vertices
            if(N_Mat.size() == 2) {
//                 it->SetValue(IS_SLIP,0);
                it->FastGetSolutionStepValue(IS_SLIP)=20.0;
            } else if(N_Mat.size() == 3) {
                it->FastGetSolutionStepValue(IS_SLIP)=30.0;
            } else if(N_Mat.size() == 1) {
                it->FastGetSolutionStepValue(IS_SLIP)=10.0;
            }
        }
 
        // For MPI: correct values on partition boundaries
        rModelPart.GetCommunicator().AssembleCurrentData(rNormalVariable);
        rModelPart.GetCommunicator().AssembleCurrentData(NODAL_PAUX);
 
        KRATOS_CATCH("");
    }
 
    template< class TValueType >
    void CalculateOnSimplexLowMemory(
        ModelPart& rModelPart,
        int Dimension,
        const Variable<TValueType>& rVariable,
        const TValueType Zero,const double rAlpha,
        const NormalVariableType& rNormalVariable = NORMAL
        )
    {
        KRATOS_TRY;
 
        // Reset normals
        //TODO: This can be parallel
        const array_1d<double,3> ZeroNormal(3,0.0);
        for(ModelPart::NodesContainerType::iterator it =  rModelPart.NodesBegin(); it !=rModelPart.NodesEnd(); it++) {
            noalias(it->GetValue(rNormalVariable)) = ZeroNormal;
            it->GetValue(NODAL_PAUX) = 0.0;
        }
 
        // Calculate new condition normals, using only conditions with rVariable == rValue
        array_1d<double,3> An(3,0.0);
 
        if ( Dimension == 2 ) {
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal2D(*itCond,An,rNormalVariable);
            }
        } else if ( Dimension == 3 ) {
            array_1d<double,3> v1(3,0.0);
            array_1d<double,3> v2(3,0.0);
 
            for ( ModelPart::ConditionIterator itCond = rModelPart.ConditionsBegin(); itCond != rModelPart.ConditionsEnd(); ++itCond ) {
                if ( itCond->GetValue(rVariable) != Zero )
                    CalculateNormal3D(*itCond,An,v1,v2,rNormalVariable);
            }
        }
 
        // Loop over nodes to set normals
        for(ModelPart::NodesContainerType::iterator it =  rModelPart.NodesBegin(); it !=rModelPart.NodesEnd(); it++) {
            std::vector< array_1d<double,3> > N_Mat;
            N_Mat.reserve(10);
            double nodal_area = 0.0;
 
            GlobalPointersVector<Condition >& ng_cond = it->GetValue(NEIGHBOUR_CONDITIONS);
 
            if(ng_cond.size() != 0){
                for(GlobalPointersVector<Condition >::iterator ic = ng_cond.begin(); ic!=ng_cond.end(); ic++) {
                Condition::GeometryType& pGeom = ic->GetGeometry();
                const auto& rNormal = ic->GetValue(rNormalVariable);
                const double Coef = 1.0 / pGeom.PointsNumber();
                double norm_normal = norm_2( rNormal );
 
                if(norm_normal != 0.0) {
                    nodal_area += Coef * norm_normal;
 
                    if(N_Mat.size() == 0.0) {
                        N_Mat.push_back( rNormal * Coef );
                    } else{
                        int added = 0;
                        for(unsigned int ii=0; ii<N_Mat.size();++ii) {
                            const array_1d<double,3>& temp_normal = N_Mat[ii];
                            double norm_temp = norm_2( temp_normal );
 
                            double cos_alpha=temp_normal[0]*rNormal[0] + temp_normal[1]*rNormal[1] +temp_normal[2]*rNormal[2];
                            cos_alpha /= (norm_temp*norm_normal);
 
                            if( cos_alpha > cos(0.017453293*rAlpha) ){
                                N_Mat[ii] += rNormal * Coef;
                                added = 1;}
                            }
                            if(!added)
                                N_Mat.push_back( rNormal*Coef );
 
                        }
                    }
                }
            }
 
            // Compute NORMAL and mark
            array_1d<double,3> sum_Normal(3,0.0);
 
            for(unsigned int ii=0; ii<N_Mat.size(); ++ii){
                sum_Normal += N_Mat[ii];
            }
 
            noalias( it->FastGetSolutionStepValue(rNormalVariable) ) = sum_Normal;
            it->FastGetSolutionStepValue(NODAL_PAUX) = nodal_area;
            //assign IS_SLIP = 0 for vertices
            if(N_Mat.size() == 2){
//                 it->SetValue(IS_SLIP,0);
                it->FastGetSolutionStepValue(IS_SLIP)=20.0;
            } else if(N_Mat.size() == 3) {
                it->FastGetSolutionStepValue(IS_SLIP)=30.0;
            } else if(N_Mat.size() == 1) {
                it->FastGetSolutionStepValue(IS_SLIP)=10.0;
 
            }
        }
 
        // For MPI: correct values on partition boundaries
        rModelPart.GetCommunicator().AssembleCurrentData(rNormalVariable);
        rModelPart.GetCommunicator().AssembleCurrentData(NODAL_PAUX);
 
        KRATOS_CATCH("");
    }
 
private:
 
 
 
 
    template<class TContainerType, bool TIsHistorical>
    void InitializeNormals(
        ModelPart& rModelPart,
        const NormalVariableType& rNormalVariable
        );
 
    template<class TContainerType, bool TIsHistorical>
    void ComputeUnitNormalsFromAreaNormals(
        ModelPart& rModelPart,
        const NormalVariableType& rNormalVariable
        );
 
    static void CalculateNormal2D(
        Condition& rCondition,
        array_1d<double,3>& rAn,
        const NormalVariableType& rNormalVariable
        );
 
    static void CalculateNormalShapeDerivative2D(
        ConditionType& rCondition
        );
 
    static void CalculateNormal3D(
        Condition& rCondition,
        array_1d<double,3>& rAn,
        array_1d<double,3>& rv1,
        array_1d<double,3>& rv2,
        const NormalVariableType& rNormalVariable
        );
 
    static void CalculateNormalShapeDerivative3D(
        ConditionType& rCondition
        );
 
    template<class TContainerType>
    TContainerType& GetContainer(ModelPart& rModelPart);
 
    template<class TContainerType, bool TIsHistorical>
    void CalculateNormalsUsingGenericAlgorithm(
        ModelPart& rModelPart,
        const bool ConsiderUnitNormal,
        const NormalVariableType& rNormalVariable
        );
 
    template<bool TIsHistorical>
    array_1d<double,3>& GetNormalValue(
        NodeType& rNode,
        const NormalVariableType& rNormalVariable
        );
 
    template<bool TIsHistorical>
    void SetNormalValue(
        NodeType& rNode,
        const NormalVariableType& rNormalVariable,
        const array_1d<double,3>& rNormalValue
        );
 
    bool CheckUseSimplex(
        const ModelPart& rModelPart,
        const bool EnforceGenericGeometryAlgorithm
        );
 
    template<bool TIsHistorical>
    void AuxiliaryCalculateOnSimplex(
        ConditionsArrayType& rConditions,
        const std::size_t Dimension,
        const NormalVariableType& rNormalVariable
        );
 
 
 
 
 
}; // Class NormalCalculationUtils
 
} // namespace Kratos
 
#endif /* KRATOS_NORMAL_CALCULATION_UTILS  defined */