distance_calculation_element_simplex.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
//
 
#if !defined(KRATOS_DISTANCE_CALCULATION_ELEMENT_H_INCLUDED )
#define  KRATOS_DISTANCE_CALCULATION_ELEMENT_H_INCLUDED
 
// System includes
#include <string>
#include <iostream>
 
 
// External includes
 
 
// Project includes
#include "containers/array_1d.h"
#include "includes/define.h"
#include "includes/element.h"
#include "includes/serializer.h"
#include "geometries/geometry.h"
#include "utilities/math_utils.h"
#include "utilities/geometry_utilities.h"
 
// Application includes
#include "includes/variables.h"
#include "includes/kratos_flags.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
template< unsigned int TDim >
class DistanceCalculationElementSimplex : public Element
{
public:
 
    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(DistanceCalculationElementSimplex);
 
    typedef Node NodeType;
 
    typedef Geometry<NodeType> GeometryType;
 
    typedef Geometry<NodeType>::PointsArrayType NodesArrayType;
 
    typedef Vector VectorType;
 
    typedef Matrix MatrixType;
 
    typedef std::size_t IndexType;
 
    typedef std::size_t SizeType;
 
    typedef Dof<double> DofType;
 
    typedef std::vector<std::size_t> EquationIdVectorType;
 
    typedef std::vector<DofType::Pointer> DofsVectorType;
 
    typedef PointerVectorSet<DofType> DofsArrayType;
 
    typedef Kratos::Vector ShapeFunctionsType;
 
    typedef Kratos::Matrix ShapeFunctionDerivativesType;
 
    typedef GeometryType::ShapeFunctionsGradientsType ShapeFunctionDerivativesArrayType;
 
 
    //Constructors.
 
 
    DistanceCalculationElementSimplex(IndexType NewId = 0) :
        Element(NewId)
    {}
 
 
    DistanceCalculationElementSimplex(IndexType NewId, const NodesArrayType& ThisNodes) :
        Element(NewId, ThisNodes)
    {}
 
 
    DistanceCalculationElementSimplex(IndexType NewId, GeometryType::Pointer pGeometry) :
        Element(NewId, pGeometry)
    {}
 
 
    DistanceCalculationElementSimplex(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) :
        Element(NewId, pGeometry, pProperties)
    {}
 
