fluid_calculation_utilities.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Suneth Warnakulasuriya
//
 
#if !defined(KRATOS_FLUID_CALCULATION_UTILITIES_H)
#define KRATOS_FLUID_CALCULATION_UTILITIES_H
 
// system includes
#include <tuple>
#include <type_traits>
 
// External includes
 
// Project includes
#include "geometries/geometry.h"
#include "includes/node.h"
#include "includes/ublas_interface.h"
 
namespace Kratos
{
 
 
class KRATOS_API(FLUID_DYNAMICS_APPLICATION) FluidCalculationUtilities
{
public:
 
    using IndexType = std::size_t;
 
    using NodeType = Node;
 
    using GeometryType = Geometry<NodeType>;
 
    template<class TDataType>
    using ComponentDataType = std::tuple<TDataType&, const TDataType&>;
 
 
    template <class... TRefVariableValuePairArgs>
    static void EvaluateInPoint(
        const GeometryType& rGeometry,
        const Vector& rShapeFunction,
        const int Step,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        KRATOS_TRY
 
        const auto& r_node = rGeometry[0];
        const double shape_function_value = rShapeFunction[0];
 
        (std::apply([&](auto&&... args) {((std::get<0>(args) = std::get<1>(args) * shape_function_value), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.FastGetSolutionStepValue(std::get<1>(rValueVariablePairs), Step))), ...);
 
        for (IndexType c = 1; c < rGeometry.PointsNumber(); ++c) {
            const auto& r_node = rGeometry[c];
            const double shape_function_value = rShapeFunction[c];
 
            (std::apply([&](auto&&... args) {((std::get<0>(args) += std::get<1>(args) * shape_function_value), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.FastGetSolutionStepValue(std::get<1>(rValueVariablePairs), Step))), ...);
        }
 
        KRATOS_CATCH("");
    }
 
    template <class... TRefVariableValuePairArgs>
    static void inline EvaluateInPoint(
        const GeometryType& rGeometry,
        const Vector& rShapeFunction,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        EvaluateInPoint<TRefVariableValuePairArgs...>(
            rGeometry, rShapeFunction, 0, rValueVariablePairs...);
    }
 
    template <class... TRefVariableValuePairArgs>
    static void EvaluateNonHistoricalInPoint(
        const GeometryType& rGeometry,
        const Vector& rShapeFunction,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        KRATOS_TRY
 
        const auto& r_node = rGeometry[0];
        const double shape_function_value = rShapeFunction[0];
 
        (std::apply([&](auto&&... args) {((std::get<0>(args) = std::get<1>(args) * shape_function_value), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.GetValue(std::get<1>(rValueVariablePairs)))), ...);
 
        for (IndexType c = 1; c < rGeometry.PointsNumber(); ++c) {
            const auto& r_node = rGeometry[c];
            const double shape_function_value = rShapeFunction[c];
 
            (std::apply([&](auto&&... args) {((std::get<0>(args) += std::get<1>(args) * shape_function_value), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.GetValue(std::get<1>(rValueVariablePairs)))), ...);
        }
 
        KRATOS_CATCH("");
    }
 
    template <class... TRefVariableValuePairArgs>
    static void EvaluateGradientInPoint(
        const GeometryType& rGeometry,
        const Matrix& rShapeFunctionDerivatives,
        const int Step,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        KRATOS_TRY
 
        const auto& r_node = rGeometry[0];
        const Vector& shape_function_derivative = row(rShapeFunctionDerivatives, 0);
 
        for (IndexType i = 0; i < rShapeFunctionDerivatives.size2(); ++i) {
            (std::apply([&](auto&&... args) {((std::get<0>(args) = std::get<1>(args) * shape_function_derivative[i]), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.FastGetSolutionStepValue(std::get<1>(rValueVariablePairs), Step), i)), ...);
        }
 
        for (IndexType c = 1; c < rGeometry.PointsNumber(); ++c) {
            const auto& r_node = rGeometry[c];
            const Vector& shape_function_derivative = row(rShapeFunctionDerivatives, c);
 
            for (IndexType i = 0; i < rShapeFunctionDerivatives.size2(); ++i) {
                (std::apply([&](auto&&... args) {((std::get<0>(args) += std::get<1>(args) * shape_function_derivative[i]), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.FastGetSolutionStepValue(std::get<1>(rValueVariablePairs), Step), i)), ...);
            }
        }
 
        KRATOS_CATCH("");
    }
 
    template <class... TRefVariableValuePairArgs>
    static void inline EvaluateGradientInPoint(
        const GeometryType& rGeometry,
        const Matrix& rShapeFunctionDerivatives,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        EvaluateGradientInPoint<TRefVariableValuePairArgs...>(
            rGeometry, rShapeFunctionDerivatives, 0, rValueVariablePairs...);
    }
 
    template <class... TRefVariableValuePairArgs>
    static void EvaluateNonHistoricalGradientInPoint(
        const GeometryType& rGeometry,
        const Matrix& rShapeFunctionDerivatives,
        const TRefVariableValuePairArgs&... rValueVariablePairs)
    {
        KRATOS_TRY
 
        const auto& r_node = rGeometry[0];
        const Vector& shape_function_derivative = row(rShapeFunctionDerivatives, 0);
 
        for (IndexType i = 0; i < rShapeFunctionDerivatives.size2(); ++i) {
            (std::apply([&](auto&&... args) {((std::get<0>(args) = std::get<1>(args) * shape_function_derivative[i]), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.GetValue(std::get<1>(rValueVariablePairs)), i)), ...);
        }
 
        for (IndexType c = 1; c < rGeometry.PointsNumber(); ++c) {
            const auto& r_node = rGeometry[c];
            const Vector& shape_function_derivative = row(rShapeFunctionDerivatives, c);
 
            for (IndexType i = 0; i < rShapeFunctionDerivatives.size2(); ++i) {
                (std::apply([&](auto&&... args) {((std::get<0>(args) += std::get<1>(args) * shape_function_derivative[i]), ...);}, GetComponentsArray(std::get<0>(rValueVariablePairs), r_node.GetValue(std::get<1>(rValueVariablePairs)), i)), ...);
            }
        }
 
        KRATOS_CATCH("");
    }
 
    template <unsigned int TSize>
    static void ReadSubVector(
        BoundedVector<double, TSize>& rOutput,
        const Vector& rInput,
        const IndexType Position)
    {
        KRATOS_TRY
 
        KRATOS_DEBUG_ERROR_IF(rInput.size() < Position + TSize)
            << "rInput size is not enough for sub vector retireval. [ "
               "rInput.size() = "
            << rInput.size() << ", Requested Position = " << Position
            << ", Requested vector size = " << TSize << " ].\n";
 
        for (IndexType i = 0; i < TSize; ++i) {
            rOutput[i] = rInput[Position + i];
        }
 
        KRATOS_CATCH("");
    }
 
    template<unsigned int TSize>
    static void AddSubVector(
        Matrix& rOutput,
        const BoundedVector<double, TSize>& rInput,
        const unsigned int RowIndex,
        const unsigned int ColumnIndex)
    {
        KRATOS_TRY
 
        KRATOS_DEBUG_ERROR_IF(RowIndex >= rOutput.size1())
            << "RowIndex is larger than or equal to rOutput number of rows. [ "
               "rOutput.size1() = "
            << rOutput.size1() << ", RowIndex = " << RowIndex << " ].\n";
 
        KRATOS_DEBUG_ERROR_IF(ColumnIndex + TSize > rOutput.size2())
            << "rOutput matrix does not have sufficient number of columns [ "
               "rOutput.size2() = "
            << rOutput.size2() << ", ColumnIndex = " << ColumnIndex
            << ", TSize = " << TSize << " ].\n";
 
        for (unsigned int i = 0; i < TSize; ++i) {
            rOutput(RowIndex, ColumnIndex + i) += rInput[i];
        }
 
        KRATOS_CATCH("");
    }
 
    static double CalculateLogarithmicYPlusLimit(
        const double Kappa,
        const double Beta,
        const int MaxIterations = 20,
        const double Tolerance = 1e-6);
 
    static double CalculateLogarithmicYPlus(
        const double WallVelocityMagnitude,
        const double WallHeight,
        const double KinematicViscosity,
        const double Kappa,
        const double Beta,
        const double YPlusLimit,
        const int MaxIterations = 20,
        const double Tolerance = 1e-6);
 
 
private:
 
    // gauss point evaluation templates
    static auto inline GetComponentsArray(double& rOutput, const double& rInput) { return std::array<ComponentDataType<double>, 1>{std::tie(rOutput, rInput)};}
    static auto inline GetComponentsArray(array_1d<double, 2>& rOutput, const array_1d<double, 3>& rInput) { return std::array<ComponentDataType<double>, 2>{std::tie(rOutput[0], rInput[0]), std::tie(rOutput[1], rInput[1])};}
    static auto inline GetComponentsArray(array_1d<double, 3>& rOutput, const array_1d<double, 3>& rInput) { return std::array<ComponentDataType<double>, 3>{std::tie(rOutput[0], rInput[0]), std::tie(rOutput[1], rInput[1]), std::tie(rOutput[2], rInput[2])};}
    static auto inline GetComponentsArray(Vector& rOutput, const Vector& rInput) { return std::array<ComponentDataType<Vector>, 1>{std::tie(rOutput, rInput)};}
    static auto inline GetComponentsArray(Matrix& rOutput, const Matrix& rInput) { return std::array<ComponentDataType<Matrix>, 1>{std::tie(rOutput, rInput)};}
 
    // gauss point gradient evaluation templates
    static auto inline GetComponentsArray(array_1d<double, 2>& rOutput, const double& rInput, const IndexType DerivativeIndex) { return std::array<ComponentDataType<double>, 1>{std::tie(rOutput[DerivativeIndex], rInput)};}
    static auto inline GetComponentsArray(array_1d<double, 3>& rOutput, const double& rInput, const IndexType DerivativeIndex) { return std::array<ComponentDataType<double>, 1>{std::tie(rOutput[DerivativeIndex], rInput)};}
    static auto inline GetComponentsArray(BoundedMatrix<double, 2, 2>& rOutput, const array_1d<double, 3>& rInput, const IndexType DerivativeIndex) { return std::array<ComponentDataType<double>, 2>{std::tie(rOutput(0, DerivativeIndex), rInput[0]), std::tie(rOutput(1, DerivativeIndex), rInput[1])};}
    static auto inline GetComponentsArray(BoundedMatrix<double, 3, 3>& rOutput, const array_1d<double, 3>& rInput, const IndexType DerivativeIndex) { return std::array<ComponentDataType<double>, 3>{std::tie(rOutput(0, DerivativeIndex), rInput[0]), std::tie(rOutput(1, DerivativeIndex), rInput[1]), std::tie(rOutput(2, DerivativeIndex), rInput[2])};}
 
};
 
 
} // namespace Kratos
 
#endif // KRATOS_FLUID_CALCULATION_UTILITIES_H