mesh_suite_modeler.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics 
//
//  License:         BSD License 
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Pooyan Dadvand
//                    
//
 
 
#if !defined(KRATOS_MESH_SUITE_MODELER_H_INCLUDED )
#define  KRATOS_MESH_SUITE_MODELER_H_INCLUDED
 
 
 
// System includes
#include <string>
#include <iostream>
#include <cstdlib>
#include <iomanip>
#include <fstream>
 
 
// External includes
#include "malla.h"
#include "voronoi.h"
 
 
// Project includes
#include "includes/define.h"
#include "includes/model_part.h"
//#include "geometries/triangle_2d.h"
#include "punto.h"
#include "nodo.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
class MeshSuiteModeler  : public malla
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(MeshSuiteModeler);
 
 
    MeshSuiteModeler() {}
 
    virtual ~MeshSuiteModeler() {}
 
 
 
 
 
    //**************************************************************************************
    //**************************************************************************************
 
    void GenerateMesh(ModelPart& ThisModelPart, Element const& rReferenceElement)
    {
        KRATOS_TRY
 
        mk_puntos(true);
        //delaunay(true,false);
 
        //erasing the elements
        ThisModelPart.Nodes().clear();
        ThisModelPart.Elements().clear();
 
        //generating in Kratos the new nodes
 
        ThisModelPart.Nodes().reserve(n.len);
        for(int i = 0 ; i < n.len ; i++)
        {
            Node::Pointer p_node = Kratos::Node::Pointer(new Kratos::Node(0,n[i][0], n[i][1], n[i][2]));
            p_node->SetId(i+1);
            ThisModelPart.Nodes().push_back(p_node);
        }
 
        Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
 
 
        ThisModelPart.Elements().reserve(e.len);
        for(int i_element = 0 ; i_element < e.len ; i_element++)
        {
            Element::NodesArrayType temp;
 
            for(int i_node = 0 ; i_node < e[i_element].nv() ; i_node++)
            {
                temp.push_back(ThisModelPart.Nodes()(e[i_element][i_node] + 1));
            }
            Element::Pointer p_element = rReferenceElement.Create(i_element + 1, temp, properties);
            ThisModelPart.Elements().push_back(p_element);
        }
        KRATOS_CATCH("")
    }
 
    //**************************************************************************************
    //**************************************************************************************
    void FixNodes(ModelPart& ThisModelPart)
    {
 
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
        {
            n[i_node->Id()-1].f.set(n_permanente);
        }
    }
 
    //**************************************************************************************
    //**************************************************************************************
    int UpdateMesh(ModelPart& ThisModelPart)
    {
        // lee los nodos
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
        {
            std::size_t id = i_node->Id() - 1;
            n[id][0] = i_node->X();
            n[id][1] = i_node->Y();
            n[id][2] = i_node->Z();
        }
 
 
        return true;
    }
 
 
    //**************************************************************************************
    //**************************************************************************************
    int SetMesh(ModelPart& ThisModelPart)
    {
        KRATOS_TRY
        ini();
 
        // lee de acuerdo a filext
        array1<nodo> nl;
        array1<elemento> el;
 
        nl.clean();
        el.clean();
 
        // lee los nodos
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
            nl+=nodo(i_node->X(), i_node->Y() , i_node->Z());
 
 
 
        // lee los elementos
        for(ModelPart::ElementIterator i_element = ThisModelPart.ElementsBegin() ; i_element != ThisModelPart.ElementsEnd() ; ++i_element)
        {
            if(i_element->GetGeometry().PointsNumber() == 2)
            {
                elemento temp(e_segmento);
                temp[0] = i_element->GetGeometry()[0].Id() - 1;
                temp[1] = i_element->GetGeometry()[1].Id() - 1;
                //      KRATOS_WATCH(temp[0])
                //      KRATOS_WATCH(temp[1])
                el+=temp;
            }
            else if(i_element->GetGeometry().PointsNumber() == 3)
            {
                elemento temp(e_triangulo);
                temp[0] = i_element->GetGeometry()[0].Id() - 1;
                temp[1] = i_element->GetGeometry()[1].Id() - 1;
                temp[2] = i_element->GetGeometry()[2].Id() - 1;
                el+=temp;
            }
        }
        //   // nombre y extensión
        //   if (!inombre[0]) {add_error(No_Mesh); return false;}
        //   renombra(inombre);
 
 
        if (!nl)
        {
            add_error(No_Mesh);
            return false;
        }
 
        int i,j,in,nc=0;
        int &nlen=nl.len;
        n.resize(nlen);
 
        // bounding box
        pmin=pmax=nl[0];
        for (i=1; i<nlen; i++)
        {
            if (nl[i].f.es(n_borrado)) continue;
            pmin.set_min(nl[i]);
            pmax.set_max(nl[i]);
        }
 
        // hay z?
        bool hayz=pmax[2]-pmin[2]>ERRADM;
        if (!hayz) tipo.set(m_planaxy);
        else tipo.reset(m_planaxy);
        if(!o)
            o=new octree(n,pmin,pmax,hayz ? 3 : 2);
        bool puesto,hayrep=false;
        double d,dmax=0;
        pline map(nlen);
        map.len=nlen;
        // los nodos frontera van primero y los
        // de h (si hay) despues
        bool haynh=false;
        for (i=0; i<nlen; i++)
        {
            if (nl[i].f.es(n_borrado)) continue;
            if (nl[i].f.es(n_h))
            {
                haynh=true;
                continue;
            }
            in=map[i]=n+=nl[i]; // presupongo que es n.len-1
            o->add_no_rep(in,nc,puesto,epsilon);
            if (puesto) continue;
            if (nc<0||nc>=nlen)
            {
                add_error(Bad_Format);
                return false;
            }
            hayrep=true;
            d=n[nc].distancia(nl[i]);
            if (d>dmax) dmax=d;
            n.len--;
            map[i]=nc;// si no es n.len-1 usar remove
            n[nc].f.set(nl[i].f); // oreo flags
            if (n[nc].h>ERRADM&&nl[i].h>ERRADM)
                set_min(n[nc].h,nl[i].h); // h es el menor
            else
                set_max(n[nc].h,nl[i].h); // h > 0
            n[nc].v=(n[nc].v+nl[i].v)/2;
        }
        // los nodos de h
        if (haynh) for (i=0; i<nlen; i++)
            {
                if (nl[i].f.es(n_borrado)) continue;
                if (nl[i].f.noes(n_h)) continue;
                in=map[i]=n+=nl[i]; // presupongo que es n.len-1
                //    o->add_no_rep(in,nc,puesto,epsilon);
                if (puesto)
                {
                    nodosh++;
                    if (n[in].h<ERRADM) n[in].h=MAXREAL;
                    continue;
                }
                n.len--;
                map[i]=nc;// si no es n.len-1 usar remove
                // no oreo flags
                if (n[nc].h>ERRADM&&nl[i].h>ERRADM)
                    set_min(n[nc].h,nl[i].h); // h es el menor
                else
                    set_max(n[nc].h,nl[i].h); // h > 0
                if (n[nc].h<ERRADM) n[nc].h=MAXREAL;
                n[nc].v=(n[nc].v+nl[i].v)/2;
            }
 
        // modifica los elementos y hace eldenod
        if (el)
        {
            tipo.reset(m_nodos);
            int nv,dim,ie;
            flagtype fdim;
            e.resize(el);
            for (i=0; i<el.len; i++)
            {
                elemento &ei=el[i];
                nv=ei.nv();
                dim=ei.dim();
                // verifica nodos h
                for (haynh=false,j=0; j<nv; j++)
                {
                    if (n[ei[j]=map[ei[j]]].f.es(n_h)) haynh=true;
                }
                // verifica si el elm esta repetido
                cpline &en=n[ei[0]].e;
                for (j=0; j<en.len; j++) if (ei==e[en[j]]) break;
                if (j<en.len)   // repetido
                {
                    add_warning(Repeated_Elements);
                    continue;
                }
                // verifica si el elm tiene nodos repetidos
                for (j=0; j<nv; j++)
                    if (ei[j]==ei.npos(j))
                        break;
                if (j<nv) continue; // no lo agrega
                // nodos h
                if (haynh)  // no agrega el elemento, solo saca el h de cada nodo
                {
                    for (j=0; j<nv; j++)
                    {
                        nodo &ni=n[ei[j]];
                        if (ni.f.noes(n_h)) continue;
                        if (ni.h!=MAXREAL) continue; // ya tiene h
                        d=ni.distancia(n[ei.nant(j)]);
                        if (d>ERRADM) ni.h=d;
                        d=ni.distancia(n[ei.npos(j)]);
                        if (d>ERRADM&&d<ni.h) ni.h=d;
                    }
                    continue;
                }
                // agrega el elemento
                ie=e+=ei;
                for (j=0; j<nv; j++) n[ei[j]].e+=ie;
                eltipos.set(ei.ftipo());
                if (dim==0) fdim.set(m_nodos);
                else if (dim==1) fdim.set(m_lin);
                else if (dim==2) fdim.set(m_sup);
                else if (dim==3) fdim.set(m_vol);
            }
            if (fdim==m_nodos) tipo.set(m_nodos);
            if (fdim==m_lin) tipo.set(m_lin);
            if (fdim==m_sup) tipo.set(m_sup);
            if (fdim==m_vol) tipo.set(m_vol);
        }
        else tipo.set(m_nodos);
 
        return true;
        KRATOS_CATCH("")
    }
 
    //*******************************************************************************************
    //*******************************************************************************************
    //void UpdateNodePosition(ModelPart& ThisModelPart)
    //{
    //    KRATOS_TRY
    //
    //        for(int i = 0 ; i < n.len ; i++)
    //        {
    //            Node& r_node = ThisModelPart.Nodes()[i+1];
    //            n[i][0] = r_node.X();
    //            n[i][1] = r_node.Y();
    //            n[i][2] = r_node.Z();
    //        }
    //        KRATOS_CATCH("")
    //}
 
    //*******************************************************************************************
    //*******************************************************************************************
    //void SetNodalH(ModelPart& ThisModelPart)
    //{
    //    KRATOS_TRY
    //
    //      ModelPart::NodesContainerType::iterator in = ThisModelPart.NodesBegin();
    //        for(int i = 0 ; i < n.len ; i++)
    //        {
    //            Node& r_node = *(in);
    //          double KratosH = r_node.FastGetSolutionStepValue(NODAL_H);
 
    //          if(r_node.FastGetSolutionStepValue(IS_BOUNDARY) == 0) //only on internal nodes
    //          {
    //              n[i].h = KratosH;
    //          }
 
    //          //KRATOS_WATCH(n[i].h);
    //          in++;
    //        }
    //        KRATOS_CATCH("")
    //}
 
    //*******************************************************************************************
    //*******************************************************************************************
    void ReGenerateMesh(
        ModelPart& ThisModelPart ,
        Element const& rReferenceElement,
        Condition const& rReferenceBoundaryCondition,
        double my_alpha = 1.4)
    {
        KRATOS_TRY
 
        delaunay(true,false);
        alpha_shape(my_alpha);
 
        //mk_vecinos_y_frontera();
        mk_nn();
        //mk_frontera();
 
        ThisModelPart.Elements().clear();
 
        for(int i = 0 ; i < n.len ; i++)
        {
            Node& r_node = ThisModelPart.Nodes()[i+1];
            r_node.X() = n[i][0];
            r_node.Y() = n[i][1];
            r_node.Z() = n[i][2];
        }
 
        Properties::Pointer properties = ThisModelPart.GetMesh().pGetProperties(1);
 
        ModelPart::NodesContainerType& r_model_nodes = ThisModelPart.Nodes();
 
        //adding domain elements
        for(int i_element = 0 ; i_element < e.len ; i_element++)
        {
            Element::GeometryType temp;
 
            Element::GeometryType::Pointer p_geometry;
 
            int nv = e[i_element].nv();
            //temp.reserve(nv);
            for(int i_node = 0 ; i_node < nv ; i_node++)
            {
                int index = e[i_element][i_node] + 1;
                temp.push_back(r_model_nodes(index));
            }
 
            p_geometry = Element::GeometryType::Pointer(new Element::GeometryType(temp));
 
            Element::Pointer p_element = rReferenceElement.Create(i_element + 1, *p_geometry, properties);
 
            ThisModelPart.Elements().push_back(p_element);
        }
 
 
        //adding boundary elements
        //ThisModelPart.Conditions().clear();
        //for(int i_bound = 0; i_bound < frontera.len ; i_bound++)
        //{
        //  for(int i_element = 0 ; i_element < frontera[i_bound].e.len ; i_element++)
        //  {
        //      Condition::NodesArrayType temp;
 
        //      for(int i_node = 0 ; i_node < frontera[i_bound].e[i_element].nv() ; i_node++)
        //      {
        //          int index = frontera[i_bound].e[i_element][i_node] + 1;
        //          temp.push_back(r_model_nodes(index));
        //      }
 
        //      int pos = ThisModelPart.Conditions().size();
        //      Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(pos, temp, properties);
 
        //      ThisModelPart.Conditions().push_back(p_cond);
 
        //  }
        //}
 
        //resetting the boundary flags needed on the frontier
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
        {
            i_node->GetSolutionStepValue(IS_BOUNDARY) = 0.00;
            i_node->GetSolutionStepValue(IS_FREE_SURFACE) = 0.00;
        }
 
        //setting the flags needed on the frontier
        for(int i_bound = 0; i_bound < frontera.len ; i_bound++)
        {
            for(int ii = 0 ; ii < frontera[i_bound].n.len ; ii++)
            {
                //attention: nodes identified with 1 are nodes in a larger domain
                //nodes identified with 2 identify nodes which are boundaries of a single isolated element
 
                if(frontera[i_bound].e.len > 3) //in this case the frontier is "big" as it has more than three edges
                {
                    ThisModelPart.Nodes()[frontera[i_bound].n[ii]+1].GetSolutionStepValue(IS_BOUNDARY) = 1;
                }
                else
                {
                    ThisModelPart.Nodes()[frontera[i_bound].n[ii]+1].GetSolutionStepValue(IS_BOUNDARY) = 2;
                }
 
                //attention: nodes which
            }
        }
 
        //free surface nodes are nodes of the boundary which are not nodes of the structure
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
        {
            if(i_node->GetSolutionStepValue(IS_BOUNDARY) == 1.00 && i_node->GetSolutionStepValue(IS_STRUCTURE) != 1)
                i_node->GetSolutionStepValue(IS_FREE_SURFACE) = 1.00;
        }
 
        for(int i = 0 ; i < nn.len ; i++)
        {
            int numb_of_neighb = nn[i].len;
            GlobalPointersVector< Node >& neighbours = (r_model_nodes[i+1].GetValues(NEIGHBOUR_NODES));
            //PointerVector< Node >& neighbours = (r_model_nodes[i+1].GetValues(NEIGHBOUR_NODES));
 
            neighbours.clear();
            neighbours.reserve(numb_of_neighb);
            for(int j = 0; j<numb_of_neighb; j++)
            {
                int ii = nn[i][j]+1;
                //neighbours.push_back(Kratos::weak_ptr< Node >( r_model_nodes(ii) ) );
                neighbours.push_back( r_model_nodes(ii)  );
            }
        }
 
        //perform an update to ensure the database is allocated correctly
        for(ModelPart::NodeIterator i_node = ThisModelPart.NodesBegin() ; i_node != ThisModelPart.NodesEnd() ; ++i_node)
        {
            (i_node->SolutionStepData()).Update();
        }
 
 
        KRATOS_CATCH("")
    }
 
    //*******************************************************************************************
    //*******************************************************************************************
    //void RefineMesh(ModelPart& ThisModelPart)
    //{
    //  KRATOS_TRY
    //      //      voronoi v(this);
    //      //      v.parcial=false;
 
    //      //      if (!v.delaunay()) // triangulacion (del para que no meta ni saque ptos)
    //      //        return ; //false
 
    //      //      double alpha=1.2;
    //      //      v.alpha_shape(alpha);// alpha shape
 
    //      //misma cantidad de nodos porque hice delaunay sin refinar
    //      for(int i=0;i<n.len;i++)  n[i].h=fabs(n[i][0]/10.00 + n[i][1]/10.00) + 0.2;
 
    //  //       bool calpha = false; // Don't do alpha_shape
    //  // //       bool calpha = true; // Do alpha_shape
    //  //      if (!v.refina_esferas(calpha)) // triangulacion
    //  //        return ; //false
    //  // //     if (slivers) v.rm_slivers();
    //  //      v.s2e();
    //  //      mk_nn();
 
    //  delaunay_refinado(true);
    //  ThisModelPart.Nodes().clear();
    //  ThisModelPart.Elements().clear();
 
    //  for(int i = 0 ; i < n.len ; i++)
    //  {
    //      Node::Pointer p_node = Kratos::Node::Pointer(new Kratos::Node(0,n[i][0], n[i][1], n[i][2]));
    //      //Node::Pointer p_node = n[i].GetKratosNode();
    //      p_node->SetId(i+1);
    //      ThisModelPart.Nodes().push_back(p_node);
    //  }
 
    //  //Node zero(0,0,0,0);
    //  Properties::Pointer properties(new Properties);
 
    //  for(int i_element = 0 ; i_element < e.len ; i_element++)
    //  {
    //      Element::GeometryType temp;
    //      Element::GeometryType::Pointer p_geometry;
 
    //      for(int i_node = 0 ; i_node < e[i_element].nv() ; i_node++)
    //          temp.push_back(ThisModelPart.Nodes()[e[i_element][i_node] + 1]);
 
 
    //      //if(e[i_element].tipo() == e_triangulo)
    //      //  p_geometry =  Element::GeometryType::Pointer(new Triangle2D<Node >(temp.Points()));
    //      //else
    //      p_geometry = Element::GeometryType::Pointer(new Element::GeometryType(temp));
 
    //      Element::Pointer p_element(new Element(i_element + 1, p_geometry, properties));
    //      ThisModelPart.Elements().push_back(p_element);
    //  }
    //  KRATOS_CATCH("")
    //}
 
 
 
 
 
 
 
 
    virtual std::string Info() const
    {
        return "";
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const {}
 
    virtual void PrintData(std::ostream& rOStream) const {}
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
 
 
 
 
 
    MeshSuiteModeler& operator=(MeshSuiteModeler const& rOther);
 
    MeshSuiteModeler(MeshSuiteModeler const& rOther);
 
 
 
}; // Class MeshSuiteModeler
 
 
 
 
 
 
inline std::istream& operator >> (std::istream& rIStream,
                                  MeshSuiteModeler& rThis);
 
inline std::ostream& operator << (std::ostream& rOStream,
                                  const MeshSuiteModeler& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_MESH_SUITE_MODELER_H_INCLUDED  defined