Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pointsettomesh.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-11-25
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __PointSetToMesh_H
30 #define __PointSetToMesh_H 1
31 
32 #include <boost/shared_ptr.hpp>
33 #include <boost/optional.hpp>
34 
35 #include <feel/feelcore/feel.hpp>
36 
37 #include <stdlib.h>
38 
39 #if defined(FEELPP_HAS_VTK)
40 
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>
49 
50 #endif /* FEELPP_HAS_VTK */
51 
54 
55 namespace Feel
56 {
57 
67 template<typename Convex, typename T>
69  :
70 public VisitorBase,
71 public Visitor<PointSet<Convex, T> >
72 {
73  typedef VisitorBase super1;
75 
76 public:
80 
81  static const uint16_type nDim = Convex::nDim;
82 
84 
88 
89  typedef T value_type;
90  typedef Convex convex_type;
91  typedef typename convex_type::template shape<convex_type::nDim, 1, convex_type::nDim>::type mesh_convex_type;
92 
94 
96  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
97 
98  typedef typename matrix_node<value_type>::type points_type;
99 
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;
105 
107 
111 
113  :
114  super1(),
115  super2(),
116  M_mesh( new mesh_type( Environment::worldComm().subWorldCommSeq() ) ),
117  M_vertices()
118  {}
119  PointSetToMesh( PointSetToMesh const & p )
120  :
121  super1( p ),
122  super2( p ),
123  M_mesh( p.M_mesh ),
124  M_vertices( p.M_vertices )
125  {}
126 
127  ~PointSetToMesh()
128  {}
129 
131 
135 
136 
138 
142 
143  mesh_ptrtype mesh()
144  {
145  return M_mesh;
146  }
147 
149 
153 
154  void addBoundaryPoints( points_type const& v )
155  {
156  M_vertices = v;
157  }
158 
160 
164 
169  void visit( pointset_type* pset )
170  {
171  visit( pset, mpl::int_<nDim>() );
172  }
173 
174 
176 
177 
178 
179 protected:
180 
181 private:
182 
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> );
186 private:
187 
188  mesh_ptrtype M_mesh;
189  boost::optional<points_type> M_vertices;
190 };
191 
192 template<typename Convex, typename T>
193 void
194 PointSetToMesh<Convex, T>::visit( pointset_type* pset, mpl::int_<1> )
195 {
196  DVLOG(2) << "[PointSetToMesh::visit(<1>)] pointset to mesh\n";
197  M_mesh = mesh_ptrtype( new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
198 
199  size_type __npts = pset->nPoints();
200 
201  for ( uint __i = 0; __i < __npts; ++__i )
202  {
203  node_type __nd( 1 );
204  __nd = pset->point( __i );
205 
206  point_type __pt( __i,__nd, false );
207  __pt.marker() = 0;
208 
209  if ( __nd[0] == -1 || __nd[0] == 1 )
210  {
211  __pt.setOnBoundary( true );
212  }
213 
214  else
215  {
216  __pt.setOnBoundary( false );
217 
218  }
219 
220  M_mesh->addPoint( __pt );
221  }
222 
223  size_type n_faces = 0;
224 
225 
226  // Add Boundary faces
227 
228  face_type* pf0 = new face_type;
229 
230  pf0->marker() = 0;
231  pf0->setPoint( 0, M_mesh->point( 0 ) );
232 
233  pf0->setId( n_faces++ );
234  pf0->setOnBoundary( true );
235  M_mesh->addFace( *pf0 );
236 
237  delete pf0;
238 
239  face_type* pf1 = new face_type;
240 
241  pf1->marker() = 0;
242  pf1->setPoint( 0, M_mesh->point( 1 ) );
243 
244  pf1->setId( n_faces++ );
245  pf1->setOnBoundary( true );
246  M_mesh->addFace( *pf1 );
247 
248  delete pf1;
249 
250  size_type __nele = __npts-1;
251 
252  for ( uint __i = 0; __i < __nele; ++__i )
253  {
254  element_type pf;
255 
256  pf.setId( __i );
257  pf.marker() = 0;
258 
259  if ( __nele == 1 )
260  {
261  pf.setPoint( 0, M_mesh->point( __i ) );
262  pf.setPoint( 1, M_mesh->point( __i+1 ) );
263  }
264 
265  else if ( __i == 0 )
266  {
267  pf.setPoint( 0, M_mesh->point( __i ) );
268  pf.setPoint( 1, M_mesh->point( __i+2 ) );
269  }
270 
271  else if ( __i == __nele-1 )
272  {
273  pf.setPoint( 0, M_mesh->point( __i+1 ) );
274  pf.setPoint( 1, M_mesh->point( 1 ) );
275  }
276 
277  else
278  {
279  pf.setPoint( 0, M_mesh->point( __i+1 ) );
280  pf.setPoint( 1, M_mesh->point( __i+2 ) );
281  }
282 
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";
287  }
288 
289 
290  FEELPP_ASSERT( n_faces == M_mesh->numFaces() )( n_faces )( M_mesh->numFaces() ).error( "invalid face container size" );
291 
292  DVLOG(2) <<"[PointSetToMesh<1>] done with element accumulation !\n";
293 
294  M_mesh->setNumVertices( __npts );
295 
296  DVLOG(2) <<"[PointSetToMesh<1>] Face Update !\n";
297 
299  //M_mesh->updateForUse( MESH_ALL_COMPONENTS & (~MESH_RENUMBER) );
300 
301  DVLOG(2) <<"[PointSetToMesh<1>] Face Update Successful !\n";
302 
303 
304  DVLOG(2) << "[PointSetToMesh::visit(<1>)] done\n";
305 }
306 template<typename Convex, typename T>
307 void
308 PointSetToMesh<Convex, T>::visit( pointset_type* pset, mpl::int_<2> )
309 {
310 #if defined(FEELPP_HAS_VTK)
311  // reinitialize mesh
312  M_mesh = mesh_ptrtype( new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
313 
314  vtkPoints *newPoints = vtkPoints::New();
315 
316  if ( M_vertices )
317  {
318  DVLOG(2) << "adding vertices\n" << M_vertices.get() << "\n";
319 
320  for ( size_type i = 0; i < M_vertices->size2(); ++i )
321  {
322  newPoints->InsertNextPoint( M_vertices.get()( 0, i ), M_vertices.get()( 1, i ), 0 );
323  }
324  }
325 
326  std::vector<double> perturbation( pset->nPoints() );
327 
328  for ( size_type i=0 ; i< pset->nPoints() ; ++i )
329  {
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";
334  }
335 
336  vtkPolyData *polyData = vtkPolyData::New();
337  polyData->SetPoints( newPoints );
338 
339  vtkDelaunay2D *delaunay2D = vtkDelaunay2D::New();
340 
341  delaunay2D->SetInput( polyData );
342 
348  //delaunay2D->SetOffset( 8 ) ;
349  delaunay2D->Update();
350 
351  vtkPolyData* outMesh = delaunay2D->GetOutput( );
352 
355  outMesh->BuildLinks();
356 
357  int Nelem = outMesh->GetNumberOfPolys();
358  int Npts = outMesh->GetNumberOfPoints();
359 
360  DVLOG(2) << "Number of cells = " << Nelem << "\n";
361  DVLOG(2) << "Number of points = " << Npts << "\n";
362 
363  for ( int i=0; i< Npts; ++i )
364  {
365  // revert back the perturbation
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] );
370  }
371 
372  for ( int i=0; i< Nelem; ++i )
373  {
374  // std::cout << "\nLa cellule num�ro : " << i << " compte " << outMesh->GetCell(i)->GetNumberOfPoints() << " points." << "\n";
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";
385 
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";
389  }
390 
391 
392 
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";
397 
398 #else
399  std::cerr << "The library was not compiled with vtk support\n";
400 #endif /* FEELPP_HAS_VTK */
401 
402 }
403 
404 template<typename Convex, typename T>
405 void
406 PointSetToMesh<Convex, T>::visit( pointset_type* pset, mpl::int_<3> )
407 {
408 #if defined(FEELPP_HAS_VTK)
409  // reinitialize mesh
410  M_mesh = mesh_ptrtype( new mesh_type( Environment::worldComm().subWorldCommSeq() ) );
411 
412  vtkPoints *newPoints = vtkPoints::New();
413 
414  //if ( M_add_bdy_pts )
415  if ( 0 )
416  {
417  newPoints->InsertNextPoint( -1, -1, -1 );
418  newPoints->InsertNextPoint( 1, -1, -1 );
419  newPoints->InsertNextPoint( -1, 1, -1 );
420  newPoints->InsertNextPoint( -1, -1, 1 );
421  }
422 
423  for ( size_type i=0 ; i< pset->nPoints() ; ++i )
424  {
425  newPoints->InsertNextPoint( pset->points()( 0,i ), pset->points()( 1,i ), pset->points()( 2,i ) );
426  }
427 
428  // create more points
429 
430  vtkPolyData *polyData = vtkPolyData::New();
431  polyData->SetPoints( newPoints );
432 
433  vtkDelaunay3D *delaunay3D = vtkDelaunay3D::New();
434  delaunay3D->SetInput( polyData );
435  delaunay3D->SetOffset( 5 );
436 
437  DVLOG(2) <<"[PointSetToMesh::visit(<3>)] Offset = " << delaunay3D->GetOffset() << "\n";
438  delaunay3D->Update();
439 
440  vtkUnstructuredGrid* outMesh = delaunay3D->GetOutput( );
441 
442 
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";
447 
448 
449 #else
450  std::cerr << "The library was not compiled with vtk support\n";
451 #endif /* FEELPP_HAS_VTK */
452 }
453 } // Feel
454 #endif /* __PointSetToMesh_H */
The base class of any Acyclic Visitor.
Definition: visitor.hpp:52
Definition: visitor.hpp:89
Definition: glas.hpp:157
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

Generated on Sun Dec 22 2013 13:11:11 for Feel++ by doxygen 1.8.5