kutta_condition_process.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
 
#ifndef KRATOS_KUTTA_CONDITION_PROCESS_H
#define KRATOS_KUTTA_CONDITION_PROCESS_H
 
 
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/kratos_flags.h"
#include "processes/process.h"
#include "geometries/geometry.h"
#include "utilities/geometry_utilities.h"
#include "utilities/math_utils.h"
#include "includes/kratos_parameters.h"
 
#include <string>
#include <iostream>
#include <sstream>
 
#include <boost/functional/hash.hpp> //TODO: remove this dependence when Kratos has en internal one
#include <unordered_map> //TODO: remove this dependence when Kratos has en internal one
#include <utility>
 
namespace Kratos
{
 
class KuttaConditionProcess: public Process
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(KuttaConditionProcess);
 
    typedef ModelPart::ElementType ElementType;
    typedef ModelPart::ConditionType ConditionType;
 
 
//     KuttaConditionProcess(ModelPart& rModelPart,
//                      KratosParameters& parameters
//                     ):
//         Process(),
//         mrModelPart(rModelPart),
//         mrOptions(Flags()),
//         mrParameters(parameters)
//     {
//     }
    KuttaConditionProcess(ModelPart& rModelPart
                    ):
        Process(),
        mrModelPart(rModelPart)
    {
    }
 
    ~KuttaConditionProcess() override {}
 
 
 
    void operator()()
    {
        Execute();
    }
 
 
 
 
    void Execute() override
    {
        KRATOS_TRY;
 
        if(mrModelPart.Elements().size() == 0)
            KRATOS_THROW_ERROR(std::invalid_argument, "the number of elements in the domain is zero. kutta condition can not be applied","")
 
        using normal_map_type = std::unordered_map<std::size_t, array_1d<double,3> >;
        
        for(auto it=mrModelPart.NodesBegin(); it!=mrModelPart.NodesEnd(); ++it)
            if(it->Is(STRUCTURE)) it->Set(VISITED,false);
 
        PointerVector<Element> kutta_elements;
        for(auto it=mrModelPart.ElementsBegin(); it!=mrModelPart.ElementsEnd(); ++it)
        {
            auto geom = it->GetGeometry();
            
            for(unsigned int i=0; i<geom.size(); ++i)
            {
                if(geom[i].Is(STRUCTURE))
                {
                    kutta_elements.push_back(*(it.base()));
                    break;
                }
            }
        }
        
        //compute one normal per each kutta node
        normal_map_type normal_map;
        for(auto it=kutta_elements.begin(); it!=kutta_elements.end(); ++it)
        {
            if(it->Is(MARKER))
            {
                //compute normal
                const Vector& elemental_distances = it->GetValue(ELEMENTAL_DISTANCES);
                
                double vol;
                array_1d<double,4> N;
                bounded_matrix<double,4,3> DN_DX;
                GeometryUtils::CalculateGeometryData(it->GetGeometry(), DN_DX, N, vol);
                
                //get a relevant norm
                array_1d<double,3> n = prod(trans(DN_DX),elemental_distances);
                n/=(norm_2(n)+1e-30);
 
                //
                auto geom = it->GetGeometry();
                for(unsigned int i=0; i<geom.size(); ++i)
                {
                    if(geom[i].Is(STRUCTURE))
                    {
                        geom[i].Set(VISITED,true);
                        normal_map[geom[i].Id()] = n;;
                    }
                }
            }
        }
        
        //now compute one elemntal distance per each element in the kutta list
        for(auto it=kutta_elements.begin(); it!=kutta_elements.end(); ++it)
        {
            it->Set(MARKER,false);
            Vector& elemental_distances = it->GetValue(ELEMENTAL_DISTANCES);
            
            array_1d<double,3> n;
            auto geom = it->GetGeometry();
            int kutta_id = -1;
            for(unsigned int i=0; i<geom.size(); ++i)
            {
                if(geom[i].Is(STRUCTURE) && geom[i].Is(VISITED))
                {
                    kutta_id = i;
                    n = normal_map[geom[i].Id()];
                    break;
                }
            }
            
            if(kutta_id == -1) //no node was found with the normal correctly calculated
            {
                it->Set(ACTIVE,false);
            }
            else
            {
                   
                
                array_1d<double,3> x0 = geom[kutta_id].Coordinates();
                double dkutta = -1e-3;
                for(int i=0; i<static_cast<int>(geom.size()); ++i)
                {
                    if(kutta_id != i)
                    {
                        
                        array_1d<double,3> v = geom[i].Coordinates() - x0;
                        elemental_distances[i] = inner_prod(n,v);
                    }
                    else
                        elemental_distances[i] = dkutta;
                }
                
                
                unsigned int npos = 0;
                unsigned int nneg = 0;
                for(unsigned int i=0; i<geom.size(); ++i)
                    if(elemental_distances[i] >= 0)
                        npos++;
                    else
                        nneg++;
                    
                if(nneg > 0 && npos >0)
                    it->Set(MARKER,true);
            }
            
        }       
        
        KRATOS_CATCH("");
    }
 
 
 
 
 
 
 
 
    std::string Info() const override
    {
        return "KuttaConditionProcess";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "KuttaConditionProcess";
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        this->PrintInfo(rOStream);
    }
 
 
 
 
 
 
private:
 
 
 
    ModelPart& mrModelPart;
    Flags mrOptions;
 
 
 
 
 
    KuttaConditionProcess& operator=(KuttaConditionProcess const& rOther);
 
    KuttaConditionProcess(KuttaConditionProcess const& rOther);
 
 
 
}; // Class Process
 
 
 
 
 
 
inline std::istream& operator >> (std::istream& rIStream,
                                  KuttaConditionProcess& rThis);
 
inline std::ostream& operator << (std::ostream& rOStream,
                                  const KuttaConditionProcess& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
 
} // namespace Kratos
 
 
#endif // KRATOS_KUTTA_CONDITION_PROCESS_H