internal_variables_interpolation_process.h

Go to the documentation of this file.

// KRATOS  __  __ _____ ____  _   _ ___ _   _  ____
//        |  \/  | ____/ ___|| | | |_ _| \ | |/ ___|
//        | |\/| |  _| \___ \| |_| || ||  \| | |  _
//        | |  | | |___ ___) |  _  || || |\  | |_| |
//        |_|  |_|_____|____/|_| |_|___|_| \_|\____| APPLICATION
//
//  License:         BSD License
//                       license: MeshingApplication/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
#if !defined(KRATOS_INTERNAL_VARIABLES_INTERPOLATION_PROCESS )
#define  KRATOS_INTERNAL_VARIABLES_INTERPOLATION_PROCESS
 
// System includes
 
// External includes
 
// Project includes
#include "meshing_application_variables.h"
#include "processes/process.h"
#include "includes/model_part.h"
#include "includes/kratos_parameters.h"
#include "includes/kratos_components.h"
#include "custom_includes/gauss_point_item.h"
#include "utilities/atomic_utilities.h"
#include "utilities/parallel_utilities.h"
// Include the point locator
#include "utilities/binbased_fast_point_locator.h"
// Include the trees
// #include "spatial_containers/bounding_volume_tree.h" // k-DOP
#include "spatial_containers/spatial_containers.h" // kd-tree
 
namespace Kratos
{
 
 
    typedef Node NodeType;
 
    typedef Geometry<NodeType> GeometryType;
 
