monolithic_dem_coupled_weak.h

Go to the documentation of this file.

/*
==============================================================================
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.
 
==============================================================================
 */
 
//
//   Project Name:        Kratos
//   Last Modified by:    $Author: gcasas $
//   Date:                $Date: 2010-10-09 10:34:00 $
//   Revision:            $Revision: 0.1 $
//
//
 
 
#ifndef MONOLITHIC_DEM_COUPLED_WEAK_H
#define MONOLITHIC_DEM_COUPLED_WEAK_H
 
// 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 "utilities/geometry_utilities.h"
 
// Application includes
#include "includes/cfd_variables.h"  //TODO: must be removed eventually
#include "fluid_dynamics_application_variables.h"
#include "swimming_dem_application_variables.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
template< unsigned int TDim,
          unsigned int TNumNodes = TDim + 1 >
class KRATOS_API(SWIMMING_DEM_APPLICATION) MonolithicDEMCoupledWeak : public Element
{
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(MonolithicDEMCoupledWeak);
 
    typedef IndexedObject BaseType;
 
    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;
 
//G
    typedef Kratos::Vector ShapeFunctionsType;
 
    typedef Kratos::Matrix ShapeFunctionDerivativesType;
 
    typedef GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType;
//Z
 
    //Constructors.
 
 
    MonolithicDEMCoupledWeak(IndexType NewId = 0) :
        Element(NewId)
    {}
 
 
    MonolithicDEMCoupledWeak(IndexType NewId, const NodesArrayType& ThisNodes) :
        Element(NewId, ThisNodes)
    {}
 
 
    MonolithicDEMCoupledWeak(IndexType NewId, GeometryType::Pointer pGeometry) :
        Element(NewId, pGeometry)
    {}
 
 
    MonolithicDEMCoupledWeak(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) :
        Element(NewId, pGeometry, pProperties)
    {}
 
    virtual ~MonolithicDEMCoupledWeak()
    {}
 
 
 
 
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,
                            PropertiesType::Pointer pProperties) const override
    {
 
        return Element::Pointer(new MonolithicDEMCoupledWeak(NewId, GetGeometry().Create(ThisNodes), pProperties));
 
    }
 
 
    virtual void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
                                      VectorType& rRightHandSideVector,
                                      const ProcessInfo& rCurrentProcessInfo) override
    {
        const unsigned int LocalSize = (TDim + 1) * TNumNodes;
 
        // Check sizes and initialize matrix
        if (rLeftHandSideMatrix.size1() != LocalSize)
            rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
 
        noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
 
        // Calculate RHS
        this->CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
    }
 
 
    virtual void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rCurrentProcessInfo) override
    {
        const unsigned int LocalSize = (TDim + 1) * TNumNodes;
 
        if (rLeftHandSideMatrix.size1() != LocalSize)
            rLeftHandSideMatrix.resize(LocalSize, LocalSize, false);
 
        noalias(rLeftHandSideMatrix) = ZeroMatrix(LocalSize, LocalSize);
    }
 
 
    virtual void CalculateRightHandSide(VectorType& rRightHandSideVector,
                                        const ProcessInfo& rCurrentProcessInfo) override
    {
        const unsigned int LocalSize = (TDim + 1) * TNumNodes;
 
        // Check sizes and initialize
        if (rRightHandSideVector.size() != LocalSize)
            rRightHandSideVector.resize(LocalSize, false);
 
        noalias(rRightHandSideVector) = ZeroVector(LocalSize);
 
        // Calculate this element's geometric parameters
        double Area;
        array_1d<double, TNumNodes> N;
        BoundedMatrix<double, TNumNodes, TDim> DN_DX;
        GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
        // Calculate this element's fluid properties
        double Density;
        this->EvaluateInPoint(Density, DENSITY, N);
 
        // Calculate Momentum RHS contribution
//G
        // Shape functions and integration points
        MatrixType NContainer;
        ShapeFunctionDerivativesArrayType DN_DXContainer;
        VectorType GaussWeights;
        this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
        const SizeType NumGauss = NContainer.size1();
 
        for (SizeType g = 0; g < NumGauss; g++)
        {
            const double GaussWeight = GaussWeights[g];
            const ShapeFunctionsType& Ng = row(NContainer, g);
            this->AddMomentumRHS(rRightHandSideVector, Density, Ng, GaussWeight);
        }
 
        const double& DeltaTime = rCurrentProcessInfo[DELTA_TIME];
        static const double arr[] = {0.5,0.5};
        std::vector<double> SchemeWeights (arr, arr + sizeof(arr) / sizeof(arr[0]));
 
        this->AddOtherContributionsRHS(rRightHandSideVector, N, SchemeWeights, DeltaTime);
//Z
 
        // For OSS: Add projection of residuals to RHS
        if (rCurrentProcessInfo[OSS_SWITCH] == 1)
        {
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            double KinViscosity;
            this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
            double Viscosity;
            this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
            // Calculate stabilization parameters
            double TauOne, TauTwo;
            this->CalculateTau(TauOne, TauTwo, AdvVel, Area,Density, Viscosity, rCurrentProcessInfo);
 
            this->AddProjectionToRHS(rRightHandSideVector, AdvVel, Density, TauOne, TauTwo, N, DN_DX, Area,rCurrentProcessInfo[DELTA_TIME]);
        }
        /*
//             else if (this->GetValue(TRACK_SUBSCALES) == 1)// Experimental: Dynamic tracking of ASGS subscales (see Codina 2002 Stabilized finite element ... using orthogonal subscales)
//             {
//                  * We want to evaluate v * d(u_subscale)/dt. This term is zero in OSS, due the orthogonality of the two terms. in ASGS, we approximate it as
//                  * d(u_s)/dt = u_s(last iteration) - u_s(previous time step) / DeltaTime where u_s(last iteration) = TauOne(MomResidual(last_iteration) + u_s(previous step)/DeltaTime)
//                  *
//                 array_1d<double, 3 > AdvVel;
//                 this->GetAdvectiveVel(AdvVel, N);
//
//                 double KinViscosity;
//                 this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
//
//                 double Viscosity;
//                 this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
//
//                 double StaticTauOne;
//                 this->CalculateStaticTau(StaticTauOne,AdvVel,Area,Viscosity);
//
//                 double TauOne,TauTwo;
//                 this->CalculateTau(TauOne,TauTwo,AdvVel,Area,Density,Viscosity,rCurrentProcessInfo);
//
//                 array_1d<double,3> ElemMomRes(3,0.0);
//                 this->ASGSMomResidual(AdvVel,Density,ElemMomRes,N,DN_DX,1.0);
//
//                 const array_1d<double,3>& rOldSubscale = this->GetValue(SUBSCALE_VELOCITY);
//                 const double DeltaTime = rCurrentProcessInfo[DELTA_TIME];
//
// //                array_1d<double,3> Subscale = TauOne * (ElemMomRes + Density * rOldSubscale / DeltaTime);
//                 const double C1 = 1.0 - TauOne/StaticTauOne;
//                 const double C2 = TauOne / (StaticTauOne*DeltaTime);
//
//                 unsigned int LocalIndex = 0;
//                 for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
//                 {
//                     for(unsigned int d = 0; d < TDim; d++)
//                         rRightHandSideVector[LocalIndex++] -= N[iNode] * Area * ( C1 * ElemMomRes[d] - C2 * rOldSubscale[d] );
// //                        rRightHandSideVector[LocalIndex++] -= N[iNode] * Area * ( Subscale[d] - rOldSubscale[d] ) / DeltaTime;
//                     LocalIndex++; // Pressure Dof
//                 }
//             }
*/
    }
 
 
    virtual void CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo) override
    {
//        rMassMatrix.resize(0,0,false);
        const unsigned int LocalSize = (TDim + 1) * TNumNodes;
 
        // Resize and set to zero
        if (rMassMatrix.size1() != LocalSize)
            rMassMatrix.resize(LocalSize, LocalSize, false);
 
        rMassMatrix = ZeroMatrix(LocalSize, LocalSize);
 
        // Get the element's geometric parameters
        double Area;
        array_1d<double, TNumNodes> N;
        BoundedMatrix<double, TNumNodes, TDim> DN_DX;
        GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
        // Calculate this element's fluid properties
        double Density;
        this->EvaluateInPoint(Density, DENSITY, N);
 
        // Add 'classical' mass matrix (lumped)
        //double Coeff = Density * Area / TNumNodes; //Optimize!
//GG
        // Shape functions and integration points
        MatrixType NContainer;
        ShapeFunctionDerivativesArrayType DN_DXContainer;
        VectorType GaussWeights;
        this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
        const SizeType NumGauss = NContainer.size1();
 
        for (SizeType g = 0; g < NumGauss; g++)
        {
            const double GaussWeight = GaussWeights[g];
            const ShapeFunctionsType& Ng = row(NContainer, g);
            const ShapeFunctionDerivativesType& DN_DXg = DN_DXContainer[g];
            AddConsistentMassMatrixContribution(rMassMatrix, Ng, Density, GaussWeight);
 
            if (rCurrentProcessInfo[OSS_SWITCH] != 1)
            {
                double KinViscosity;
                this->EvaluateInPoint(KinViscosity, VISCOSITY, Ng);
 
                double Viscosity;
                this->GetEffectiveViscosity(Density,KinViscosity, Ng, DN_DXg, Viscosity, rCurrentProcessInfo);
 
                // Get Advective velocity
                array_1d<double, 3 > AdvVel;
 
                this->GetAdvectiveVel(AdvVel, Ng);
 
                // Calculate stabilization parameters
                double TauOne, TauTwo;
                this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
 
                // Add dynamic stabilization terms ( all terms involving a delta(u) )
                this->AddMassStabTerms(rMassMatrix, Density, AdvVel, TauOne, Ng, DN_DXg, GaussWeight);
            }
        }
/*
        //this->CalculateLumpedMassMatrix(rMassMatrix, Coeff);
 
 
         For ASGS: add dynamic stabilization terms.
         These terms are not used in OSS, as they belong to the finite element
         space and cancel out with their projections.
         //
          if (rCurrentProcessInfo[OSS_SWITCH] != 1)
          {
              double KinViscosity;
              this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
              double Viscosity;
              this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
              // Get Advective velocity
              array_1d<double, 3 > AdvVel;
 
              this->GetAdvectiveVel(AdvVel, N);
 
              // Calculate stabilization parameters
              double TauOne, TauTwo;
              this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
 
              // Add dynamic stabilization terms ( all terms involving a delta(u) )
              this->AddMassStabTerms(rMassMatrix, Density, AdvVel, TauOne, N, DN_DX, Area);
          }
 
*/
//ZZ
    }
 
 
    virtual void CalculateLocalVelocityContribution(MatrixType& rDampingMatrix,
            VectorType& rRightHandSideVector,
            const ProcessInfo& rCurrentProcessInfo) override
    {
        const unsigned int LocalSize = (TDim + 1) * TNumNodes;
 
        // Resize and set to zero the matrix
        // Note that we don't clean the RHS because it will already contain body force (and stabilization) contributions
        if (rDampingMatrix.size1() != LocalSize)
            rDampingMatrix.resize(LocalSize, LocalSize, false);
 
        noalias(rDampingMatrix) = ZeroMatrix(LocalSize, LocalSize);
 
        // Get this element's geometric properties
        double Area;
        array_1d<double, TNumNodes> N;
        BoundedMatrix<double, TNumNodes, TDim> DN_DX;
        GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
        // Calculate this element's fluid properties
        double Density, KinViscosity;
        this->EvaluateInPoint(Density, DENSITY, N);
        this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
        double Viscosity;
        this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
        // Get Advective velocity
        array_1d<double, 3 > AdvVel;
        this->GetAdvectiveVel(AdvVel, N);
 
        // Calculate stabilization parameters
        double TauOne, TauTwo;
        this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
 
        /*            if(this->GetValue(TRACK_SUBSCALES)==1)
                    {
                        const double DeltaTime = rCurrentProcessInfo[DELTA_TIME];
                        this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, N, DN_DX, Area,DeltaTime);
                    }
                    else
                    {*/
//GG
        // Shape functions and integration points
        MatrixType NContainer;
        ShapeFunctionDerivativesArrayType DN_DXContainer;
        VectorType GaussWeights;
        this->CalculateWeights(DN_DXContainer, NContainer, GaussWeights);
        const SizeType NumGauss = NContainer.size1();
 
        for (SizeType g = 0; g < NumGauss; g++)
        {
            const double GaussWeight = GaussWeights[g];
            const ShapeFunctionsType& Ng = row(NContainer, g);
            const ShapeFunctionDerivativesType& DN_DXg = DN_DXContainer[g];
            this->GetAdvectiveVel(AdvVel, Ng);
            this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, Ng, DN_DXg, GaussWeight);
        }
 
        //this->AddIntegrationPointVelocityContribution(rDampingMatrix, rRightHandSideVector, Density, Viscosity, AdvVel, TauOne, TauTwo, N, DN_DX, Area);
