CGAL users discussion list

[Lukasz Kaczmarczyk <>]

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);
  }
}