    typedef GaussPointItem                                        PointType;
    typedef PointType::Pointer                             PointTypePointer;
    typedef std::vector<PointTypePointer>                       PointVector;
    typedef PointVector::iterator                             PointIterator;
    typedef std::vector<double>                              DistanceVector;
    typedef DistanceVector::iterator                       DistanceIterator;
 
 
 
 
class KRATOS_API(MESHING_APPLICATION) InternalVariablesInterpolationProcess
    : public Process
{
public:
 
    typedef Bucket< 3ul, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator > BucketType;
    typedef Tree< KDTreePartition<BucketType> > KDTree;
 
    typedef Variable<double>             DoubleVarType;
    typedef Variable<array_1d<double, 3>> ArrayVarType;
    typedef Variable<Vector>             VectorVarType;
    typedef Variable<Matrix>             MatrixVarType;
 
    // General type definitions
    typedef ModelPart::NodesContainerType                    NodesArrayType;
    typedef ModelPart::ElementsContainerType              ElementsArrayType;
    typedef ModelPart::ConditionsContainerType          ConditionsArrayType;
    typedef Node                                                NodeType;
    typedef Geometry<NodeType>                                 GeometryType;
 
    KRATOS_CLASS_POINTER_DEFINITION( InternalVariablesInterpolationProcess );
 
 
    enum class InterpolationTypes {
        CLOSEST_POINT_TRANSFER = 0, 
        LEAST_SQUARE_TRANSFER = 1, 
        SHAPE_FUNCTION_TRANSFER = 2  
        };
 
 
    // Class Constructor
 
    InternalVariablesInterpolationProcess(
        ModelPart& rOriginMainModelPart,
        ModelPart& rDestinationMainModelPart,
        Parameters ThisParameters =  Parameters(R"({})")
        );
 
    ~InternalVariablesInterpolationProcess() override= default;
 
 
    void operator()()
    {
        Execute();
    }
 
 
    void Execute() override;
 
    const Parameters GetDefaultParameters() const override;
 
 
 
 
    /************************************ GET INFO *************************************/
    /***********************************************************************************/
 
    std::string Info() const override
    {
        return "InternalVariablesInterpolationProcess";
    }
 
    /************************************ PRINT INFO ***********************************/
    /***********************************************************************************/
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
private:
 
    // Auxiliar search struct
    template<std::size_t TDim>
    struct auxiliar_search
    {
        auxiliar_search(ModelPart& rModelPart)
            : point_locator(rModelPart)
        {
            // We create the locator
            point_locator.UpdateSearchDatabase();
        }
 
        BinBasedFastPointLocator<TDim> point_locator;
        Vector N;
        Element::Pointer p_element;
    };
 
 
 
    // The model parts
    ModelPart& mrOriginMainModelPart;                    
    ModelPart& mrDestinationMainModelPart;               
    const std::size_t mDimension;                        
 
    // The allocation parameters
    std::size_t mAllocationSize;                    
    std::size_t mBucketSize;                        
 
    // The seatch variables
    double mSearchFactor;                           
    PointVector mPointListOrigin;                   
 
    // Variables to interpolate
    std::vector<std::string> mInternalVariableList; 
    InterpolationTypes mThisInterpolationType;      
 
 
 
    PointVector CreateGaussPointList(ModelPart& ThisModelPart);
 
    void InterpolateGaussPointsClosestPointTransfer();
 
    void InterpolateGaussPointsLeastSquareTransfer();
 
    void InterpolateGaussPointsShapeFunctionTransfer();
 
    std::size_t ComputeTotalNumberOfVariables();
 
    template<class TVarType>
    static inline void SaveValuesOnGaussPoint(
        const Variable<TVarType>& rThisVar,
        PointTypePointer pPointOrigin,
        ElementsArrayType::iterator itElemOrigin,
        const IndexType GaussPointId,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        std::vector<TVarType> values;
        itElemOrigin->CalculateOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
        pPointOrigin->SetValue(rThisVar, values[GaussPointId]);
    }
 
    template<class TVarType>
    static inline void GetAndSetDirectVariableOnConstitutiveLaw(
        const Variable<TVarType>& rThisVar,
        ConstitutiveLaw::Pointer pOriginConstitutiveLaw,
        ConstitutiveLaw::Pointer pDestinationConstitutiveLaw,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        TVarType origin_value;
        origin_value = pOriginConstitutiveLaw->GetValue(rThisVar, origin_value);
 
        pDestinationConstitutiveLaw->SetValue(rThisVar, origin_value, rCurrentProcessInfo);
    }
 
    template<class TVarType>
    static inline void GetAndSetDirectVariableOnElements(
        const Variable<TVarType>& rThisVar,
        PointTypePointer pPointOrigin,
        ElementsArrayType::iterator itElemDestination,
        const IndexType GaussPointId,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        std::vector<TVarType> values;
        itElemDestination->CalculateOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
        TVarType aux_value;
        values[GaussPointId] = pPointOrigin->GetValue(rThisVar, aux_value);
        itElemDestination->SetValuesOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
    }
 
    template<class TVarType>
    static inline void GetAndSetWeightedVariableOnConstitutiveLaw(
        const Variable<TVarType>& rThisVar,
        const std::size_t NumberOfPointsFound,
        PointVector& PointsFound,
        const std::vector<double>& PointDistances,
        const double CharacteristicLenght,
        ConstitutiveLaw::Pointer pDestinationConstitutiveLaw,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        TVarType weighting_function_numerator = rThisVar.Zero();
        double weighting_function_denominator = 0.0;
        TVarType origin_value;
 
        for (std::size_t i_point_found = 0; i_point_found < NumberOfPointsFound; ++i_point_found) {
            PointTypePointer p_gp_origin = PointsFound[i_point_found];
 
            const double distance = PointDistances[i_point_found];
 
            origin_value = (p_gp_origin->GetConstitutiveLaw())->GetValue(rThisVar, origin_value);
 
            const double ponderated_weight = p_gp_origin->GetWeight() * std::exp( -4.0 * distance * distance /std::pow(CharacteristicLenght, 2));
 
            weighting_function_numerator   += ponderated_weight * origin_value;
            weighting_function_denominator += ponderated_weight;
        }
 
        const TVarType destination_value = weighting_function_numerator/weighting_function_denominator;
 
        pDestinationConstitutiveLaw->SetValue(rThisVar, destination_value, rCurrentProcessInfo);
    }
 
    template<class TVarType>
    static inline void GetAndSetWeightedVariableOnElements(
        const Variable<TVarType>& rThisVar,
        const std::size_t NumberOfPointsFound,
        PointVector& PointsFound,
        const std::vector<double>& PointDistances,
        const double CharacteristicLenght,
        ElementsArrayType::iterator itElemDestination,
        const IndexType GaussPointId,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        TVarType weighting_function_numerator = rThisVar.Zero();
        double weighting_function_denominator = 0.0;
        TVarType origin_value;
 
        for (std::size_t i_point_found = 0; i_point_found < NumberOfPointsFound; ++i_point_found) {
            PointTypePointer p_gp_origin = PointsFound[i_point_found];
 
            const double distance = PointDistances[i_point_found];
 
            origin_value = p_gp_origin->GetValue(rThisVar, origin_value);
 
            const double ponderated_weight = p_gp_origin->GetWeight() * std::exp( -4.0 * distance * distance /std::pow(CharacteristicLenght, 2));
 
            weighting_function_numerator   += ponderated_weight * origin_value;
            weighting_function_denominator += ponderated_weight;
        }
 
        const TVarType destination_value = weighting_function_numerator/weighting_function_denominator;
 
        std::vector<TVarType> values;
        itElemDestination->CalculateOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
        values[GaussPointId] = destination_value;
        itElemDestination->SetValuesOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
    }
 
    template<class TVarType>
    static inline void InterpolateAddVariableOnConstitutiveLaw(
        GeometryType& rThisGeometry,
        const Variable<TVarType>& rThisVar,
        const Vector& N,
        ConstitutiveLaw::Pointer& pConstitutiveLaw,
        const double Weight
        );
 
    template<class TVarType>
    static inline void InterpolateAddVariableOnElement(
        GeometryType& rThisGeometry,
        const Variable<TVarType>& rThisVar,
        const Vector& N,
        Element& rElement,
        const IndexType GaussPointId,
        const double Weight,
        const ProcessInfo& rCurrentProcessInfo
        );
 
    template<class TVarType>
    static inline void PonderateVariable(
        GeometryType& rThisGeometry,
        const Variable<TVarType>& rThisVar,
        const double TotalWeight
        );
 
    template<class TVarType>
    static inline void InterpolateToNode(
        const Variable<TVarType>& rThisVar,
        const Vector& N,
        NodeType& rNode,
        Element::Pointer pElement
        )
    {
        // An auxiliar value
        TVarType aux_value = rThisVar.Zero();
 
        // Interpolate with shape function
        const std::size_t number_nodes = pElement->GetGeometry().size();
        for (std::size_t i_node = 0; i_node < number_nodes; ++i_node)
            aux_value += N[i_node] * pElement->GetGeometry()[i_node].GetValue(rThisVar);
 
        rNode.SetValue(rThisVar, aux_value);
    }
 
    template<class TVarType>
    static inline void SetInterpolatedValueOnConstitutiveLaw(
        GeometryType& rThisGeometry,
        const Variable<TVarType>& rThisVar,
        const Vector& N,
        ConstitutiveLaw::Pointer pDestinationConstitutiveLaw,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        // An auxiliar value
        TVarType destination_value = rThisVar.Zero();
 
        // Interpolate with shape function
        const std::size_t number_nodes = rThisGeometry.size();
        for (std::size_t i_node = 0; i_node < number_nodes; ++i_node)
            destination_value += N[i_node] * rThisGeometry[i_node].GetValue(rThisVar);
 
        pDestinationConstitutiveLaw->SetValue(rThisVar, destination_value, rCurrentProcessInfo);
    }
 
    template<class TVarType>
    static inline void SetInterpolatedValueOnElement(
        GeometryType& rThisGeometry,
        const Variable<TVarType>& rThisVar,
        const Vector& N,
        Element& rElement,
        const IndexType GaussPointId,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        // An auxiliar value
        TVarType destination_value = rThisVar.Zero();
 
        // Interpolate with shape function
        const std::size_t number_nodes = rThisGeometry.size();
        for (std::size_t i_node = 0; i_node < number_nodes; ++i_node)
            destination_value += N[i_node] * rThisGeometry[i_node].GetValue(rThisVar);
 
        std::vector<TVarType> values;
        rElement.CalculateOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
        values[GaussPointId] = destination_value;
        rElement.SetValuesOnIntegrationPoints(rThisVar, values, rCurrentProcessInfo);
    }
 
    template<std::size_t TDim>
    void InterpolateToNodes()
    {
        // Iterate over nodes
        NodesArrayType& r_nodes_array = mrDestinationMainModelPart.Nodes();
 
        /* Nodes */
        block_for_each(r_nodes_array, auxiliar_search<TDim>(mrOriginMainModelPart),
        [this](Node& rNode, auxiliar_search<TDim>& aux) {
 
            const bool old_entity = rNode.IsDefined(OLD_ENTITY) ? rNode.Is(OLD_ENTITY) : false;
            if (!old_entity) {
                const bool found = aux.point_locator.FindPointOnMeshSimplified(rNode.Coordinates(), aux.N, aux.p_element, mAllocationSize);
 
                if (!found) {
                    KRATOS_WARNING("InternalVariablesInterpolationProcess") << "WARNING: Node "<< rNode.Id() << " not found (interpolation not posible)" <<  "\t X:"<< rNode.X() << "\t Y:"<< rNode.Y() << "\t Z:"<< rNode.Z() << std::endl;
                } else {
                    for (auto& variable_name : mInternalVariableList) {
                        if (KratosComponents<DoubleVarType>::Has(variable_name)) {
                            const auto& r_variable = KratosComponents<DoubleVarType>::Get(variable_name);
                            InterpolateToNode(r_variable, aux.N, rNode, aux.p_element);
                        } else if (KratosComponents<ArrayVarType>::Has(variable_name)) {
                            const auto& r_variable = KratosComponents<ArrayVarType>::Get(variable_name);
                            InterpolateToNode(r_variable, aux.N, rNode, aux.p_element);
                        } else if (KratosComponents<VectorVarType>::Has(variable_name)) {
                            const auto& r_variable = KratosComponents<VectorVarType>::Get(variable_name);
                            InterpolateToNode(r_variable, aux.N, rNode, aux.p_element);
                        } else if (KratosComponents<MatrixVarType>::Has(variable_name)) {
                            const auto& r_variable = KratosComponents<MatrixVarType>::Get(variable_name);
                            InterpolateToNode(r_variable, aux.N, rNode, aux.p_element);
                        } else {
                            KRATOS_WARNING("InternalVariablesInterpolationProcess") << "WARNING:: " << variable_name << " is not registered as any type of compatible variable: DOUBLE or ARRAY_1D or VECTOR or Matrix" << std::endl;
                        }
                    }
                }
            }
        });
    }
 
    InterpolationTypes ConvertInter(const std::string& Str);
 
 
 
 
 
}; // Class InternalVariablesInterpolationProcess
 
 
template<> 
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnConstitutiveLaw<double>(
    GeometryType& rThisGeometry,
    const Variable<double>& rThisVar,
    const Vector& N,
    ConstitutiveLaw::Pointer& pConstitutiveLaw,
    const double Weight
    );
template<> 
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnConstitutiveLaw<array_1d<double, 3>>(
    GeometryType& rThisGeometry,
    const Variable<array_1d<double, 3>>& rThisVar,
    const Vector& N,
    ConstitutiveLaw::Pointer& pConstitutiveLaw,
    const double Weight
    );
template<> 
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnConstitutiveLaw<Vector>(
    GeometryType& rThisGeometry,
    const Variable<Vector>& rThisVar,
    const Vector& N,
    ConstitutiveLaw::Pointer& pConstitutiveLaw,
    const double Weight
    );
template<> 
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnConstitutiveLaw<Matrix>(
    GeometryType& rThisGeometry,
    const Variable<Matrix>& rThisVar,
    const Vector& N,
    ConstitutiveLaw::Pointer& pConstitutiveLaw,
    const double Weight
    );
 
template<>
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnElement(
    GeometryType& rThisGeometry,
    const Variable<double>& rThisVar,
    const Vector& N,
    Element& rElement,
    const IndexType GaussPointId,
    const double Weight,
    const ProcessInfo& rCurrentProcessInfo
    );
 
template<>
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnElement(
    GeometryType& rThisGeometry,
    const Variable<array_1d<double, 3>>& rThisVar,
    const Vector& N,
    Element& rElement,
    const IndexType GaussPointId,
    const double Weight,
    const ProcessInfo& rCurrentProcessInfo
    );
 
template<>
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnElement(
    GeometryType& rThisGeometry,
    const Variable<Vector>& rThisVar,
    const Vector& N,
    Element& rElement,
    const IndexType GaussPointId,
    const double Weight,
    const ProcessInfo& rCurrentProcessInfo
    );
 
template<>
void InternalVariablesInterpolationProcess::InterpolateAddVariableOnElement(
    GeometryType& rThisGeometry,
    const Variable<Matrix>& rThisVar,
    const Vector& N,
    Element& rElement,
    const IndexType GaussPointId,
    const double Weight,
    const ProcessInfo& rCurrentProcessInfo
    );
 
template<>
void InternalVariablesInterpolationProcess::PonderateVariable(
    GeometryType& rThisGeometry,
    const Variable<double>& rThisVar,
    const double TotalWeight
    );
 
template<>
void InternalVariablesInterpolationProcess::PonderateVariable(
    GeometryType& rThisGeometry,
    const Variable<array_1d<double, 3>>& rThisVar,
    const double TotalWeight
    );
 
template<>
void InternalVariablesInterpolationProcess::PonderateVariable(
    GeometryType& rThisGeometry,
    const Variable<Vector>& rThisVar,
    const double TotalWeight
    );
 
template<>
void InternalVariablesInterpolationProcess::PonderateVariable(
    GeometryType& rThisGeometry,
    const Variable<Matrix>& rThisVar,
    const double TotalWeight
    );
 
 
/****************************** INPUT STREAM FUNCTION ******************************/
/***********************************************************************************/
 
template<class TPointType, class TPointerType>
inline std::istream& operator >> (std::istream& rIStream,
                                  InternalVariablesInterpolationProcess& rThis);
 
/***************************** OUTPUT STREAM FUNCTION ******************************/
/***********************************************************************************/
 
template<class TPointType, class TPointerType>
inline std::ostream& operator << (std::ostream& rOStream,
                                  const InternalVariablesInterpolationProcess& rThis)
{
    return rOStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_INTERNAL_VARIABLES_INTERPOLATION_PROCESS  defined