flow_rate_slip_utility.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_FLOW_RATE_SLIP_UTILITY_H_INCLUDED
#define KRATOS_FLOW_RATE_SLIP_UTILITY_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "utilities/parallel_utilities.h"
#include "utilities/coordinate_transformation_utilities.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
class FlowRateSlipUtility
    : public CoordinateTransformationUtils<TLocalMatrixType,TLocalVectorType,TValueType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(FlowRateSlipUtility);
 
    typedef CoordinateTransformationUtils<TLocalMatrixType,TLocalVectorType,TValueType> BaseType;
 
    typedef std::size_t SizeType;
 
    typedef Node NodeType;
 
    typedef Geometry<NodeType> GeometryType;
 
 
    FlowRateSlipUtility() : BaseType(2,3,SLIP) {}
 
    virtual ~FlowRateSlipUtility() {}
 
 
 
    virtual void ApplySlipCondition(
        TLocalMatrixType& rLocalMatrix,
        TLocalVectorType& rLocalVector,
        GeometryType& rGeometry) const override
    {
        const SizeType LocalSize = rLocalVector.size(); // We expect this to work both with elements and conditions
 
        if (LocalSize > 0)
        {
            for (SizeType it_node = 0; it_node < rGeometry.PointsNumber(); ++it_node)
            {
                if (this->IsSlip(rGeometry[it_node]))
                {
                    // We fix the first dof (normal velocity) for each rotated block
                    SizeType j = it_node * BaseType::GetBlockSize();
 
                    array_1d<double,3> vel = rGeometry[it_node].FastGetSolutionStepValue(MOMENTUM);
                    array_1d<double,3> n = rGeometry[it_node].FastGetSolutionStepValue(NORMAL);
                    this->Normalize(n);
 
                    for (SizeType i = 0; i < j; ++i) // Skip term (i,i)
                    {
                        rLocalMatrix(i,j) = 0.0;
                        rLocalMatrix(j,i) = 0.0;
                    }
                    for (SizeType i = j+1; i < LocalSize; ++i)
                    {
                        rLocalMatrix(i,j) = 0.0;
                        rLocalMatrix(j,i) = 0.0;
                    }
 
                    rLocalVector(j) = - inner_prod(n, vel);
                    rLocalMatrix(j,j) = 1.0;
                }
            }
        }
    }
 
    virtual void ApplySlipCondition(
        TLocalVectorType& rLocalVector,
        GeometryType& rGeometry) const override
    {
        if (rLocalVector.size() > 0)
        {
            for (SizeType it_node = 0; it_node < rGeometry.PointsNumber(); ++it_node)
            {
                if (this->IsSlip(rGeometry[it_node]))
                {
                    // We fix the first dof (normal velocity) for each rotated block
                    SizeType j = it_node * BaseType::GetBlockSize();
 
                    array_1d<double,3> vel = rGeometry[it_node].FastGetSolutionStepValue(MOMENTUM);
                    array_1d<double,3> n = rGeometry[it_node].FastGetSolutionStepValue(NORMAL);
                    this->Normalize(n);
 
                    rLocalVector[j] = inner_prod(n, vel);
                }
            }
        }
    }
 
    virtual void RotateVelocities(ModelPart& rModelPart) const override
    {
        struct TLS {
            TLocalVectorType vel;
            TLocalVectorType tmp;
        };
        TLS tls;
        tls.vel.resize(BaseType::GetDomainSize());
        tls.tmp.resize(BaseType::GetDomainSize());
 
        block_for_each(rModelPart.Nodes(), tls, [&](NodeType& rNode, TLS& rTLS){
            if (this->IsSlip(rNode))
            {
                // For shallow water problems, domain size is always 2
                BoundedMatrix<double,2,2> rot;
                BaseType::LocalRotationOperatorPure(rot, rNode);
 
                array_1d<double,3>& r_velocity = rNode.FastGetSolutionStepValue(MOMENTUM);
                for(SizeType i = 0; i < 2; i++) {
                    rTLS.vel[i] = r_velocity[i];
                }
                noalias(rTLS.tmp) = prod(rot, rTLS.vel);
                for(SizeType i = 0; i < 2; i++) {
                    r_velocity[i] = rTLS.tmp[i];
                }
            }
        });
    }
 
    virtual void RecoverVelocities(ModelPart& rModelPart) const override
    {
        struct TLS {
            TLocalVectorType vel;
            TLocalVectorType tmp;
        };
        TLS tls;
        tls.vel.resize(BaseType::GetDomainSize());
        tls.tmp.resize(BaseType::GetDomainSize());
 
        block_for_each(rModelPart.Nodes(), tls, [&](NodeType& rNode, TLS& rTLS){
            if( this->IsSlip(rNode) )
            {
                // For shallow water problems, domain size is always 2
                BoundedMatrix<double,2,2> rot;
                BaseType::LocalRotationOperatorPure(rot, rNode);
 
                array_1d<double,3>& r_velocity = rNode.FastGetSolutionStepValue(MOMENTUM);
                for(SizeType i = 0; i < 2; i++) {
                    rTLS.vel[i] = r_velocity[i];
                }
                noalias(rTLS.tmp) = prod(trans(rot), rTLS.vel);
                for(SizeType i = 0; i < 2; i++) {
                    r_velocity[i] = rTLS.tmp[i];
                }
            }
        });
    }
 
 
 
 
    virtual std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "FlowRateSlipUtility";
        return buffer.str();
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "FlowRateSlipUtility";
    }
 
 
 
 
private:
 
    // FlowRateSlipUtility& operator=(FlowRateSlipUtility const& rOther) {}
 
    // FlowRateSlipUtility(FlowRateSlipUtility const& rOther) {}
 
 
}; // Class FlowRateSlipUtility
 
 
 
 
template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
inline std::istream& operator >> (
    std::istream& rIStream,
    FlowRateSlipUtility<TLocalMatrixType, TLocalVectorType,TValueType>& rThis)
{
    return rIStream;
}
 
template<class TLocalMatrixType, class TLocalVectorType, class TValueType>
inline std::ostream& operator << (
    std::ostream& rOStream,
    const FlowRateSlipUtility<TLocalMatrixType, TLocalVectorType,TValueType>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_FLOW_RATE_SLIP_UTILITY_H_INCLUDED  defined