fractional_step.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ \.
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Jordi Cotela
//
 
#if !defined(KRATOS_FRACTIONAL_STEP_H_INCLUDED )
#define  KRATOS_FRACTIONAL_STEP_H_INCLUDED
 
// System includes
#include <string>
#include <iostream>
 
 
// External includes
 
 
// Project includes
#include "containers/array_1d.h"
#include "includes/define.h"
#include "includes/element.h"
#include "includes/serializer.h"
#include "geometries/geometry.h"
#include "utilities/math_utils.h"
 
// Application includes
#include "fluid_dynamics_application_variables.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
    template< unsigned int TDim >
    class FractionalStep : public Element
    {
    public:
 
        KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(FractionalStep);
 
        typedef Node NodeType;
 
        typedef Geometry<NodeType> GeometryType;
 
        typedef Geometry<NodeType>::PointsArrayType NodesArrayType;
 
        typedef Vector VectorType;
 
        typedef Matrix MatrixType;
 
        typedef std::size_t IndexType;
 
        typedef std::size_t SizeType;
 
        typedef std::vector<std::size_t> EquationIdVectorType;
 
        typedef std::vector< Dof<double>::Pointer > DofsVectorType;
 
        typedef PointerVectorSet<Dof<double>, IndexedObject> DofsArrayType;
 
        typedef Kratos::Vector ShapeFunctionsType;
 
        typedef Kratos::Matrix ShapeFunctionDerivativesType;
 
        typedef GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType;
 
 
        //Constructors.
 
 
        FractionalStep(IndexType NewId = 0) :
            Element(NewId)
        {}
 
 
        FractionalStep(IndexType NewId, const NodesArrayType& ThisNodes) :
            Element(NewId, ThisNodes)
        {}
 
 
        FractionalStep(IndexType NewId, GeometryType::Pointer pGeometry) :
            Element(NewId, pGeometry)
        {}
 
 
        FractionalStep(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) :
            Element(NewId, pGeometry, pProperties)
        {}
 
        ~FractionalStep() override
        {}
 
 
 
 
 
 
        Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,
                                Element::PropertiesType::Pointer pProperties) const override
        {
            return Kratos::make_intrusive< FractionalStep<TDim> >(NewId, this->GetGeometry().Create(ThisNodes), pProperties);
        }
 
        Element::Pointer Create(IndexType NewId, Element::GeometryType::Pointer pGeom, Element::PropertiesType::Pointer pProperties) const override
        {
            return Kratos::make_intrusive< FractionalStep<TDim> >(NewId, pGeom, pProperties);
        }
 
        Element::Pointer Clone(IndexType NewId, NodesArrayType const& rThisNodes) const override
        {
            Element::Pointer pNewElement = Create(NewId, GetGeometry().Create( rThisNodes ), pGetProperties() );
 
            pNewElement->SetData(this->GetData());
            pNewElement->SetFlags(this->GetFlags());
 
            return pNewElement;
        }
 
        void Initialize(const ProcessInfo &rCurrentProcessInfo) override;
 
        void InitializeSolutionStep(const ProcessInfo &rCurrentProcessInfo) override;
 
 
        void InitializeNonLinearIteration(const ProcessInfo &rCurrentProcessInfo) override;
 
        void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
                                          VectorType& rRightHandSideVector,
                                          const ProcessInfo& rCurrentProcessInfo) override;
 
 
        void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix,
                                           const ProcessInfo& rCurrentProcessInfo) override
        {
            KRATOS_TRY;
            KRATOS_THROW_ERROR(std::logic_error,"FractionalStep::CalculateLeftHandSide not implemented","");
            KRATOS_CATCH("");
        }
 
        void CalculateRightHandSide(VectorType& rRightHandSideVector,
                                    const ProcessInfo& rCurrentProcessInfo) override
        {
            KRATOS_TRY;
            KRATOS_THROW_ERROR(std::logic_error,"FractionalStep::CalculateRightHandSide not implemented","");
            KRATOS_CATCH("");
        }
 
        void Calculate(const Variable<double>& rVariable,
                               double& rOutput,
                               const ProcessInfo& rCurrentProcessInfo) override;
 
 
        void Calculate(const Variable<array_1d<double, 3 > >& rVariable,
                               array_1d<double, 3 > & rOutput,
                               const ProcessInfo& rCurrentProcessInfo) override;
 
 
        // The following methods have different implementations depending on TDim
 
        void EquationIdVector(EquationIdVectorType& rResult,
                              const ProcessInfo& rCurrentProcessInfo) const override;
 
 
        void GetDofList(DofsVectorType& rElementalDofList,
                        const ProcessInfo& rCurrentProcessInfo) const override;
 
 
        GeometryData::IntegrationMethod GetIntegrationMethod() const override;
 
 
        void CalculateOnIntegrationPoints(
            const Variable<array_1d<double, 3 > >& rVariable,
            std::vector<array_1d<double, 3 > >& rValues,
            const ProcessInfo& rCurrentProcessInfo) override;
 
 
        void CalculateOnIntegrationPoints(
            const Variable<double>& rVariable,
            std::vector<double>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override;
 
 
        void CalculateOnIntegrationPoints(
            const Variable<array_1d<double, 6 > >& rVariable,
            std::vector<array_1d<double, 6 > >& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
        {
            this->GetElementalValueForOutput< array_1d<double,6> >(rVariable,rValues);
        }
 
 
        void CalculateOnIntegrationPoints(
            const Variable<Vector>& rVariable,
            std::vector<Vector>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
        {
            this->GetElementalValueForOutput< Vector >(rVariable,rValues);
        }
 
 
        void CalculateOnIntegrationPoints(
            const Variable<Matrix>& rVariable,
            std::vector<Matrix>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
        {
            this->GetElementalValueForOutput< Matrix >(rVariable,rValues);
        }
 
 
 
 
        int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 
 
 
        const Parameters GetSpecifications() const override;
 
        std::string Info() const override
        {
            std::stringstream buffer;
            buffer << "FractionalStep #" << Id();
            return buffer.str();
        }
 
        void PrintInfo(std::ostream& rOStream) const override
        {
            rOStream << "FractionalStep" << TDim << "D";
        }
 
//        /// Print object's data.
//        virtual void PrintData(std::ostream& rOStream) const;
 
 
 
 
    protected:
 
 
 
 
 
 
 
        virtual void CalculateLocalFractionalVelocitySystem(
            MatrixType& rLeftHandSideMatrix,
            VectorType& rRightHandSideVector,
            const ProcessInfo& rCurrentProcessInfo);
 
        virtual void CalculateLocalPressureSystem(MatrixType& rLeftHandSideMatrix,
                                          VectorType& rRightHandSideVector,
                                          const ProcessInfo& rCurrentProcessInfo);
 
        void CalculateEndOfStepSystem(MatrixType& rLeftHandSideMatrix,
                                      VectorType& rRightHandSideVector,
                                      const ProcessInfo& rCurrentProcessInfo);
 
        void VelocityEquationIdVector(EquationIdVectorType& rResult,
                                      const ProcessInfo& rCurrentProcessInfo) const ;
 
        void PressureEquationIdVector(EquationIdVectorType& rResult,
                                      const ProcessInfo& rCurrentProcessInfo) const ;
 
        void GetVelocityDofList(DofsVectorType& rElementalDofList,
                                const ProcessInfo& rCurrentProcessInfo) const ;
 
        void GetPressureDofList(DofsVectorType& rElementalDofList,
                                const ProcessInfo& rCurrentProcessInfo) const ;
 
        void GetPressureValues(Vector& rValues,
                               const int Step = 0) const ;
 
        void GetVelocityValues(Vector& rValues,
                               const int Step = 0) const ;
 
        virtual void CalculateGeometryData(ShapeFunctionDerivativesArrayType& rDN_DX,
                                           Matrix& rNContainer,
                                           Vector& rGaussWeights);
 
        double ElementSize(/*ShapeFunctionDerivativesType& rDN_DX*/);
 
 
        virtual double EffectiveViscosity(double Density,
                                          const ShapeFunctionsType &rN,
                                          const ShapeFunctionDerivativesType &rDN_DX,
                                          double ElemSize,
                                          const ProcessInfo &rProcessInfo);
 
        double EquivalentStrainRate(const ShapeFunctionDerivativesType &rDN_DX) const;
 
 
 
        void AddMomentumMassTerm(Matrix& rMassMatrix,
                                 const ShapeFunctionsType& rN,
                                 const double Weight);
 
        virtual void AddMomentumSystemTerms(Matrix& rLHSMatrix,
                                    Vector& rRHSVector,
                                    const double Density,
                                    const Vector& rConvOperator,
                                    const array_1d<double,3>& rBodyForce,
                                    const double OldPressure,
                                    const double TauOne,
                                    const double TauTwo,
                                    const array_1d<double,3>& rMomentumProjection,
                                    const double MassProjection,
                                    const ShapeFunctionsType& rN,
                                    const ShapeFunctionDerivativesType& rDN_DX,
                                    const double Weight);
 
        void AddViscousTerm(MatrixType& rDampingMatrix,
                            const ShapeFunctionDerivativesType& rShapeDeriv,
                            const double Weight);
 
 
//        virtual void CalculateTau(double& TauOne,
//                                  double& TauTwo,
//                                  double ElemSize,
//                                  const ProcessInfo& rCurrentProcessInfo);
 
 
        virtual void CalculateTau(double& TauOne,
                                  double& TauTwo,
                                  double ElemSize,
                                  const array_1d< double, 3 > & rAdvVel,
                                  const double Density,
                                  const double Viscosity,
                                  const ProcessInfo& rCurrentProcessInfo);
 
        virtual void CalculateProjectionRHS(VectorType& rMomentumRHS,
                                            VectorType& rMassRHS,
                                            const ShapeFunctionsType& rN,
                                            const ShapeFunctionDerivativesType& rDN_DX,
                                            const double Weight);
 
        virtual void CalculateProjectionRHS(VectorType& rConvTerm,
                                            VectorType& rPresTerm,
                                            VectorType& rDivTerm,
                                            const ShapeFunctionsType& rN,
                                            const ShapeFunctionDerivativesType& rDN_DX,
                                            const double Weight);
 
 
        void ModulatedGradientDiffusion(MatrixType& rDampingMatrix,
                                        const ShapeFunctionDerivativesType& rDN_DX,
                                        const double Weight);
 
 
 
        void ConvectionOperator(Vector& rResult,
                                const array_1d<double,3>& rConvVel,
                                const ShapeFunctionDerivativesType& DN_DX);
 
 
        virtual void EvaluateConvVelocity(array_1d<double,3>& rConvVel,
                                          const ShapeFunctionsType& N);
 
 
        template< class TVariableType >
        void EvaluateInPoint(TVariableType& rResult,
                             const Kratos::Variable<TVariableType>& Var,
                             const ShapeFunctionsType& rShapeFunc)
        {
            GeometryType& rGeom = this->GetGeometry();
            const SizeType NumNodes = rGeom.PointsNumber();
 
            rResult = rShapeFunc[0] * rGeom[0].FastGetSolutionStepValue(Var);
 
            for(SizeType i = 1; i < NumNodes; i++)
            {
                rResult += rShapeFunc[i] * rGeom[i].FastGetSolutionStepValue(Var);
            }
        }
 
 
        template< class TVariableType >
        void EvaluateInPoint(TVariableType& rResult,
                             const Kratos::Variable<TVariableType>& Var,
                             const ShapeFunctionsType& rShapeFunc,
                             const IndexType Step)
        {
            GeometryType& rGeom = this->GetGeometry();
            const SizeType NumNodes = rGeom.PointsNumber();
 
            rResult = rShapeFunc[0] * rGeom[0].FastGetSolutionStepValue(Var,Step);
 
            for(SizeType i = 1; i < NumNodes; i++)
            {
                rResult += rShapeFunc[i] * rGeom[i].FastGetSolutionStepValue(Var,Step);
            }
        }
 
        void EvaluateGradientInPoint(array_1d<double,TDim>& rResult,
                                     const Kratos::Variable<double>& Var,
                                     const ShapeFunctionDerivativesType& rDN_DX)
        {
            GeometryType& rGeom = this->GetGeometry();
            const SizeType NumNodes = rGeom.PointsNumber();
 
            const double& var = rGeom[0].FastGetSolutionStepValue(Var);
            for (SizeType d = 0; d < TDim; ++d)
                rResult[d] = rDN_DX(0,d) * var;
 
            for(SizeType i = 1; i < NumNodes; i++)
            {
              const double& var = rGeom[i].FastGetSolutionStepValue(Var);
              for (SizeType d = 0; d < TDim; ++d)
                        rResult[d] += rDN_DX(i,d) * var;
 
            }
        }
 
        void EvaluateDivergenceInPoint(double& rResult,
                                       const Kratos::Variable< array_1d<double,3> >& Var,
                                       const ShapeFunctionDerivativesType& rDN_DX)
        {
            GeometryType& rGeom = this->GetGeometry();
            const SizeType NumNodes = rGeom.PointsNumber();
 
            rResult = 0.0;
            for(SizeType i = 0; i < NumNodes; i++)
            {
                const array_1d<double,3>& var = rGeom[i].FastGetSolutionStepValue(Var);
                for (SizeType d = 0; d < TDim; ++d)
                {
                    rResult += rDN_DX(i,d) * var[d];
                }
            }
        }
 
 
        template<class TValueType>
        void GetElementalValueForOutput(const Kratos::Variable<TValueType>& rVariable,
                                        std::vector<TValueType>& rOutput)
        {
            unsigned int NumValues = this->GetGeometry().IntegrationPointsNumber(GeometryData::IntegrationMethod::GI_GAUSS_2);
            rOutput.resize(NumValues);
            /*
             The cast is done to avoid modification of the element's data. Data modification
             would happen if rVariable is not stored now (would initialize a pointer to &rVariable
             with associated value of 0.0). This is catastrophic if the variable referenced
             goes out of scope.
             */
            const FractionalStep<TDim>* const_this = static_cast<const FractionalStep<TDim>*> (this);
            const TValueType& Val = const_this->GetValue(rVariable);
 
            for (unsigned int i = 0; i < NumValues; i++)
                rOutput[i] = Val;
        }
 
 
 
 
 
 
 
 
    private:
 
 
 
        friend class Serializer;
 
        void save(Serializer& rSerializer) const override
        {
            KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Element );
        }
 
        void load(Serializer& rSerializer) override
        {
            KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element);
        }
 
 
 
 
 
 
 
 
 
 
        FractionalStep & operator=(FractionalStep const& rOther);
 
        FractionalStep(FractionalStep const& rOther);
 
 
    }; // Class FractionalStep
 
 
 
 
 
 
    template< unsigned int TDim >
    inline std::istream& operator >>(std::istream& rIStream,
                                     FractionalStep<TDim>& rThis)
    {
        return rIStream;
    }
 
    template< unsigned int TDim >
    inline std::ostream& operator <<(std::ostream& rOStream,
                                     const FractionalStep<TDim>& rThis)
    {
        rThis.PrintInfo(rOStream);
        rOStream << std::endl;
        rThis.PrintData(rOStream);
 
        return rOStream;
    }
 
 
} // namespace Kratos.
 
#endif // KRATOS_FRACTIONAL_STEP_H_INCLUDED  defined