shallow_water_utilities.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Miguel Maso Sotomayor
//
 
#ifndef KRATOS_SHALLOW_WATER_UTILITIES_H_INCLUDED
#define KRATOS_SHALLOW_WATER_UTILITIES_H_INCLUDED
 
 
// System includes
 
 
// External includes
 
 
// Project includes
#include "includes/model_part.h"
#include "utilities/parallel_utilities.h"
#include "utilities/reduction_utilities.h"
#include "shallow_water_application_variables.h"
 
 
namespace Kratos
{
 
 
 
 
 
 
class KRATOS_API(SHALLOW_WATER_APPLICATION) ShallowWaterUtilities
{
public:
 
    typedef Node NodeType;
 
    typedef Geometry<NodeType> GeometryType;
 
    typedef ModelPart::NodesContainerType NodesContainerType;
 
 
    KRATOS_CLASS_POINTER_DEFINITION(ShallowWaterUtilities);
 
 
 
    void ComputeFreeSurfaceElevation(ModelPart& rModelPart);
 
    void ComputeHeightFromFreeSurface(ModelPart& rModelPart);
 
    void ComputeVelocity(ModelPart& rModelPart, bool PerformProjection = false);
 
    void ComputeSmoothVelocity(ModelPart& rModelPart);
 
    void ComputeMomentum(ModelPart& rModelPart);
 
    void ComputeLinearizedMomentum(ModelPart& rModelPart);
 
    template<bool THistorical>
    void ComputeFroude(ModelPart& rModelPart, const double Epsilon);
 
    template<bool THistorical>
    void ComputeEnergy(ModelPart& rModelPart);
 
    void FlipScalarVariable(Variable<double>& rOriginVariable, Variable<double>& rDestinationVariable, ModelPart& rModelPart);
 
    void IdentifySolidBoundary(ModelPart& rModelPart, double SeaWaterLevel, Flags SolidBoundaryFlag);
 
    void FlagWetElements(ModelPart& rModelPart, Flags WetFlag, double RelativeDryHeight = -1.0);
 
    void ExtrapolateElementalFlagToNodes(ModelPart& rModelPart, Flags Flag);
 
    void NormalizeVector(ModelPart& rModelPart, const Variable<array_1d<double,3>>& rVariable);
 
    template<class TDataType, class TVarType = Variable<TDataType>>
    void SmoothHistoricalVariable(
        const TVarType& rVariable,
        NodesContainerType& rNodes,
        const double ElapsedTime,
        const double SemiPeriod)
    {
        const double smooth = -std::expm1(-ElapsedTime / SemiPeriod);
        block_for_each(rNodes, [&](NodeType& rNode){
            TDataType& initial = rNode.FastGetSolutionStepValue(rVariable, 1);
            TDataType& current = rNode.FastGetSolutionStepValue(rVariable);
            TDataType increment = current - initial;
            current = initial + smooth * increment;
        });
    }
 
    template<class TVarType>
    void CopyVariableToPreviousTimeStep(ModelPart& rModelPart, const TVarType& rVariable)
    {
        block_for_each(rModelPart.Nodes(), [&](NodeType& rNode){
            rNode.FastGetSolutionStepValue(rVariable, 1) = rNode.FastGetSolutionStepValue(rVariable);
        });
    }
 
    void SetMinimumValue(ModelPart& rModelPart, const Variable<double>& rVariable, double MinValue);
 
    void SetMeshZCoordinateToZero(ModelPart& rModelPart);
 
    void SetMeshZ0CoordinateToZero(ModelPart& rModelPart);
 
    void SetMeshZCoordinate(ModelPart& rModelPart, const Variable<double>& rVariable);
 
    void OffsetMeshZCoordinate(ModelPart& rModelPart, const double Increment);
 
    void SwapYZCoordinates(ModelPart& rModelPart);
 
    void SwapY0Z0Coordinates(ModelPart& rModelPart);
 
    void StoreNonHistoricalGiDNoDataIfDry(ModelPart& rModelPart, const Variable<double>& rVariable);
 
    void SwapYZComponents(const Variable<array_1d<double,3>>& rVariable, NodesContainerType& rNodes)
    {
        block_for_each(rNodes, [&](NodeType& rNode){
            array_1d<double,3>& r_value = rNode.FastGetSolutionStepValue(rVariable);
            std::swap(r_value[1], r_value[2]);
        });
    }
 
    template<class TContainerType>
    void SwapYZComponentsNonHistorical(const Variable<array_1d<double,3>>& rVariable, TContainerType& rContainer)
    {
        block_for_each(rContainer, [&](typename TContainerType::value_type& rEntity){
            array_1d<double,3>& r_value = rEntity.GetValue(rVariable);
            std::swap(r_value[1], r_value[2]);
        });
    }
 
    template<class TContainerType>
    void OffsetIds(TContainerType& rContainer, const double Offset)
    {
        block_for_each(rContainer, [&](typename TContainerType::value_type& rEntity){
            rEntity.SetId(rEntity.Id() + Offset);
        });
    }
 
    template<class TContainerType>
    void OffsetIds(TContainerType& rContainer)
    {
        const std::size_t offset = rContainer.size();
        OffsetIds(rContainer, offset);
    }
 
    template<bool THistorical>
    double ComputeL2Norm(ModelPart& rModelPart, const Variable<double>& rVariable);
 
    template<bool THistorical>
    double ComputeL2NormAABB(
        ModelPart& rModelPart,
        const Variable<double>& rVariable,
        Point& rLow,
        Point& rHigh);
 
    template<class TContainerType>
    array_1d<double,3> ComputeHydrostaticForces(
        TContainerType& rContainer,
        const ProcessInfo& rProcessInfo,
        const double RelativeDryHeight = -1.0)
    {
        KRATOS_ERROR_IF_NOT(rProcessInfo.Has(GRAVITY)) << "ShallowWaterUtilities::ComputeHydrostaticForces : GRAVITY is not defined in the ProcessInfo" << std::endl;
        if (rContainer.size() > 0) {
            const auto& r_prop = rContainer.begin()->GetProperties();
            KRATOS_ERROR_IF_NOT(r_prop.Has(DENSITY)) << "ShallowWaterUtilities::ComputeHydrostaticForces : DENSITY is not defined in the Properties" << std::endl;
        }
 
        array_1d<double,3> forces = ZeroVector(3);
        forces = block_for_each<SumReduction<array_1d<double,3>>>(
            rContainer, [&](typename TContainerType::value_type& rEntity){
                array_1d<double,3> local_force = ZeroVector(3);
                if (RelativeDryHeight >= 0.0) {
                    if (IsWet(rEntity.GetGeometry(), RelativeDryHeight)) {
                        rEntity.Calculate(FORCE, local_force, rProcessInfo);
                    }
                } else {
                    rEntity.Calculate(FORCE, local_force, rProcessInfo);
                }
                return local_force;
            }
        );
        return forces;
    }
 
 
private:
 
 
    void CalculateMassMatrix(Matrix& rMassMatrix, const GeometryType& rGeometry);
 
    template<bool THistorical>
    double& GetValue(NodeType& rNode, const Variable<double>& rVariable);
 
    bool IsWet(const GeometryType& rGeometry, const double RelativeDryHeight);
 
    bool IsWet(const GeometryType& rGeometry, const double Height, const double RelativeDryHeight);
 
    bool IsWet(const double Height, const double DryHeight);
 
 
}; // Class ShallowWaterUtilities
 
 
 
 
 
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_SHALLOW_WATER_UTILITIES_H_INCLUDED  defined