//ZZ
 
//        this->ModulatedGradientDiffusion(rDampingMatrix,DN_DX,Density*Area);
//             }
 
        // Now calculate an additional contribution to the residual: r -= rDampingMatrix * (u,p)
        VectorType U = ZeroVector(LocalSize);
        int LocalIndex = 0;
 
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
        {
            array_1d< double, 3 > & rVel = this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY);
            for (unsigned int d = 0; d < TDim; ++d) // Velocity Dofs
            {
                U[LocalIndex] = rVel[d];
                ++LocalIndex;
            }
            U[LocalIndex] = this->GetGeometry()[iNode].FastGetSolutionStepValue(PRESSURE); // Pressure Dof
            ++LocalIndex;
        }
 
        noalias(rRightHandSideVector) -= prod(rDampingMatrix, U);
    }
 
    virtual void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override
    {
//             if(this->GetValue(TRACK_SUBSCALES) == 1)
//             {
//                 // Get this element's geometric properties
//                 double Area;
//                 array_1d<double, TNumNodes> N;
//                 BoundedMatrix<double, TNumNodes, TDim> DN_DX;
//                 GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
//
//                 // Calculate this element's fluid properties
//                 double Density, KinViscosity;
//                 this->EvaluateInPoint(Density, DENSITY, N);
//                 this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
//
//                 double Viscosity;
//                 this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
//
//                 // Get Advective velocity
//                 array_1d<double, 3 > AdvVel;
//                 this->GetAdvectiveVel(AdvVel, N);
//
//                 // Calculate stabilization parameters
//                 double StaticTauOne;
//                 this->CalculateStaticTau(StaticTauOne,AdvVel,Area,Viscosity);
//
//                 array_1d<double,3> ElementalMomRes(3,0.0);
//
//                 if ( rCurrentProcessInfo[OSS_SWITCH] != 1 ) // ASGS
//                 {
//                     this->ASGSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);
//                 }
//                 else // OSS
//                 {
//                     this->OSSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);;
//                 }
//
//                 // Update subscale term
//                 const double DeltaTime = rCurrentProcessInfo.GetValue(DELTA_TIME);
//                 array_1d<double,3>& OldSubscaleVel = this->GetValue(SUBSCALE_VELOCITY);
//                 array_1d<double,3> Tmp = ( 1.0/( 1.0/DeltaTime + 1.0/StaticTauOne)) * (Density*OldSubscaleVel/DeltaTime + ElementalMomRes);
//                 OldSubscaleVel = Tmp;
//             }
    }
 
 
    virtual void Calculate(const Variable<double>& rVariable,
                           double& rOutput,
                           const ProcessInfo& rCurrentProcessInfo) override
    {
        if (rVariable == ERROR_RATIO)
        {
            // Get the element's geometric parameters
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            // Calculate this element's fluid properties
            double Density, KinViscosity;
            this->EvaluateInPoint(Density, DENSITY, N);
            this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
            double Viscosity;
            this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
            // Get Advective velocity
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            // Output container
            array_1d< double, 3 > ElementalMomRes(3, 0.0);
 
            // Calculate stabilization parameter. Note that to estimate the subscale velocity, the dynamic coefficient in TauOne is assumed zero.
            double TauOne;
            this->CalculateStaticTau(TauOne, AdvVel, Area,Density, Viscosity);
 
            if ( rCurrentProcessInfo[OSS_SWITCH] != 1 ) // ASGS
            {
                this->ASGSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);
                ElementalMomRes *= TauOne;
            }
            else // OSS
            {
                this->OSSMomResidual(AdvVel,Density,ElementalMomRes,N,DN_DX,1.0);;
                ElementalMomRes *= TauOne;
            }
 
            // Error estimation ( ||U'|| / ||Uh_gauss|| ), taking ||U'|| = TauOne ||MomRes||
            double ErrorRatio(0.0);//, UNorm(0.0);
//                array_1d< double, 3 > UGauss(3, 0.0);
//                this->AddPointContribution(UGauss, VELOCITY, N);
 
            for (unsigned int i = 0; i < TDim; ++i)
            {
                ErrorRatio += ElementalMomRes[i] * ElementalMomRes[i];
//                    UNorm += UGauss[i] * UGauss[i];
            }
            ErrorRatio = sqrt(ErrorRatio); // / UNorm);
            ErrorRatio /= Density;
            this->SetValue(ERROR_RATIO, ErrorRatio);
            rOutput = ErrorRatio;
        }
        else if (rVariable == NODAL_AREA)
        {
            // Get the element's geometric parameters
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            // Carefully write results to nodal variables, to avoid parallelism problems
            for (unsigned int i = 0; i < TNumNodes; ++i)
            {
                this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
                this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
                this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
            }
        }
    }
 
 
    virtual void Calculate(const Variable<array_1d<double, 3 > >& rVariable,
                           array_1d<double, 3 > & rOutput,
                           const ProcessInfo& rCurrentProcessInfo) override
    {
        if (rVariable == ADVPROJ) // Compute residual projections for OSS
        {
            // Get the element's geometric parameters
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            // Calculate this element's fluid properties
            double Density;
            this->EvaluateInPoint(Density, DENSITY, N);
 
            // Get Advective velocity
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            // Output containers
            array_1d< double, 3 > ElementalMomRes(3, 0.0);
            double ElementalMassRes(0);
 
            this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes, N, DN_DX, Area);
 
            if (rCurrentProcessInfo[OSS_SWITCH] == 1)
            {
                // Carefully write results to nodal variables, to avoid parallelism problems
                for (unsigned int i = 0; i < TNumNodes; ++i)
                {
                    this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
                    array_1d< double, 3 > & rAdvProj = this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ);
                    for (unsigned int d = 0; d < TDim; ++d)
                        rAdvProj[d] += N[i] * ElementalMomRes[d];
 
                    this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ) += N[i] * ElementalMassRes;
                    this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
                    this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
                }
            }
 
            rOutput = ElementalMomRes;
        }
        else if (rVariable == SUBSCALE_VELOCITY)
        {
            // Get the element's geometric parameters
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            // Calculate this element's fluid properties
            double Density;
            this->EvaluateInPoint(Density, DENSITY, N);
 
            // Get Advective velocity
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            // Output containers
            array_1d< double, 3 > ElementalMomRes(3,0.0);
            double ElementalMassRes(0.0);
 
            this->AddProjectionResidualContribution(AdvVel, Density, ElementalMomRes, ElementalMassRes, N, DN_DX, Area);
 
            if (rCurrentProcessInfo[OSS_SWITCH] == 1)
            {
                /* Projections of the elemental residual are computed with
                 * Newton-Raphson iterations of type M(lumped) dx = ElemRes - M(consistent) * x
                 */
                const double Weight = ConsistentMassCoef(Area); // Consistent mass matrix is Weigth * ( Ones(TNumNodes,TNumNodes) + Identity(TNumNodes,TNumNodes) )
                // Carefully write results to nodal variables, to avoid parallelism problems
                for (unsigned int i = 0; i < TNumNodes; ++i)
                {
                    this->GetGeometry()[i].SetLock(); // So it is safe to write in the node in OpenMP
 
                    // Add elemental residual to RHS
                    array_1d< double, 3 > & rMomRHS = this->GetGeometry()[i].GetValue(ADVPROJ);
                    double& rMassRHS = this->GetGeometry()[i].GetValue(DIVPROJ);
                    for (unsigned int d = 0; d < TDim; ++d)
                        rMomRHS[d] += N[i] * ElementalMomRes[d];
 
                    rMassRHS += N[i] * ElementalMassRes;
 
                    // Write nodal area
                    this->GetGeometry()[i].FastGetSolutionStepValue(NODAL_AREA) += Area * N[i];
 
                    // Substract M(consistent)*x(i-1) from RHS
                    for(unsigned int j = 0; j < TNumNodes; ++j) // RHS -= Weigth * Ones(TNumNodes,TNumNodes) * x(i-1)
                    {
                        for(unsigned int d = 0; d < TDim; ++d)
                            rMomRHS[d] -= Weight * this->GetGeometry()[j].FastGetSolutionStepValue(ADVPROJ)[d];
                        rMassRHS -= Weight * this->GetGeometry()[j].FastGetSolutionStepValue(DIVPROJ);
                    }
                    for(unsigned int d = 0; d < TDim; ++d) // RHS -= Weigth * Identity(TNumNodes,TNumNodes) * x(i-1)
                        rMomRHS[d] -= Weight * this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ)[d];
                    rMassRHS -= Weight * this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ);
 
                    this->GetGeometry()[i].UnSetLock(); // Free the node for other threads
                }
            }
 
            rOutput = ElementalMomRes;
        }
    }
 
    // The following methods have different implementations depending on TDim
 
    virtual void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo& rCurrentProcessInfo) const override;
 
 
    virtual void GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo& rCurrentProcessInfo) const override;
 
 
    virtual void GetFirstDerivativesVector(Vector& Values, int Step = 0) const override;
 
 
    virtual void GetSecondDerivativesVector(Vector& Values, int Step = 0) const override;
 
 
    virtual void CalculateOnIntegrationPoints(const Variable<array_1d<double, 3 > >& rVariable,
            std::vector<array_1d<double, 3 > >& rOutput,
            const ProcessInfo& rCurrentProcessInfo) override;
 
 
    virtual void CalculateOnIntegrationPoints(const Variable<double>& rVariable,
            std::vector<double>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
    {
        if (rVariable == TAUONE || rVariable == TAUTWO || rVariable == MU)
        {
            double TauOne, TauTwo;
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            double Density,KinViscosity;
            this->EvaluateInPoint(Density, DENSITY, N);
            this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
            double Viscosity;
            this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
            this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
 
            rValues.resize(1, false);
            if (rVariable == TAUONE)
            {
                rValues[0] = TauOne;
            }
            else if (rVariable == TAUTWO)
            {
                rValues[0] = TauTwo;
            }
            else if (rVariable == MU)
            {
                rValues[0] = Density * Viscosity;
            }
        }
        else if(rVariable == SUBSCALE_PRESSURE)
        {
            double TauOne, TauTwo;
            double Area;
            array_1d<double, TNumNodes> N;
            BoundedMatrix<double, TNumNodes, TDim> DN_DX;
            GeometryUtils::CalculateGeometryData(this->GetGeometry(), DN_DX, N, Area);
 
            array_1d<double, 3 > AdvVel;
            this->GetAdvectiveVel(AdvVel, N);
 
            double Density,KinViscosity;
            this->EvaluateInPoint(Density, DENSITY, N);
            this->EvaluateInPoint(KinViscosity, VISCOSITY, N);
 
            double Viscosity;
            this->GetEffectiveViscosity(Density,KinViscosity, N, DN_DX, Viscosity, rCurrentProcessInfo);
 
            this->CalculateTau(TauOne, TauTwo, AdvVel, Area, Density, Viscosity, rCurrentProcessInfo);
 
            double DivU = 0.0;
            for(unsigned int i=0; i < TNumNodes; i++)
            {
                for(unsigned int d = 0; d < TDim; d++)
                    DivU -= DN_DX(i,d) * this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY)[d];
            }
 
            rValues.resize(1, false);
 
            rValues[0] = TauTwo * DivU;// *Density?? decide on criteria and use the same for SUBSCALE_VELOCITY
 
            if(rCurrentProcessInfo[OSS_SWITCH]==1)
            {
                double Proj = 0.0;
                for(unsigned int i=0; i < TNumNodes; i++)
                {
                    Proj += N[i]*this->GetGeometry()[i].FastGetSolutionStepValue(DIVPROJ);
                }
                rValues[0] -= TauTwo*Proj;
            }
        }
        else if (rVariable == NODAL_AREA && TDim == 3)
        {
            MatrixType J = ZeroMatrix(3,3);
            const array_1d<double,3>& X0 = this->GetGeometry()[0].Coordinates();
            const array_1d<double,3>& X1 = this->GetGeometry()[1].Coordinates();
            const array_1d<double,3>& X2 = this->GetGeometry()[2].Coordinates();
            const array_1d<double,3>& X3 = this->GetGeometry()[3].Coordinates();
 
            J(0,0) = X1[0]-X0[0];
            J(0,1) = X2[0]-X0[0];
            J(0,2) = X3[0]-X0[0];
            J(1,0) = X1[1]-X0[1];
            J(1,1) = X2[1]-X0[1];
            J(1,2) = X3[1]-X0[1];
            J(2,0) = X1[2]-X0[2];
            J(2,1) = X2[2]-X0[2];
            J(2,2) = X3[2]-X0[2];
 
            double DetJ = J(0,0)*( J(1,1)*J(2,2) - J(1,2)*J(2,1) ) + J(0,1)*( J(1,2)*J(2,0) - J(1,0)*J(2,2) ) + J(0,2)*( J(1,0)*J(2,1) - J(1,1)*J(2,0) );
            rValues.resize(1, false);
            rValues[0] = DetJ;
        }
        else // Default behaviour (returns elemental data)
        {
            rValues.resize(1, false);
            /*
             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 MonolithicDEMCoupledWeak<TDim, TNumNodes>* const_this = static_cast<const MonolithicDEMCoupledWeak<TDim, TNumNodes>*> (this);
            rValues[0] = const_this->GetValue(rVariable);
        }
    }
 
    virtual void CalculateOnIntegrationPoints(const Variable<array_1d<double, 6 > >& rVariable,
            std::vector<array_1d<double, 6 > >& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
    {}
 
    virtual void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
            std::vector<Vector>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
    {}
 
    virtual void CalculateOnIntegrationPoints(const Variable<Matrix>& rVariable,
            std::vector<Matrix>& rValues,
            const ProcessInfo& rCurrentProcessInfo) override
    {}
 
 
 
 
    virtual int Check(const ProcessInfo& rCurrentProcessInfo) const override
    {
        KRATOS_TRY
 
        // Perform basic element checks
        int ErrorCode = Kratos::Element::Check(rCurrentProcessInfo);
        if(ErrorCode != 0) return ErrorCode;
 
        // Additional variables, only required to print results:
        // SUBSCALE_VELOCITY, SUBSCALE_PRESSURE, TAUONE, TAUTWO, MU, VORTICITY.
 
        // 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) {
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(VELOCITY, this->GetGeometry()[i])
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(PRESSURE, this->GetGeometry()[i])
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(MESH_VELOCITY, this->GetGeometry()[i])
            KRATOS_CHECK_VARIABLE_IN_NODAL_DATA(ACCELERATION, this->GetGeometry()[i])
 
            KRATOS_CHECK_DOF_IN_NODE(VELOCITY_X, this->GetGeometry()[i])
            KRATOS_CHECK_DOF_IN_NODE(VELOCITY_Y, this->GetGeometry()[i])
            KRATOS_CHECK_DOF_IN_NODE(VELOCITY_Z, this->GetGeometry()[i])
            KRATOS_CHECK_DOF_IN_NODE(PRESSURE, this->GetGeometry()[i])
        }
        // Not checking OSS related variables NODAL_AREA, ADVPROJ, DIVPROJ, which are only required as SolutionStepData if OSS_SWITCH == 1
 
        // If this is a 2D problem, check that nodes are in XY plane
        if (this->GetGeometry().WorkingSpaceDimension() == 2) {
            for (unsigned int i=0; i<this->GetGeometry().size(); ++i) {
                KRATOS_ERROR_IF(this->GetGeometry()[i].Z() != 0.0) << "Node with non-zero Z coordinate found. Id: "<< this->GetGeometry()[i].Id() << std::endl;
            }
        }
 
        return 0;
 
        KRATOS_CATCH("");
    }
 
 
 
 
    virtual std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "MonolithicDEMCoupledWeak #" << Id();
        return buffer.str();
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "MonolithicDEMCoupledWeak" << TDim << "D";
    }
 
//        /// Print object's data.
//        virtual void PrintData(std::ostream& rOStream) const;
 
 
 
 
protected:
 
//GG
    void CalculateWeights(ShapeFunctionDerivativesArrayType& rDN_DX, Matrix& rNContainer, Vector& rGaussWeights);
//ZZ
 
 
 
 
 
 
 
    virtual void CalculateTau(double& TauOne,
                              double& TauTwo,
                              const array_1d< double, 3 > & rAdvVel,
                              const double Area,
                              const double Density,
                              const double KinViscosity,
                              const ProcessInfo& rCurrentProcessInfo)
    {
        // Compute mean advective velocity norm
        double AdvVelNorm = 0.0;
        for (unsigned int d = 0; d < TDim; ++d)
            AdvVelNorm += rAdvVel[d] * rAdvVel[d];
 
        AdvVelNorm = sqrt(AdvVelNorm);
 
//GG
       // const double Element_Diam = sqrt(this->FilterWidth());
//ZZ
        const double Element_Size = this->ElementSize(Area);
 
//G
    //  TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 4 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
//GG
    //  TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 5.6666666666 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
        double L0 = 0.015; // 0.1 * (volumen probeta) ^(1/3)
        double lp = sqrt(L0 * Element_Size);
        double lu = lp;
        double c2u = 2;
        double c2p = c2u;
        double sigma = 0.0;
        //this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
        if (sigma > 0.00001){
            TauOne = Element_Size * Element_Size / (c2u * sigma * lu * lu);
            TauTwo = lp * lp * sigma * c2p;
          }
        else{
           TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 5.6666666666 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
           TauTwo = Density * (KinViscosity + 0.5 * Element_Size * AdvVelNorm);
          }
//ZZ
//Z
 
//GG
    //   TauTwo = Density * (KinViscosity + 0.5 * Element_Size * AdvVelNorm);
 
//ZZ
        //TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 12.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
        //TauTwo = Density * (KinViscosity + Element_Size * AdvVelNorm / 6.0);
 
 
    }
 
 
    virtual void CalculateStaticTau(double& TauOne,
                                    const array_1d< double, 3 > & rAdvVel,
                                    const double Area,
                                    const double Density,
                                    const double KinViscosity)
    {
        // Compute mean advective velocity norm
        double AdvVelNorm = 0.0;
        for (unsigned int d = 0; d < TDim; ++d)
            AdvVelNorm += rAdvVel[d] * rAdvVel[d];
 
        AdvVelNorm = sqrt(AdvVelNorm);
 
//GG
       // const double Element_Diam = sqrt(this->FilterWidth());
//ZZ
        const double Element_Size = this->ElementSize(Area);
 
        TauOne = 1.0 / (Density*(4.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size));
//G
    //  TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 4 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
//GG
    //  TauOne = 1.0 / (Density * ( rCurrentProcessInfo[DYNAMIC_TAU] / rCurrentProcessInfo[DELTA_TIME] + 5.6666666666 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size) );
        double L0 = 0.015;
        double lp = sqrt(L0 * Element_Size);
        double lu = lp;
        double c2u = 2;
        double sigma = 0.0;
        //this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
 
        if (sigma > 0.00001){
            TauOne = Element_Size * Element_Size / (c2u * sigma * lu * lu);
          }
        else{
            TauOne = 1.0 / (Density*(4.0 * KinViscosity / (Element_Size * Element_Size) + 2.0 * AdvVelNorm / Element_Size));
          }
//ZZ
    }
 
    virtual void AddMomentumRHS(VectorType& F,
                                const double Density,
                                const array_1d<double, TNumNodes>& rShapeFunc,
                                const double Weight)
    {
        double Coef = Density * Weight;
 
        array_1d<double, 3 > BodyForce(3, 0.0);
        this->EvaluateInPoint(BodyForce, BODY_FORCE, rShapeFunc);
 
        // Add the results to the velocity components (Local Dofs are vx, vy, [vz,] p for each node)
        int LocalIndex = 0;
 
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
        {
            for (unsigned int d = 0; d < TDim; ++d)
            {
                F[LocalIndex++] += Coef * rShapeFunc[iNode] * BodyForce[d];
            }
            ++LocalIndex; // Skip pressure Dof
        }
    }
//G
    virtual void AddOtherContributionsRHS(VectorType& F,
    const array_1d<double, TNumNodes>& rShapeFunc,
    const std::vector<double>& TimeSchemeWeights,
    const double& DeltaTime)
    {
 
      double FluidFractionRate;
      this->EvaluateTimeDerivativeInPoint(FluidFractionRate, FLUID_FRACTION_RATE, rShapeFunc, DeltaTime, TimeSchemeWeights);
 
      // Add the results to the velocity components (Local Dofs are vx, vy, [vz,] p for each node)
      int LocalIndex = 0;
 
      for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
 
          for (unsigned int d = 0; d < TDim; ++d){
              F[LocalIndex++] -= FluidFractionRate;
          }
 
          ++LocalIndex; // Skip pressure Dof
      }
 
    }
//Z
 
    virtual void AddProjectionToRHS(VectorType& RHS,
                                    const array_1d<double, 3 > & rAdvVel,
                                    const double Density,
                                    const double TauOne,
                                    const double TauTwo,
                                    const array_1d<double, TNumNodes>& rShapeFunc,
                                    const BoundedMatrix<double, TNumNodes, TDim>& rShapeDeriv,
                                    const double Weight,
                                    const double DeltaTime = 1.0)
    {
        const unsigned int BlockSize = TDim + 1;
 
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
 
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
        array_1d<double,3> MomProj(3,0.0);
        double DivProj = 0.0;
        this->EvaluateInPoint(MomProj,ADVPROJ,rShapeFunc);
        this->EvaluateInPoint(DivProj,DIVPROJ,rShapeFunc);
 
        MomProj *= TauOne;
        DivProj *= TauTwo;
 
        unsigned int FirstRow = 0;
 
        for (unsigned int i = 0; i < TNumNodes; i++)
        {
//G
            double FluidFraction = this->GetGeometry()[i].FastGetSolutionStepValue(FLUID_FRACTION);
            array_1d<double,3> FluidFractionGradient(3,0.0);
            FluidFractionGradient[0] += rShapeDeriv(i, 0) * FluidFraction;
            FluidFractionGradient[1] += rShapeDeriv(i, 1) * FluidFraction;
            FluidFractionGradient[2] += rShapeDeriv(i, 2) * FluidFraction;
//Z
 
            for (unsigned int d = 0; d < TDim; d++)
            {
//G
//              RHS[FirstRow+d] -= Weight * (Density * AGradN[i] * MomProj[d] + rShapeDeriv(i,d) * DivProj); // TauOne * ( a * Grad(v) ) * MomProjection + TauTwo * Div(v) * MassProjection
//GG
    //          RHS[FirstRow+d] -= Weight * (Density * AGradN[i] * MomProj[d] + (FluidFraction * rShapeDeriv(i,d) + FluidFractionGradient[d] * rShapeFunc[i]) * DivProj); // TauOne * (a * Grad(v)) * MomProjection + TauTwo * Div(v) * MassProjection
                RHS[FirstRow+d] -= Weight * ((Density * AGradN[i] - sigma * rShapeFunc[i]) * MomProj[d] + (FluidFraction * rShapeDeriv(i,d) + FluidFractionGradient[d] * rShapeFunc[i]) * DivProj); // TauOne * (a * Grad(v)) * MomProjection + TauTwo * Div(v) * MassProjection
//ZZ
//Z
                RHS[FirstRow+TDim] -= Weight * rShapeDeriv(i,d) * MomProj[d]; // TauOne * Grad(q) * MomProjection
            }
            FirstRow += BlockSize;
        }
    }
 
 
    void CalculateLumpedMassMatrix(MatrixType& rLHSMatrix,
                                   const double Mass)
    {
        unsigned int DofIndex = 0;
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
        {
            for (unsigned int d = 0; d < TDim; ++d)
            {
                rLHSMatrix(DofIndex, DofIndex) += Mass;
                ++DofIndex;
            }
            ++DofIndex; // Skip pressure Dof
        }
    }
 
    void AddConsistentMassMatrixContribution(MatrixType& rLHSMatrix,
            const array_1d<double,TNumNodes>& rShapeFunc,
            const double Density,
            const double Weight)
    {
        const unsigned int BlockSize = TDim + 1;
 
        double Coef = Density * Weight;
        unsigned int FirstRow(0), FirstCol(0);
        double K; // Temporary results
 
        // Note: Dof order is (vx,vy,[vz,]p) for each node
        for (unsigned int i = 0; i < TNumNodes; ++i)
        {
            // Loop over columns
            for (unsigned int j = 0; j < TNumNodes; ++j)
            {
                // Delta(u) * TauOne * [ AdvVel * Grad(v) ] in velocity block
                K = Coef * rShapeFunc[i] * rShapeFunc[j];
 
                for (unsigned int d = 0; d < TDim; ++d) // iterate over dimensions for velocity Dofs in this node combination
                {
                    rLHSMatrix(FirstRow + d, FirstCol + d) += K;
                }
                // Update column index
                FirstCol += BlockSize;
            }
            // Update matrix indices
            FirstRow += BlockSize;
            FirstCol = 0;
        }
    }
 
 
 
    void AddMassStabTerms(MatrixType& rLHSMatrix,
                          const double Density,
                          const array_1d<double, 3 > & rAdvVel,
                          const double TauOne,
                          const array_1d<double, TNumNodes>& rShapeFunc,
                          const BoundedMatrix<double, TNumNodes, TDim>& rShapeDeriv,
                          const double Weight)
    {
        const unsigned int BlockSize = TDim + 1;
 
        double Coef = Weight * TauOne;
        unsigned int FirstRow(0), FirstCol(0);
        double K; // Temporary results
 
        // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//G
        double FluidFraction;
        array_1d<double,3> FluidFractionGradient(3,0.0);
        this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
        this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
 
        for (unsigned int i = 0; i < TNumNodes; ++i) {
            this->GetGeometry()[i].FastGetSolutionStepValue(FLUID_FRACTION_GRADIENT) = FluidFractionGradient;
        }
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
//Z
 
        // Note: Dof order is (vx,vy,[vz,]p) for each node
        for (unsigned int i = 0; i < TNumNodes; ++i)
        {
            // Loop over columns
            for (unsigned int j = 0; j < TNumNodes; ++j)
            {
                // Delta(u) * TauOne * [ AdvVel * Grad(v) ] in velocity block
                K = Coef * Density * AGradN[i] * Density * rShapeFunc[j];
                // K -= Coef * sigma * rShapeFunc[i] * rShapeFunc[j];
 
                for (unsigned int d = 0; d < TDim; ++d) // iterate over dimensions for velocity Dofs in this node combination
                {
                    rLHSMatrix(FirstRow + d, FirstCol + d) += K;
                    // Delta(u) * TauOne * Grad(q) in q * Div(u) block
//G
    //              rLHSMatrix(FirstRow + TDim, FirstCol + d) += Coef * Density * rShapeDeriv(i, d) * rShapeFunc[j];
                    rLHSMatrix(FirstRow + TDim, FirstCol + d) += Coef * Density * (FluidFraction * rShapeDeriv(i, d) + FluidFractionGradient[d] * rShapeFunc[i]) * rShapeFunc[j];
//Z
                }
                // Update column index
                FirstCol += BlockSize;
            }
            // Update matrix indices
            FirstRow += BlockSize;
            FirstCol = 0;
        }
    }
 
    void AddIntegrationPointVelocityContribution(MatrixType& rDampingMatrix,
            VectorType& rDampRHS,
            const double Density,
            const double Viscosity,
            const array_1d< double, 3 > & rAdvVel,
            const double TauOne,
            const double TauTwo,
            const array_1d< double, TNumNodes >& rShapeFunc,
            const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
            const double Weight)
    {
        const unsigned int BlockSize = TDim + 1;
 
        // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
        // Build the local matrix and RHS
        unsigned int FirstRow(0), FirstCol(0); // position of the first term of the local matrix that corresponds to each node combination
        double K, G, PDivV, L, qF; // Temporary results
 
        array_1d<double,3> BodyForce(3,0.0);
        this->EvaluateInPoint(BodyForce,BODY_FORCE,rShapeFunc);
        BodyForce *= Density;
//G
        double DivPEpsilon, GEps, FluidFraction;
        array_1d<double,3> FluidFractionGradient(3,0.0);
        this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
        this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
//GG
        // Lumped version for 1 integration point of the sigma u*v terms
//        unsigned int DofIndex = 0;
//        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
//        {
//            for (unsigned int d = 0; d < TDim; ++d)
//            {
//                rDampingMatrix(DofIndex, DofIndex) += Weight * (sigma - TauOne * sigma * sigma) * rShapeFunc[iNode]; // lumped version of sigma * u * v + Stabilization (- sigma * v * TauOne * sigma * u)
//                ++DofIndex;
//            }
//            ++DofIndex; // Skip pressure Dof
//        }
//ZZ
//Z
 
        for (unsigned int i = 0; i < TNumNodes; ++i) // iterate over rows
        {
            for (unsigned int j = 0; j < TNumNodes; ++j) // iterate over columns
            {
                // Calculate the part of the contributions that is constant for each node combination
 
                // Velocity block
                K = Density * rShapeFunc[i] * AGradN[j]; // Convective term: v * ( a * Grad(u) )
 
                //K = 0.5 * Density * (rShapeFunc[i] * AGradN[j] - AGradN[i] * rShapeFunc[j]); // Skew-symmetric convective term 1/2( v*grad(u)*u - grad(v) uu )
                K += TauOne * Density * AGradN[i] * Density * AGradN[j]; // Stabilization: (a * Grad(v)) * TauOne * (a * Grad(u))
//GG
                K += (sigma - TauOne * sigma * sigma) * rShapeFunc[i] * rShapeFunc[j]; // sigma * u * v + Stabilization (- sigma * v * TauOne * sigma * u)
//ZZ
                K *= Weight;
                // q-p stabilization block (reset result)
                L = 0;
 
                for (unsigned int m = 0; m < TDim; ++m) // iterate over v components (vx,vy[,vz])
                {
                    // Velocity block
                    //K += Weight * Density * Viscosity * rShapeDeriv(i, m) * rShapeDeriv(j, m); // Diffusive term: Viscosity * Grad(v) * Grad(u)
 
                    // v * Grad(p) block
 
                    G = TauOne * Density * AGradN[i] * rShapeDeriv(j, m); // Stabilization: (a * Grad(v)) * TauOne * Grad(p)
//GG
                    G -= TauOne * sigma * rShapeFunc[i] * rShapeDeriv(j, m); // Stabilization: (- sigma * v * TauOne * Grad(p)
//ZZ
//G
                    GEps =  TauOne * Density   * AGradN[i] * (FluidFraction * rShapeDeriv(j, m) + FluidFractionGradient[m] * rShapeFunc[j]); // Stabilization: (a * Grad(u)) * TauOne * (eps * Grad(q) - q * Grad(eps))
//GG
                    GEps += TauOne * sigma * rShapeFunc[i] * (FluidFraction * rShapeDeriv(j, m) + FluidFractionGradient[m] * rShapeFunc[j]); // Stabilization: sigma * u * TauOne * (eps * Grad(q) - q * Grad(eps))
//ZZ
                    DivPEpsilon = (FluidFraction * rShapeDeriv(i, m) + FluidFractionGradient[m] * rShapeFunc[i]) * rShapeFunc[j]; // eps * q * Div(v) + q * Grad(eps) * u
//Z
                    PDivV = rShapeDeriv(i, m) * rShapeFunc[j]; // Div(v) * p
 
                    // Write v * Grad(p) component
                    rDampingMatrix(FirstRow + m, FirstCol + TDim) += Weight * (G - PDivV);
                    // Use symmetry to write the q * Div(u) component
//G
    //              rDampingMatrix(FirstCol + TDim, FirstRow + m) += Weight * (G + PDivV);
                    rDampingMatrix(FirstCol + TDim, FirstRow + m) += Weight * (GEps + DivPEpsilon);
 
    //              q-p stabilization block
    //              L += rShapeDeriv(i, m) * rShapeDeriv(j, m); // Stabilization: Grad(q) * TauOne * Grad(p)
                    L += (FluidFraction * rShapeDeriv(i, m) + FluidFractionGradient[m] * rShapeFunc[i]) * rShapeDeriv(j, m); // Stabilization: (eps * Grad(q) + q * Grad(eps)) * TauOne * Grad(p)
//Z
                    for (unsigned int n = 0; n < TDim; ++n) // iterate over u components (ux,uy[,uz])
                    {
                        // Velocity block
//G
    //                  rDampingMatrix(FirstRow + m, FirstCol + n) += Weight * TauTwo * rShapeDeriv(i, m) * rShapeDeriv(j, n); // Stabilization: Div(v) * TauTwo * Div(u)
                        rDampingMatrix(FirstRow + m, FirstCol + n) += Weight * TauTwo * rShapeDeriv(i, m) * (FluidFraction * rShapeDeriv(j, n) + FluidFractionGradient[n] * rShapeFunc[j]); // Stabilization: Div(v) * TauTwo * (eps * Div(u) + Grad(eps) * u)
//Z
                    }
 
                }
 
                // Write remaining terms to velocity block
                for (unsigned int d = 0; d < TDim; ++d)
                    rDampingMatrix(FirstRow + d, FirstCol + d) += K;
 
                // Write q-p stabilization block
                rDampingMatrix(FirstRow + TDim, FirstCol + TDim) += Weight * TauOne * L;
 
 
                // Update reference column index for next iteration
                FirstCol += BlockSize;
            }
 
            // Operate on RHS
            qF = 0.0;
            for (unsigned int d = 0; d < TDim; ++d)
            {
//G
//GG
    //          rDampRHS[FirstRow + d] += Weight * TauOne * Density * AGradN[i] * BodyForce[d]; // ( a * Grad(v) ) * TauOne * (Density * BodyForce)
                rDampRHS[FirstRow + d] += Weight * TauOne * (Density * AGradN[i] - sigma * rShapeFunc[i]) * BodyForce[d];
//ZZ
    //          qF += rShapeDeriv(i, d) * BodyForce[d];
                qF += (FluidFraction * rShapeDeriv(i, d) + FluidFractionGradient[d] * rShapeFunc[i]) * BodyForce[d];
 
//Z
            }
 
            rDampRHS[FirstRow + TDim] += Weight * TauOne * qF; // Grad(q) * TauOne * (Density * BodyForce)
 
            // Update reference indices
            FirstRow += BlockSize;
            FirstCol = 0;
        }
 
//            this->AddBTransCB(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
 
        this->AddViscousTerm(rDampingMatrix,rShapeDeriv,Viscosity*Density*Weight);
    }
/*
//         void AddIntegrationPointVelocityContribution(MatrixType& rDampingMatrix,
//                                                      VectorType& rDampRHS,
//                                                      const double Density,
//                                                      const double Viscosity,
//                                                      const array_1d< double, 3 > & rAdvVel,
//                                                      const double TauOne,
//                                                      const double TauTwo,
//                                                      const array_1d< double, TNumNodes >& rShapeFunc,
//                                                      const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
//                                                      const double Weight,
//                                                      const double DeltaTime)
//         {
//             const unsigned int BlockSize = TDim + 1;
//
//             const double TauCoef = TauOne*TauTwo / DeltaTime;
//
//             // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
//             array_1d<double, TNumNodes> AGradN;
//             this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//
//             // Build the local matrix and RHS
//             unsigned int FirstRow(0), FirstCol(0); // position of the first term of the local matrix that corresponds to each node combination
//             double K, G, PDivV, L; // Temporary results
//             double Coef = Density * Weight;
//
//             // Note that we iterate first over columns, then over rows to read the Body Force only once per node
//             for (unsigned int j = 0; j < TNumNodes; ++j) // iterate over colums
//             {
//                 // Get Body Force
//                 const array_1d<double, 3 > & rBodyForce = this->GetGeometry()[j].FastGetSolutionStepValue(BODY_FORCE);
//
//                 const array_1d<double,3>& OldVelocity = this->GetGeometry()[j].FastGetSolutionStepValue(VELOCITY,1);
//
//                 for (unsigned int i = 0; i < TNumNodes; ++i) // iterate over rows
//                 {
//                     // Calculate the part of the contributions that is constant for each node combination
//
//                     // Velocity block
//                     K = rShapeFunc[i] * AGradN[j]; // Convective term: v * ( a * Grad(u) )
//                     K += TauOne * AGradN[i] * AGradN[j]; // Stabilization: (a * Grad(v)) * TauOne * (a * Grad(u))
//
//                     // q-p stabilization block (reset result)
//                     L = 0;
//
//                     for (unsigned int m = 0; m < TDim; ++m) // iterate over v components (vx,vy[,vz])
//                     {
//                         // Velocity block
// //                        K += Viscosity * rShapeDeriv(i, m) * rShapeDeriv(j, m); // Diffusive term: Viscosity * Grad(v) * Grad(u)
//                         // Note that we are usig kinematic viscosity, as we will multiply it by density later
//
//                         // v * Grad(p) block
//                         G = TauOne * AGradN[i] * rShapeDeriv(j, m); // Stabilization: (a * Grad(v)) * TauOne * Grad(p)
//                         PDivV = rShapeDeriv(i, m) * rShapeFunc[j]; // Div(v) * p
//
//                         // Write v * Grad(p) component
//                         rDampingMatrix(FirstRow + m, FirstCol + TDim) += Weight * (G - PDivV);
//                         // Use symmetry to write the q * rho * Div(u) component
//                         rDampingMatrix(FirstCol + TDim, FirstRow + m) += Coef * (G + PDivV);
//
//                         // q-p stabilization block
//                         L += rShapeDeriv(i, m) * rShapeDeriv(j, m); // Stabilization: Grad(q) * TauOne * Grad(p)
//
//                         for (unsigned int n = 0; n < TDim; ++n) // iterate over u components (ux,uy[,uz])
//                         {
//                             // Velocity block
//                             rDampingMatrix(FirstRow + m, FirstCol + n) += Coef * (TauTwo + TauCoef) * rShapeDeriv(i, m) * rShapeDeriv(j, n); // Stabilization: Div(v) * TauTwo *( 1+TauOne/Dt) * Div(u)
//                             rDampRHS[FirstRow + m] -= Coef * TauCoef * rShapeDeriv(i, m) * rShapeDeriv(j, n) * OldVelocity[n]; // Stabilization: Div(v) * TauTwo*TauOne/Dt * Div(u_old)
//                         }
//
//                     }
//
//                     // Write remaining terms to velocity block
//                     K *= Coef; // Weight by nodal area and density
//                     for (unsigned int d = 0; d < TDim; ++d)
//                         rDampingMatrix(FirstRow + d, FirstCol + d) += K;
//
//                     // Write q-p stabilization block
//                     rDampingMatrix(FirstRow + TDim, FirstCol + TDim) += Weight * TauOne * L;
//
//                     // Operate on RHS
//                     L = 0; // We reuse one of the temporary variables for the pressure RHS
//
//                     for (unsigned int d = 0; d < TDim; ++d)
//                     {
//                         rDampRHS[FirstRow + d] += Coef * TauOne * AGradN[i] * rShapeFunc[j] * rBodyForce[d]; // ( a * Grad(v) ) * TauOne * (Density * BodyForce)
//                         L += rShapeDeriv(i, d) * rShapeFunc[j] * rBodyForce[d];
//                     }
//                     rDampRHS[FirstRow + TDim] += Coef * TauOne * L; // Grad(q) * TauOne * (Density * BodyForce)
//
//                     // Update reference row index for next iteration
//                     FirstRow += BlockSize;
//                 }
//
//                 // Update reference indices
//                 FirstRow = 0;
//                 FirstCol += BlockSize;
//             }
//
// //            this->AddBTransCB(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
//             this->AddViscousTerm(rDampingMatrix,rShapeDeriv,Viscosity*Coef);
//         }
 */
 
    void AddProjectionResidualContribution(const array_1d< double, 3 > & rAdvVel,
                                           const double Density,
                                           array_1d< double, 3 > & rElementalMomRes,
                                           double& rElementalMassRes,
                                           const array_1d< double, TNumNodes >& rShapeFunc,
                                           const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                                           const double Weight)
    {
        // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
//G
        double FluidFraction;
        array_1d<double,3> FluidFractionGradient(3,0.0);
        this->EvaluateInPoint(FluidFraction, FLUID_FRACTION, rShapeFunc);
        this->EvaluateGradientOfScalarInPoint(FluidFractionGradient, FLUID_FRACTION, rShapeDeriv);
//Z
        // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
        for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
        {
 
            // Variable references
            const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
            const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
            const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
 
            // Compute this node's contribution to the residual (evaluated at inegration point)
            for (unsigned int d = 0; d < TDim; ++d)
            {
//GG
    //          rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
                rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - sigma * rShapeFunc[i] * rVelocity[d] - rShapeDeriv(i, d) * rPressure);
//ZZ
//G
    //          rElementalMassRes -= Weight * rShapeDeriv(i, d) * rVelocity[d];
                rElementalMassRes -= Weight * (FluidFraction * rShapeDeriv(i, d) * rVelocity[d] + FluidFractionGradient[d] * rVelocity[d]);
            }
 
            rElementalMassRes +=  Weight * this->GetGeometry()[i].FastGetSolutionStepValue(FLUID_FRACTION_RATE);
//Z
 
        }
    }
 
 
    void ASGSMomResidual(const array_1d< double, 3 > & rAdvVel,
                         const double Density,
                         array_1d< double, 3 > & rElementalMomRes,
                         const array_1d< double, TNumNodes >& rShapeFunc,
                         const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                         const double Weight)
    {
        // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
        // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
        for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
        {
 
            // Variable references
            const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
            const array_1d< double, 3 > & rAcceleration = this->GetGeometry()[i].FastGetSolutionStepValue(ACCELERATION);
            const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
            const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
 
            // Compute this node's contribution to the residual (evaluated at inegration point)
            for (unsigned int d = 0; d < TDim; ++d)
            {
//GG
    //          rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * (rBodyForce[d] - rAcceleration[d]) - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
                rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * (rBodyForce[d] - rAcceleration[d]) - AGradN[i] * rVelocity[d]) - sigma * rShapeFunc[i] * rVelocity[d] - rShapeDeriv(i, d) * rPressure);
//ZZ
            }
        }
    }
 
 
    void OSSMomResidual(const array_1d< double, 3 > & rAdvVel,
                        const double Density,
                        array_1d< double, 3 > & rElementalMomRes,
                        const array_1d< double, TNumNodes >& rShapeFunc,
                        const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                        const double Weight)
    {
        // If we want to use more than one Gauss point to integrate the convective term, this has to be evaluated once per integration point
        array_1d<double, TNumNodes> AGradN;
        this->GetConvectionOperator(AGradN, rAdvVel, rShapeDeriv); // Get a * grad(Ni)
//GG
        for (int i = 0; i != TNumNodes; ++i){
            AGradN[i] = 0.0;
          }
        double sigma = 0.0;
        this->EvaluateInPoint(sigma,PERMEABILITY_1_DAY,rShapeFunc);
//ZZ
        // Compute contribution to Kij * Uj, with Kij = Ni * Residual(Nj); Uj = (v,p)Node_j (column vector)
        for (unsigned int i = 0; i < TNumNodes; ++i) // Iterate over element nodes
        {
 
            // Variable references
            const array_1d< double, 3 > & rVelocity = this->GetGeometry()[i].FastGetSolutionStepValue(VELOCITY);
            const array_1d< double, 3 > & rBodyForce = this->GetGeometry()[i].FastGetSolutionStepValue(BODY_FORCE);
            const array_1d< double, 3 > & rProjection = this->GetGeometry()[i].FastGetSolutionStepValue(ADVPROJ);
            const double& rPressure = this->GetGeometry()[i].FastGetSolutionStepValue(PRESSURE);
 
            // Compute this node's contribution to the residual (evaluated at inegration point)
            for (unsigned int d = 0; d < TDim; ++d)
            {
//GG
    //          rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - rShapeDeriv(i, d) * rPressure);
                rElementalMomRes[d] += Weight * (Density * (rShapeFunc[i] * rBodyForce[d] - AGradN[i] * rVelocity[d]) - sigma * rShapeFunc[i] * rVelocity[d] - rShapeDeriv(i, d) * rPressure);
//ZZ
                rElementalMomRes[d] -= Weight * rShapeFunc[i] * rProjection[d];
 
            }
        }
    }
 
 
 
    virtual void GetEffectiveViscosity(const double Density,
                                       const double MolecularViscosity,
                                       const array_1d<double, TNumNodes>& rShapeFunc,
                                       const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                                       double& TotalViscosity,
                                       const ProcessInfo& rCurrentProcessInfo)
    {
        const double C = this->GetValue(C_SMAGORINSKY);
 
        TotalViscosity = MolecularViscosity;
        if (C != 0.0 )
        {
            // The filter width in Smagorinsky is typically the element size h. We will store the square of h, as the final formula involves the squared filter width
            const double FilterWidth = this->FilterWidth(rShapeDeriv);
 
            const double NormS = this->SymmetricGradientNorm(rShapeDeriv);
 
            // Total Viscosity
            TotalViscosity += 2.0 * C * C * FilterWidth * NormS;
        }
    }
 
 
    virtual void GetAdvectiveVel(array_1d< double, 3 > & rAdvVel,
                                 const array_1d< double, TNumNodes >& rShapeFunc)
    {
        // Compute the weighted value of the advective velocity in the (Gauss) Point
        rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY));
        for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
            rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY));
    }
 
 
    virtual void GetAdvectiveVel(array_1d< double, 3 > & rAdvVel,
                                 const array_1d< double, TNumNodes >& rShapeFunc,
                                 const std::size_t Step)
    {
        // Compute the weighted value of the advective velocity in the (Gauss) Point
        rAdvVel = rShapeFunc[0] * (this->GetGeometry()[0].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY, Step));
        for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
            rAdvVel += rShapeFunc[iNode] * (this->GetGeometry()[iNode].FastGetSolutionStepValue(VELOCITY, Step) - this->GetGeometry()[iNode].FastGetSolutionStepValue(MESH_VELOCITY, Step));
    }
 
 
    void GetConvectionOperator(array_1d< double, TNumNodes >& rResult,
                               const array_1d< double, 3 > & rVelocity,
                               const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
    {
        // Evaluate (and weight) the a * Grad(Ni) operator in the integration point, for each node i
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode) // Loop over nodes
        {
            // Initialize result
            rResult[iNode] = rVelocity[0] * rShapeDeriv(iNode, 0);
            for (unsigned int d = 1; d < TDim; ++d) // loop over components
                rResult[iNode] += rVelocity[d] * rShapeDeriv(iNode, d);
        }
    }
 
 
    virtual void AddPointContribution(double& rResult,
                                      const Variable< double >& rVariable,
                                      const array_1d< double, TNumNodes >& rShapeFunc,
                                      const double Weight = 1.0)
    {
        // Compute the weighted value of the nodal variable in the (Gauss) Point
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
            rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
    }
 
 
    virtual void EvaluateInPoint(double& rResult,
                                 const Variable< double >& rVariable,
                                 const array_1d< double, TNumNodes >& rShapeFunc)
    {
        // Compute the weighted value of the nodal variable in the (Gauss) Point
        rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
        for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
            rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
    }
 
