MC_connectivity_index.cpp

Go to the documentation of this file.

#include <MC_int_pair.hpp>

#include <MC_connectivity_index.hpp>
#include <MC_double_vector.hpp>

namespace mesh_conv
{

    MC_connectivity_index::MC_connectivity_index() :MC_int_vector_vector(){}
    MC_connectivity_index::MC_connectivity_index(const MC_int_vector& v0) :MC_int_vector_vector(v0){}
    MC_connectivity_index::MC_connectivity_index(const MC_int_vector& v0,const MC_int_vector& v1) :MC_int_vector_vector(v0,v1){}
    MC_connectivity_index::MC_connectivity_index(const MC_int_vector& v0,const MC_int_vector& v1,const MC_int_vector& v2) :MC_int_vector_vector(v0,v1,v2) {}
    MC_connectivity_index::MC_connectivity_index(const MC_int_vector& v0,const MC_int_vector& v1,const MC_int_vector& v2,const MC_int_vector& v3) :MC_int_vector_vector(v0,v1,v2,v3) {}
    MC_connectivity_index::MC_connectivity_index(const MC_int_vector_vector& v0)
            :MC_int_vector_vector(v0)
    {update_neighbors();}
    MC_connectivity_index::MC_connectivity_index(const MC_connectivity_index& v0) :MC_int_vector_vector(v0)
    {
        if(v0.ring_1.size()==0 || v0.star_1.size()==0)
            update_neighbors();
        else
        {ring_1=v0.ring_1;star_1=v0.star_1;}
    }

    MC_connectivity_index& MC_connectivity_index::build_ring()
    {
        ring_1.clear();

        int N_polygon=size();
        for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
        {
            MC_int_vector poly_index=(*this)[k_polygon];
            int N_poly=poly_index.size();
            for(int k=0;k<N_poly;++k)
            {
                int u0=poly_index[k];
                int u1=poly_index[(k+1)%N_poly];

                // set u0->u1
                std::map<int,std::set <int> > :: iterator it_1=ring_1.find(u0);
                if(it_1==ring_1.end())//create the new entry
                    it_1=ring_1.insert(std::pair<int,std::set<int> >(u0,std::set<int>())).first;
                it_1->second.insert(u1);

                // set u1->u0
                std::map<int,std::set <int> > ::iterator it_2=ring_1.find(u1);
                if(it_2==ring_1.end())//create the new entry
                    it_2=ring_1.insert(std::pair<int,std::set<int> >(u1,std::set <int>() )).first;
                it_2->second.insert(u0);

            }
        }
        return *this;
    }
    const std::map<int,std::set<int> >& MC_connectivity_index::ring() const
    {return ring_1;}

    MC_connectivity_index& MC_connectivity_index::build_star()
    {
        star_1.clear();
        int N_polygon=size();
        for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
        {
            MC_int_vector poly_index=(*this)[k_polygon];
            int N_poly=poly_index.size();
            for(int k=0;k<N_poly;++k)
            {
                std::map<int,std::set<int> > :: iterator it=star_1.find(poly_index[k]);
                if(it==star_1.end())
                    it=star_1.insert(std::pair<int,std::set<int> > (poly_index[k],std::set<int>() )).first;
                it->second.insert(k_polygon);
            }
        }
        return *this;

    }

    const std::map<int,std::set<int> >& MC_connectivity_index::star() const
    {return star_1;}

    MC_connectivity_index& MC_connectivity_index::update_neighbors()
    {
        if(size()<20000)
        {
            build_ring();
            build_star();
        }
       // else
       // {std::cout<<"Warning, large mesh do not update automatically the neighbooring, do it manually"<<std::endl;}
        return *this;
    }

