wall_condition_discontinuous.h

Go to the documentation of this file.

/*
==============================================================================
Kratos Fluid Dynamics Application
Kratos
A General Purpose Software for Multi-Physics Finite Element Analysis
Version 1.0 (Released on march 05, 2007).
 
Copyright 2007
Pooyan Dadvand, Riccardo Rossi
pooyan@cimne.upc.edu
rrossi@cimne.upc.edu
CIMNE (International Center for Numerical Methods in Engineering),
Gran Capita' s/n, 08034 Barcelona, Spain
 
Permission is hereby granted, free  of charge, to any person obtaining
a  copy  of this  software  and  associated  documentation files  (the
"Software"), to  deal in  the Software without  restriction, including
without limitation  the rights to  use, copy, modify,  merge, publish,
distribute,  sublicense and/or  sell copies  of the  Software,  and to
permit persons to whom the Software  is furnished to do so, subject to
the following condition:
 
Distribution of this code for  any  commercial purpose  is permissible
ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNER.
 
The  above  copyright  notice  and  this permission  notice  shall  be
included in all copies or substantial portions of the Software.
 
THE  SOFTWARE IS  PROVIDED  "AS  IS", WITHOUT  WARRANTY  OF ANY  KIND,
EXPRESS OR  IMPLIED, INCLUDING  BUT NOT LIMITED  TO THE  WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT  SHALL THE AUTHORS OR COPYRIGHT HOLDERS  BE LIABLE FOR ANY
CLAIM, DAMAGES OR  OTHER LIABILITY, WHETHER IN AN  ACTION OF CONTRACT,
TORT  OR OTHERWISE, ARISING  FROM, OUT  OF OR  IN CONNECTION  WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
==============================================================================
 */
 
#ifndef KRATOS_WALL_CONDITION_DISCONTINUOUS_H
#define KRATOS_WALL_CONDITION_DISCONTINUOUS_H
 
// System includes
#include <string>
#include <iostream>
 
 
// External includes
 
 
// Project includes
#include "includes/define.h"
#include "includes/serializer.h"
#include "wall_condition.h"
#include "includes/process_info.h"
#include "includes/deprecated_variables.h"
 
// Application includes
#include "fluid_dynamics_application_variables.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
template< unsigned int TDim, unsigned int TNumNodes = TDim >
class WallConditionDiscontinuous : public WallCondition<TDim,TNumNodes>
{
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(WallConditionDiscontinuous);
 
    typedef Node NodeType;
 
    typedef Properties PropertiesType;
 
    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;
 
 
 
    WallConditionDiscontinuous(IndexType NewId = 0):
        WallCondition<TDim,TNumNodes>(NewId)
    {
    }
 
 
    WallConditionDiscontinuous(IndexType NewId,
                               const NodesArrayType& ThisNodes):
        WallCondition<TDim,TNumNodes>(NewId,ThisNodes)
    {
    }
 
 
    WallConditionDiscontinuous(IndexType NewId,
                               GeometryType::Pointer pGeometry):
        WallCondition<TDim,TNumNodes>(NewId,pGeometry)
    {
    }
 
 
    WallConditionDiscontinuous(IndexType NewId,
                               GeometryType::Pointer pGeometry,
                               PropertiesType::Pointer pProperties):
        WallCondition<TDim,TNumNodes>(NewId,pGeometry,pProperties)
    {
    }
 
    WallConditionDiscontinuous(WallConditionDiscontinuous const& rOther):
        WallCondition<TDim,TNumNodes>(rOther)
    {
    }
 
    ~WallConditionDiscontinuous() override {}
 
 
 
    WallConditionDiscontinuous & operator=(WallConditionDiscontinuous const& rOther)
    {
        Condition::operator=(rOther);
 
        return *this;
    }
 
 
 
    Condition::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override
    {
        return Kratos::make_intrusive<WallConditionDiscontinuous>(NewId, this->GetGeometry().Create(ThisNodes), pProperties);
    }
 
 
 
    Condition::Pointer Create(IndexType NewId, Condition::GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
    {
        return Kratos::make_intrusive<WallConditionDiscontinuous>(NewId, pGeom, pProperties);
    }
 
 
 
    void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
                                      VectorType& rRightHandSideVector,
                                      const ProcessInfo& rCurrentProcessInfo) override
    {
        const ProcessInfo& r_process_info = rCurrentProcessInfo;
        unsigned int step = r_process_info[FRACTIONAL_STEP];
        if ( step == 1)
        {
            // Initialize local contributions
            const SizeType LocalSize = TDim * TNumNodes;
 
            if (rLeftHandSideMatrix.size1() != LocalSize)
                rLeftHandSideMatrix.resize(LocalSize,LocalSize);
            if (rRightHandSideVector.size() != LocalSize)
                rRightHandSideVector.resize(LocalSize);
 
            noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize,LocalSize);
            noalias(rRightHandSideVector) = ZeroVector(LocalSize);
 
            this->ApplyInflowCondition(rLeftHandSideMatrix,rRightHandSideVector);
 
            this->ApplyWallLaw(rLeftHandSideMatrix,rRightHandSideVector);
        }
        else if(step == 5)
        {
            // Initialize local contributions
            const SizeType LocalSize = TNumNodes;
 
            if (rLeftHandSideMatrix.size1() != LocalSize)
                rLeftHandSideMatrix.resize(LocalSize,LocalSize);
            if (rRightHandSideVector.size() != LocalSize)
                rRightHandSideVector.resize(LocalSize);
 
            noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize,LocalSize);
            noalias(rRightHandSideVector) = ZeroVector(LocalSize);
 