//G
 
    virtual void EvaluateGradientOfScalarInPoint(array_1d< double, 3 >& rResult,
                                                 const Variable< double >& rVariable,
                                                 const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
    {
 
        for (unsigned int i = 0; i < TNumNodes; ++i) {
            double& scalar = this->GetGeometry()[i].FastGetSolutionStepValue(rVariable);
 
            for (unsigned int d = 0; d < TDim; ++d){
                rResult[d] += rShapeDeriv(i, d) * scalar;
            }
        }
    }
 
    virtual void EvaluateTimeDerivativeInPoint(double& rResult,
                                               const Variable< double >& rVariable,
                                               const array_1d< double, TNumNodes >& rShapeFunc,
                                               const double& DeltaTime,
                                               const std::vector<double>& rSchemeWeigths)
    {
        // Compute the time derivative of a nodal variable as a liner contribution of weighted value of the nodal variable in the (Gauss) Point
 
      if (rVariable == FLUID_FRACTION_RATE){
          rResult = 0.0;
 
          for (unsigned int iWeight = 0; iWeight < rSchemeWeigths.size(); ++iWeight){
 
              for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode){
                  rResult += rSchemeWeigths[iWeight] * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable, iWeight);
              }
 
          }
 
          rResult /= DeltaTime;
        }
    }
