wave_element.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_WAVE_ELEMENT_H_INCLUDED
#define KRATOS_WAVE_ELEMENT_H_INCLUDED
 
// System includes
 
 
// External includes
 
 
// Project includes
#include "includes/define.h"
#include "includes/element.h"
#include "includes/serializer.h"
#include "custom_friction_laws/friction_law.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<std::size_t TNumNodes>
class WaveElement : public Element
{
public:
 
    typedef std::size_t IndexType;
 
    typedef Node NodeType;
 
    typedef array_1d<double, 3*TNumNodes> LocalVectorType;
 
    typedef BoundedMatrix<double, 3*TNumNodes, 3*TNumNodes> LocalMatrixType;
 
    typedef GeometryType::ShapeFunctionsGradientsType ShapeFunctionsGradientsType;
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(WaveElement);
 
 
 
    WaveElement() : Element(){}
 
    WaveElement(IndexType NewId, const NodesArrayType& ThisNodes) : Element(NewId, ThisNodes){}
 
    WaveElement(IndexType NewId, GeometryType::Pointer pGeometry) : Element(NewId, pGeometry){}
 
    WaveElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) : Element(NewId, pGeometry, pProperties){}
 
    ~ WaveElement() override {};
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override
    {
        return Kratos::make_intrusive<WaveElement<TNumNodes>>(NewId, this->GetGeometry().Create(ThisNodes), pProperties);
    }
 
    Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override
    {
        return Kratos::make_intrusive<WaveElement<TNumNodes>>(NewId, pGeom, pProperties);
    }
 
    Element::Pointer Clone(IndexType NewId, NodesArrayType const& ThisNodes) const override
    {
        Element::Pointer p_new_elem = Create(NewId, this->GetGeometry().Create(ThisNodes), this->pGetProperties());
        p_new_elem->SetData(this->GetData());
        p_new_elem->Set(Flags(*this));
        return p_new_elem;
    }
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
    void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const override;
 
    void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo& rCurrentProcessInfo) const override;
 
    void GetValuesVector(Vector& rValues, int Step = 0) const override;
 
    void GetFirstDerivativesVector(Vector& rValues, int Step = 0) const override;
 
    void GetSecondDerivativesVector(Vector& rValues, int Step = 0) const override;
 
    void CalculateOnIntegrationPoints(
        const Variable<double>& rVariable,
        std::vector<double>& rValues,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void Calculate(
        const Variable<array_1d<double,3>>& rVariable,
        array_1d<double,3>& rOutput,
        const ProcessInfo& rCurrentProcessInfo) override;
 
    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo) override;
 
    void CalculateDampingMatrix(MatrixType& rDampingMatrix, const ProcessInfo& rCurrentProcessInfo) override;
 
 
    GeometryData::IntegrationMethod GetIntegrationMethod() const override;
 
    const Parameters GetSpecifications() const override;
 
 
    std::string Info() const override
    {
        return "WaveElement";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info() << " : " << Id();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << GetGeometry();
    }
 
 
 
 
protected:
 
    static constexpr IndexType mLocalSize = 3 * TNumNodes;
 
 
 
 
    struct ElementData
    {
        bool integrate_by_parts;
        double stab_factor;
        double shock_stab_factor;
        double relative_dry_height;
        double gravity;
        double length;
        double absorbing_distance;
        double absorbing_damping;
 
        double depth;
        double height;
        array_1d<double,3> velocity;
 
        BoundedMatrix<double,3,3> A1;
        BoundedMatrix<double,3,3> A2;
        array_1d<double,3> b1;
        array_1d<double,3> b2;
 
        array_1d<double,TNumNodes> nodal_f;
        array_1d<double,TNumNodes> nodal_h;
        array_1d<double,TNumNodes> nodal_z;
        array_1d<double,TNumNodes> nodal_w;
        array_1d<array_1d<double,3>,TNumNodes> nodal_v;
        array_1d<array_1d<double,3>,TNumNodes> nodal_q;
        array_1d<array_1d<double,3>,TNumNodes> nodal_a;
        array_1d<array_1d<double,3>,TNumNodes> nodal_Jh;
        array_1d<array_1d<double,3>,TNumNodes> nodal_Ju;
        array_1d<array_1d<double,3>,TNumNodes> nodal_v_mesh;
 
        FrictionLaw::Pointer p_bottom_friction;
    };
 
 
    virtual const Variable<double>& GetUnknownComponent(int Index) const;
 
    virtual LocalVectorType GetUnknownVector(const ElementData& rData) const;
 
    virtual void InitializeData(ElementData& rData, const ProcessInfo& rCurrentProcessInfo);
 
    virtual void GetNodalData(ElementData& rData, const GeometryType& rGeometry, int Step = 0);
 
    virtual void UpdateGaussPointData(ElementData& rData, const array_1d<double,TNumNodes>& rN);
 
    void CalculateGeometryData(
        const GeometryType& rGeometry,
        Vector &rGaussWeights,
        Matrix &rNContainer,
        ShapeFunctionsGradientsType &rDN_DX);
 
    static double ShapeFunctionProduct(
        const array_1d<double,TNumNodes>& rN,
        const std::size_t I,
        const std::size_t J);
 
    static array_1d<double,3> VectorProduct(
        const array_1d<array_1d<double,3>,TNumNodes>& rV,
        const array_1d<double,TNumNodes>& rN);
 
    static array_1d<double,3> ScalarGradient(
        const array_1d<double,TNumNodes>& rS,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX);
 
    static double VectorDivergence(
        const array_1d<array_1d<double,3>,TNumNodes>& rV,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX);
 
    static BoundedMatrix<double,3,3> VectorGradient(
        const array_1d<array_1d<double,3>,TNumNodes>& rV,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX);
 
    static double InverseHeight(const ElementData& rData);
 
    static double WetFraction(const ElementData& rData);
 
    virtual void CalculateArtificialViscosity(
        BoundedMatrix<double,3,3>& rViscosity,
        BoundedMatrix<double,2,2>& rDiffusion,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX);
 
    virtual void CalculateArtificialDamping(
        BoundedMatrix<double,3,3>& rDamping,
        const ElementData& rData);
 
    void AddWaveTerms(
        LocalMatrixType& rMatrix,
        LocalVectorType& rVector,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX,
        const double Weight = 1.0);
 
    void AddFrictionTerms(
        LocalMatrixType& rMatrix,
        LocalVectorType& rVector,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX,
        const double Weight = 1.0);
 
    void AddArtificialViscosityTerms(
        LocalMatrixType& rMatrix,
        LocalVectorType& rVector,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX,
        const double Weight = 1.0);
 
    virtual void AddDispersiveTerms(
        LocalVectorType& rVector,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX,
        const double Weight = 1.0);
 
    virtual void AddMassTerms(
        LocalMatrixType& rMatrix,
        const ElementData& rData,
        const array_1d<double,TNumNodes>& rN,
        const BoundedMatrix<double,TNumNodes,2>& rDN_DX,
        const double Weight = 1.0);
 
    virtual double StabilizationParameter(const ElementData& rData) const;
 
 
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);
    }
 
 
 
 
}; // Class WaveElement
 
 
 
 
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_WAVE_ELEMENT_H_INCLUDED  defined