00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <iostream>
00022
00023
00024 #include "vtkPoints.h"
00025 #include "vtkPolyData.h"
00026 #include "vtkDelaunay3D.h"
00027 #include "vtkShrinkFilter.h"
00028 #include "vtkDataSetMapper.h"
00029 #include "vtkOutlineFilter.h"
00030 #include "vtkProperty.h"
00031 #include "vtkPolyDataMapper.h"
00032 #include "vtkGeometryFilter.h"
00033 #include "vtkUnstructuredGridGeometryFilter.h"
00034 #include "vtkUnstructuredGridWriter.h"
00035 #include "vtkPolyDataWriter.h"
00036
00037
00038 #include "util/exception.h"
00039 #include "numerical/meshgrid.h"
00040 #include "vtkwrap/vtkwin.h"
00041
00042 using std::cout;
00043 using std::endl;
00044
00045 int main(int argc, char **argv) try
00046 {
00047 if (argc!=2)
00048 {
00049 cout << "Usage:\n";
00050 cout << " " << argv[0] << " NPTS\n";
00051 cout << endl;
00052 cout << "ex.: " << argv[0] << " 6\n";
00053 cout << endl;
00054 return 1;
00055 }
00056
00057
00058 int N = atoi(argv[1]); if (N<2) throw new Fatal(_("NPTS must be greater than 1"));
00059 REAL L = 10.0;
00060 MeshGrid mg(0.0,L,atoi(argv[1]),
00061 0.0,L,atoi(argv[1]),
00062 0.0,L,atoi(argv[1]));
00063
00064
00065 vtkPoints * points = vtkPoints::New();
00066 points->Allocate(mg.Length());
00067 for (int i=0; i<mg.Length(); ++i)
00068 {
00069 double P[3] = {mg.X(i), mg.Y(i), mg.Z(i)};
00070 points->InsertPoint(i,P);
00071 }
00072
00073
00074
00075
00076
00077
00078 vtkPolyData * vertices = vtkPolyData::New();
00079 vtkDelaunay3D * delaunay = vtkDelaunay3D::New();
00080 vertices -> SetPoints(points);
00081 delaunay -> SetInput(vertices);
00082 delaunay -> SetTolerance(0.01);
00083
00084 delaunay -> SetAlpha(0);
00085 delaunay -> BoundingTriangulationOff();
00086
00087
00088 vtkGeometryFilter * geofil = vtkGeometryFilter::New();
00089
00090 geofil -> SetInputConnection(delaunay->GetOutputPort());
00091
00092
00093 vtkShrinkFilter * shrink = vtkShrinkFilter::New();
00094
00095 shrink -> SetInputConnection(geofil->GetOutputPort());
00096 shrink -> SetShrinkFactor(0.85);
00097
00098
00099 vtkDataSetMapper * mapper = vtkDataSetMapper::New();
00100 vtkActor * actor = vtkActor::New();
00101 mapper -> SetInputConnection (shrink->GetOutputPort());
00102 actor -> SetMapper (mapper);
00103 actor -> GetProperty () -> SetColor(0,1,0);
00104
00105
00106 vtkOutlineFilter * outline = vtkOutlineFilter::New();
00107 vtkPolyDataMapper * outline_mapper = vtkPolyDataMapper::New();
00108 vtkActor * outline_actor = vtkActor::New();
00109 outline -> SetInputConnection(delaunay->GetOutputPort());
00110 outline_mapper -> SetInput (outline->GetOutput());
00111 outline_actor -> SetMapper (outline_mapper);
00112 outline_actor -> GetProperty () -> SetColor(0,0,0);
00113
00114
00115 VTKWin win;
00116 win.AddActor(actor);
00117 win.AddActor(outline_actor);
00118 win.Show();
00119
00120
00121 vtkUnstructuredGridWriter * writer = vtkUnstructuredGridWriter::New();
00122 writer -> SetInputConnection(delaunay->GetOutputPort());
00123 writer -> SetFileName("test10.elements.vtk");
00124 writer -> Write();
00125
00126
00127 vtkPolyDataWriter * fwriter = vtkPolyDataWriter::New();
00128
00129 fwriter -> SetInputConnection(geofil->GetOutputPort());
00130 fwriter -> SetFileName("test10.faces.vtk");
00131 fwriter -> Write();
00132
00133
00134 points -> Delete();
00135 vertices -> Delete();
00136 delaunay -> Delete();
00137 shrink -> Delete();
00138 mapper -> Delete();
00139 actor -> Delete();
00140 outline -> Delete();
00141 outline_mapper -> Delete();
00142 outline_actor -> Delete();
00143 writer -> Delete();
00144 fwriter -> Delete();
00145
00146 return 0;
00147 }
00148 catch (Exception * e)
00149 {
00150 e->Cout();
00151 if (e->IsFatal()) {delete e; exit(1);}
00152 delete e;
00153 }
00154 catch (char const * m)
00155 {
00156 std::cout << "Fatal: " << m << std::endl;
00157 exit (1);
00158 }
00159 catch (...)
00160 {
00161 std::cout << "Some exception (...) ocurred\n";
00162 }
00163
00164