//Z
 
 
    virtual void AddPointContribution(array_1d< double, 3 > & rResult,
                                      const Variable< array_1d< double, 3 > >& rVariable,
                                      const array_1d< double, TNumNodes>& rShapeFunc,
                                      const double Weight = 1.0)
    {
        // Compute the weighted value of the nodal variable in the (Gauss) Point
        for (unsigned int iNode = 0; iNode < TNumNodes; ++iNode)
            rResult += Weight * rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
    }
 
 
    virtual void EvaluateInPoint(array_1d< double, 3 > & rResult,
                                 const Variable< array_1d< double, 3 > >& rVariable,
                                 const array_1d< double, TNumNodes >& rShapeFunc)
    {
        // Compute the weighted value of the nodal variable in the (Gauss) Point
        rResult = rShapeFunc[0] * this->GetGeometry()[0].FastGetSolutionStepValue(rVariable);
        for (unsigned int iNode = 1; iNode < TNumNodes; ++iNode)
            rResult += rShapeFunc[iNode] * this->GetGeometry()[iNode].FastGetSolutionStepValue(rVariable);
    }
 
 
    double ElementSize(const double);
 
    double FilterWidth();
 
    double FilterWidth(const BoundedMatrix<double, TNumNodes, TDim >& DN_DX);
 
 
    double SymmetricGradientNorm(const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv)
    {
        const unsigned int GradientSize = (TDim*(TDim+1))/2; // Number of different terms in the symmetric gradient matrix
        array_1d<double,GradientSize> GradientVector( GradientSize, 0.0 );
        unsigned int Index;
 
        // Compute Symmetric Grad(u). Note that only the lower half of the matrix is calculated
        for (unsigned int k = 0; k < TNumNodes; ++k)
        {
            const array_1d< double, 3 > & rNodeVel = this->GetGeometry()[k].FastGetSolutionStepValue(VELOCITY);
            Index = 0;
            for (unsigned int i = 0; i < TDim; ++i)
            {
                for (unsigned int j = 0; j < i; ++j) // Off-diagonal
                    GradientVector[Index++] += 0.5 * (rShapeDeriv(k, j) * rNodeVel[i] + rShapeDeriv(k, i) * rNodeVel[j]);
                GradientVector[Index++] += rShapeDeriv(k, i) * rNodeVel[i]; // Diagonal
            }
        }
 
        // Norm[ Symmetric Grad(u) ] = ( 2 * Sij * Sij )^(1/2)
        Index = 0;
        double NormS(0.0);
        for (unsigned int i = 0; i < TDim; ++i)
        {
            for (unsigned int j = 0; j < i; ++j)
            {
                NormS += 2.0 * GradientVector[Index] * GradientVector[Index]; // Using symmetry, lower half terms of the matrix are added twice
                ++Index;
            }
            NormS += GradientVector[Index] * GradientVector[Index]; // Diagonal terms
            ++Index; // Diagonal terms
        }
 
        NormS = sqrt( 2.0 * NormS );
        return NormS;
    }
 
 
    virtual void AddViscousTerm(MatrixType& rDampingMatrix,
                                const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                                const double Weight);
 
 
    void AddBTransCB(MatrixType& rDampingMatrix,
                     const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv,
                     const double Weight)
    {
        BoundedMatrix<double, (TDim * TNumNodes)/2, TDim*TNumNodes > B;
        BoundedMatrix<double, (TDim * TNumNodes)/2, (TDim*TNumNodes)/2 > C;
        this->CalculateB(B,rShapeDeriv);
        this->CalculateC(C,Weight);
 
        const unsigned int BlockSize = TDim + 1;
        const unsigned int StrainSize = (TDim*TNumNodes)/2;
 
        DenseVector<unsigned int> aux(TDim*TNumNodes);
        for(unsigned int i=0; i<TNumNodes; i++)
        {
            int base_index = TDim*i;
            int aux_index = BlockSize*i;
            for(unsigned int j=0; j<TDim; j++)
            {
                aux[base_index+j] = aux_index+j;
            }
        }
 
        for(unsigned int k=0; k< StrainSize; k++)
        {
            for(unsigned int l=0; l< StrainSize; l++)
            {
                const double Ckl = C(k,l);
                for (unsigned int i = 0; i < TDim*TNumNodes; ++i) // iterate over v components (vx,vy[,vz])
                {
                    const double Bki=B(k,i);
                    for (unsigned int j = 0; j < TDim*TNumNodes; ++j) // iterate over u components (ux,uy[,uz])
                    {
                        rDampingMatrix(aux[i],aux[j]) += Bki*Ckl*B(l,j);
                    }
 
                }
 
            }
        }
    }
 
    void ModulatedGradientDiffusion(MatrixType& rDampingMatrix,
            const BoundedMatrix<double, TNumNodes, TDim >& rDN_DX,
            const double Weight)
    {
        const GeometryType& rGeom = this->GetGeometry();
 
        // Velocity gradient
        MatrixType GradU = ZeroMatrix(TDim,TDim);
        for (unsigned int n = 0; n < TNumNodes; n++)
        {
            const array_1d<double,3>& rVel = this->GetGeometry()[n].FastGetSolutionStepValue(VELOCITY);
            for (unsigned int i = 0; i < TDim; i++)
                for (unsigned int j = 0; j < TDim; j++)
                    GradU(i,j) += rDN_DX(n,j)*rVel[i];
        }
 
        // Element lengths
        array_1d<double,3> Delta(3,0.0);
        Delta[0] = fabs(rGeom[TNumNodes-1].X()-rGeom[0].X());
        Delta[1] = fabs(rGeom[TNumNodes-1].Y()-rGeom[0].Y());
        Delta[2] = fabs(rGeom[TNumNodes-1].Z()-rGeom[0].Z());
 
        for (unsigned int n = 1; n < TNumNodes; n++)
        {
            double hx = fabs(rGeom[n].X()-rGeom[n-1].X());
            if (hx > Delta[0]) Delta[0] = hx;
            double hy = fabs(rGeom[n].Y()-rGeom[n-1].Y());
            if (hy > Delta[1]) Delta[1] = hy;
            double hz = fabs(rGeom[n].Z()-rGeom[n-1].Z());
            if (hz > Delta[2]) Delta[2] = hz;
        }
 
        double AvgDeltaSq = Delta[0];
        for (unsigned int d = 1; d < TDim; d++)
            AvgDeltaSq *= Delta[d];
        AvgDeltaSq = std::pow(AvgDeltaSq,2./TDim);
 
        Delta[0] = Delta[0]*Delta[0]/12.0;
        Delta[1] = Delta[1]*Delta[1]/12.0;
        Delta[2] = Delta[2]*Delta[2]/12.0;
 
        // Gij
        MatrixType G = ZeroMatrix(TDim,TDim);
        for (unsigned int i = 0; i < TDim; i++)
            for (unsigned int j = 0; j < TDim; j++)
                for (unsigned int d = 0; d < TDim; d++)
                    G(i,j) += Delta[d]*GradU(i,d)*GradU(j,d);
 
        // Gij:Sij
        double GijSij = 0.0;
        for (unsigned int i = 0; i < TDim; i++)
            for (unsigned int j = 0; j < TDim; j++)
                GijSij += 0.5*G(i,j)*( GradU(i,j) + GradU(j,i) );
 
        if (GijSij < 0.0) // Otherwise model term is clipped
        {
            // Gkk
            double Gkk = G(0,0);
            for (unsigned int d = 1; d < TDim; d++)
                Gkk += G(d,d);
 
            // C_epsilon
            const double Ce = 1.0;
 
            // ksgs
            double ksgs = -4*AvgDeltaSq*GijSij/(Ce*Ce*Gkk);
 
            // Assembly of model term
            unsigned int RowIndex = 0;
            unsigned int ColIndex = 0;
 
            for (unsigned int i = 0; i < TNumNodes; i++)
            {
                for (unsigned int j = 0; j < TNumNodes; j++)
                {
                    for (unsigned int d = 0; d < TDim; d++)
                    {
                        double Aux = rDN_DX(i,d) * Delta[0] * G(d,0)*rDN_DX(j,0);
                        for (unsigned int k = 1; k < TDim; k++)
                            Aux += rDN_DX(i,d) *Delta[k] * G(d,k)*rDN_DX(j,k);
                        rDampingMatrix(RowIndex+d,ColIndex+d) += Weight * 2.0*ksgs *  Aux;
                    }
 
                    ColIndex += TDim;
                }
                RowIndex += TDim;
                ColIndex = 0;
            }
        }
 
    }
 
 
    void CalculateB( BoundedMatrix<double, (TDim * TNumNodes) / 2, TDim * TNumNodes >& rB,
                     const BoundedMatrix<double, TNumNodes, TDim >& rShapeDeriv);
 
 
    virtual void CalculateC( BoundedMatrix<double, (TDim * TNumNodes) / 2, (TDim * TNumNodes) / 2 >& rC,
                             const double Viscosity);
 
    double ConsistentMassCoef(const double Area);
 
 
 
 
 
 
 
 
private:
 
 
 
 
    friend class Serializer;
 
    virtual void save(Serializer& rSerializer) const override
    {
        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Element );
    }
 
    virtual void load(Serializer& rSerializer) override
    {
        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element);
    }
 
 
 
 
 
 
 
 
 
 
    MonolithicDEMCoupledWeak & operator=(MonolithicDEMCoupledWeak const& rOther);
 
    MonolithicDEMCoupledWeak(MonolithicDEMCoupledWeak const& rOther);
 
 
}; // Class MonolithicDEMCoupledWeak
 
 
 
 
 
 
template< unsigned int TDim,
          unsigned int TNumNodes >
inline std::istream& operator >>(std::istream& rIStream,
                                 MonolithicDEMCoupledWeak<TDim, TNumNodes>& rThis)
{
    return rIStream;
}
 
template< unsigned int TDim,
          unsigned int TNumNodes >
inline std::ostream& operator <<(std::ostream& rOStream,
                                 const MonolithicDEMCoupledWeak<TDim, TNumNodes>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
} // namespace Kratos.
 
#endif // MONOLITHIC_DEM_COUPLED_WEAK_H