    std::set <MC_int_pair,MC_int_pair_less> MC_connectivity_index::boundary() const
    {
        std::vector <std::set <int> > index_edge;
        int N_polygon=size();
        for(int k_polygon=0;k_polygon<N_polygon;++k_polygon)
        {
            MC_int_vector index_polygon=(*this)[k_polygon];
            int size_polygon=index_polygon.size();
            for(int k_edge=0;k_edge<size_polygon;++k_edge)
            {
                int v0=index_polygon[k_edge];
                int v1=index_polygon[(k_edge+1)%size_polygon];

                int u0=std::min(v0,v1),u1=std::max(v0,v1);
                if(u0>=static_cast <int> (index_edge.size()))
                {
                    index_edge.resize(u0+1);
                    index_edge[u0].insert(u1);
                }
                else
                {
                    std::set <int> :: iterator it;
                    it=index_edge[u0].find(u1);
                    if(it!=index_edge[u0].end())
                        index_edge[u0].erase(it);
                    else
                        index_edge[u0].insert(u1);
                }

            }
        }


        std::set <MC_int_pair,MC_int_pair_less> boundary_edge;
        int N_vertex=index_edge.size();
        for(int k_vertex=0;k_vertex<N_vertex;++k_vertex)
            if(index_edge[k_vertex].size()>0)
            {
              std::set <int> :: const_iterator it;
              std::set <int> :: const_iterator it_end=index_edge[k_vertex].end();
              for(it=index_edge[k_vertex].begin();it!=it_end;++it)
                  boundary_edge.insert(MC_int_pair(k_vertex,*it));
            }

        return boundary_edge;
    }

    std::set<int> MC_connectivity_index::used_vertex() const
    {
        std::set <int> vertices;
        int N=size();
        for(int k_poly=0;k_poly<N;++k_poly)
        {
            MC_int_vector pol=(*this)[k_poly];
            int N_vertex=pol.size();
            for(int k_vertex=0;k_vertex<N_vertex;++k_vertex)
                vertices.insert(pol[k_vertex]);
        }
        return vertices;

    }

    bool MC_connectivity_index::is_triangle_only() const
    {
        int N=size();
        for(int k_poly=0;k_poly<N;++k_poly)
            if((*this)[k_poly].size()!=3)
                return false;
        return true;
    }

    std::set<int> MC_connectivity_index::boundary_vertex() const
    {
        std::set <MC_int_pair,MC_int_pair_less> bnd=boundary();
        std::set <MC_int_pair,MC_int_pair_less> :: const_iterator it_bnd=bnd.begin();
        std::set <MC_int_pair,MC_int_pair_less> :: const_iterator it_bnd_end=bnd.end();

        std::set<int> index_boundary;
        for(;it_bnd!=it_bnd_end;++it_bnd)
            for(int k_dim=0;k_dim<2;++k_dim)
                index_boundary.insert((*it_bnd)[k_dim]);

        return index_boundary;
    }
    std::set<int> MC_connectivity_index::boundary_polygon() const
    {
        std::map<int,std::set<int> > ::const_iterator it_star;
        std::map<int,std::set<int> > ::const_iterator it_star_end=star_1.end();

        std::set<int> set_bnd_polygon;

        MC_int_vector vertex_bnd=boundary_vertex();
        int N_vertex_bnd=vertex_bnd.size();
        for(int k_vertex=0;k_vertex<N_vertex_bnd;++k_vertex)
        {
            it_star=star_1.find(vertex_bnd(k_vertex));
            if(it_star==it_star_end)
            {std::cout<<"Error in MC_connectivity_index::boundary_polygon(), cannot find vertex"<<vertex_bnd(k_vertex)<<", at k="<<k_vertex<<" for star_1.size()="<<star_1.size()<<std::endl;exit(-1);}

           std::set<int> current_neigh=it_star->second;
           set_bnd_polygon.insert(current_neigh.begin(),current_neigh.end());
        }

        return set_bnd_polygon;
    }

    MC_connectivity_index MC_connectivity_index::triangulated() const
    {
        MC_connectivity_index connectivity_triangle;
        for(int k=0,N=size();k<N;++k)
        {
            MC_int_vector poly=(*this)(k);
            if(poly.size()>3)
            {
                for(int k_v=0,N_tri=poly.size()-2;k_v<N_tri;++k_v)
                    connectivity_triangle.add(MC_int_vector(poly[0],poly[k_v+1],poly[k_v+2]));
            }
            else
                connectivity_triangle.add(poly);
        }

        return connectivity_triangle;
    }

    MC_connectivity_index MC_connectivity_index::fliped_polygon() const
    {
        MC_connectivity_index new_connectivity;new_connectivity.resize(size());
        for(int k=0,N=size();k<N;++k)
        {
            MC_int_vector current=(*this)[k];
            new_connectivity[k]=current(MC_int_vector::linspace(current.size()-1,0,-1));
        }
        return new_connectivity;
    }

