RINGMesh  Version 5.0.0
A programming library for geological model meshes
geomodel_mesh.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2012-2017, Association Scientifique pour la Geologie et ses
3  * Applications (ASGA). All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  * * Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * * Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * * Neither the name of ASGA nor the
13  * names of its contributors may be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
18  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ASGA BE LIABLE FOR ANY DIRECT,
20  * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
25  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * http://www.ring-team.org
28  *
29  * RING Project
30  * Ecole Nationale Superieure de Geologie - GeoRessources
31  * 2 Rue du Doyen Marcel Roubault - TSA 70605
32  * 54518 VANDOEUVRE-LES-NANCY
33  * FRANCE
34  */
35 
37 
38 #include <numeric>
39 #include <stack>
40 
41 #include <geogram/basic/algorithm.h>
42 
43 #include <geogram/basic/permutation.h>
44 #include <geogram/mesh/mesh_geometry.h>
45 
47 
50 
52 #include <ringmesh/mesh/well.h>
53 
58 namespace
59 {
60  using namespace RINGMesh;
61 
62  template < index_t DIMENSION >
63  class GeoModelMeshPolygonsBaseSort
64  {
65  public:
66  GeoModelMeshPolygonsBaseSort( const SurfaceMesh< DIMENSION >& mesh,
67  const std::vector< index_t >& surface_id )
68  : mesh_( mesh ), surface_id_( surface_id )
69  {
70  }
71 
72  bool operator()( index_t i, index_t j ) const
73  {
74  if( surface_id_[i] != surface_id_[j] )
75  {
76  return surface_id_[i] < surface_id_[j];
77  }
78  else
79  {
80  return mesh_.nb_polygon_vertices( i )
81  < mesh_.nb_polygon_vertices( j );
82  }
83  }
84 
85  private:
86  const SurfaceMesh< DIMENSION >& mesh_;
87  const std::vector< index_t >& surface_id_;
88  };
89 
90  template < index_t DIMENSION >
91  class GeoModelMeshCellsSort
92  {
93  public:
94  GeoModelMeshCellsSort( const VolumeMesh< DIMENSION >& mesh,
95  const std::vector< index_t >& region_id )
96  : mesh_( mesh ), region_id_( region_id )
97  {
98  }
99 
100  bool operator()( index_t i, index_t j ) const
101  {
102  if( region_id_[i] != region_id_[j] )
103  {
104  return region_id_[i] < region_id_[j];
105  }
106  return mesh_.cell_type( i ) < mesh_.cell_type( j );
107  }
108 
109  private:
110  const VolumeMesh< DIMENSION >& mesh_;
111  const std::vector< index_t >& region_id_;
112  };
113 
114  std::string vertex_map_name()
115  {
116  return "model_vertex_map";
117  }
118 
119  template < index_t DIMENSION >
120  std::vector< index_t > cell_facets_around_vertex(
121  const VolumeMesh< DIMENSION >& mesh, index_t cell, index_t vertex_id )
122  {
123  std::vector< index_t > facets;
124  facets.reserve( mesh.nb_cell_facets( cell ) );
125  for( auto f : range( mesh.nb_cell_facets( cell ) ) )
126  {
127  for( auto v : range( mesh.nb_cell_facet_vertices(
128  CellLocalFacet( cell, f ) ) ) )
129  {
130  if( mesh.cell_facet_vertex( CellLocalFacet( cell, f ), v )
131  == vertex_id )
132  {
133  facets.push_back( f );
134  break;
135  }
136  }
137  }
138  return facets;
139  }
140 
141  template < index_t DIMENSION >
142  void copy_vertices( MeshBaseBuilder< DIMENSION >* builder,
143  const MeshBase< DIMENSION >& mesh )
144  {
145  builder->clear_vertices( true, false );
146  builder->create_vertices( mesh.nb_vertices() );
147  for( auto v : range( mesh.nb_vertices() ) )
148  {
149  builder->set_vertex( v, mesh.vertex( v ) );
150  }
151  }
152 } // namespace
153 
154 namespace RINGMesh
155 {
156  template < index_t DIMENSION >
159  : gmm_( gmm ), geomodel_( gm ), mesh_base_( nullptr )
160  {
161  }
162 
163  template < index_t DIMENSION >
166  const GeoModel< DIMENSION >& geomodel )
167  : geomodel_vertices_( geomodel_vertices ), geomodel_( geomodel )
168  {
175  }
176 
177  template < index_t DIMENSION >
179  {
180  gme_vertices_.clear();
182  }
183 
184  template < index_t DIMENSION >
186  geomodel_vertex_index( const gmme_id& mesh_entity_id,
187  index_t mesh_entity_vertex_index ) const
188  {
190  mesh_entity_vertex_index
191  < geomodel_.mesh_entity( mesh_entity_id ).nb_vertices() );
192 
193  return vertex_map( mesh_entity_id )[mesh_entity_vertex_index];
194  }
195 
196  template < index_t DIMENSION >
197  const std::vector< GMEVertex >& GeoModelMeshVerticesBase< DIMENSION >::
199  {
200  ringmesh_assert( v < gme_vertices_.size() );
201  return gme_vertices_[v];
202  }
203 
204  template < index_t DIMENSION >
205  std::vector< GMEVertex > GeoModelMeshVerticesBase< DIMENSION >::
207  index_t v, const MeshEntityType& mesh_entity_type ) const
208  {
209  const auto& all_gmes = mesh_entity_vertex_indices( v );
210  std::vector< GMEVertex > result;
211  result.reserve( all_gmes.size() );
212  for( const auto& vertex : all_gmes )
213  {
214  if( vertex.gmme.type() == mesh_entity_type )
215  {
216  result.push_back( vertex );
217  }
218  }
219  return result;
220  }
221 
222  template < index_t DIMENSION >
223  std::vector< index_t > GeoModelMeshVerticesBase< DIMENSION >::
225  index_t v, const gmme_id& mesh_entity_id ) const
226  {
227  std::vector< index_t > result;
228  auto all_gmes = mesh_entity_vertex_indices( v );
229  for( const auto& vertex : all_gmes )
230  {
231  if( vertex.gmme == mesh_entity_id )
232  {
233  result.push_back( vertex.v_index );
234  }
235  }
236  return result;
237  }
238 
239  template < index_t DIMENSION >
240  const std::vector< index_t >&
242  const gmme_id& mesh_entity_id ) const
243  {
244  return (
245  *vertex_maps_.at( mesh_entity_id.type() ) )[mesh_entity_id.index()];
246  }
247 
248  template < index_t DIMENSION >
249  std::vector< index_t >&
251  const gmme_id& mesh_entity_id )
252  {
253  return ( *vertex_maps_[mesh_entity_id.type()] )[mesh_entity_id.index()];
254  }
255 
256  template < index_t DIMENSION >
258  set_vertex_map_value( const gmme_id& mesh_entity_id,
259  index_t mesh_entity_vertex_index,
260  index_t geomodel_entity_vertex_index )
261  {
263  vertex_map( mesh_entity_id )[mesh_entity_vertex_index] =
264  geomodel_entity_vertex_index;
265  }
266 
267  template < index_t DIMENSION >
270  const GMEVertex& gme_vertex, index_t geomodel_vertex_index )
271  {
272  gme_vertices_[geomodel_vertex_index].push_back( gme_vertex );
273  }
274 
275  template < index_t DIMENSION >
278  {
279  const auto& all_mesh_entity_types =
280  geomodel_.entity_type_manager()
281  .mesh_entity_manager.mesh_entity_types();
282  for( const auto& cur_entity_type : all_mesh_entity_types )
283  {
284  auto nb_cur_type_entities =
285  geomodel_.nb_mesh_entities( cur_entity_type );
286  vertex_maps_.at( cur_entity_type )->clear();
287  vertex_maps_.at( cur_entity_type )->resize( nb_cur_type_entities );
288  for( auto e : range( nb_cur_type_entities ) )
289  {
290  resize_vertex_map( { cur_entity_type, e } );
291  }
292  }
293  }
294 
295  template < index_t DIMENSION >
296  std::vector< index_t >& GeoModelMeshVerticesBase< DIMENSION >::
298  {
299  ringmesh_assert( mesh_entity_id.index()
300  < vertex_maps_[mesh_entity_id.type()]->size() );
301  if( geomodel_vertices_.is_initialized() )
302  {
303  const auto& mesh_entity = geomodel_.mesh_entity( mesh_entity_id );
304  vertex_maps_.at( mesh_entity_id.type() )
305  ->at( mesh_entity_id.index() )
306  .resize( mesh_entity.nb_vertices(), NO_ID );
307  }
308  return vertex_map( mesh_entity_id );
309  }
310 
311  template < index_t DIMENSION >
314  const std::vector< index_t >& old2new )
315  {
316  const auto& all_mesh_entity_types =
317  geomodel_.entity_type_manager()
318  .mesh_entity_manager.mesh_entity_types();
319  for( const auto& cur_entity_type : all_mesh_entity_types )
320  {
321  for( auto e :
322  range( geomodel_.nb_mesh_entities( cur_entity_type ) ) )
323  {
324  const auto& E = geomodel_.mesh_entity( cur_entity_type, e );
325  auto id = E.gmme();
326  for( auto v : range( E.nb_vertices() ) )
327  {
328  auto old_m_id = geomodel_vertex_index( id, v );
329  auto new_m_id = old2new[old_m_id];
330  set_vertex_map_value( id, v, new_m_id );
331 
332  // Merge gme_vertices information
333  if( find( gme_vertices_[new_m_id], GMEVertex( id, v ) )
334  == NO_ID )
335  {
336  gme_vertices_[new_m_id].emplace_back( id, v );
337  }
338  }
339  }
340  }
341  }
342 
343  template < index_t DIMENSION >
345  clear_vertex_map( const gmme_id& mesh_entity_id )
346  {
347  resize_all_mesh_entity_vertex_maps( mesh_entity_id.type() );
348  if( !vertex_maps_.at( mesh_entity_id.type() )
349  ->at( mesh_entity_id.index() )
350  .empty() )
351  {
352  vertex_maps_.at( mesh_entity_id.type() )
353  ->at( mesh_entity_id.index() )
354  .clear();
355  }
356  }
357 
358  template < index_t DIMENSION >
361  {
362  auto& mesh_entity_vertex_map = resize_vertex_map( mesh_entity_id );
363  const auto& E = geomodel_.mesh_entity( mesh_entity_id );
364  for( auto v : range( E.nb_vertices() ) )
365  {
366  mesh_entity_vertex_map[v] =
367  geomodel_vertices_.nn_search().get_closest_neighbor(
368  E.vertex( v ) );
369  }
370  }
371 
372  template < index_t DIMENSION >
375  const gmme_id& mesh_entity_id )
376  {
377  resize_all_mesh_entity_vertex_maps( mesh_entity_id.type() );
378  if( !is_mesh_entity_vertex_map_initialized( mesh_entity_id ) )
379  {
380  initialize_mesh_entity_vertex_map( mesh_entity_id );
381  return false;
382  }
383  return true;
384  }
385 
386  template < index_t DIMENSION >
389  const gmme_id& mesh_entity_id ) const
390  {
391  return !vertex_maps_.find( mesh_entity_id.type() )
392  ->second->at( mesh_entity_id.index() )
393  .empty();
394  }
395 
396  template < index_t DIMENSION >
399  {
400  for( auto& vertex_map : vertex_maps_ )
401  {
402  for( auto e : range( vertex_map.second->size() ) )
403  {
404  vertex_map.second->at( e ).clear();
405  }
406  vertex_map.second->clear();
407  }
408  }
409 
410  template < index_t DIMENSION >
413  {
414  vertex_maps_.at( type )->resize( geomodel_.nb_mesh_entities( type ) );
415  }
416 
417  template < index_t DIMENSION >
418  GEO::AttributesManager& GeoModelMeshVerticesBase< DIMENSION >::
420  const gmme_id& mesh_entity_id ) const
421  {
422  const auto& mesh_entity = geomodel_.mesh_entity( mesh_entity_id );
423  return mesh_entity.vertex_attribute_manager();
424  }
425 
426  template < index_t DIMENSION >
430  std::unique_ptr< PointSetMesh< DIMENSION > >& mesh )
431  : GeoModelMeshCommon< DIMENSION >( gmm, gm ),
432  mesh_( mesh ),
433  vertex_mapper_( *this, gmm.geomodel() )
434  {
435  this->set_mesh( mesh_.get() );
436  }
437 
438  template < index_t DIMENSION >
440  {
441  return mesh_->nb_vertices() > 0;
442  }
443 
444  template < index_t DIMENSION >
446  {
447  if( !is_initialized() )
448  {
449  const_cast< GeoModelMeshVerticesBase* >( this )->initialize();
450  }
451  }
452 
453  template < index_t DIMENSION >
454  index_t nb_entity_vertices( const GeoModel< DIMENSION >& geomodel,
455  const MeshEntityType& entity_type )
456  {
457  index_t count{ 0 };
458  for( auto i : range( geomodel.nb_mesh_entities( entity_type ) ) )
459  {
460  count += geomodel.mesh_entity( entity_type, i ).nb_vertices();
461  }
462  return count;
463  }
464 
465  template < index_t DIMENSION >
467  const GeoModel< DIMENSION >& geomodel,
468  const MeshEntityType& entity_type,
469  index_t& count )
470  {
471  auto mesh_builder =
473  for( auto i : range( geomodel.nb_mesh_entities( entity_type ) ) )
474  {
475  const auto& E = geomodel.mesh_entity( entity_type, i );
476  if( E.nb_vertices() == 0 )
477  {
478  continue;
479  }
480 
481  // Map and vertex
482  auto id = E.gmme();
483  for( auto v : range( E.nb_vertices() ) )
484  {
485  auto local_count = count + v;
486  mesh_builder->set_vertex( local_count, E.vertex( v ) );
487  // Map from vertices of MeshEntities to GeoModelMeshVerticesBase
488  vertex_mapper_.set_vertex_map_value( id, v, local_count );
489  vertex_mapper_.add_to_gme_vertices(
490  GMEVertex( id, v ), local_count );
491  }
492  // Global vertex index increment
493  count += E.nb_vertices();
494  }
495  }
496 
497  template < index_t DIMENSION >
499  {
500  index_t nb{ 0 };
501  nb += nb_entity_vertices(
502  this->geomodel_, Corner< DIMENSION >::type_name_static() );
503  nb += nb_entity_vertices(
504  this->geomodel_, Line< DIMENSION >::type_name_static() );
505  nb += nb_entity_vertices(
506  this->geomodel_, Surface< DIMENSION >::type_name_static() );
507  return nb;
508  }
509 
510  template < index_t DIMENSION >
512  {
513  // Total number of vertices in the
514  // Corners, Lines, Surfaces and Regions of the GeoModel
515  auto nb = nb_total_vertices();
516 
517  // Get out if no vertices
518  if( nb == 0 )
519  {
520  return;
521  }
522 
523  // Fill the vertices
524  auto builder =
526  builder->create_vertices( nb );
527  vertex_mapper_.clear_and_resize_geomodel_vertex_gmes( nb );
528  vertex_mapper_.bind_all_mesh_entity_vertex_maps();
529 
530  fill_vertices();
531 
532  // Remove colocated vertices
533  remove_colocated();
534  }
535 
536  template < index_t DIMENSION >
538  {
539  index_t count{ 0 };
540  fill_vertices_for_entity_type(
541  this->geomodel_, Corner< DIMENSION >::type_name_static(), count );
542  fill_vertices_for_entity_type(
543  this->geomodel_, Line< DIMENSION >::type_name_static(), count );
544  fill_vertices_for_entity_type(
545  this->geomodel_, Surface< DIMENSION >::type_name_static(), count );
546  return count;
547  }
548 
549  template < index_t DIMENSION >
551  {
552  this->gmm_.polygons.clear();
553  this->gmm_.wells.clear();
554  vertex_mapper_.clear();
555 
556  auto builder =
558  builder->clear( true, false );
559  }
560 
561  template < index_t DIMENSION >
563  const gmme_id& mesh_entity_id )
564  {
565  vertex_mapper_.clear_vertex_map( mesh_entity_id );
566  }
567 
568  template < index_t DIMENSION >
570  const gmme_id& mesh_entity_id )
571  {
572  vertex_mapper_.resize_vertex_map( mesh_entity_id );
573  }
574 
575  template < index_t DIMENSION >
577  {
578  test_and_initialize();
579  return mesh_->nb_vertices();
580  }
581 
582  template < index_t DIMENSION >
584  index_t v ) const
585  {
586  test_and_initialize();
587  ringmesh_assert( v < nb() );
588  return mesh_->vertex( v );
589  }
590 
591  template < index_t DIMENSION >
593  const vecn< DIMENSION >& p ) const
594  {
595  test_and_initialize();
596  const auto& colocator = mesh_->vertex_nn_search();
597  auto vertices = colocator.get_neighbors( p, this->geomodel_.epsilon() );
598  if( vertices.empty() )
599  {
600  return NO_ID;
601  }
602  ringmesh_assert( vertices.size() == 1 );
603  return vertices[0];
604  }
605 
606  template < index_t DIMENSION >
608  const gmme_id& mesh_entity, index_t entity_vertex_index ) const
609  {
610  test_and_initialize();
611  return vertex_mapper_.geomodel_vertex_index(
612  mesh_entity, entity_vertex_index );
613  }
614 
615  template < index_t DIMENSION >
617  const gmme_id& mesh_entity,
618  const ElementLocalVertex& element_local_vertex ) const
619  {
620  auto entity_vertex_index =
621  this->geomodel_.mesh_entity( mesh_entity )
622  .mesh_element_vertex_index( element_local_vertex );
623  return geomodel_vertex_id( mesh_entity, entity_vertex_index );
624  }
625 
626  template < index_t DIMENSION >
627  std::vector< index_t >
629  const gmme_id& mesh_entity, index_t geomodel_vertex_index ) const
630  {
631  test_and_initialize();
632  ringmesh_assert( geomodel_vertex_index < nb() );
633  return vertex_mapper_.mesh_entity_vertex_indices(
634  geomodel_vertex_index, mesh_entity );
635  }
636 
637  template < index_t DIMENSION >
638  const std::vector< GMEVertex >&
640  {
641  test_and_initialize();
642  return vertex_mapper_.mesh_entity_vertex_indices( v );
643  }
644 
645  template < index_t DIMENSION >
646  std::vector< GMEVertex >
648  const MeshEntityType& entity_type, index_t v ) const
649  {
650  test_and_initialize();
651  return vertex_mapper_.mesh_entity_vertex_indices( v, entity_type );
652  }
653 
654  template < index_t DIMENSION >
656  const vecn< DIMENSION >& point )
657  {
658  auto builder =
660  const auto index = builder->create_vertex( point );
661  vertex_mapper_.resize_geomodel_vertex_gmes( nb() );
662  return index;
663  }
664 
665  template < index_t DIMENSION >
667  const std::vector< vecn< DIMENSION > >& points )
668  {
669  ringmesh_assert( !points.empty() );
670  auto builder =
672  const auto start_index = builder->create_vertex( points[0] );
673  for( auto i : range( 1, points.size() ) )
674  {
675  builder->create_vertex( points[i] );
676  }
677  vertex_mapper_.resize_geomodel_vertex_gmes( nb() );
678  return start_index;
679  }
680 
681  template < index_t DIMENSION >
683  index_t v, const vecn< DIMENSION >& point )
684  {
685  test_and_initialize();
686  ringmesh_assert( v < nb() );
687  // Change the position of the unique_vertex
688  auto mesh_builder =
690  mesh_builder->set_vertex( v, point );
691 
692  GeoModelBuilder< DIMENSION > builder( this->geomodel_ );
693 
694  const auto& gme_v = gme_vertices( v );
695  for( const auto& info : gme_v )
696  {
697  builder.geometry.set_mesh_entity_vertex(
698  info.gmme, info.v_index, point, false );
699  }
700  }
701 
702  template < index_t DIMENSION >
704  const gmme_id& entity_id,
705  index_t entity_vertex_index,
706  index_t geomodel_vertex_index )
707  {
708  vertex_mapper_.set_vertex_map_value(
709  entity_id, entity_vertex_index, geomodel_vertex_index );
710  vertex_mapper_.add_to_gme_vertices(
711  GMEVertex( entity_id, entity_vertex_index ),
712  geomodel_vertex_index );
713  }
714 
715  template < index_t DIMENSION >
717  {
718  // Get out if nothing to do
719  // and compute the points if they are not initialized yet
720  if( nb() == 0 )
721  {
722  return;
723  }
724  // Identify and invalidate colocated vertices
725  index_t nb_colocalised_vertices{ NO_ID };
726  std::vector< index_t > old2new;
727  std::tie( nb_colocalised_vertices, old2new ) =
728  mesh_->vertex_nn_search().get_colocated_index_mapping(
729  this->geomodel_.epsilon() );
730  if( nb_colocalised_vertices > 0 )
731  {
732  erase_vertices( old2new );
733  }
734  }
735 
736  template < index_t DIMENSION >
738  std::vector< index_t >& to_delete )
739  {
740  ringmesh_assert( to_delete.size() == nb() );
741 
742  // For mesh vertices deletion
743  std::vector< bool > to_delete_bool( nb(), false );
744 
745  // Fill the delete information for geogram
746  // Recycle the to_delete vertex to get the mapping between
747  // new and old points. This is implemented to be the same
748  // as what is done in the delete_elements function in geogram
749  index_t nb_todelete{ 0 };
750  index_t cur{ 0 };
751  for( auto v : range( nb() ) )
752  {
753  if( to_delete[v] != v )
754  {
755  to_delete_bool[v] = true;
756  nb_todelete++;
757  if( to_delete[v] != NO_ID )
758  {
759  ringmesh_assert( to_delete[v] < v );
760  to_delete[v] = to_delete[to_delete[v]];
761  }
762  }
763  else
764  {
765  to_delete[v] = cur;
766  ++cur;
767  }
768  }
769  if( nb_todelete == 0 )
770  {
771  return;
772  }
773  if( nb_todelete == nb() )
774  {
775  // Clear everything
776  clear();
777  return;
778  }
779 
780  // Empty the gme_vertices_ of the deleted vertices and erase them
781  for( auto v : range( nb() ) )
782  {
783  vertex_mapper_.clear_geomodel_vertex_gmes( v );
784  }
785 
786  // Delete the vertices - false is to not remove
787  // isolated vertices (here all the vertices)
789  ->delete_vertices( to_delete_bool );
790 
791  vertex_mapper_.update_mesh_entity_maps_and_gmes( to_delete );
792  }
793 
794  template < index_t DIMENSION >
798  std::unique_ptr< PointSetMesh< DIMENSION > >& mesh )
799  : GeoModelMeshVerticesBase< DIMENSION >( gmm, gm, mesh )
800  {
801  }
802 
804  GeoModel3D& gm,
805  std::unique_ptr< PointSetMesh3D >& mesh )
806  : GeoModelMeshVerticesBase< 3 >( gmm, gm, mesh )
807  {
808  }
809 
811  {
812  this->gmm_.cells.clear();
813  GeoModelMeshVerticesBase3D::clear();
814  }
815 
817  {
818  auto nb = GeoModelMeshVerticesBase3D::nb_total_vertices();
819  nb +=
820  nb_entity_vertices( this->geomodel_, Region3D::type_name_static() );
821  return nb;
822  }
823 
825  {
826  auto count = GeoModelMeshVerticesBase3D::fill_vertices();
828  this->geomodel_, Region3D::type_name_static(), count );
829  return count;
830  }
831 
832  template <>
834  GeoModelMeshVerticesBase& geomodel_vertices,
835  const GeoModel3D& geomodel )
836  : geomodel_vertices_( geomodel_vertices ), geomodel_( geomodel )
837  {
838  vertex_maps_[Corner3D::type_name_static()] = &corner_vertex_maps_;
839  vertex_maps_[Line3D::type_name_static()] = &line_vertex_maps_;
840  vertex_maps_[Surface3D::type_name_static()] = &surface_vertex_maps_;
841  vertex_maps_[Region3D::type_name_static()] = &region_vertex_maps_;
842  }
843  /*******************************************************************************/
844 
845  template < index_t DIMENSION >
849  std::unique_ptr< VolumeMesh< DIMENSION > >& mesh )
850  : GeoModelMeshCommon< DIMENSION >( gmm, gm ), mesh_( mesh )
851  {
852  this->set_mesh( mesh_.get() );
853  }
854 
855  template < index_t DIMENSION >
857  {
858  return mesh_->nb_cells() > 0;
859  }
860 
861  template < index_t DIMENSION >
863  {
864  if( !is_initialized() )
865  {
866  const_cast< GeoModelMeshCells* >( this )->initialize();
867  }
868  }
869 
870  template < index_t DIMENSION >
872  {
873  this->gmm_.vertices.test_and_initialize();
874  auto mesh_builder =
876  if( mesh_->nb_vertices() != this->gmm_.vertices.nb() )
877  {
878  copy_vertices( mesh_builder.get(), *this->gmm_.vertices.mesh_ );
879  }
880 
881  region_cell_ptr_.resize(
882  this->geomodel_.nb_regions()
884  + 1,
885  0 );
886 
887  // Total number of cells
888  std::vector< index_t > nb_cells_per_type(
890  index_t nb{ 0 };
891 
892  for( const auto& region : this->geomodel_.regions() )
893  {
894  nb += region.nb_mesh_elements();
895  }
896 
897  // Get out if no cells
898  if( nb == 0 )
899  {
900  return;
901  }
902 
903  // Compute the number of cell per type and per region
904  for( const auto& region : this->geomodel_.regions() )
905  {
906  auto r = region.index();
907  for( auto c : range( region.nb_mesh_elements() ) )
908  {
909  CellType cur_cell_type = region.cell_type( c );
910  switch( cur_cell_type )
911  {
913  nb_cells_per_type[to_underlying_type(
916  * r
919  + 1]++;
920  break;
922  nb_cells_per_type[to_underlying_type(
927  break;
928  case CellType::PRISM:
929  nb_cells_per_type[to_underlying_type( CellType::PRISM )]++;
931  * r
933  + 1]++;
934  break;
935  case CellType::PYRAMID:
936  nb_cells_per_type[to_underlying_type(
937  CellType::PYRAMID )]++;
939  * r
941  + 1]++;
942  break;
944  nb_cells_per_type[to_underlying_type(
947  * r
950  + 1]++;
951  break;
952  default:
954  break;
955  }
956  }
957  }
958 
959  // Compute the cell offsets
960  std::vector< index_t > cells_offset_per_type(
962  for( auto t : range( to_underlying_type( CellType::TETRAHEDRON ) + 1,
964  {
965  cells_offset_per_type[t] += cells_offset_per_type[t - 1];
966  cells_offset_per_type[t] += nb_cells_per_type[t - 1];
967  }
968 
969  for( auto i : range( 1, region_cell_ptr_.size() - 1 ) )
970  {
971  region_cell_ptr_[i + 1] += region_cell_ptr_[i];
972  }
973 
974  // Create "empty" tet, hex, pyr and prism
975  for( auto i : range( to_underlying_type( CellType::UNDEFINED ) ) )
976  {
977  mesh_builder->create_cells( nb_cells_per_type[i], CellType( i ) );
978  }
979 
980  // Fill the cells with vertices
982  std::vector< index_t > cur_cell_per_type(
984  const auto& geomodel_vertices = this->gmm_.vertices;
985  for( const auto& region : this->geomodel_.regions() )
986  {
987  for( auto c : range( region.nb_mesh_elements() ) )
988  {
989  CellType cur_cell_type = region.cell_type( c );
990  index_t cur_cell =
991  cells_offset_per_type[to_underlying_type( cur_cell_type )]
992  + cur_cell_per_type[to_underlying_type( cur_cell_type )]++;
993  for( auto v : range( mesh_->nb_cell_vertices( cur_cell ) ) )
994  {
995  auto region_vertex_index =
996  region.mesh_element_vertex_index( { c, v } );
997  auto global_vertex_id =
998  geomodel_vertices.geomodel_vertex_id(
999  region.gmme(), region_vertex_index );
1000  mesh_builder->set_cell_vertex(
1001  ElementLocalVertex( cur_cell, v ), global_vertex_id );
1002  }
1003  region_id_[cur_cell] = region.index();
1004  cell_id_[cur_cell] = c;
1005  }
1006  }
1007 
1008  sort_cells();
1009 
1010  // Retrieve the adjacencies
1011  mesh_builder->connect_cells();
1012 
1013  // Cache some values
1014  nb_tets_ =
1015  nb_cells_per_type[to_underlying_type( CellType::TETRAHEDRON )];
1016  nb_hexs_ =
1017  nb_cells_per_type[to_underlying_type( CellType::HEXAHEDRON )];
1018  nb_prisms_ = nb_cells_per_type[to_underlying_type( CellType::PRISM )];
1019  nb_pyramids_ =
1020  nb_cells_per_type[to_underlying_type( CellType::PYRAMID )];
1021  nb_connectors_ =
1022  nb_cells_per_type[to_underlying_type( CellType::UNCLASSIFIED )];
1023  }
1024 
1025  template < index_t DIMENSION >
1027  {
1028  auto mesh_builder =
1030  std::vector< index_t > sorted_indices( mesh_->nb_cells() );
1031  std::iota( sorted_indices.begin(), sorted_indices.end(), 0 );
1032  GeoModelMeshCellsSort< DIMENSION > action( *mesh_, region_id_ );
1033  std::sort( sorted_indices.begin(), sorted_indices.end(), action );
1034  mesh_builder->permute_cells( sorted_indices );
1035 
1036  auto sorted_indices_geo =
1037  copy_std_vector_to_geo_vector( sorted_indices );
1038  GEO::Permutation::apply(
1039  region_id_.data(), sorted_indices_geo, sizeof( index_t ) );
1040  GEO::Permutation::apply(
1041  cell_id_.data(), sorted_indices_geo, sizeof( index_t ) );
1042  }
1043 
1044  template < index_t DIMENSION >
1046  {
1047  region_id_.resize( mesh_->nb_cells(), NO_ID );
1048  cell_id_.resize( mesh_->nb_cells(), NO_ID );
1049  }
1050 
1051  template < index_t DIMENSION >
1053  {
1054  region_id_.clear();
1055  cell_id_.clear();
1056  polygon_id_.clear();
1057  }
1058 
1059  template < index_t DIMENSION >
1061  {
1063  return mesh_->nb_cells();
1064  }
1065 
1066  template < index_t DIMENSION >
1068  {
1070  ringmesh_assert( c < mesh_->nb_cells() );
1071  return mesh_->nb_cell_vertices( c );
1072  }
1073 
1074  template < index_t DIMENSION >
1076  const ElementLocalVertex& cell_local_vertex ) const
1077  {
1079  ringmesh_assert( cell_local_vertex.element_id_ < mesh_->nb_cells() );
1081  cell_local_vertex.local_vertex_id_
1082  < mesh_->nb_cell_vertices( cell_local_vertex.element_id_ ) );
1083  return mesh_->cell_vertex( cell_local_vertex );
1084  }
1085 
1086  template < index_t DIMENSION >
1087  index_t GeoModelMeshCells< DIMENSION >::nb_edges( index_t c ) const
1088  {
1090  ringmesh_assert( c < mesh_->nb_cells() );
1091  return mesh_->nb_cell_edges( c );
1092  }
1093 
1094  template < index_t DIMENSION >
1096  {
1098  ringmesh_assert( c < mesh_->nb_cells() );
1099  return mesh_->nb_cell_facets( c );
1100  }
1101 
1102  template < index_t DIMENSION >
1104  const CellLocalFacet& cell_local_facet ) const
1105  {
1107  ringmesh_assert( cell_local_facet.cell_id_ < mesh_->nb_cells() );
1108  ringmesh_assert( cell_local_facet.local_facet_id_
1109  < mesh_->nb_cell_facets( cell_local_facet.cell_id_ ) );
1110  return mesh_->nb_cell_facet_vertices( cell_local_facet );
1111  }
1112 
1113  template < index_t DIMENSION >
1115  const CellLocalFacet& cell_local_facet, index_t lv ) const
1116  {
1118  ringmesh_assert( cell_local_facet.cell_id_ < mesh_->nb_cells() );
1119  ringmesh_assert( cell_local_facet.local_facet_id_
1120  < mesh_->nb_cell_facets( cell_local_facet.cell_id_ ) );
1121  return mesh_->cell_facet_vertex( cell_local_facet, lv );
1122  }
1123 
1124  template < index_t DIMENSION >
1126  index_t c, index_t le, index_t lv ) const
1127  {
1128  geo_debug_assert( le < nb_edges( c ) );
1129  geo_debug_assert( lv < 2 );
1130  return mesh_->cell_edge_vertex( c, le, lv );
1131  }
1132 
1133  template < index_t DIMENSION >
1135  index_t c, index_t f ) const
1136  {
1138  ringmesh_assert( c < mesh_->nb_cells() );
1139  ringmesh_assert( f < mesh_->nb_cell_facets( c ) );
1140  return mesh_->cell_adjacent( CellLocalFacet( c, f ) );
1141  }
1142 
1143  template < index_t DIMENSION >
1144  index_t GeoModelMeshCells< DIMENSION >::region( index_t c ) const
1145  {
1147  ringmesh_assert( c < mesh_->nb_cells() );
1148  return region_id_[c];
1149  }
1150 
1151  template < index_t DIMENSION >
1153  {
1155  ringmesh_assert( c < mesh_->nb_cells() );
1156  return cell_id_[c];
1157  }
1158 
1159  template < index_t DIMENSION >
1161  {
1163  ringmesh_assert( c < mesh_->nb_cells() );
1164  return mesh_->cell_type( c );
1165  }
1166 
1167  template < index_t DIMENSION >
1169  {
1171  switch( type )
1172  {
1173  case CellType::TETRAHEDRON:
1174  return nb_tet();
1175  case CellType::HEXAHEDRON:
1176  return nb_hex();
1177  case CellType::PRISM:
1178  return nb_prism();
1179  case CellType::PYRAMID:
1180  return nb_pyramid();
1182  return nb_connector();
1183  case CellType::UNDEFINED:
1184  return nb();
1185  default:
1187  return 0;
1188  }
1189  }
1190 
1191  template < index_t DIMENSION >
1193  index_t r, CellType type ) const
1194  {
1196  switch( type )
1197  {
1198  case CellType::TETRAHEDRON:
1199  return nb_tet( r );
1200  case CellType::HEXAHEDRON:
1201  return nb_hex( r );
1202  case CellType::PRISM:
1203  return nb_prism( r );
1204  case CellType::PYRAMID:
1205  return nb_pyramid( r );
1207  return nb_connector( r );
1208  case CellType::UNDEFINED:
1211  * ( r + 1 )]
1213  * r]
1214  == this->geomodel_.region( r ).nb_mesh_elements() );
1216  * ( r + 1 )]
1218  * r];
1219  default:
1221  return 0;
1222  }
1223  }
1224 
1225  template < index_t DIMENSION >
1227  index_t r, index_t c, CellType type ) const
1228  {
1230  switch( type )
1231  {
1232  case CellType::TETRAHEDRON:
1233  return tet( r, c );
1234  case CellType::HEXAHEDRON:
1235  return hex( r, c );
1236  case CellType::PRISM:
1237  return prism( r, c );
1238  case CellType::PYRAMID:
1239  return pyramid( r, c );
1241  return connector( r, c );
1242  case CellType::UNDEFINED:
1244  * r]
1245  + c;
1246  default:
1248  return 0;
1249  }
1250  }
1251 
1252  template < index_t DIMENSION >
1254  {
1256  return nb_tets_;
1257  }
1258 
1259  template < index_t DIMENSION >
1260  index_t GeoModelMeshCells< DIMENSION >::nb_tet( index_t r ) const
1261  {
1263  ringmesh_assert( r < this->geomodel_.nb_regions() );
1266  + 1 )]
1270  }
1271 
1272  template < index_t DIMENSION >
1273  index_t GeoModelMeshCells< DIMENSION >::tet( index_t r, index_t t ) const
1274  {
1276  ringmesh_assert( r < this->geomodel_.nb_regions() );
1279  + t;
1280  }
1281 
1282  template < index_t DIMENSION >
1284  {
1286  return nb_hexs_;
1287  }
1288 
1289  template < index_t DIMENSION >
1290  index_t GeoModelMeshCells< DIMENSION >::nb_hex( index_t r ) const
1291  {
1293  ringmesh_assert( r < this->geomodel_.nb_regions() );
1296  + 1 )]
1299  }
1300 
1301  template < index_t DIMENSION >
1302  index_t GeoModelMeshCells< DIMENSION >::hex( index_t r, index_t h ) const
1303  {
1305  ringmesh_assert( r < this->geomodel_.nb_regions() );
1308  + h;
1309  }
1310 
1311  template < index_t DIMENSION >
1313  {
1315  return nb_prisms_;
1316  }
1317 
1318  template < index_t DIMENSION >
1319  index_t GeoModelMeshCells< DIMENSION >::nb_prism( index_t r ) const
1320  {
1322  ringmesh_assert( r < this->geomodel_.nb_regions() );
1324  + ( to_underlying_type( CellType::PRISM ) + 1 )]
1327  }
1328 
1329  template < index_t DIMENSION >
1330  index_t GeoModelMeshCells< DIMENSION >::prism( index_t r, index_t p ) const
1331  {
1333  ringmesh_assert( r < this->geomodel_.nb_regions() );
1336  + p;
1337  }
1338 
1339  template < index_t DIMENSION >
1341  {
1343  return nb_pyramids_;
1344  }
1345 
1346  template < index_t DIMENSION >
1348  {
1350  ringmesh_assert( r < this->geomodel_.nb_regions() );
1353  + 1 )]
1356  }
1357 
1358  template < index_t DIMENSION >
1360  index_t r, index_t p ) const
1361  {
1363  ringmesh_assert( r < this->geomodel_.nb_regions() );
1366  + p;
1367  }
1368 
1369  template < index_t DIMENSION >
1371  {
1373  return nb_connectors_;
1374  }
1375 
1376  template < index_t DIMENSION >
1378  {
1380  ringmesh_assert( r < this->geomodel_.nb_regions() );
1383  + 1 )]
1387  }
1388 
1389  template < index_t DIMENSION >
1391  index_t r, index_t c ) const
1392  {
1394  ringmesh_assert( r < this->geomodel_.nb_regions() );
1397  + c;
1398  }
1399 
1400  template < index_t DIMENSION >
1402  {
1403  return mode_ == this->gmm_.duplicate_mode();
1404  }
1405 
1406  template < index_t DIMENSION >
1408  {
1410  {
1411  const_cast< GeoModelMeshCells< DIMENSION >* >( this )
1413  }
1414  }
1415 
1416  template < index_t DIMENSION >
1418  {
1420 
1422  std::vector< vec3 > corner_vertices(
1423  mesh_->cell_end( mesh_->nb_cells() - 1 ) );
1424  for( auto c : range( mesh_->nb_cells() ) )
1425  {
1426  index_t begin = mesh_->cell_begin( c );
1427  for( auto v : range( mesh_->nb_cell_vertices( c ) ) )
1428  {
1429  corner_vertices[begin + v] = mesh_->vertex(
1430  mesh_->cell_vertex( ElementLocalVertex( c, v ) ) );
1431  }
1432  }
1433 
1435  std::vector< ActionOnSurface > actions_on_surfaces(
1436  this->geomodel_.nb_surfaces(), SKIP );
1437  std::vector< bool > is_vertex_to_duplicate(
1438  corner_vertices.size(), false );
1439  {
1440  NNSearch< DIMENSION > nn_search( corner_vertices, false );
1441  for( const auto& surface : this->geomodel_.surfaces() )
1442  {
1443  if( !is_surface_to_duplicate( surface.index() ) )
1444  {
1445  continue;
1446  }
1447  actions_on_surfaces[surface.index()] = TO_PROCESS;
1448  for( auto v : range( surface.nb_vertices() ) )
1449  {
1450  auto colocated_corners = nn_search.get_neighbors(
1451  surface.vertex( v ), this->geomodel_.epsilon() );
1452  for( auto co : colocated_corners )
1453  {
1454  is_vertex_to_duplicate[co] = true;
1455  }
1456  }
1457  }
1458  }
1459  // Free some memory
1460  corner_vertices.clear();
1461 
1463  /* The goal is to visit the corners of the GeoModelMesh
1464  * that are on one side of a Surface. We propagate through the cells
1465  * that have one vertex on a Surface without crossing the Surface.
1466  * All the corners visited during this propagation around the vertex
1467  * are duplicated if needed.
1468  */
1469  this->gmm_.polygons.test_and_initialize();
1470  for( auto c : range( mesh_->nb_cells() ) )
1471  {
1472  for( auto v : range( mesh_->nb_cell_vertices( c ) ) )
1473  {
1474  // get the index of the corner inside cell_corners_
1475  index_t co = mesh_->cell_begin( c ) + v;
1476 
1477  if( !is_vertex_to_duplicate[co] )
1478  continue;
1479  // The vertex is on a surface to duplicate
1480 
1481  // Propagate on the cells around the corresponding vertex.
1482  // The propagation process cannot cross any surface.
1483  auto vertex_id =
1484  mesh_->cell_vertex( ElementLocalVertex( c, v ) );
1485 
1486  // all the cell corners resulting of the propagation
1487  std::vector< index_t > corner_used;
1488 
1489  // all cells used during the propagation, used to provide
1490  // adding the same cell several times into the stack
1491  std::vector< index_t > cell_added;
1492 
1493  // all the surfaces encountered during the propagation
1494  // and which side stopped the propagation
1495  std::vector< action_on_surface > surfaces;
1496 
1497  // stack of the front of cells
1498  std::stack< index_t > S;
1499  S.push( c );
1500  cell_added.push_back( c );
1501  do
1502  {
1503  auto cur_c = S.top();
1504  S.pop();
1505  // Find which corner of the current cell matches vertex_id
1506  auto cur_co = mesh_->find_cell_corner( cur_c, vertex_id );
1507  ringmesh_assert( cur_co != NO_ID );
1508  is_vertex_to_duplicate[cur_co] = false;
1509  corner_used.push_back( cur_co );
1510 
1511  // Find the cell facets including the vertex
1512  auto facets =
1513  cell_facets_around_vertex( *mesh_, cur_c, vertex_id );
1514  for( auto cur_f : facets )
1515  {
1516  // Find if the facet is on a surface or inside the
1517  // domain
1518  index_t polygon{ NO_ID };
1519  bool side;
1521  cur_c, cur_f, polygon, side ) )
1522  {
1523  auto surface_id =
1524  this->gmm_.polygons.surface( polygon );
1525  surfaces.emplace_back(
1526  surface_id, ActionOnSurface( side ) );
1527  }
1528  else
1529  {
1530  // The cell facet is not on a surface.
1531  // Add the adjacent cell to the stack if it exists
1532  // and has not already been processed or added into
1533  // the stack
1534  auto cur_adj =
1535  mesh_->cell_adjacent( { cur_c, cur_f } );
1536  if( cur_adj != NO_ID
1537  && !contains( cell_added, cur_adj ) )
1538  {
1539  cell_added.push_back( cur_adj );
1540  S.push( cur_adj );
1541  }
1542  }
1543  }
1544  } while( !S.empty() );
1545 
1546  // Remove redundant occurrences and sort the remaining ones
1547  sort_unique( surfaces );
1548 
1549  // Determine if the corners should be duplicated or not because
1550  // we need to duplicate only one side of the surface
1551  if( are_corners_to_duplicate( surfaces, actions_on_surfaces ) )
1552  {
1553  // Add a new duplicated vertex and its associated vertex
1554 
1555  /* @todo Review : Use the total_nb_vertices function [JP]
1556  * why mm_.vertices.nb_vertices() and not nb_vertices() ?
1557  * Please help the reader !! same thing 2 lines below [JP]
1558  */
1559  auto duplicated_vertex_id =
1560  this->gmm_.vertices.nb()
1561  + static_cast< index_t >(
1562  duplicated_vertex_indices_.size() );
1563  duplicated_vertex_indices_.push_back( vertex_id );
1564 
1565  // Update all the cell corners on this side of the surface
1566  // to the new duplicated vertex index
1567  auto mesh_builder =
1569  *mesh_ );
1570  for( auto cur_co : corner_used )
1571  {
1572  mesh_builder->set_cell_corner_vertex_index(
1573  cur_co, duplicated_vertex_id );
1574  }
1575  }
1576  }
1577  }
1578  mode_ = this->gmm_.duplicate_mode();
1579  }
1580 
1581  template < index_t DIMENSION >
1583  index_t c, index_t f, index_t& polygon, bool& side ) const
1584  {
1586  polygon = polygon_id_[mesh_->cell_facet( { c, f } )];
1587  if( polygon != NO_ID )
1588  {
1589  auto facet_normal = this->gmm_.polygons.normal( polygon );
1590  auto cell_facet_normal = mesh_->cell_facet_normal( { c, f } );
1591  side = dot( facet_normal, cell_facet_normal ) > 0;
1592  }
1593  return polygon != NO_ID;
1594  }
1595 
1596  template < index_t DIMENSION >
1598  const std::vector< action_on_surface >& surfaces,
1599  std::vector< ActionOnSurface >& info )
1600  {
1601  // Temporary vector, it is equal to surfaces
1602  std::vector< action_on_surface > temp_surfaces;
1603  temp_surfaces.reserve( temp_surfaces.size() );
1604 
1605  // Find if a free border was found, if so we encountered the
1606  // two sides of the surface during the propagation around a vertex
1607  for( auto i : range( 1, surfaces.size() ) )
1608  {
1609  if( surfaces[i - 1].first == surfaces[i].first )
1610  {
1611  // Found free border -> skip
1612  ringmesh_assert( surfaces[i - 1].second != surfaces[i].second );
1613  i++; // skip the action_on_surface (i-1) and i
1614  }
1615  else
1616  {
1617  temp_surfaces.push_back( surfaces[i - 1] );
1618  }
1619  }
1620  temp_surfaces.push_back( surfaces.back() );
1621 
1622  for( const auto& action : temp_surfaces )
1623  {
1624  auto s = action.first;
1625  switch( info[s] )
1626  {
1627  case SKIP:
1628  break;
1629  case TO_PROCESS:
1630  // First time we encounter this surface, do not duplicate
1631  // this side but wait to see if we encounter the other.
1632  // In the case of surfaces in the VOI, it is encountered only
1633  // once
1634  ringmesh_assert( action.second > TO_PROCESS );
1635  info[s] = ActionOnSurface( !action.second );
1636  break;
1637  default:
1638  // If the side matches -> duplicate
1639  if( info[s] == action.second )
1640  {
1641  return true;
1642  }
1643  }
1644  }
1645 
1646  return false;
1647  }
1648 
1649  template < index_t DIMENSION >
1651  index_t surface_id ) const
1652  {
1653  const auto& cur_surface = this->geomodel_.surface( surface_id );
1654  if( cur_surface.is_on_voi() )
1655  {
1656  return false;
1657  }
1658  switch( this->gmm_.duplicate_mode() )
1659  {
1660  case ALL:
1661  return true;
1662  case FAULT:
1663  {
1664  auto parent_interface = cur_surface.parent_gmge(
1666  if( parent_interface.is_defined() )
1667  {
1668  auto feature =
1669  this->geomodel_.geological_entity( parent_interface )
1670  .geological_feature();
1672  feature );
1673  }
1674  return false;
1675  }
1676  default:
1677  return false;
1678  }
1679  return false;
1680  }
1681 
1682  template < index_t DIMENSION >
1684  {
1686  return static_cast< index_t >( duplicated_vertex_indices_.size() );
1687  }
1688 
1689  template < index_t DIMENSION >
1691  {
1692  return nb_duplicated_vertices() + mesh_->nb_vertices();
1693  }
1694 
1695  template < index_t DIMENSION >
1697  const ElementLocalVertex& cell_local_vertex ) const
1698  {
1700  ringmesh_assert( cell_local_vertex.element_id_ < mesh_->nb_cells() );
1702  cell_local_vertex.local_vertex_id_
1703  < mesh_->nb_cell_vertices( cell_local_vertex.element_id_ ) );
1704  auto corner_value = mesh_->cell_vertex( cell_local_vertex );
1705  if( corner_value < mesh_->nb_vertices() )
1706  {
1707  return NO_ID;
1708  }
1709  return corner_value - mesh_->nb_vertices();
1710  }
1711 
1712  template < index_t DIMENSION >
1714  index_t duplicate_vertex_index ) const
1715  {
1718  duplicate_vertex_index < duplicated_vertex_indices_.size() );
1719  return duplicated_vertex_indices_[duplicate_vertex_index];
1720  }
1721 
1722  template < index_t DIMENSION >
1724  {
1725  auto mesh_builder =
1727  mesh_builder->clear( true, false );
1728  region_cell_ptr_.clear();
1729  nb_tets_ = 0;
1730  nb_hexs_ = 0;
1731  nb_prisms_ = 0;
1732  nb_pyramids_ = 0;
1733  nb_connectors_ = 0;
1734 
1735  mode_ = NONE;
1737  }
1738 
1739  template < index_t DIMENSION >
1741  {
1742  auto mesh_builder =
1744  for( auto c : range( mesh_->nb_cells() ) )
1745  {
1746  for( auto v : range( mesh_->nb_cell_vertices( c ) ) )
1747  {
1748  auto index = duplicated_corner_index( { c, v } );
1749  if( index != NO_ID )
1750  {
1751  mesh_builder->set_cell_corner_vertex_index(
1752  c, duplicated_vertex( index ) );
1753  }
1754  }
1755  }
1756 
1757  mode_ = NONE;
1759  }
1760 
1761  template < index_t DIMENSION >
1763  {
1764  if( polygon_id_.empty() )
1765  {
1766  const_cast< GeoModelMeshCells* >( this )->initialize_cell_facet();
1767  }
1768  }
1769 
1770  template < index_t DIMENSION >
1772  {
1773  this->gmm_.polygons.test_and_initialize();
1774 
1775  polygon_id_.resize( mesh_->nb_cell_facets(), NO_ID );
1776  const auto& nn_search = this->gmm_.polygons.nn_search();
1777  for( auto c : range( mesh_->nb_cells() ) )
1778  {
1779  for( auto f : range( mesh_->nb_cell_facets( c ) ) )
1780  {
1781  auto result = nn_search.get_neighbors(
1782  mesh_->cell_facet_barycenter( { c, f } ),
1783  this->geomodel_.epsilon() );
1784  if( !result.empty() )
1785  {
1786  polygon_id_[mesh_->cell_facet( { c, f } )] = result[0];
1787  // If there are more than 1 matching facet, this is WRONG
1788  // and the vertex indices should be checked too [Jeanne]
1789  ringmesh_assert( result.size() == 1 );
1790  }
1791  }
1792  }
1793  }
1794 
1795  template < index_t DIMENSION >
1797  {
1799  return mesh_->cell_barycenter( c );
1800  }
1801 
1802  template < index_t DIMENSION >
1803  double GeoModelMeshCells< DIMENSION >::volume( index_t c ) const
1804  {
1806  return mesh_->cell_volume( c );
1807  }
1808 
1809  template < index_t DIMENSION >
1812  {
1814  return mesh_->cell_aabb();
1815  }
1816 
1817  /*******************************************************************************/
1818  template < index_t DIMENSION >
1819  const std::string GeoModelMeshEdges< DIMENSION >::line_att_name = "line";
1820  template < index_t DIMENSION >
1822  "edge_line";
1823 
1824  template < index_t DIMENSION >
1828  std::unique_ptr< LineMesh< DIMENSION > >& mesh )
1829  : GeoModelMeshCommon< DIMENSION >( gmm, gm ), mesh_( mesh )
1830  {
1831  this->set_mesh( mesh_.get() );
1832  }
1833 
1834  template < index_t DIMENSION >
1836  {
1837  clear_edge_data();
1838  }
1839 
1840  template < index_t DIMENSION >
1842  {
1843  line_id_.resize( mesh_->nb_edges(), NO_ID );
1844  edge_id_.resize( mesh_->nb_edges(), NO_ID );
1845  }
1846 
1847  template < index_t DIMENSION >
1849  {
1850  line_id_.clear();
1851  edge_id_.clear();
1852  }
1853 
1854  template < index_t DIMENSION >
1856  {
1857  return mesh_->nb_edges() > 0;
1858  }
1859 
1860  template < index_t DIMENSION >
1862  {
1864  return mesh_->nb_edges();
1865  }
1866 
1867  template < index_t DIMENSION >
1869  const ElementLocalVertex& edge_local_vertex ) const
1870  {
1872  ringmesh_assert( edge_local_vertex.element_id_ < mesh_->nb_edges() );
1873  ringmesh_assert( edge_local_vertex.local_vertex_id_ < 2 );
1874  return mesh_->edge_vertex( edge_local_vertex );
1875  }
1876 
1877  template < index_t DIMENSION >
1878  index_t GeoModelMeshEdges< DIMENSION >::line( index_t e ) const
1879  {
1881  ringmesh_assert( e < mesh_->nb_edges() );
1882  return line_id_[e];
1883  }
1884 
1885  template < index_t DIMENSION >
1887  {
1889  ringmesh_assert( e < mesh_->nb_edges() );
1890  return edge_id_[e];
1891  }
1892 
1893  template < index_t DIMENSION >
1894  index_t GeoModelMeshEdges< DIMENSION >::nb_edges( index_t l ) const
1895  {
1897  ringmesh_assert( l < this->geomodel_.nb_lines() );
1898  return line_edge_ptr_[l + 1];
1899  }
1900 
1901  template < index_t DIMENSION >
1902  index_t GeoModelMeshEdges< DIMENSION >::edge( index_t l, index_t e ) const
1903  {
1905  ringmesh_assert( l < this->geomodel_.nb_lines() );
1906  return line_edge_ptr_[l] + e;
1907  }
1908 
1909  template < index_t DIMENSION >
1911  {
1912  line_edge_ptr_.clear();
1913  nb_edges_ = 0;
1914  auto mesh_builder =
1916  mesh_builder->clear( true, false );
1917  }
1918 
1919  template < index_t DIMENSION >
1921  {
1922  if( !is_initialized() )
1923  {
1924  const_cast< GeoModelMeshEdges* >( this )->initialize();
1925  }
1926  }
1927 
1928  template < index_t DIMENSION >
1930  {
1931  this->gmm_.vertices.test_and_initialize();
1932  line_edge_ptr_.resize( this->geomodel_.nb_lines() + 1, 0 );
1933  auto mesh_builder =
1935  if( mesh_->nb_vertices() != this->gmm_.vertices.nb() )
1936  {
1937  copy_vertices( mesh_builder.get(), *this->gmm_.vertices.mesh_ );
1938  }
1939 
1940  // Compute the total number of edges per line
1941  for( auto l : range( this->geomodel_.nb_lines() ) )
1942  {
1943  const Line< DIMENSION >& line = this->geomodel_.line( l );
1944  line_edge_ptr_[l + 1] = line_edge_ptr_[l] + line.nb_mesh_elements();
1945  nb_edges_ += line.nb_mesh_elements();
1946  }
1947 
1948  // Create edges
1949  mesh_builder->create_edges( nb_edges_ );
1950  resize_edge_data();
1951  const auto& geomodel_vertices = this->gmm_.vertices;
1952  index_t cur_edge{ 0 };
1953  for( auto l : range( this->geomodel_.nb_lines() ) )
1954  {
1955  const auto& line = this->geomodel_.line( l );
1956  auto line_id = line.gmme();
1957  for( auto e : range( line.nb_mesh_elements() ) )
1958  {
1959  for( auto v : range( 2 ) )
1960  {
1961  auto v_id = geomodel_vertices.geomodel_vertex_id(
1962  line_id, ElementLocalVertex( e, v ) );
1963  ringmesh_assert( v_id != NO_ID );
1964  mesh_builder->set_edge_vertex(
1965  ElementLocalVertex( cur_edge, v ), v_id );
1966  }
1967  line_id_[cur_edge] = l;
1968  edge_id_[cur_edge] = e;
1969  cur_edge++;
1970  }
1971  }
1972  }
1973 
1974  template < index_t DIMENSION >
1976  {
1978  return mesh_->edge_barycenter( e );
1979  }
1980 
1981  template < index_t DIMENSION >
1982  double GeoModelMeshEdges< DIMENSION >::length( index_t e ) const
1983  {
1985  return mesh_->edge_length( e );
1986  }
1987 
1988  template < index_t DIMENSION >
1991  {
1993  return mesh_->edge_aabb();
1994  }
1995 
1996  /*******************************************************************************/
1997  template < index_t DIMENSION >
1999  "surface";
2000  template < index_t DIMENSION >
2001  const std::string
2003  "polygon_surface";
2004 
2005  template < index_t DIMENSION >
2009  std::unique_ptr< SurfaceMesh< DIMENSION > >& mesh )
2010  : GeoModelMeshCommon< DIMENSION >( gmm, gm ), mesh_( mesh )
2011  {
2012  this->set_mesh( mesh_.get() );
2013  }
2014 
2015  template < index_t DIMENSION >
2017  {
2018  clear_polygones_data();
2019  }
2020 
2021  template < index_t DIMENSION >
2023  {
2024  surface_id_.resize( mesh_->nb_polygons(), NO_ID );
2025  polygon_id_.resize( mesh_->nb_polygons(), NO_ID );
2026  }
2027 
2028  template < index_t DIMENSION >
2030  {
2031  surface_id_.clear();
2032  polygon_id_.clear();
2033  }
2034 
2035  template < index_t DIMENSION >
2037  {
2038  return mesh_->nb_polygons() > 0;
2039  }
2040 
2041  template < index_t DIMENSION >
2043  {
2044  test_and_initialize();
2045  return mesh_->nb_polygons();
2046  }
2047 
2048  template < index_t DIMENSION >
2050  index_t p ) const
2051  {
2052  test_and_initialize();
2053  ringmesh_assert( p < mesh_->nb_polygons() );
2054  return mesh_->nb_polygon_vertices( p );
2055  }
2056 
2057  template < index_t DIMENSION >
2059  const ElementLocalVertex& polygon_local_vertex ) const
2060  {
2061  test_and_initialize();
2063  polygon_local_vertex.element_id_ < mesh_->nb_polygons() );
2065  polygon_local_vertex.local_vertex_id_
2066  < mesh_->nb_polygon_vertices( polygon_local_vertex.element_id_ ) );
2067  return mesh_->polygon_vertex( polygon_local_vertex );
2068  }
2069 
2070  template < index_t DIMENSION >
2072  const PolygonLocalEdge& polygon_local_edge ) const
2073  {
2074  test_and_initialize();
2076  polygon_local_edge.polygon_id_ < mesh_->nb_polygons() );
2078  polygon_local_edge.local_edge_id_
2079  < mesh_->nb_polygon_vertices( polygon_local_edge.polygon_id_ ) );
2080  return mesh_->polygon_adjacent( polygon_local_edge );
2081  }
2082 
2083  template < index_t DIMENSION >
2085  {
2086  test_and_initialize();
2087  ringmesh_assert( p < mesh_->nb_polygons() );
2088  return surface_id_[p];
2089  }
2090 
2091  template < index_t DIMENSION >
2093  index_t p ) const
2094  {
2095  test_and_initialize();
2096  ringmesh_assert( p < mesh_->nb_polygons() );
2097  return polygon_id_[p];
2098  }
2099 
2100  template < index_t DIMENSION >
2101  std::tuple< PolygonType, index_t >
2103  {
2104  test_and_initialize();
2105  ringmesh_assert( p < mesh_->nb_polygons() );
2106  auto polygon = index_in_surface( p );
2107  auto s = surface( p );
2108  for( auto t : range( to_underlying_type( PolygonType::TRIANGLE ),
2109  to_underlying_type( PolygonType::UNDEFINED ) ) )
2110  {
2111  auto T = static_cast< PolygonType >( t );
2112  if( polygon < nb_polygons( s, T ) )
2113  {
2114  return std::make_tuple( T, polygon );
2115  }
2116  polygon -= nb_polygons( s, T );
2117  }
2119  return std::make_tuple( PolygonType::UNDEFINED, NO_ID );
2120  }
2121 
2122  template < index_t DIMENSION >
2124  PolygonType type ) const
2125  {
2126  test_and_initialize();
2127  switch( type )
2128  {
2129  case PolygonType::TRIANGLE:
2130  return nb_triangle();
2131  case PolygonType::QUAD:
2132  return nb_quad();
2133  case PolygonType::UNCLASSIFIED:
2134  return nb_unclassified_polygon();
2135  case PolygonType::UNDEFINED:
2136  return nb();
2137  default:
2139  return 0;
2140  }
2141  }
2142 
2143  template < index_t DIMENSION >
2145  index_t s, PolygonType type ) const
2146  {
2147  test_and_initialize();
2148  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2149  switch( type )
2150  {
2151  case PolygonType::TRIANGLE:
2152  return nb_triangle( s );
2153  case PolygonType::QUAD:
2154  return nb_quad( s );
2155  case PolygonType::UNCLASSIFIED:
2156  return nb_unclassified_polygon( s );
2157  case PolygonType::UNDEFINED:
2158  return surface_polygon_ptr_[to_underlying_type(
2159  PolygonType::UNDEFINED )
2160  * ( s + 1 )]
2161  - surface_polygon_ptr_[to_underlying_type(
2162  PolygonType::UNDEFINED )
2163  * s];
2164  default:
2166  return 0;
2167  }
2168  }
2169 
2170  template < index_t DIMENSION >
2172  index_t s, index_t p, PolygonType type ) const
2173  {
2174  test_and_initialize();
2175  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2176  switch( type )
2177  {
2178  case PolygonType::TRIANGLE:
2179  return triangle( s, p );
2180  case PolygonType::QUAD:
2181  return quad( s, p );
2182  case PolygonType::UNCLASSIFIED:
2183  return unclassified_polygon( s, p );
2184  case PolygonType::UNDEFINED:
2185  return surface_polygon_ptr_[to_underlying_type(
2186  PolygonType::UNDEFINED )
2187  * s]
2188  + p;
2189  default:
2191  return 0;
2192  }
2193  }
2194 
2195  template < index_t DIMENSION >
2197  {
2198  test_and_initialize();
2199  return nb_triangles_;
2200  }
2201 
2202  template < index_t DIMENSION >
2204  index_t s ) const
2205  {
2206  test_and_initialize();
2207  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2208  return surface_polygon_ptr_
2209  [to_underlying_type( PolygonType::UNDEFINED ) * s
2210  + ( to_underlying_type( PolygonType::TRIANGLE ) + 1 )]
2211  - surface_polygon_ptr_
2212  [to_underlying_type( PolygonType::UNDEFINED ) * s
2213  + to_underlying_type( PolygonType::TRIANGLE )];
2214  }
2215 
2216  template < index_t DIMENSION >
2218  index_t s, index_t t ) const
2219  {
2220  test_and_initialize();
2221  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2222  return surface_polygon_ptr_[to_underlying_type( PolygonType::UNDEFINED )
2223  * s
2225  PolygonType::TRIANGLE )]
2226  + t;
2227  }
2228 
2229  template < index_t DIMENSION >
2231  {
2232  test_and_initialize();
2233  return nb_quads_;
2234  }
2235 
2236  template < index_t DIMENSION >
2238  {
2239  test_and_initialize();
2240  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2241  return surface_polygon_ptr_[to_underlying_type( PolygonType::UNDEFINED )
2242  * s
2243  + ( to_underlying_type( PolygonType::QUAD )
2244  + 1 )]
2245  - surface_polygon_ptr_
2246  [to_underlying_type( PolygonType::UNDEFINED ) * s
2247  + to_underlying_type( PolygonType::QUAD )];
2248  }
2249 
2250  template < index_t DIMENSION >
2252  index_t s, index_t q ) const
2253  {
2254  test_and_initialize();
2255  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2256  return surface_polygon_ptr_[to_underlying_type( PolygonType::UNDEFINED )
2257  * s
2258  + to_underlying_type( PolygonType::QUAD )]
2259  + q;
2260  }
2261 
2262  template < index_t DIMENSION >
2263  index_t
2265  {
2266  test_and_initialize();
2267  return nb_unclassified_polygons_;
2268  }
2269 
2270  template < index_t DIMENSION >
2272  index_t s ) const
2273  {
2274  test_and_initialize();
2275  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2276  return surface_polygon_ptr_
2277  [to_underlying_type( PolygonType::UNDEFINED ) * s
2278  + ( to_underlying_type( PolygonType::UNCLASSIFIED )
2279  + 1 )]
2280  - surface_polygon_ptr_
2281  [to_underlying_type( PolygonType::UNDEFINED ) * s
2282  + to_underlying_type( PolygonType::UNCLASSIFIED )];
2283  }
2284 
2285  template < index_t DIMENSION >
2287  index_t s, index_t p ) const
2288  {
2289  test_and_initialize();
2290  ringmesh_assert( s < this->geomodel_.nb_surfaces() );
2291  return surface_polygon_ptr_[to_underlying_type( PolygonType::UNDEFINED )
2292  * s
2294  PolygonType::UNCLASSIFIED )]
2295  + p;
2296  }
2297 
2298  template < index_t DIMENSION >
2300  {
2301  surface_polygon_ptr_.clear();
2302  nb_triangles_ = 0;
2303  nb_quads_ = 0;
2304  auto mesh_builder =
2306  mesh_builder->clear( true, false );
2307  }
2308 
2309  template < index_t DIMENSION >
2311  {
2312  if( !is_initialized() )
2313  {
2314  const_cast< GeoModelMeshPolygonsBase* >( this )->initialize();
2315  }
2316  }
2317 
2318  template < index_t DIMENSION >
2320  {
2321  this->gmm_.vertices.test_and_initialize();
2322  surface_polygon_ptr_.resize(
2323  this->geomodel_.nb_surfaces()
2324  * to_underlying_type( PolygonType::UNDEFINED )
2325  + 1,
2326  0 );
2327  auto mesh_builder =
2329  if( mesh_->nb_vertices() != this->gmm_.vertices.nb() )
2330  {
2331  copy_vertices( mesh_builder.get(), *this->gmm_.vertices.mesh_ );
2332  }
2333 
2334  // Compute the total number of polygons per type and per surface
2335  std::map< PolygonType, index_t > nb_polygon_per_type = {
2336  { PolygonType::TRIANGLE, 0 }, { PolygonType::QUAD, 0 },
2337  { PolygonType::UNCLASSIFIED, 0 }
2338  };
2339  for( auto s : range( this->geomodel_.nb_surfaces() ) )
2340  {
2341  const auto& surface = this->geomodel_.surface( s );
2342  if( surface.is_simplicial() )
2343  {
2344  nb_polygon_per_type[PolygonType::TRIANGLE] +=
2345  surface.nb_mesh_elements();
2346  surface_polygon_ptr_
2347  [to_underlying_type( PolygonType::UNDEFINED ) * s
2348  + to_underlying_type( PolygonType::TRIANGLE ) + 1] +=
2349  surface.nb_mesh_elements();
2350  }
2351  else
2352  {
2353  for( auto p : range( surface.nb_mesh_elements() ) )
2354  {
2355  switch( surface.nb_mesh_element_vertices( p ) )
2356  {
2357  case 3:
2358  nb_polygon_per_type[PolygonType::TRIANGLE]++;
2359  surface_polygon_ptr_
2360  [to_underlying_type( PolygonType::UNDEFINED ) * s
2361  + to_underlying_type( PolygonType::TRIANGLE )
2362  + 1]++;
2363  break;
2364  case 4:
2365  nb_polygon_per_type[PolygonType::QUAD]++;
2366  surface_polygon_ptr_
2367  [to_underlying_type( PolygonType::UNDEFINED ) * s
2368  + to_underlying_type( PolygonType::QUAD )
2369  + 1]++;
2370  break;
2371  default:
2372  nb_polygon_per_type[PolygonType::UNCLASSIFIED]++;
2373  surface_polygon_ptr_[to_underlying_type(
2374  PolygonType::UNDEFINED )
2375  * s
2377  PolygonType::UNCLASSIFIED )
2378  + 1]++;
2379  break;
2380  }
2381  }
2382  }
2383  }
2384 
2385  // Get out if no polygons
2386  auto nb_total_polygons =
2387  nb_polygon_per_type[PolygonType::TRIANGLE]
2388 
2389  + nb_polygon_per_type[PolygonType::QUAD]
2390  + nb_polygon_per_type[PolygonType::UNCLASSIFIED];
2391 
2392  if( nb_total_polygons == 0 )
2393  {
2394  return;
2395  }
2396 
2397  // Create triangles and quads, the polygons will be handle later
2398  if( nb_polygon_per_type[PolygonType::TRIANGLE] )
2399  {
2400  mesh_builder->create_triangles(
2401  nb_polygon_per_type[PolygonType::TRIANGLE] );
2402  }
2403  if( nb_polygon_per_type[PolygonType::QUAD] )
2404  {
2405  mesh_builder->create_quads(
2406  nb_polygon_per_type[PolygonType::QUAD] );
2407  }
2408 
2409  // Compute the polygon offset
2410  std::map< PolygonType, index_t > polygon_offset_per_type = {
2411  { PolygonType::TRIANGLE, 0 }, { PolygonType::QUAD, 0 },
2412  { PolygonType::UNCLASSIFIED, 0 }
2413  };
2414  polygon_offset_per_type[PolygonType::QUAD] =
2415  nb_polygon_per_type[PolygonType::TRIANGLE];
2416  polygon_offset_per_type[PolygonType::UNCLASSIFIED] =
2417  nb_polygon_per_type[PolygonType::TRIANGLE]
2418  + nb_polygon_per_type[PolygonType::QUAD];
2419 
2420  for( auto i : range( 1, surface_polygon_ptr_.size() - 1 ) )
2421  {
2422  surface_polygon_ptr_[i + 1] += surface_polygon_ptr_[i];
2423  }
2424 
2425  // Fill the triangles and quads created above
2426  // Create and fill polygons
2427  resize_polygones_data();
2428  const auto& geomodel_vertices = this->gmm_.vertices;
2429  std::vector< index_t > cur_polygon_per_type(
2430  to_underlying_type( PolygonType::UNDEFINED ), 0 );
2431  for( auto s : range( this->geomodel_.nb_surfaces() ) )
2432  {
2433  const auto& surface = this->geomodel_.surface( s );
2434  auto surface_id = surface.gmme();
2435  for( auto p : range( surface.nb_mesh_elements() ) )
2436  {
2437  auto nb_vertices = surface.nb_mesh_element_vertices( p );
2438  index_t cur_polygon{ NO_ID };
2439  if( nb_vertices < 5 )
2440  {
2441  auto T = static_cast< PolygonType >( nb_vertices - 3 );
2442  cur_polygon =
2443  polygon_offset_per_type[T]
2444  + cur_polygon_per_type[to_underlying_type( T )]++;
2445  for( auto v : range( nb_vertices ) )
2446  {
2447  auto v_id = geomodel_vertices.geomodel_vertex_id(
2448  surface_id, ElementLocalVertex( p, v ) );
2449  ringmesh_assert( v_id != NO_ID );
2450  mesh_builder->set_polygon_vertex(
2451  ElementLocalVertex( cur_polygon, v ), v_id );
2452  }
2453  }
2454  else
2455  {
2456  std::vector< index_t > vertices( nb_vertices );
2457  for( auto v : range( nb_vertices ) )
2458  {
2459  vertices[v] = geomodel_vertices.geomodel_vertex_id(
2460  surface_id, ElementLocalVertex( p, v ) );
2461  }
2462  cur_polygon = mesh_builder->create_polygon( vertices );
2463  }
2464  surface_id_[cur_polygon] = s;
2465  polygon_id_[cur_polygon] = p;
2466  }
2467  }
2468 
2469  sort_polygons();
2470 
2471  // Compute polygon adjacencies
2472  mesh_builder->connect_polygons();
2473  disconnect_along_lines();
2474 
2475  // Cache some values
2476  nb_triangles_ = nb_polygon_per_type[PolygonType::TRIANGLE];
2477  nb_quads_ = nb_polygon_per_type[PolygonType::QUAD];
2478  nb_unclassified_polygons_ =
2479  nb_polygon_per_type[PolygonType::UNCLASSIFIED];
2480  }
2481 
2482  template < index_t DIMENSION >
2484  {
2485  auto mesh_builder =
2487  std::vector< index_t > sorted_indices( mesh_->nb_polygons() );
2488  std::iota( sorted_indices.begin(), sorted_indices.end(), 0 );
2489  GeoModelMeshPolygonsBaseSort< DIMENSION > action( *mesh_, surface_id_ );
2490  std::sort( sorted_indices.begin(), sorted_indices.end(), action );
2491  mesh_builder->permute_polygons( sorted_indices );
2492 
2493  auto sorted_indices_geo =
2494  copy_std_vector_to_geo_vector( sorted_indices );
2495  GEO::Permutation::apply(
2496  surface_id_.data(), sorted_indices_geo, sizeof( index_t ) );
2497  GEO::Permutation::apply(
2498  polygon_id_.data(), sorted_indices_geo, sizeof( index_t ) );
2499  }
2500 
2501  template < index_t DIMENSION >
2503  {
2504  auto mesh_builder =
2506  for( auto s : range( this->geomodel_.nb_surfaces() ) )
2507  {
2508  const auto& surface = this->geomodel_.surface( s );
2509  for( auto p : range( nb_polygons( s ) ) )
2510  {
2511  auto polygon_id = polygon( s, p );
2512  auto surface_polygon_id = index_in_surface( polygon_id );
2513  for( auto v : range( nb_vertices( polygon_id ) ) )
2514  {
2515  auto adj = surface.polygon_adjacent_index(
2516  PolygonLocalEdge( surface_polygon_id, v ) );
2517  if( adj == NO_ID )
2518  {
2519  mesh_builder->set_polygon_adjacent(
2520  ElementLocalVertex( polygon_id, v ), NO_ID );
2521  }
2522  }
2523  }
2524  }
2525  }
2526 
2527  template < index_t DIMENSION >
2529  index_t p ) const
2530  {
2531  test_and_initialize();
2532  return mesh_->polygon_barycenter( p );
2533  }
2534 
2535  template < index_t DIMENSION >
2537  {
2538  test_and_initialize();
2539  return mesh_->polygon_area( p );
2540  }
2541 
2542  template < index_t DIMENSION >
2545  {
2546  test_and_initialize();
2547  return mesh_->polygon_aabb();
2548  }
2549 
2550  template < index_t DIMENSION >
2554  std::unique_ptr< SurfaceMesh< DIMENSION > >& mesh )
2555  : GeoModelMeshPolygonsBase< DIMENSION >( gmm, gm, mesh )
2556  {
2557  }
2558 
2560  GeoModel3D& gm,
2561  std::unique_ptr< SurfaceMesh3D >& mesh )
2562  : GeoModelMeshPolygonsBase< 3 >( gmm, gm, mesh )
2563  {
2564  }
2565 
2567  {
2569  return mesh_->polygon_normal( p );
2570  }
2571 
2572  /*******************************************************************************/
2573 
2574  template < index_t DIMENSION >
2578  std::unique_ptr< LineMesh< DIMENSION > >& mesh )
2579  : GeoModelMeshCommon< DIMENSION >( gmm, gm ), mesh_( mesh )
2580  {
2581  this->set_mesh( mesh_.get() );
2582  }
2583 
2584  template < index_t DIMENSION >
2586  {
2588  return this->geomodel_.wells() ? this->geomodel_.wells()->nb_wells()
2589  : 0;
2590  }
2591 
2592  template < index_t DIMENSION >
2594  {
2596  return mesh_->nb_edges();
2597  }
2598 
2599  template < index_t DIMENSION >
2600  index_t GeoModelMeshWells< DIMENSION >::nb_edges( index_t w ) const
2601  {
2603  return well_ptr_[w + 1] - well_ptr_[w];
2604  }
2605 
2606  template < index_t DIMENSION >
2608  index_t w, index_t e, index_t v ) const
2609  {
2611  return mesh_->edge_vertex( ElementLocalVertex( well_ptr_[w] + e, v ) );
2612  }
2613 
2614  template < index_t DIMENSION >
2616  {
2617  auto mesh_builder =
2619  mesh_builder->clear( true, false );
2620  well_ptr_.clear();
2621  }
2622 
2623  template < index_t DIMENSION >
2625  {
2626  return mesh_->nb_edges() > 0;
2627  }
2628 
2629  template < index_t DIMENSION >
2631  {
2632  if( !is_initialized() )
2633  {
2634  const_cast< GeoModelMeshWells* >( this )->initialize();
2635  }
2636  }
2637 
2638  template < index_t DIMENSION >
2640  {
2641  if( !this->geomodel_.wells() )
2642  {
2643  return;
2644  }
2645  this->gmm_.vertices.test_and_initialize();
2646  auto mesh_builder =
2648  if( mesh_->nb_vertices() != this->gmm_.vertices.nb() )
2649  {
2650  copy_vertices( mesh_builder.get(), *this->gmm_.vertices.mesh_ );
2651  }
2652 
2653  // Compute the total number of edge per well
2654  const auto& wells = *this->geomodel_.wells();
2655  well_ptr_.resize( wells.nb_wells() + 1, 0 );
2656  index_t nb_edges{ 0 };
2657  for( auto w : range( wells.nb_wells() ) )
2658  {
2659  nb_edges += wells.well( w ).nb_edges();
2660  well_ptr_[w + 1] = nb_edges;
2661  }
2662 
2663  // Compute the edge offset
2664  for( auto i : range( 1, well_ptr_.size() - 1 ) )
2665  {
2666  well_ptr_[i + 1] += well_ptr_[i];
2667  }
2668 
2669  // Create edges
2670  mesh_builder->create_edges( well_ptr_.back() );
2671 
2672  // Fill edges
2673  index_t cur_edge{ 0 };
2674  for( auto w : range( 0, wells.nb_wells() ) )
2675  {
2676  const Well< DIMENSION >& well = wells.well( w );
2677  for( auto p : range( well.nb_parts() ) )
2678  {
2679  for( auto e : range( well.part( p ).nb_edges() ) )
2680  {
2681  const auto& e0 = well.part( p ).edge_vertex( { e, 0 } );
2682  mesh_builder->set_edge_vertex(
2683  ElementLocalVertex( cur_edge, 0 ),
2684  this->gmm_.vertices.index( e0 ) );
2685  const auto& e1 = well.part( p ).edge_vertex( { e, 1 } );
2686  mesh_builder->set_edge_vertex(
2687  ElementLocalVertex( cur_edge, 1 ),
2688  this->gmm_.vertices.index( e1 ) );
2689  cur_edge++;
2690  }
2691  }
2692  }
2693  }
2694 
2695  template < index_t DIMENSION >
2698  {
2700  return mesh_->edge_aabb();
2701  }
2702 
2703  /*******************************************************************************/
2704 
2705  template < index_t DIMENSION >
2708  : geomodel_( geomodel ),
2709  vertices( gmm, geomodel, mesh_set_.point_set_mesh ),
2710  edges( gmm, geomodel, mesh_set_.line_mesh ),
2711  wells( gmm, geomodel, mesh_set_.well_mesh ),
2712  polygons( gmm, geomodel, mesh_set_.surface_mesh )
2713  {
2714  }
2715 
2716  template < index_t DIMENSION >
2718  {
2719  polygons.clear_polygones_data();
2720  }
2721 
2722  template < index_t DIMENSION >
2724  {
2725  vertices.remove_colocated();
2726  }
2727 
2728  template < index_t DIMENSION >
2730  std::vector< index_t >& to_delete )
2731  {
2732  vertices.erase_vertices( to_delete );
2733  }
2734 
2735  template < index_t DIMENSION >
2737  const MeshType& type )
2738  {
2739  if( mesh_set_.point_set_mesh->type_name() != type )
2740  {
2741  vertices.clear();
2742  mesh_set_.create_point_set_mesh( type );
2743  }
2744  }
2745 
2746  template < index_t DIMENSION >
2748  const MeshType& type )
2749  {
2750  if( mesh_set_.line_mesh->type_name() != type )
2751  {
2752  wells.clear();
2753  mesh_set_.create_line_mesh( type );
2754  }
2755  }
2756 
2757  template < index_t DIMENSION >
2759  const MeshType& type )
2760  {
2761  if( mesh_set_.surface_mesh->type_name() != type )
2762  {
2763  polygons.clear();
2764  mesh_set_.create_surface_mesh( type );
2765  }
2766  }
2767 
2768  template < index_t DIMENSION >
2770  : GeoModelMeshBase< DIMENSION >( *this, geomodel )
2771  {
2772  }
2773 
2775  : GeoModelMeshBase< 3 >( *this, geomodel ),
2776  cells( *this, geomodel, mesh_set_.volume_mesh )
2777  {
2778  }
2779 
2781  {
2782  }
2783 
2785  const MeshType& type )
2786  {
2787  if( mesh_set_.volume_mesh->type_name() != type )
2788  {
2789  cells.clear();
2790  mesh_set_.create_volume_mesh( type );
2791  }
2792  }
2793 
2795  {
2796  transfer_vertex_attributes_from_gmm_to_gm_regions();
2797  transfer_cell_attributes_from_gmm_to_gm_regions();
2798  }
2799 
2801  {
2802  transfer_vertex_attributes_from_gm_regions_to_gmm();
2803  transfer_cell_attributes_from_gm_regions_to_gmm();
2804  }
2805 
2807  const
2808  {
2809  auto& gmm_v_attr_mgr = vertices.attribute_manager();
2810  GEO::vector< std::string > att_v_names;
2811  gmm_v_attr_mgr.list_attribute_names( att_v_names );
2812  for( const auto& cur_attr_name : att_v_names )
2813  {
2814  // It is not necessary to copy the coordinates. There are already
2815  // there.
2816  if( cur_attr_name == "point" )
2817  {
2818  continue;
2819  }
2820 
2821  auto* cur_v_att_store_in_gmm =
2822  gmm_v_attr_mgr.find_attribute_store( cur_attr_name );
2823  ringmesh_assert( cur_v_att_store_in_gmm != nullptr );
2824  auto dim = cur_v_att_store_in_gmm->dimension();
2825 
2826  for( auto v : range( vertices.nb() ) )
2827  {
2828  auto vertices_on_geomodel_region = vertices.gme_type_vertices(
2829  Region3D::type_name_static(), v );
2830  for( const auto& cur_vertex_on_geomodel :
2831  vertices_on_geomodel_region )
2832  {
2833  const auto& cur_region =
2834  geomodel_.region( cur_vertex_on_geomodel.gmme.index() );
2835  auto& reg_v_attr_mgr =
2836  cur_region.vertex_attribute_manager();
2837  GEO::AttributeStore* cur_v_att_store_in_reg{ nullptr };
2838 
2839  if( !reg_v_attr_mgr.is_defined( cur_attr_name ) )
2840  {
2841  const auto cur_type_name = GEO::AttributeStore::
2842  element_type_name_by_element_typeid_name(
2843  cur_v_att_store_in_gmm->element_typeid_name() );
2845  GEO::AttributeStore::element_type_name_is_known(
2846  cur_type_name ) );
2847  cur_v_att_store_in_reg = GEO::AttributeStore::
2848  create_attribute_store_by_element_type_name(
2849  cur_type_name, dim );
2850  reg_v_attr_mgr.bind_attribute_store(
2851  cur_attr_name, cur_v_att_store_in_reg );
2852  }
2853  else
2854  {
2855  cur_v_att_store_in_reg =
2856  reg_v_attr_mgr.find_attribute_store(
2857  cur_attr_name );
2858  }
2859  ringmesh_assert( cur_v_att_store_in_reg != nullptr );
2861  cur_v_att_store_in_reg->element_size()
2862  == cur_v_att_store_in_gmm->element_size() );
2863 
2864  GEO::Memory::copy(
2865  static_cast< GEO::Memory::pointer >(
2866  cur_v_att_store_in_reg->data() )
2867  + cur_vertex_on_geomodel.v_index * dim
2868  * cur_v_att_store_in_reg->element_size(),
2869  static_cast< GEO::Memory::pointer >(
2870  cur_v_att_store_in_gmm->data() )
2871  + v * dim * cur_v_att_store_in_gmm->element_size(),
2872  dim * cur_v_att_store_in_reg->element_size() );
2873  }
2874  }
2875  }
2876  }
2877 
2879  const
2880  {
2881  for( const auto& cur_reg : geomodel().regions() )
2882  {
2883  GEO::vector< std::string > att_v_names;
2884  GEO::AttributesManager& reg_vertex_attr_mgr =
2885  cur_reg.vertex_attribute_manager();
2886  reg_vertex_attr_mgr.list_attribute_names( att_v_names );
2887  for( const auto& cur_attr_name : att_v_names )
2888  {
2889  // It is not necessary to copy the coordinates. There are
2890  // already there.
2891  if( cur_attr_name == "point" )
2892  {
2893  continue;
2894  }
2895 
2896  auto* cur_v_att_store_in_reg =
2897  reg_vertex_attr_mgr.find_attribute_store( cur_attr_name );
2898  ringmesh_assert( cur_v_att_store_in_reg != nullptr );
2899  index_t dim = cur_v_att_store_in_reg->dimension();
2900  GEO::AttributeStore* cur_v_att_store{ nullptr };
2901  if( !vertices.attribute_manager().is_defined( cur_attr_name ) )
2902  {
2903  const auto cur_type_name = GEO::AttributeStore::
2904  element_type_name_by_element_typeid_name(
2905  cur_v_att_store_in_reg->element_typeid_name() );
2907  GEO::AttributeStore::element_type_name_is_known(
2908  cur_type_name ) );
2909  cur_v_att_store = GEO::AttributeStore::
2910  create_attribute_store_by_element_type_name(
2911  cur_type_name, dim );
2912  vertices.attribute_manager().bind_attribute_store(
2913  cur_attr_name, cur_v_att_store );
2914  }
2915  else
2916  {
2917  cur_v_att_store =
2918  vertices.attribute_manager().find_attribute_store(
2919  cur_attr_name );
2920  }
2921  ringmesh_assert( cur_v_att_store != nullptr );
2922  ringmesh_assert( cur_v_att_store->element_size()
2923  == cur_v_att_store_in_reg->element_size() );
2924 
2925  for( auto v_in_reg_itr : range( cur_reg.nb_vertices() ) )
2926  {
2927  auto v_id_in_gmm = vertices.geomodel_vertex_id(
2928  cur_reg.gmme(), v_in_reg_itr );
2929  GEO::Memory::copy(
2930  static_cast< GEO::Memory::pointer >(
2931  cur_v_att_store->data() )
2932  + v_id_in_gmm * dim
2933  * cur_v_att_store->element_size(),
2934  static_cast< GEO::Memory::pointer >(
2935  cur_v_att_store_in_reg->data() )
2936  + v_in_reg_itr * dim
2937  * cur_v_att_store_in_reg->element_size(),
2938  dim * cur_v_att_store->element_size() );
2939  }
2940  }
2941  }
2942  }
2943 
2945  const
2946  {
2947  auto& gmm_c_attr_mgr = cells.attribute_manager();
2948  GEO::vector< std::string > att_c_names;
2949  gmm_c_attr_mgr.list_attribute_names( att_c_names );
2950  const auto& nn_search = cells.cell_nn_search();
2951 
2952  for( const auto& cur_attr_name : att_c_names )
2953  {
2954  auto* cur_c_att_store_in_gmm =
2955  gmm_c_attr_mgr.find_attribute_store( cur_attr_name );
2956  ringmesh_assert( cur_c_att_store_in_gmm != nullptr );
2957  auto dim = cur_c_att_store_in_gmm->dimension();
2958 
2959  for( const auto& cur_region : geomodel_.regions() )
2960  {
2961  auto& reg_c_attr_mgr = cur_region.cell_attribute_manager();
2962  GEO::AttributeStore* cur_c_att_store_in_reg{ nullptr };
2963 
2964  if( !reg_c_attr_mgr.is_defined( cur_attr_name ) )
2965  {
2966  const auto cur_type_name = GEO::AttributeStore::
2967  element_type_name_by_element_typeid_name(
2968  cur_c_att_store_in_gmm->element_typeid_name() );
2970  GEO::AttributeStore::element_type_name_is_known(
2971  cur_type_name ) );
2972  cur_c_att_store_in_reg = GEO::AttributeStore::
2973  create_attribute_store_by_element_type_name(
2974  cur_type_name, dim );
2975  reg_c_attr_mgr.bind_attribute_store(
2976  cur_attr_name, cur_c_att_store_in_reg );
2977  }
2978  else
2979  {
2980  cur_c_att_store_in_reg =
2981  reg_c_attr_mgr.find_attribute_store( cur_attr_name );
2982  }
2983  ringmesh_assert( cur_c_att_store_in_reg != nullptr );
2984  ringmesh_assert( cur_c_att_store_in_reg->element_size()
2985  == cur_c_att_store_in_gmm->element_size() );
2986 
2987  for( auto c : range( cur_region.nb_mesh_elements() ) )
2988  {
2989  auto center = cur_region.mesh_element_barycenter( c );
2990  auto c_in_geom_model_mesh =
2991  nn_search.get_neighbors( center, geomodel_.epsilon() );
2992  ringmesh_assert( c_in_geom_model_mesh.size() == 1 );
2993  GEO::Memory::copy(
2994  static_cast< GEO::Memory::pointer >(
2995  cur_c_att_store_in_reg->data() )
2996  + c * dim * cur_c_att_store_in_reg->element_size(),
2997  static_cast< GEO::Memory::pointer >(
2998  cur_c_att_store_in_gmm->data() )
2999  + c_in_geom_model_mesh[0] * dim
3000  * cur_c_att_store_in_gmm->element_size(),
3001  dim * cur_c_att_store_in_reg->element_size() );
3002  }
3003  }
3004  }
3005  }
3006 
3008  const
3009  {
3010  const auto& nn_search = cells.cell_nn_search();
3011  for( const auto& cur_reg : geomodel().regions() )
3012  {
3013  GEO::vector< std::string > att_c_names;
3014  auto& reg_cell_attr_mgr = cur_reg.cell_attribute_manager();
3015  reg_cell_attr_mgr.list_attribute_names( att_c_names );
3016  for( const auto& cur_attr_name : att_c_names )
3017  {
3018  auto* cur_c_att_store_in_reg =
3019  reg_cell_attr_mgr.find_attribute_store( cur_attr_name );
3020  ringmesh_assert( cur_c_att_store_in_reg != nullptr );
3021  index_t dim = cur_c_att_store_in_reg->dimension();
3022  GEO::AttributeStore* cur_c_att_store{ nullptr };
3023  if( !cells.attribute_manager().is_defined( cur_attr_name ) )
3024  {
3025  const std::string cur_type_name = GEO::AttributeStore::
3026  element_type_name_by_element_typeid_name(
3027  cur_c_att_store_in_reg->element_typeid_name() );
3029  GEO::AttributeStore::element_type_name_is_known(
3030  cur_type_name ) );
3031  cur_c_att_store = GEO::AttributeStore::
3032  create_attribute_store_by_element_type_name(
3033  cur_type_name, dim );
3034  cells.attribute_manager().bind_attribute_store(
3035  cur_attr_name, cur_c_att_store );
3036  }
3037  else
3038  {
3039  cur_c_att_store =
3040  cells.attribute_manager().find_attribute_store(
3041  cur_attr_name );
3042  }
3043  ringmesh_assert( cur_c_att_store != nullptr );
3044  ringmesh_assert( cur_c_att_store->element_size()
3045  == cur_c_att_store_in_reg->element_size() );
3046 
3047  for( auto c_in_reg_itr : range( cur_reg.nb_mesh_elements() ) )
3048  {
3049  auto center =
3050  cur_reg.mesh_element_barycenter( c_in_reg_itr );
3051  auto c_in_geom_model_mesh =
3052  nn_search.get_neighbors( center, geomodel_.epsilon() );
3053  ringmesh_assert( c_in_geom_model_mesh.size() == 1 );
3054  GEO::Memory::copy(
3055  static_cast< GEO::Memory::pointer >(
3056  cur_c_att_store->data() )
3057  + c_in_geom_model_mesh[0] * dim
3058  * cur_c_att_store->element_size(),
3059  static_cast< GEO::Memory::pointer >(
3060  cur_c_att_store_in_reg->data() )
3061  + c_in_reg_itr * dim
3062  * cur_c_att_store_in_reg->element_size(),
3063  dim * cur_c_att_store->element_size() );
3064  }
3065  }
3066  }
3067  }
3068 
3069  template class RINGMESH_API GeoModelMeshBase< 2 >;
3070  template class RINGMESH_API GeoModelMesh< 2 >;
3071  template class RINGMESH_API GeoModelMeshVerticesBase< 2 >;
3072  template class RINGMESH_API GeoModelMeshWells< 2 >;
3073  template class RINGMESH_API GeoModelMeshEdges< 2 >;
3074  template class RINGMESH_API GeoModelMeshPolygonsBase< 2 >;
3075 
3076  template class RINGMESH_API GeoModelMeshBase< 3 >;
3077  template class RINGMESH_API GeoModelMeshVerticesBase< 3 >;
3078  template class RINGMESH_API GeoModelMeshWells< 3 >;
3079  template class RINGMESH_API GeoModelMeshEdges< 3 >;
3080  template class RINGMESH_API GeoModelMeshPolygonsBase< 3 >;
3081  template class RINGMESH_API GeoModelMeshCells< 3 >;
3082 
3083 } // namespace RINGMesh
const LineAABBTree< DIMENSION > & aabb() const
return the AABB tree for the edges of the mesh
index_t tet(index_t region, index_t tet) const
index_t hex(index_t region, index_t hex) const
const std::vector< index_t > & vertex_map(const gmme_id &mesh_entity_id) const
void resize_all_mesh_entity_vertex_maps(const MeshEntityType &type)
void update_mesh_entity_maps_and_gmes(const std::vector< index_t > &old2new)
Updates all the vertex maps with regards to the global indexing changes.
GeoModelMeshVertices< DIMENSION > vertices
void clear_vertex_map(const gmme_id &mesh_entity_id)
std::unique_ptr< PointSetMesh< DIMENSION > > & mesh_
Attached Mesh.
GEO::vecng< DIMENSION, double > vecn
Definition: types.h:74
index_t create_vertex()
Creates a new vertex.
Definition: mesh_builder.h:126
static std::unique_ptr< LineMeshBuilder > create_builder(LineMesh< DIMENSION > &mesh)
void add_vertices(LineMeshBuilder< DIMENSION > *builder, index_t size)
Definition: test-aabb.cpp:69
void initialize()
Initialize the cells from the cells of the GeoModel Region cells.
GeoModelMesh< DIMENSION > & gmm_
Attached GeoModelMesh.
index_t nb_hexs_
Number of hex in the GeoModelMesh.
index_t nb_pyramids_
Number of pyramid in the GeoModelMesh.
const VolumeAABBTree< DIMENSION > & aabb() const
return the AABB tree for the cells of the mesh
bool is_cell_facet_on_surface(index_t cell, index_t facet_index, index_t &colocated_facet_index, bool &side) const
virtual const vecn< DIMENSION > & vertex(index_t v_id) const =0
Gets a point.
GeoModelBuilderGeometry< DIMENSION > geometry
GeoModelMeshPolygonsBase(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< SurfaceMesh< DIMENSION > > &mesh)
index_t vertex(index_t well, index_t edge, index_t vertex) const
void initialize_mesh_entity_vertex_map(const gmme_id &mesh_entity_id)
Initializes the given GeoModelMeshEntity vertex map.
CellType type(index_t cell) const
virtual const GeoModelMeshEntity< DIMENSION > & mesh_entity(const gmme_id &id) const
Generic access to a meshed entity.
Definition: geomodel.cpp:131
index_t nb_edges(index_t line) const
void delete_vertices(const std::vector< bool > &to_delete)
Deletes a set of vertices.
Definition: mesh_builder.h:171
vecn< 3 > vec3
Definition: types.h:76
GeoModelMeshVertices(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< PointSetMesh< DIMENSION > > &mesh)
GeoModel< DIMENSION > & geomodel_
Attached GeoModel.
std::vector< std::vector< index_t > > line_vertex_maps_
index_t prism(index_t region, index_t prism) const
GeoModelMeshVerticesBase(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< PointSetMesh< DIMENSION > > &mesh)
index_t nb_facets(index_t cell) const
index_t nb_connectors_
Number of connector in the GeoModelMesh.
index_t nb_entity_vertices(const GeoModel< DIMENSION > &geomodel, const MeshEntityType &entity_type)
index_t cell(index_t region, index_t cell, CellType type=CellType::UNDEFINED) const
void set_mesh(MeshBase< DIMENSION > *mesh)
Definition: geomodel_mesh.h:89
A GeoModelEntity of type CORNER.
index_t duplicated_corner_index(const ElementLocalVertex &cell_local_vertex) const
std::map< MeshEntityType, std::vector< std::vector< index_t > > *> vertex_maps_
void clear_all_mesh_entity_vertex_map()
Unbinds all the GeoModelMeshEntity vertex maps.
index_t vertex(const ElementLocalVertex &cell_local_vertex) const
static MeshEntityType type_name_static()
index_t nb_edges(index_t cell) const
GeoModelMeshCommon(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &geomodel)
std::unique_ptr< VolumeMesh< DIMENSION > > & mesh_
Attached Mesh.
bool are_corners_to_duplicate(const std::vector< action_on_surface > &surfaces, std::vector< ActionOnSurface > &info)
virtual index_t nb_total_vertices() const
return the number of all vertices it is computed summing all entities.nb().
double length(index_t edge) const
GeoModelMeshBase(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &geomodel)
index_t nb_facet_vertices(const CellLocalFacet &cell_local_facet) const
Entity_type_template type() const
Definition: entity_type.h:202
GeoModelMesh(GeoModel< DIMENSION > &geomodel)
GeoModelMeshCells(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< VolumeMesh< DIMENSION > > &mesh)
std::unique_ptr< LineMesh< DIMENSION > > & mesh_
Attached Mesh.
std::vector< index_t > region_cell_ptr_
GEO::AttributesManager & mesh_entity_vertex_attribute_manager(const gmme_id &mesh_entity_id) const
Returns the vertex attribute of a GeoModelMeshEntity.
index_t index_in_line(index_t edge) const
auto to_underlying_type(Enum e) -> typename std::underlying_type< Enum >::type
Definition: types.h:114
double volume(index_t cell) const
index_t nb() const
Number of cells stored.
vecn< DIMENSION > center(index_t edge) const
static MeshEntityType type_name_static()
index_t nb() const
Number of non colocated vertices stored.
static std::unique_ptr< VolumeMeshBuilder< DIMENSION > > create_builder(VolumeMesh< DIMENSION > &mesh)
GeoModelMeshPolygons(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< SurfaceMesh< DIMENSION > > &mesh)
void sort_unique(CONTAINER &container, const CMP &cmp)
Sorts a container and suppresses all duplicated entities.
Definition: algorithm.h:135
std::vector< index_t > cell_id_
Vector storing the cell index in region per cell.
index_t index_in_region(index_t cell) const
index_t adjacent(index_t cell, index_t facet) const
index_t nb_edges_
Number of edges in the GeoModelMesh.
index_t find(const container &in, const T &value)
Returns the position of the first entity matching.
Definition: algorithm.h:55
index_t nb_tets_
Number of tet in the GeoModelMesh.
index_t nb_cells(CellType type=CellType::UNDEFINED) const
DuplicateMode mode_
Current duplicate mode applied on the mesh.
index_t connector(index_t region, index_t connector) const
index_t geomodel_vertex_index(const gmme_id &mesh_entity_id, index_t mesh_entity_vertex_index) const
Returns the index of a GeoModelMeshEntity vertex in the geomodel global indexing. ...
CellType
Definition: types.h:89
std::string MeshType
Definition: mesh.h:69
void test_and_initialize_cell_facet() const
bool test_and_initialize_mesh_entity_vertex_map(const gmme_id &mesh_entity_id)
Tests if the given GeoModelMeshEntity vertex map is initialized. If not, initializes it...
index_t create_vertices(index_t nb)
Creates a contiguous chunk of vertices.
Definition: mesh_builder.h:148
index_t edge(index_t line, index_t edge) const
std::vector< index_t > duplicated_vertex_indices_
Vector of duplicated vertices.
GeoModelVertexMapper(GeoModelMeshVerticesBase &geomodel_vertices, const GeoModel< DIMENSION > &geomodel)
static MeshEntityType type_name_static()
index_t nb_vertices(index_t cell) const
index_t vertex(const ElementLocalVertex &edge_local_vertex) const
void set_vertex(index_t v_id, const vecn< DIMENSION > &vertex)
Sets a point.
Definition: mesh_builder.h:117
PolygonType
Definition: types.h:105
const vecn< DIMENSION > & vertex(index_t v) const
Coordinates of a vertex of the GeoModel.
index_t nb_polygons(PolygonType type=PolygonType::UNDEFINED) const
encapsulate adimensional mesh functionalities in order to provide an API on which we base the RINGMes...
Definition: mesh.h:138
std::vector< std::vector< index_t > > corner_vertex_maps_
Vertex maps.
index_t facet_vertex(const CellLocalFacet &cell_local_facet, index_t local_vertex) const
Gets a cell vertex by local facet index and local vertex index in the edge.
virtual index_t nb_cell_facets(index_t cell_id) const =0
Gets the number of facet in a cell.
std::vector< index_t > line_edge_ptr_
std::vector< index_t > region_id_
Vector storing the region index per cell.
need to be duplicated (don&#39;t know which side yet)
virtual index_t nb_cell_facet_vertices(const CellLocalFacet &cell_local_facet) const =0
Gets the number of vertices of a facet in a cell.
void test_and_initialize_duplication() const
index_t edge_vertex(index_t cell, index_t local_edge, index_t local_vertex) const
Gets a cell vertex by local edge index and local vertex index in the edge.
std::vector< index_t > edge_id_
Vector storing the edge index in line per edge.
GeoModelVertexMapper vertex_mapper_
Mapper from/to GeoModelMeshEntity vertices.
virtual index_t nb_vertices() const =0
std::vector< index_t > & resize_vertex_map(const gmme_id &mesh_entity_id)
static bool is_fault(GEOL_FEATURE feature)
std::vector< std::vector< GMEVertex > > gme_vertices_
GeoModelEntity Vertices for each geomodel vertex.
GeoModelMeshVerticesBase< DIMENSION > & geomodel_vertices_
std::vector< std::vector< index_t > > region_vertex_maps_
Vertex in a GeoModelEntity.
#define ringmesh_assert(x)
vec3 barycenter(index_t cell) const
const WellPart< DIMENSION > & part(index_t part_id) const
Definition: well.h:325
bool contains(const container &in, const T &value)
Definition: algorithm.h:87
index_t nb() const
Number of edges stored.
index_t pyramid(index_t region, index_t pyramid) const
void clear()
Clears all the information about vertex mapping (vector maps and vectors of GME_Vertices.
GEO::vector< U > copy_std_vector_to_geo_vector(const std::vector< T > &in, index_t from, index_t to)
void fill_vertices_for_entity_type(const GeoModel< DIMENSION > &M, const MeshEntityType &entity_type, index_t &count)
index_t nb_parts() const
Definition: well.h:360
std::vector< std::vector< index_t > > surface_vertex_maps_
The MeshEntityType described the type of the meshed entities There are 4 MeshEntityTypes correspondin...
Definition: entity_type.h:117
bool is_mesh_entity_vertex_map_initialized(const gmme_id &mesh_entity_id) const
Tests if the given GeoModelMeshEntity vertex map exists.
void add_to_gme_vertices(const GMEVertex &gme_vertex, index_t geomodel_vertex_index)
void clear_vertices(bool keep_attributes, bool keep_memory)
Removes all the vertices and attributes.
Definition: mesh_builder.h:184
std::vector< index_t > polygon_id_
Vector storing the colocalised polygon index per cell facet If a cell facet is on a surface...
Classes to build GeoModel from various inputs.
Definition: algorithm.h:48
std::vector< index_t > line_id_
Vector storing the line index per edge.
const GeoModel< DIMENSION > & geomodel() const
std::unique_ptr< SurfaceMesh< DIMENSION > > & mesh_
Attached Mesh.
A GeoModelEntity of type LINE.
void erase_vertices(std::vector< index_t > &to_delete)
Delete vertices for which to_delete[i] != i The global vertices are deleted, gme_vertices_ is update...
GeoModelMeshEdges(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< LineMesh< DIMENSION > > &mesh)
std::vector< index_t > get_neighbors(const vecn< DIMENSION > &v, double threshold_distance) const
Definition: nn_search.cpp:205
index_t line(index_t edge) const
index_t nb_mesh_elements() const final
index_t region(index_t cell) const
virtual index_t nb_mesh_entities(const MeshEntityType &type) const
Returns the number of mesh entities of the given type.
Definition: geomodel.cpp:108
index_t index() const
Definition: entity_type.h:197
virtual index_t cell_facet_vertex(const CellLocalFacet &cell_local_facet, index_t vertex_id) const =0
Gets a vertex by cell facet and local vertex index.
MeshSet< DIMENSION > mesh_set_
Mesh storing all the elements.
virtual void clear()
Clear the vertices - clear the gme_vertices_ - clear global vertex information in the all BMME...
const LineAABBTree< DIMENSION > & aabb() const
return the AABB tree for the edges of the mesh
index_t nb_duplicated_vertices() const
const GeoModel< DIMENSION > & geomodel_
This template is a specialization of a gme_id to the GeoModelMeshEntity.
Definition: entity_type.h:285
bool is_surface_to_duplicate(index_t surface) const
GeoModelMeshWells(GeoModelMesh< DIMENSION > &gmm, GeoModel< DIMENSION > &gm, std::unique_ptr< LineMesh< DIMENSION > > &mesh)
std::unique_ptr< LineMesh< DIMENSION > > & mesh_
Attached Mesh.
index_t duplicated_vertex(index_t duplicate_vertex_index) const
index_t local_facet_id_
Definition: mesh.h:127
std::vector< index_t > well_ptr_
index_t nb_prisms_
Number of prism in the GeoModelMesh.
void clear(bool keep_attributes, bool keep_memory)
Removes all the entities and attributes of this mesh.
Definition: mesh_builder.h:90
const std::vector< GMEVertex > & mesh_entity_vertex_indices(index_t v) const
Returns all the corresponding vertices in GeoModelMeshEntities to a given geomodel vertex...
#define ringmesh_assert_not_reached
void set_vertex_map_value(const gmme_id &mesh_entity_id, index_t mesh_entity_vertex_index, index_t geomodel_entity_vertex_index)
Sets the geomodel vertex mapping value of a given vertex in a GeoModelMeshEntity. ...