            if(!this->Is(SLIP) )
            {
                const GeometryType& rGeom = this->GetGeometry();
                const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(GeometryData::IntegrationMethod::GI_GAUSS_2);
                const unsigned int NumGauss = IntegrationPoints.size();
                Vector GaussWeights = ZeroVector(NumGauss);
 
                MatrixType N = rGeom.ShapeFunctionsValues(GeometryData::IntegrationMethod::GI_GAUSS_2);
 
                array_1d<double,3> Normal;
                this->CalculateNormal(Normal); //this already contains the area
                double A = norm_2(Normal);
                Normal /= A;
 
                for (unsigned int g = 0; g < NumGauss; g++)
                    GaussWeights[g] = 2.0*A * IntegrationPoints[g].Weight();
 
                for (unsigned int g = 0; g < NumGauss; g++)
                {
                    double Weight = GaussWeights[g];
 
                    array_1d<double,3> vgauss = N(0,g)*rGeom[0].FastGetSolutionStepValue(VELOCITY);
                    for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
                    {
                        vgauss += N(iNode,g)*rGeom[iNode].FastGetSolutionStepValue(VELOCITY);
                    }
 
                    double aux =inner_prod(Normal,vgauss);
 
                    for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
                    {
                        rRightHandSideVector[iNode] -= N(iNode,g)*Weight*aux;
                    }
                }
            }
 
 
 