    std::map<int,std::set<int> > MC_connectivity_index::polygon_neighbors() const
    {
        std::map<int,std::set<int> > neigh;
        std::map<int,std::set<int> > s=star();

        for(int k=0,N=size();k<N;++k)
        {
            MC_int_vector index_polygon=(*this)(k);

            for(int k_v=0,N_v=index_polygon.size();k_v<N_v;++k_v)
            {
                std::set <int> current_neighbors=s[index_polygon[k_v]];
                current_neighbors.erase(k);

                std::map<int,std::set<int> > ::iterator it=neigh.find(k);
                if(it==neigh.end())
                    neigh.insert(std::make_pair(k,current_neighbors));
                else
                {
                    for(std::set<int>::iterator it_2=current_neighbors.begin(),it_2_end=current_neighbors.end();it_2!=it_2_end;++it_2)
                    {
                        it->second.insert(*it_2);
                    }
                }
            }

        }
        return neigh;
    }

    MC_int_vector MC_connectivity_index::triangle_propagation_using_vertex_field(const MC_int_vector& seed_vertices,const MC_double_vector& field_per_vertex,const double& threshold,const std::set<int>& invalidated_triangle) const
    {
        if(static_cast<int>(star_1.size())!=field_per_vertex.size())
        {std::cout<<"Error in MC_connectivity_index::triangle_propagation_using_vertex_field(), size problem with star_1 size: "<<star_1.size()<<std::endl;exit(-1);}

        if(static_cast<int>(ring_1.size())!=field_per_vertex.size())
        {std::cout<<"Error in MC_connectivity_index::triangle_propagation_using_vertex_field(), size problem with ring_1 size: "<<ring_1.size()<<std::endl;exit(-1);}


        std::set<int> selected_triangle;

        std::set<int> selected_vertices;
        std::set<int> vertices_to_treat=seed_vertices.to_set();


        do
        {
            int current_vertex=*vertices_to_treat.begin();
            vertices_to_treat.erase(vertices_to_treat.begin());


            bool is_forbidden_triangle=false;
            const std::set<int>& current_star=star_1.find(current_vertex)->second;
            for(std::set<int>::const_iterator it_star=current_star.begin(),it_star_end=current_star.end();it_star!=it_star_end;++it_star)
            {
                if(invalidated_triangle.find(*it_star)!=invalidated_triangle.end())
                {
                    is_forbidden_triangle=true;break;
                }
            }

            if(is_forbidden_triangle==false)
            {
                if(field_per_vertex[current_vertex]<threshold)
                    selected_vertices.insert(current_vertex);

                const std::set<int>& current_neighbors=ring_1.find(current_vertex)->second;
                for(std::set<int>::const_iterator it=current_neighbors.begin(),it_end=current_neighbors.end();it!=it_end;++it)
                {
                    if(field_per_vertex[*it]<threshold && selected_vertices.find(*it)==selected_vertices.end() && invalidated_triangle.find(*it)==invalidated_triangle.end())
                        vertices_to_treat.insert(*it);
                }
            }

        }while(vertices_to_treat.size()>0);

        for(std::set<int>::const_iterator it=selected_vertices.begin(),it_end=selected_vertices.end();it!=it_end;++it)
        {
            const std::set<int>& tri_neighbors=star_1.find(*it)->second;
            for(std::set<int>::const_iterator it_tri=tri_neighbors.begin(),it_tri_end=tri_neighbors.end();it_tri!=it_tri_end;++it_tri)
                selected_triangle.insert(*it_tri);
        }

        return MC_int_vector(selected_triangle);
    }

    std::set<MC_int_pair,MC_int_pair_less> MC_connectivity_index::edges() const
    {
        std::set<MC_int_pair,MC_int_pair_less> e;
        for(int k=0,N=size();k<N;++k)
        {
            const MC_int_vector& poly=(*this)(k);
            int N_poly=poly.size();

            for(int k_v=0;k_v<N_poly;++k_v)
                e.insert(MC_int_pair(poly[k_v],poly[(k_v+1)%N_poly]));
        }
        return e;
    }


}