CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

A Complex_3_in_Triangulation_3 is already based on a Regular
triangulation which is based on a Triangulation_data_structure_3.
You can use the copy_tds() [1] member function.


Sebastien.

[1] 
http://www.cgal.org/Manual/latest/doc_html/cgal_manual/TriangulationDS_3_ref/Concept_TriangulationDataStructure_3.html

On 06/10/2012 10:07 PM, Lukasz Kaczmarczyk wrote:
> On 9 Jun 2012, at 23:39, Lukasz Kaczmarczyk wrote:
>
> > Hello,
> >
> > I try to copy mesh for MOAB (MeshOrientedDatabase) into 
> > TriangulationDataStructure_3. I hope that I could re-triangulate some ports 
> > of it and add/move some nodes in volume. However I have some problem with 
> > external tets and its neighbours, could you tell me how I should construct 
> > neighbours for boundary faces?
>
>
> Hello again,
>
> I made some improvements, code is still far from prefect, after I add 
> additional nodes I get,
>
> Linux,
> sum of the indices != 6
> sum of the indices != 6
>
> MacOS X (the same code)
> terminate called after throwing an instance of 'CGAL::Precondition_exception'
>    what():  CGAL ERROR: precondition violation!
> Expr: orientation(p0, p1, p2, p3) == POSITIVE
> File: 
> /opt/build_for_gcc-mp-4.4/CGAL-3.9/include/CGAL/Regular_triangulation_3.h
>
> I suspect that Its something wrong with my ordering. However if I not insert 
> any nodes there are any complains form CGAL.
>
> Regards,
> Lukasz
>
>
> void make_triangulation(Interface&  moab,const EntityHandle meshset_in,const 
> EntityHandle meshset_out,Range&add_nodes) {
>    Tag th_weight;
>    rval = moab.tag_get_handle(RT_TRIANGULATION_WEIGHT,th_weight);  
> CHKERR(rval);
>    Rt T;
>    T.tds().set_dimension(3);
>    //cout<<"nb. cells"<<T.number_of_cells()<<endl;
>    //cout<<"nb. finite cells"<<T.number_of_finite_cells()<<endl;
>    std::vector<Weighted_point>  P;
>    // generate points on a 3D grid
>    Range elems;
>    moab.get_entities_by_type(meshset_in,MBTET,elems);
>    Range moab_nodes;
>    moab.get_connectivity(elems,moab_nodes);
>    map<EntityHandle,Vertex_handle>  moab_cgal_vertex_map;
>    for(Range::iterator moab_nit = 
> moab_nodes.begin();moab_nit!=moab_nodes.end();moab_nit++) {
>      double weigth;
>      rval = moab.tag_get_data(th_weight,&*moab_nit,1,&weigth); CHKERR(rval);
>      double coords[3];
>      moab.get_coords(&*moab_nit,1,coords);
>      Point_3_with_moab_handle p(coords,*moab_nit);
>      //P.push_back(Weighted_point(p,weigth));
>      moab_cgal_vertex_map[*moab_nit] = T.tds().create_vertex();
>      moab_cgal_vertex_map[*moab_nit]->set_point(Weighted_point(p,weigth));
>    }
>    //cout<<"nb 
> vertexs"<<T.tds().number_of_vertices()<<""<<moab_nodes.size()<<endl;
>    map<EntityHandle,Cell_handle>  moab_cgal_cell_map;
>    for(Range::iterator moab_eit = 
> elems.begin();moab_eit!=elems.end();moab_eit++) {
>      //cout<<"nb. cells"<<T.number_of_cells()<<endl;
>      //cout<<"nb. finite cells"<<T.number_of_finite_cells()<<endl;
>      const EntityHandle* conn;
>      int num_nodes;
>      moab.get_connectivity(*moab_eit,conn,num_nodes);
>      assert(num_nodes == 4);
>      moab_cgal_cell_map[*moab_eit] = T.tds().create_cell(
>        moab_cgal_vertex_map[conn[0]],
>        moab_cgal_vertex_map[conn[1]],
>        moab_cgal_vertex_map[conn[2]],
>        moab_cgal_vertex_map[conn[3]] );
>    }
>    //cout<<"nb. cells"<<T.number_of_cells()<<endl;
>    //cout<<"nb. finite cells"<<T.number_of_finite_cells()<<endl;
>    map<EntityHandle,Cell_handle>  moab_cgal_face_map;
>    Range skin_faces;
>    Skinner tool(&moab);
>    tool.find_skin(elems,2,skin_faces);
>    for(Range::iterator moab_fit = 
> skin_faces.begin();moab_fit!=skin_faces.end();moab_fit++) {
>      Range adj_elems;
>      moab.get_adjacencies(&*moab_fit,1,3,false,adj_elems);
>      assert(adj_elems.size() == 1);
>      const EntityHandle* conn_face;
>      int num_nodes;
>      moab.get_connectivity(*moab_fit,conn_face,num_nodes);
>      assert(num_nodes == 3);
>      const EntityHandle* conn_elem;
>      moab.get_connectivity(*adj_elems.begin(),conn_elem,num_nodes);
>      assert(num_nodes == 4);
>      moab_cgal_face_map[*moab_fit] = T.tds().create_cell(T.infinite_cell());
>      for(int ii = 0;ii<4;ii++) {
>        int jj = ii;
>        if(ii == 0) jj = 1;
>        if(ii == 1) jj = 0;
>        if(find(&conn_face[0],&conn_face[3],conn_elem[ii])!=&conn_face[3]) {
>         
> moab_cgal_face_map[*moab_fit]->set_vertex(jj,moab_cgal_vertex_map[conn_elem[ii]]);
>        } else {
>         moab_cgal_face_map[*moab_fit]->set_vertex(jj,T.infinite_vertex());
>         
> moab_cgal_face_map[*moab_fit]->set_neighbor(jj,moab_cgal_cell_map[*adj_elems.begin()]);
>        }
>      }
>    }
>    //cout<<"nb cells"<<T.tds().number_of_cells()<<""<<elems.size()<<endl;
>    //cout<<"nb. finite cells"<<T.number_of_finite_cells()<<endl;
>    for(Range::iterator moab_nit = 
> moab_nodes.begin();moab_nit!=moab_nodes.end();moab_nit++) {
>      Range adj_elems;
>      moab.get_adjacencies(&*moab_nit,1,3,false,adj_elems);
>      
> moab_cgal_vertex_map[*moab_nit]->set_cell(moab_cgal_cell_map[*adj_elems.begin()]);
>    }
>    const int moab_cgall_cannonical_map[4] = { 1,2,0,3 };
>    for(Range::iterator moab_eit = 
> elems.begin();moab_eit!=elems.end();moab_eit++) {
>      for(int ff = 0;ff<4;ff++) {
>        EntityHandle face;
>        moab.side_element(*moab_eit,2,moab_cgall_cannonical_map[ff],face);
>        Range adj_elems;
>        moab.get_adjacencies(&face,1,3,false,adj_elems);
>        assert(adj_elems.size()>= 1);
>        assert(adj_elems.size()<= 2);
>        if(adj_elems.size() == 1) {
>         
> moab_cgal_cell_map[*moab_eit]->set_neighbor(ff,moab_cgal_face_map[face]);
>        } else {
>         adj_elems.erase(*moab_eit);
>         
> moab_cgal_cell_map[*moab_eit]->set_neighbor(ff,moab_cgal_cell_map[*adj_elems.begin()]);
>        }
>      }
>    }
>    for(Range::iterator moab_fit = 
> skin_faces.begin();moab_fit!=skin_faces.end();moab_fit++) {
>      Range adj_edges;
>      moab.get_adjacencies(&*moab_fit,1,1,true,adj_edges);
>      assert(adj_edges.size() == 3);
>      for(int ee = 0;ee<3;ee++) {
>        EntityHandle edge = adj_edges[ee];
>        Range adj_faces;
>        moab.get_adjacencies(&edge,1,2,true,adj_faces);
>        adj_faces = intersect(adj_faces,skin_faces);     
>        assert(adj_faces.size()==2);
>        adj_faces.erase(*moab_fit);
>        assert(adj_faces.size()==1);
>        const EntityHandle* conn_edge;
>        int num_nodes;
>        moab.get_connectivity(edge,conn_edge,num_nodes);
>        assert(num_nodes==2);
>        int ii = 0;
>        for(;ii<4;ii++) {
>         if( T.is_infinite(moab_cgal_face_map[*moab_fit]->vertex(ii)) ) {
>           
> moab_cgal_face_map[*moab_fit]->vertex(ii)->set_cell(moab_cgal_face_map[*moab_fit]);
>           continue;
>         }
>         
> assert(moab_cgal_face_map[*moab_fit]->vertex(ii)->point().h!=(EntityHandle)-1);
>         if(
>           moab_cgal_face_map[*moab_fit]->vertex(ii)->point().h == 
> conn_edge[0] ||
>           moab_cgal_face_map[*moab_fit]->vertex(ii)->point().h == 
> conn_edge[1] ) continue;
>         
> moab_cgal_face_map[*moab_fit]->set_neighbor(ii,moab_cgal_face_map[*adj_faces.begin()]);
>         break;
>        }
>      }
>    }
>    //T.tds().delete_cell(T.all_cells_begin());
>    //Add nodes;
>    for(Range::iterator moab_nit = 
> add_nodes.begin();moab_nit!=add_nodes.end();moab_nit++) {
>      double weigth;
>      rval = moab.tag_get_data(th_weight,&*moab_nit,1,&weigth); CHKERR(rval);
>      double coords[3];
>      moab.get_coords(&*moab_nit,1,coords);
>      Point_3_with_moab_handle p(coords,*moab_nit);
>      P.push_back(Weighted_point(p,weigth));
>    }
>    T.insert(P.begin(),P.end());
>    assert( T.dimension() == 3 );
>    assert( T.is_valid(true) );
>    Rt::All_vertices_iterator cgal_vit_all = T.all_vertices_begin();
>    for(;cgal_vit_all!=T.all_vertices_end();++cgal_vit_all) {
>      //printf("aaa\n");
>      //if(T.is_infinite(cgal_eit_all)) printf("inf\n");
>      T.tds().is_valid(cgal_vit_all,true);
>    }
>    Rt::All_cells_iterator cgal_eit_all = T.all_cells_begin();
>    for(;cgal_eit_all!=T.all_cells_end();++cgal_eit_all) {
>      //printf("aaa\n");
>      //if(T.is_infinite(cgal_eit_all)) printf("inf\n");
>      T.tds().is_valid(cgal_eit_all,true);
>    }
>    // insert all points in a row (this is faster than one insert() at a time).
>    Rt::Finite_cells_iterator cgal_eit = T.finite_cells_begin();
>    for(;cgal_eit!=T.finite_cells_end();++cgal_eit) {
>      EntityHandle conn[4];
>      for(int ii = 0;ii<4;ii++) {
>        assert(cgal_eit->vertex(ii)->point().h!=(EntityHandle)-1);
>        conn[ii] = cgal_eit->vertex(ii)->point().h;
>      }
>      EntityHandle ent;
>      moab.create_element(MBTET,conn,4,ent);
>      moab.add_entities(meshset_out,&ent,1);
>    }
> }