            // Initialize local contributions
//             const SizeType LocalSize = TNumNodes;
//
//             if (rLeftHandSideMatrix.size1() != LocalSize)
//                 rLeftHandSideMatrix.resize(LocalSize,LocalSize);
//             if (rRightHandSideVector.size() != LocalSize)
//                 rRightHandSideVector.resize(LocalSize);
//
//             noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize,LocalSize);
//             noalias(rRightHandSideVector) = ZeroVector(LocalSize);
//
//             if (!this->Is(SLIP))
//             {
//                 const double N = 1.0 / static_cast<double>(TNumNodes);
//                 array_1d<double,3> rNormal;
//                 this->CalculateNormal(rNormal); //this already contains the area
//
//                 for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
//                 {
//                     const array_1d<double,3>& rVelocity = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY);
//                     for (unsigned int d = 0; d < TDim; ++d)
//                         rRightHandSideVector[iNode] -= N*rVelocity[d]*rNormal[d];
//                 }
//             }
        }
        else
        {
            if (rLeftHandSideMatrix.size1() != 0)
                rLeftHandSideMatrix.resize(0,0,false);
 
            if (rRightHandSideVector.size() != 0)
                rRightHandSideVector.resize(0,false);
        }
    }
 
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const override
    {
        KRATOS_TRY;
 
        int Check = Condition::Check(rCurrentProcessInfo); // Checks id > 0 and area > 0
 
        if (Check != 0)
        {
            return Check;
        }
        else
        {
                // Checks on nodes
                // Check that the element's nodes contain all required SolutionStepData and Degrees of freedom
                for(unsigned int i=0; i<this->GetGeometry().size(); ++i)
                {
 
                    if(this->GetGeometry()[i].SolutionStepsDataHas(VELOCITY) == false)
                        KRATOS_THROW_ERROR(std::invalid_argument,"missing VELOCITY variable on solution step data for node ",this->GetGeometry()[i].Id());
                    if(this->GetGeometry()[i].SolutionStepsDataHas(MESH_VELOCITY) == false)
                        KRATOS_THROW_ERROR(std::invalid_argument,"missing MESH_VELOCITY variable on solution step data for node ",this->GetGeometry()[i].Id());
                    if(this->GetGeometry()[i].SolutionStepsDataHas(NORMAL) == false)
                        KRATOS_THROW_ERROR(std::invalid_argument,"missing NORMAL variable on solution step data for node ",this->GetGeometry()[i].Id());
                    if(this->GetGeometry()[i].HasDofFor(VELOCITY_X) == false ||
                            this->GetGeometry()[i].HasDofFor(VELOCITY_Y) == false ||
                            this->GetGeometry()[i].HasDofFor(VELOCITY_Z) == false)
                        KRATOS_THROW_ERROR(std::invalid_argument,"missing VELOCITY component degree of freedom on node ",this->GetGeometry()[i].Id());
                }
 
            return Check;
        }
 
        KRATOS_CATCH("");
    }
 
 
    void ApplyInflowCondition(MatrixType& rLocalMatrix,
                              VectorType& rLocalVector)
    {
        if (!this->Is(SLIP))
        {
            const unsigned int LocalSize = TNumNodes;
            const GeometryType& rGeom = this->GetGeometry();
            const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(GeometryData::IntegrationMethod::GI_GAUSS_2);
            const unsigned int NumGauss = IntegrationPoints.size();
            Vector GaussWeights = ZeroVector(NumGauss);
 
            MatrixType NContainer = rGeom.ShapeFunctionsValues(GeometryData::IntegrationMethod::GI_GAUSS_2);
 
            array_1d<double,3> Normal;
            this->CalculateNormal(Normal); //this already contains the area
            double A = std::sqrt(Normal[0]*Normal[0]+Normal[1]*Normal[1]+Normal[2]*Normal[2]);
            Normal /= A;
 
            for (unsigned int g = 0; g < NumGauss; g++)
                GaussWeights[g] = 2.0*A * IntegrationPoints[g].Weight();
 
            for (unsigned int g = 0; g < NumGauss; g++)
            {
                Vector N = row(NContainer,g);
                double Weight = GaussWeights[g];
 
                // Neumann boundary condition
                //             for (unsigned int i = 0; i < TNumNodes; i++)
                //             {
                //                 //unsigned int row = i*LocalSize;
                //                 const NodeType& rConstNode = this->GetGeometry()[i];
                //                 if ( rConstNode.IsFixed(PRESSURE)==true)
                //                 {
                //                     const double pext = rConstNode.FastGetSolutionStepValue(PRESSURE);
                //                     for (unsigned int j = 0; j < TNumNodes; j++)
                //                     {
                //                         unsigned int row = j*LocalSize;
                //                         for (unsigned int d = 0; d < TDim;d++)
                //                             rLocalVector[row+d] -= Weight*N[j]*N[i]*pext*Normal[d];
                //                     }
                //                 }
                //             }
 
                // Velocity inflow correction
                array_1d<double,3> Vel = ZeroVector(3);
                double Density = 0.0;
 
                for (unsigned int i = 0; i < TNumNodes; i++)
                {
                    const NodeType& rConstNode = this->GetGeometry()[i];
                    Vel += N[i]*rConstNode.FastGetSolutionStepValue(VELOCITY);
                    Density += N[i]*rConstNode.FastGetSolutionStepValue(DENSITY);
                }
 
                double Proj = Vel[0]*Normal[0] + Vel[1]*Normal[1] + Vel[2]*Normal[2];
 
                if (Proj < 0)
                {
                    const double W = Weight*Density*Proj;
                    for (unsigned int i = 0; i < TNumNodes; i++)
                    {
                        double row = i*LocalSize;
                        for (unsigned int j = 0; j < TNumNodes; j++)
                        {
                            double col = j*LocalSize;
                            const array_1d<double,3>& rVel = this->GetGeometry()[j].FastGetSolutionStepValue(VELOCITY);
                            for (unsigned int d = 0; d < TDim; d++)
                            {
                                double Tij = W*N[i]*N[j];
                                rLocalMatrix(row+d,col+d) -= Tij;
                                rLocalVector[row+d] += Tij*rVel[d];
                            }
                        }
                    }
                }
            }
        }
    }
 
 
 
    void EquationIdVector(EquationIdVectorType& rResult,
                                  const ProcessInfo& rCurrentProcessInfo) const override;
 
 
 
    void GetDofList(DofsVectorType& ConditionDofList,
                            const ProcessInfo& CurrentProcessInfo) const override;
 
 
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        this->PrintInfo(buffer);
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "WallConditionDiscontinuous" << TDim << "D #" << this->Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        this->pGetGeometry()->PrintData(rOStream);
    }
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
    friend class Serializer;
 
    void save(Serializer& rSerializer) const override
    {
        typedef WallCondition<TDim,TNumNodes> BaseType;
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType );
    }
 
    void load(Serializer& rSerializer) override
    {
        typedef WallCondition<TDim,TNumNodes> BaseType;
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType );
    }
 
 
 
 
 
 
 
 
 
 
 
 
}; // Class WallConditionDiscontinuous
 
 
 
 
 
 
 
template< unsigned int TDim, unsigned int TNumNodes >
inline std::istream& operator >> (std::istream& rIStream,
                                  WallConditionDiscontinuous<TDim,TNumNodes>& rThis)
{
    return rIStream;
}
 
template< unsigned int TDim, unsigned int TNumNodes >
inline std::ostream& operator << (std::ostream& rOStream,
                                  const WallConditionDiscontinuous<TDim,TNumNodes>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_WALL_CONDITION_DISCONTINUOUS_H