integration_utilities.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//                   Vicente Mataix Ferrandiz
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "geometries/geometry.h"
#include "geometrical_sensitivity_utility.h"
namespace Kratos
{
 
class IntegrationUtilities
{
public:
    template<class TPointType>
    static GeometryData::IntegrationMethod GetIntegrationMethodForExactMassMatrixEvaluation(const Geometry<TPointType>& rGeometry)
    {
        GeometryData::IntegrationMethod integration_method = rGeometry.GetDefaultIntegrationMethod();
        if (integration_method == GeometryData::IntegrationMethod::GI_GAUSS_1)
            integration_method = GeometryData::IntegrationMethod::GI_GAUSS_2;
        else if(integration_method == GeometryData::IntegrationMethod::GI_GAUSS_2)
            integration_method = GeometryData::IntegrationMethod::GI_GAUSS_3;
        else if(integration_method == GeometryData::IntegrationMethod::GI_GAUSS_3)
            integration_method = GeometryData::IntegrationMethod::GI_GAUSS_4;
        else if(integration_method == GeometryData::IntegrationMethod::GI_GAUSS_4)
            integration_method = GeometryData::IntegrationMethod::GI_GAUSS_5;
        return integration_method;
    }
 
    template<class TGeometryType>
    static inline double ComputeDomainSize(const TGeometryType& rGeometry)
    {
        const auto& r_integration_points = rGeometry.IntegrationPoints();
        const auto number_gp = r_integration_points.size();
        Vector temp(number_gp);
        temp = rGeometry.DeterminantOfJacobian(temp);
        double domain_size = 0.0;
        for (unsigned int i = 0; i < number_gp; ++i) {
            domain_size += temp[i] * r_integration_points[i].Weight();
        }
        return domain_size;
    }
 
    template<class TPointType>
    static inline double ComputeDomainSize(
        const Geometry<TPointType>& rGeometry,
        const typename Geometry<TPointType>::IntegrationMethod IntegrationMethod
        )
    {
        const auto& r_integration_points = rGeometry.IntegrationPoints( IntegrationMethod );
        const auto number_gp = r_integration_points.size();
        Vector temp(number_gp);
        temp = rGeometry.DeterminantOfJacobian(temp, IntegrationMethod);
        double domain_size = 0.0;
        for (unsigned int i = 0; i < number_gp; ++i) {
            domain_size += temp[i] * r_integration_points[i].Weight();
        }
        return domain_size;
    }
 
    template<class TPointType>
    static inline double ComputeArea2DGeometry(const Geometry<TPointType>& rGeometry)
    {
        const auto integration_method = rGeometry.GetDefaultIntegrationMethod();
        const auto& r_integration_points = rGeometry.IntegrationPoints( integration_method );
        double volume = 0.0;
        Matrix J(2, 2);
        for ( unsigned int i = 0; i < r_integration_points.size(); i++ ) {
            rGeometry.Jacobian( J, i, integration_method);
            volume += MathUtils<double>::Det2(J) * r_integration_points[i].Weight();
        }
 
        return volume;
    }
 
    template<class TPointType>
    static inline double ComputeArea2DGeometryDerivative(
        const unsigned int DerivativeNodeIndex,
        const unsigned int DerivativeDirectionIndex,
        const Geometry<TPointType>& rGeometry)
    {
        const auto integration_method = rGeometry.GetDefaultIntegrationMethod();
        const auto& r_integration_points = rGeometry.IntegrationPoints( integration_method );
 
        Matrix jacobian(2, 2);
        Matrix shape_function_local_gradients(rGeometry.size(), 2), rdNdXDerivative;
        double detJ_derivative;
        ShapeParameter shape_param;
 
        shape_param.NodeIndex = DerivativeNodeIndex;
        shape_param.Direction = DerivativeDirectionIndex;
 
        double area_derivative = 0.0;
        for (unsigned int g = 0; g < r_integration_points.size(); g++) {
            noalias(shape_function_local_gradients) = rGeometry.ShapeFunctionLocalGradient(g, integration_method);
            rGeometry.Jacobian(jacobian, g, integration_method);
 
            GeometricalSensitivityUtility geometrical_sensitivity_utility(jacobian, shape_function_local_gradients);
            geometrical_sensitivity_utility.CalculateSensitivity(shape_param, detJ_derivative, rdNdXDerivative);
 
            area_derivative += detJ_derivative * r_integration_points[g].Weight();
        }
 
        return area_derivative;
    }
 
    template<class TPointType>
    static inline double ComputeVolume3DGeometry(const Geometry<TPointType>& rGeometry)
    {
        const auto integration_method = rGeometry.GetDefaultIntegrationMethod();
        const auto& r_integration_points = rGeometry.IntegrationPoints( integration_method );
        double volume = 0.0;
        Matrix J(3, 3);
        for ( unsigned int i = 0; i < r_integration_points.size(); i++ ) {
            rGeometry.Jacobian( J, i, integration_method);
            volume += MathUtils<double>::Det3(J) * r_integration_points[i].Weight();
        }
 
        return volume;
    }
 
    template<class TPointType>
    static inline double ComputeVolume3DGeometryDerivative(
        const unsigned int DerivativeNodeIndex,
        const unsigned int DerivativeDirectionIndex,
        const Geometry<TPointType>& rGeometry)
    {
        const auto integration_method = rGeometry.GetDefaultIntegrationMethod();
        const auto& r_integration_points = rGeometry.IntegrationPoints( integration_method );
 
        Matrix jacobian(3, 3);
        Matrix shape_function_local_gradients(rGeometry.size(), 3), rdNdXDerivative;
        double detJ_derivative;
        ShapeParameter shape_param;
 
        shape_param.NodeIndex = DerivativeNodeIndex;
        shape_param.Direction = DerivativeDirectionIndex;
 
        double volume_derivative = 0.0;
        for (unsigned int g = 0; g < r_integration_points.size(); g++) {
            noalias(shape_function_local_gradients) = rGeometry.ShapeFunctionLocalGradient(g, integration_method);
            rGeometry.Jacobian(jacobian, g, integration_method);
 
            GeometricalSensitivityUtility geometrical_sensitivity_utility(jacobian, shape_function_local_gradients);
            geometrical_sensitivity_utility.CalculateSensitivity(shape_param, detJ_derivative, rdNdXDerivative);
 
            volume_derivative += detJ_derivative * r_integration_points[g].Weight();
        }
 
        return volume_derivative;
    }
 
};
 
}  // namespace Kratos.