29 #ifndef __PointSetToMesh_H
30 #define __PointSetToMesh_H 1
32 #include <boost/shared_ptr.hpp>
33 #include <boost/optional.hpp>
35 #include <feel/feelcore/feel.hpp>
39 #if defined(FEELPP_HAS_VTK)
41 #include <vtkPointSet.h>
42 #include <vtkDelaunay2D.h>
43 #include <vtkDelaunay3D.h>
44 #include <vtkPolyData.h>
45 #include <vtkCellArray.h>
46 #include <vtkExtractEdges.h>
47 #include <vtkPolyDataMapper.h>
48 #include <vtkDataSetMapper.h>
67 template<
typename Convex,
typename T>
71 public Visitor<PointSet<Convex, T> >
81 static const uint16_type nDim = Convex::nDim;
91 typedef typename convex_type::template shape<convex_type::nDim, 1, convex_type::nDim>::type mesh_convex_type;
96 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
100 typedef typename mesh_type::point_type point_type;
101 typedef typename point_type::node_type node_type;
102 typedef typename mesh_type::edge_type edge_type;
103 typedef typename mesh_type::face_type face_type;
104 typedef typename mesh_type::element_type element_type;
116 M_mesh(
new mesh_type( Environment::worldComm().subWorldCommSeq() ) ),
124 M_vertices( p.M_vertices )
154 void addBoundaryPoints( points_type
const& v )
171 visit( pset, mpl::int_<nDim>() );
183 void visit( pointset_type* pset, mpl::int_<1> );
184 void visit( pointset_type* pset, mpl::int_<2> );
185 void visit( pointset_type* pset, mpl::int_<3> );
189 boost::optional<points_type> M_vertices;
192 template<
typename Convex,
typename T>
196 DVLOG(2) <<
"[PointSetToMesh::visit(<1>)] pointset to mesh\n";
197 M_mesh = mesh_ptrtype(
new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
201 for ( uint __i = 0; __i < __npts; ++__i )
204 __nd = pset->point( __i );
206 point_type __pt( __i,__nd,
false );
209 if ( __nd[0] == -1 || __nd[0] == 1 )
211 __pt.setOnBoundary(
true );
216 __pt.setOnBoundary(
false );
220 M_mesh->addPoint( __pt );
228 face_type* pf0 =
new face_type;
231 pf0->setPoint( 0, M_mesh->point( 0 ) );
233 pf0->setId( n_faces++ );
234 pf0->setOnBoundary(
true );
235 M_mesh->addFace( *pf0 );
239 face_type* pf1 =
new face_type;
242 pf1->setPoint( 0, M_mesh->point( 1 ) );
244 pf1->setId( n_faces++ );
245 pf1->setOnBoundary(
true );
246 M_mesh->addFace( *pf1 );
252 for ( uint __i = 0; __i < __nele; ++__i )
261 pf.setPoint( 0, M_mesh->point( __i ) );
262 pf.setPoint( 1, M_mesh->point( __i+1 ) );
267 pf.setPoint( 0, M_mesh->point( __i ) );
268 pf.setPoint( 1, M_mesh->point( __i+2 ) );
271 else if ( __i == __nele-1 )
273 pf.setPoint( 0, M_mesh->point( __i+1 ) );
274 pf.setPoint( 1, M_mesh->point( 1 ) );
279 pf.setPoint( 0, M_mesh->point( __i+1 ) );
280 pf.setPoint( 1, M_mesh->point( __i+2 ) );
283 element_type
const& e = M_mesh->addElement( pf );
284 DVLOG(2) <<
"o element " << e.id() <<
"\n"
285 <<
" p1 = " << e.point( 0 ).node() <<
"\n"
286 <<
" p2 = " << e.point( 1 ).node() <<
"\n";
290 FEELPP_ASSERT( n_faces == M_mesh->numFaces() )( n_faces )( M_mesh->numFaces() ).error(
"invalid face container size" );
292 DVLOG(2) <<
"[PointSetToMesh<1>] done with element accumulation !\n";
294 M_mesh->setNumVertices( __npts );
296 DVLOG(2) <<
"[PointSetToMesh<1>] Face Update !\n";
301 DVLOG(2) <<
"[PointSetToMesh<1>] Face Update Successful !\n";
304 DVLOG(2) <<
"[PointSetToMesh::visit(<1>)] done\n";
306 template<
typename Convex,
typename T>
310 #if defined(FEELPP_HAS_VTK)
312 M_mesh = mesh_ptrtype(
new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
314 vtkPoints *newPoints = vtkPoints::New();
318 DVLOG(2) <<
"adding vertices\n" << M_vertices.get() <<
"\n";
320 for (
size_type i = 0; i < M_vertices->size2(); ++i )
322 newPoints->InsertNextPoint( M_vertices.get()( 0, i ), M_vertices.get()( 1, i ), 0 );
326 std::vector<double> perturbation( pset->nPoints() );
328 for (
size_type i=0 ; i< pset->nPoints() ; ++i )
330 perturbation[i] = std::rand()*1e-6/RAND_MAX;
331 uint16_type index = newPoints->InsertNextPoint( pset->points()( 0,i )+perturbation[i], pset->points()( 1,i ), 0 );
332 DVLOG(2) <<
"Inserting point with id " << index <<
"\n";
333 DVLOG(2) <<
"pset.point( " << i <<
" )= " << pset->point( i ) <<
"\n";
336 vtkPolyData *polyData = vtkPolyData::New();
337 polyData->SetPoints( newPoints );
339 vtkDelaunay2D *delaunay2D = vtkDelaunay2D::New();
341 delaunay2D->SetInput( polyData );
349 delaunay2D->Update();
351 vtkPolyData* outMesh = delaunay2D->GetOutput( );
355 outMesh->BuildLinks();
357 int Nelem = outMesh->GetNumberOfPolys();
358 int Npts = outMesh->GetNumberOfPoints();
360 DVLOG(2) <<
"Number of cells = " << Nelem <<
"\n";
361 DVLOG(2) <<
"Number of points = " << Npts <<
"\n";
363 for (
int i=0; i< Npts; ++i )
366 outMesh->GetPoints()->SetPoint( i,
367 outMesh->GetPoints()->GetPoint( i )[0]- perturbation[i],
368 outMesh->GetPoints()->GetPoint( i )[1],
369 outMesh->GetPoints()->GetPoint( i )[2] );
372 for (
int i=0; i< Nelem; ++i )
375 DVLOG(2) <<
"Element Id = " << i <<
"\n";
376 DVLOG(2) <<
"Point 0 (" << ( int )outMesh->GetCell( i )->GetPointId( 0 ) <<
") =" ;
377 DVLOG(2) <<
"(" << outMesh->GetCell( i )->GetPoints()->GetPoint( 0 )[0] <<
" , "
378 << outMesh->GetCell( i )->GetPoints()->GetPoint( 0 )[1]<<
")" <<
"\n";
379 DVLOG(2) <<
"Point 1 (" << ( int )outMesh->GetCell( i )->GetPointId( 1 ) <<
") =" ;
380 DVLOG(2) <<
"(" << outMesh->GetCell( i )->GetPoints()->GetPoint( 1 )[0] <<
" , "
381 << outMesh->GetCell( i )->GetPoints()->GetPoint( 1 )[1]<<
")" <<
"\n";
382 DVLOG(2) <<
"Point 2 (" << ( int )outMesh->GetCell( i )->GetPointId( 2 ) <<
") =" ;
383 DVLOG(2) <<
"(" << outMesh->GetCell( i )->GetPoints()->GetPoint( 2 )[0] <<
" , "
384 << outMesh->GetCell( i )->GetPoints()->GetPoint( 2 )[1]<<
")" <<
"\n";
386 DVLOG(2) << outMesh->GetCell( i )->GetNumberOfEdges() <<
"\n";
387 DVLOG(2) << ( int )outMesh->GetCell( i )->GetEdge( 0 )->GetPointId( 0 ) <<
"\n";
388 DVLOG(2) << ( int )outMesh->GetCell( i )->GetEdge( 0 )->GetPointId( 1 ) <<
"\n";
393 DVLOG(2) <<
"[PointSetToMesh::visit(<2>)] delaunay done, now vtk to Mesh<>\n";
394 FilterFromVtk<mesh_type> meshfromvtk( outMesh );
395 meshfromvtk.visit( M_mesh.get() );
396 DVLOG(2) <<
"[PointSetToMesh::visit(<2>)] done\n";
399 std::cerr <<
"The library was not compiled with vtk support\n";
404 template<
typename Convex,
typename T>
408 #if defined(FEELPP_HAS_VTK)
410 M_mesh = mesh_ptrtype(
new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
412 vtkPoints *newPoints = vtkPoints::New();
417 newPoints->InsertNextPoint( -1, -1, -1 );
418 newPoints->InsertNextPoint( 1, -1, -1 );
419 newPoints->InsertNextPoint( -1, 1, -1 );
420 newPoints->InsertNextPoint( -1, -1, 1 );
423 for (
size_type i=0 ; i< pset->nPoints() ; ++i )
425 newPoints->InsertNextPoint( pset->points()( 0,i ), pset->points()( 1,i ), pset->points()( 2,i ) );
430 vtkPolyData *polyData = vtkPolyData::New();
431 polyData->SetPoints( newPoints );
433 vtkDelaunay3D *delaunay3D = vtkDelaunay3D::New();
434 delaunay3D->SetInput( polyData );
435 delaunay3D->SetOffset( 5 );
437 DVLOG(2) <<
"[PointSetToMesh::visit(<3>)] Offset = " << delaunay3D->GetOffset() <<
"\n";
438 delaunay3D->Update();
440 vtkUnstructuredGrid* outMesh = delaunay3D->GetOutput( );
443 DVLOG(2) <<
"[PointSetToMesh::visit(<3>)] delaunay done, now vtk to Mesh<>\n";
444 FilterFromVtk3D<mesh_type> meshfromvtk( outMesh );
445 meshfromvtk.visit( M_mesh.get() );
446 DVLOG(2) <<
"[PointSetToMesh::visit(<3>)] done\n";
450 std::cerr <<
"The library was not compiled with vtk support\n";
The base class of any Acyclic Visitor.
Definition: visitor.hpp:52
Definition: visitor.hpp:89
void visit(pointset_type *pset)
Definition: pointsettomesh.hpp:169
Convex base class.
Definition: convex.hpp:43
transform a point set to a mesh data structure using a Delaunay
Definition: pointsettomesh.hpp:68
unifying mesh class
Definition: createsubmesh.hpp:41
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Class of all PointSet on a Convex.
Definition: pointset.hpp:54