    ~DistanceCalculationElementSimplex() override
    {}
 
 
 
 
 
 
    Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,
                            PropertiesType::Pointer pProperties) const override
    {
        return Element::Pointer(new DistanceCalculationElementSimplex(NewId, GetGeometry().Create(ThisNodes), pProperties));
    }
 
 
    Element::Pointer Create(IndexType NewId,
                            GeometryType::Pointer pGeom,
                            PropertiesType::Pointer pProperties) const override
    {
        return Element::Pointer(new DistanceCalculationElementSimplex(NewId, pGeom, pProperties));
    }
 
    void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix,
                              VectorType& rRightHandSideVector,
                              const ProcessInfo& rCurrentProcessInfo) override
    {
        const unsigned int number_of_points = TDim+1;
 
        BoundedMatrix<double, TDim+1, TDim > DN_DX;
        array_1d<double, TDim+1 > N;
 
        if (rLeftHandSideMatrix.size1() != number_of_points)
            rLeftHandSideMatrix.resize(number_of_points, number_of_points, false);
 
        if (rRightHandSideVector.size() != number_of_points)
            rRightHandSideVector.resize(number_of_points, false);
 
        //getting data for the given geometry
        double Area;
        GeometryUtils::CalculateGeometryData(GetGeometry(), DN_DX, N, Area);
 
        //get distances at the nodes
        array_1d<double, TDim+1 > distances;
        for(unsigned int i=0; i<number_of_points; i++)
        {
            distances[i] = GetGeometry()[i].FastGetSolutionStepValue(DISTANCE);
//          if (distances[i] >= 0.0 && distances[i] < 1e-30) distances[i] = 1e-30;
        }
 
       const double dgauss = inner_prod(N,distances);
 
        const double coeff1 = rCurrentProcessInfo.Has(VARIATIONAL_REDISTANCE_COEFFICIENT_FIRST) ? rCurrentProcessInfo[VARIATIONAL_REDISTANCE_COEFFICIENT_FIRST] : 0.01;
        const double coeff2 = rCurrentProcessInfo.Has(VARIATIONAL_REDISTANCE_COEFFICIENT_SECOND) ? rCurrentProcessInfo[VARIATIONAL_REDISTANCE_COEFFICIENT_SECOND] : 0.1;
        const unsigned int step = rCurrentProcessInfo[FRACTIONAL_STEP];
 
        if(step == 1) //solve a poisson problem with a positive/negative heat source depending on the sign of the existing distance function
        {
            //compute distance on gauss point
 
            this->SetValue(DISTANCE,dgauss); //saving the distance, to see if it changed sign between iterations
            //compute LHS
            noalias(rLeftHandSideMatrix) = Area*prod(DN_DX,trans(DN_DX));
 
            //compute RHS
            double source = 1.0;
            if(dgauss < 0.0) source=-1.0;
            noalias(rRightHandSideVector) = source*Area*N;
 
            noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,distances);
 
 
            //impose that the normal gradient is 1 on outer faces
            unsigned int nboundary = 0;
            for(unsigned int i=0; i<TDim+1; i++)
                if(GetGeometry()[i].Is(BOUNDARY)) nboundary +=1;
 
            if(nboundary == TDim)
            {
          array_1d<double,TDim> DN_out(TDim, 0.0);
                for(unsigned int i=0; i<TDim+1; i++)
                    if(GetGeometry()[i].IsNot(BOUNDARY))
                    {
                        noalias(DN_out) = row(DN_DX,i);
                        break;
                    }
 
                double normDn = norm_2(DN_out);
                for(unsigned int i=0; i<TDim+1; i++)
                    if(GetGeometry()[i].Is(BOUNDARY))
                    {
                        rRightHandSideVector[i] += coeff1*source*normDn*Area*(TDim-1); //TODO: check this! it should be TDim*(TDim-1)*N[i] with N[i] on the face and then equal to 1/TDim
                        // using the area as weighting factor is excessive. Reduced it to get a closer to constant gradient between the regions close and far away from the interface
                    }
 
            }
 
 
        }
        else //solve an optimization problem with the goal of achievieng a gradient of one for the distance function
        {
            //for debuggin purposes:
            if( dgauss * (this->GetValue(DISTANCE)) < 0.0 )
                std::cout << "Element " << this->Id() << " changed sign while redistancing!!" << std::endl;
 
            //compute the gradient of the distance
            const array_1d<double,TDim> grad = prod(trans(DN_DX),distances);
            double grad_norm = norm_2(grad);
 
            //compute RHS ad grad N_i \cdot ( 1/norm_grad * grad - grad)
            //and multiply everything by grad_norm
            noalias(rRightHandSideVector) = Area*(1.0 - grad_norm)* prod(DN_DX,grad);
 
 
            //compute the LHS as an approximation of the tangent.
            //such approximation is taken as a laplacian, which comes from the assumption that the
            //direction of n does not change when d changes
            //
            //
            //note that the exact tangent could be computed as "P1+P2" with
            //n = grad/grad_norm
            //P1 = (1.0 - 1.0 / (grad_norm + eps) )    * DN_DX * DN_DX.transpose()
            //P2 = 1.0/(grad_norm + eps) * dot(DN_DX * outer(n,n) * DN_DX.transpose() )
            //unfortunately the numerical experiments tell that this in too unstable to be used unless a very
            //good initial approximation is used
//            noalias(rLeftHandSideMatrix) = (Area*(grad_norm - 1.0))*rod(DN_DX,trans(DN_DX) ); //RISKY!!
            noalias(rLeftHandSideMatrix) = Area*std::max(grad_norm,coeff2)*prod( DN_DX,trans(DN_DX) );
        }
    }
 
 
 
    void EquationIdVector(EquationIdVectorType& rResult,
                          const ProcessInfo& rCurrentProcessInfo) const override
    {
 
        unsigned int number_of_nodes = TDim+1;
        if (rResult.size() != number_of_nodes)
            rResult.resize(number_of_nodes, false);
 
        for (unsigned int i = 0; i < number_of_nodes; i++)
            rResult[i] = GetGeometry()[i].GetDof(DISTANCE).EquationId();
    }
 
 
    void GetDofList(DofsVectorType& rElementalDofList,
                    const ProcessInfo& rCurrentProcessInfo) const override
    {
        unsigned int number_of_nodes = TDim+1;
 
        if (rElementalDofList.size() != number_of_nodes)
            rElementalDofList.resize(number_of_nodes);
 
        for (unsigned int i = 0; i < number_of_nodes; i++)
            rElementalDofList[i] = GetGeometry()[i].pGetDof(DISTANCE);
 
    }
 
 
 
 
 
    int Check(const ProcessInfo& rCurrentProcessInfo) const override
{
        KRATOS_TRY
 
        // Perform basic element checks
        int ErrorCode = Kratos::Element::Check(rCurrentProcessInfo);
        if(ErrorCode != 0) return ErrorCode;
 
        if(this->GetGeometry().size() != TDim+1)
            KRATOS_THROW_ERROR(std::invalid_argument,"wrong number of nodes for element",this->Id());
 
        // Checks on nodes
 
        // Check that the element's nodes contain all required SolutionStepData and Degrees of freedom
        for(unsigned int i=0; i<this->GetGeometry().size(); ++i)
        {
            if(this->GetGeometry()[i].SolutionStepsDataHas(DISTANCE) == false)
                KRATOS_THROW_ERROR(std::invalid_argument,"missing DISTANCE variable on solution step data for node ",this->GetGeometry()[i].Id());
        }
 
 
 
        return 0;
 
        KRATOS_CATCH("");
    }
 
    void CalculateRightHandSide(VectorType& rRightHandSideVector,
                                const ProcessInfo& rCurrentProcessInfo) override
    {
        KRATOS_ERROR << "should not enter here" << std::endl;
        if (rRightHandSideVector.size() != 0)
      rRightHandSideVector.resize(0, false);
    }
 
 
 
    const Parameters GetSpecifications() const override
    {
        const Parameters specifications = Parameters(R"({
            "time_integration"           : ["static"],
            "framework"                  : "eulerian",
            "symmetric_lhs"              : true,
            "positive_definite_lhs"      : true,
            "output"                     : {
                "gauss_point"            : [],
                "nodal_historical"       : ["DISTANCE"],
                "nodal_non_historical"   : [],
                "entity"                 : []
            },
            "required_variables"         : ["DISTANCE"],
            "required_dofs"              : ["DISTANCE"],
            "flags_used"                 : ["BOUNDARY"],
            "compatible_geometries"      : ["Triangle2D3","Tetrahedra3D4"],
            "element_integrates_in_time" : false,
            "compatible_constitutive_laws": {
                "type"        : [],
                "dimension"   : [],
                "strain_size" : []
            },
            "required_polynomial_degree_of_geometry" : 1,
            "documentation"   :
                "This element is intended to be used in combination with the VariationalDistanceCalculationProcess. It implements a two-step resolution of an Eikonal equation in order to obtain a distance field with unit gradient norm."
        })");
        return specifications;
    }
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "DistanceCalculationElementSimplex #" << Id();
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "DistanceCalculationElementSimplex" << TDim << "D";
    }
 
//        /// Print object's data.
//        virtual void PrintData(std::ostream& rOStream) const;
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
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);
    }
 
 
 
 
 
 
 
 
 
 
    DistanceCalculationElementSimplex & operator=(DistanceCalculationElementSimplex const& rOther);
 
    DistanceCalculationElementSimplex(DistanceCalculationElementSimplex const& rOther);
 
 
}; // Class DistanceCalculationElementSimplex
 
 
 
 
 
 
template< unsigned int TDim >
inline std::istream& operator >>(std::istream& rIStream,
                                 DistanceCalculationElementSimplex<TDim>& rThis)
{
    return rIStream;
}
 
template< unsigned int TDim >
inline std::ostream& operator <<(std::ostream& rOStream,
                                 const DistanceCalculationElementSimplex<TDim>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
} // namespace Kratos.
 
#endif // KRATOS_DISTANCE_CALCULATION_ELEMENT_H_INCLUDED  defined