PRAgMaTIc  master
VTKTools.h
Go to the documentation of this file.
1 /* Copyright (C) 2010 Imperial College London and others.
2  *
3  * Please see the AUTHORS file in the main source directory for a
4  * full list of copyright holders.
5  *
6  * Gerard Gorman
7  * Applied Modelling and Computation Group
8  * Department of Earth Science and Engineering
9  * Imperial College London
10  *
11  * g.gorman@imperial.ac.uk
12  *
13  * Redistribution and use in source and binary forms, with or without
14  * modification, are permitted provided that the following conditions
15  * are met:
16  * 1. Redistributions of source code must retain the above copyright
17  * notice, this list of conditions and the following disclaimer.
18  * 2. Redistributions in binary form must reproduce the above
19  * copyright notice, this list of conditions and the following
20  * disclaimer in the documentation and/or other materials provided
21  * with the distribution.
22  *
23  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
24  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
25  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
26  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
27  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
28  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
30  * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
31  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
32  * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
33  * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
34  * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35  * SUCH DAMAGE.
36  */
37 
38 #ifndef VTK_TOOLS_H
39 #define VTK_TOOLS_H
40 
41 #include <vtkCellType.h>
42 #include <vtkXMLUnstructuredGridReader.h>
43 #include <vtkXMLPUnstructuredGridReader.h>
44 #include <vtkXMLUnstructuredGridWriter.h>
45 #include <vtkXMLPUnstructuredGridWriter.h>
46 #include <vtkUnstructuredGrid.h>
47 #include <vtkIntArray.h>
48 #include <vtkIdTypeArray.h>
49 #include <vtkUnsignedCharArray.h>
50 #include <vtkUnsignedIntArray.h>
51 #include <vtkDoubleArray.h>
52 #include <vtkCell.h>
53 #include <vtkPoints.h>
54 #include <vtkPointData.h>
55 #include <vtkDataArray.h>
56 #include <vtkCellData.h>
57 #include <vtkSmartPointer.h>
58 
59 #include <vtkJPEGReader.h>
60 #include <vtkImageData.h>
61 #include <vtkDataSetTriangleFilter.h>
62 #include <vtkImageGaussianSmooth.h>
63 #include <vtkImageAnisotropicDiffusion2D.h>
64 
65 #ifndef vtkFloatingPointType
66 #define vtkFloatingPointType vtkFloatingPointType
67 typedef float vtkFloatingPointType;
68 #endif
69 
70 #include <vector>
71 #include <string>
72 #include <cfloat>
73 #include <typeinfo>
74 
75 #include "Mesh.h"
76 #include "MetricTensor.h"
77 #include "ElementProperty.h"
78 
79 extern "C" {
80 #include "metis.h"
81 }
82 
83 #ifdef HAVE_MPI
84 #include "mpi_tools.h"
85 #endif
86 
87 #ifdef HAVE_BOOST_UNORDERED_MAP_HPP
88 #include <boost/unordered_map.hpp>
89 #endif
90 
95 template<typename real_t> class VTKTools{
96  public:
97  static Mesh<real_t>* import_vtu(std::string filename){
98  vtkSmartPointer<vtkUnstructuredGrid> ug = vtkSmartPointer<vtkUnstructuredGrid>::New();
99 
100  if(filename.substr(filename.find_last_of('.'))==".pvtu"){
101  vtkSmartPointer<vtkXMLPUnstructuredGridReader> reader = vtkSmartPointer<vtkXMLPUnstructuredGridReader>::New();
102  reader->SetFileName(filename.c_str());
103  reader->Update();
104 
105  ug->DeepCopy(reader->GetOutput());
106  }else{
107  vtkSmartPointer<vtkXMLUnstructuredGridReader> reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
108  reader->SetFileName(filename.c_str());
109  reader->Update();
110 
111  ug->DeepCopy(reader->GetOutput());
112  }
113 
114  size_t NNodes = ug->GetNumberOfPoints();
115  size_t NElements = ug->GetNumberOfCells();
116 
117  std::map<int, int> renumber;
118  std::map<int, int> gnns;
119 
120  std::vector<real_t> x,y,z;
121  int ncnt=0;
122  bool have_global_ids = ug->GetPointData()->GetArray("GlobalId")!=NULL;
123  if(have_global_ids){
124  for(size_t i=0;i<NNodes;i++){
125  int gnn = ug->GetPointData()->GetArray("GlobalId")->GetTuple1(i);
126  if(gnns.find(gnn)==gnns.end()){
127  gnns[gnn] = ncnt;;
128  renumber[i] = ncnt;
129  ncnt++;
130 
131  real_t r[3];
132  ug->GetPoints()->GetPoint(i, r);
133  x.push_back(r[0]);
134  y.push_back(r[1]);
135  z.push_back(r[2]);
136  }else{
137  renumber[i] = gnns[gnn];
138  }
139  }
140  }else{
141  for(size_t i=0;i<NNodes;i++){
142  real_t r[3];
143  ug->GetPoints()->GetPoint(i, r);
144  x.push_back(r[0]);
145  y.push_back(r[1]);
146  z.push_back(r[2]);
147  }
148  }
149 
150  NNodes = x.size();
151 
152  int cell_type = ug->GetCell(0)->GetCellType();
153 
154  int nloc = -1;
155  int ndims = -1;
156  if(cell_type==VTK_TRIANGLE){
157  nloc = 3;
158  ndims = 2;
159  }else if(cell_type==VTK_TETRA){
160  nloc = 4;
161  ndims = 3;
162  }else{
163  std::cerr<<"ERROR("<<__FILE__<<"): unsupported element type\n";
164  exit(-1);
165  }
166 
167  std::vector<int> ENList;
168  for(size_t i=0;i<NElements;i++){
169  int ghost=0;
170  if(ug->GetCellData()->GetArray("vtkGhostLevels")!=NULL)
171  ghost = ug->GetCellData()->GetArray("vtkGhostLevels")->GetTuple1(i);
172 
173  if(ghost>0)
174  continue;
175 
176  vtkCell *cell = ug->GetCell(i);
177  assert(cell->GetCellType()==cell_type);
178  for(int j=0;j<nloc;j++){
179  if(have_global_ids)
180  ENList.push_back(renumber[cell->GetPointId(j)]);
181  else
182  ENList.push_back(cell->GetPointId(j));
183  }
184  }
185  NElements = ENList.size()/nloc;
186 
187  int nparts=1;
188  Mesh<real_t> *mesh=NULL;
189 
190  // Handle mpi parallel run.
191  MPI_Comm_size(MPI_COMM_WORLD, &nparts);
192 
193  if(nparts>1){
194  std::vector<index_t> owner_range;
195  std::vector<index_t> lnn2gnn;
196  std::vector<int> node_owner;
197 
198  std::vector<int> epart(NElements, 0), npart(NNodes, 0);
199  int rank;
200  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
201 
202  if(rank==0){
203  int edgecut;
204 
205  std::vector<int> eind(NElements*nloc);
206  eind[0] = 0;
207  for(size_t i=0;i<NElements*nloc;i++)
208  eind[i] = ENList[i];
209 
210  int intNElements = NElements;
211  int intNNodes = NNodes;
212 
213 #ifdef METIS_VER_MAJOR
214  int vsize = nloc - 1;
215  std::vector<int> eptr(NElements+1);
216  for(size_t i=0;i<NElements;i++)
217  eptr[i+1] = eptr[i]+nloc;
218  METIS_PartMeshNodal(&intNElements,
219  &intNNodes,
220  &(eptr[0]),
221  &(eind[0]),
222  NULL,
223  &vsize,
224  &nparts,
225  NULL,
226  NULL,
227  &edgecut,
228  &(epart[0]),
229  &(npart[0]));
230 #else
231  std::vector<int> etype(NElements);
232  for(size_t i=0;i<NElements;i++)
233  etype[i] = ndims-1;
234  int numflag = 0;
235  METIS_PartMeshNodal(&intNElements,
236  &intNNodes,
237  &(eind[0]),
238  &(etype[0]),
239  &numflag,
240  &nparts,
241  &edgecut,
242  &(epart[0]),
243  &(npart[0]));
244 #endif
245  }
246 
247  mpi_type_wrapper<index_t> mpi_index_t_wrapper;
248  MPI_Datatype MPI_INDEX_T = mpi_index_t_wrapper.mpi_type;
249 
250  MPI_Bcast(&(epart[0]), NElements, MPI_INDEX_T, 0, MPI_COMM_WORLD);
251  MPI_Bcast(&(npart[0]), NNodes, MPI_INDEX_T, 0, MPI_COMM_WORLD);
252 
253  // Separate out owned nodes.
254  std::vector< std::vector<index_t> > node_partition(nparts);
255  for(size_t i=0;i<NNodes;i++)
256  node_partition[npart[i]].push_back(i);
257 
258 #ifdef HAVE_BOOST_UNORDERED_MAP_HPP
259  boost::unordered_map<index_t, index_t> renumber;
260 #else
261  std::map<index_t, index_t> renumber;
262 #endif
263  {
264  index_t pos=0;
265  owner_range.push_back(0);
266  for(int i=0;i<nparts;i++){
267  int pNNodes = node_partition[i].size();
268  owner_range.push_back(owner_range[i]+pNNodes);
269  for(int j=0;j<pNNodes;j++)
270  renumber[node_partition[i][j]] = pos++;
271  }
272  }
273  std::vector<index_t> element_partition;
274  std::set<index_t> halo_nodes;
275  for(size_t i=0;i<NElements;i++){
276  std::set<index_t> residency;
277  for(int j=0;j<nloc;j++)
278  residency.insert(npart[ENList[i*nloc+j]]);
279 
280  if(residency.count(rank)){
281  element_partition.push_back(i);
282 
283  for(int j=0;j<nloc;j++){
284  index_t nid = ENList[i*nloc+j];
285  if(npart[nid]!=rank)
286  halo_nodes.insert(nid);
287  }
288  }
289  }
290 
291  // Append halo nodes to local node partition.
292  for(typename std::set<index_t>::const_iterator it=halo_nodes.begin();it!=halo_nodes.end();++it){
293  node_partition[rank].push_back(*it);
294  }
295 
296  // Global numbering to partition numbering look up table.
297  NNodes = node_partition[rank].size();
298  lnn2gnn.resize(NNodes);
299  node_owner.resize(NNodes);
300  for(size_t i=0;i<NNodes;i++){
301  index_t nid = node_partition[rank][i];
302  index_t gnn = renumber[nid];
303  lnn2gnn[i] = gnn;
304  node_owner[i] = npart[nid];
305  }
306 
307  // Construct local mesh.
308  std::vector<real_t> lx(NNodes), ly(NNodes), lz(NNodes);
309  if(ndims==2){
310  for(size_t i=0;i<NNodes;i++){
311  lx[i] = x[node_partition[rank][i]];
312  ly[i] = y[node_partition[rank][i]];
313  }
314  }else{
315  for(size_t i=0;i<NNodes;i++){
316  lx[i] = x[node_partition[rank][i]];
317  ly[i] = y[node_partition[rank][i]];
318  lz[i] = z[node_partition[rank][i]];
319  }
320  }
321 
322  NElements = element_partition.size();
323  std::vector<index_t> lENList(NElements*nloc);
324  for(size_t i=0;i<NElements;i++){
325  for(int j=0;j<nloc;j++){
326  index_t nid = renumber[ENList[element_partition[i]*nloc+j]];
327  lENList[i*nloc+j] = nid;
328  }
329  }
330 
331  // Swap
332  x.swap(lx);
333  y.swap(ly);
334  if(ndims==3)
335  z.swap(lz);
336  ENList.swap(lENList);
337 
338  MPI_Comm comm = MPI_COMM_WORLD;
339 
340  if(ndims==2)
341  mesh = new Mesh<real_t>(NNodes, NElements, &(ENList[0]), &(x[0]), &(y[0]), &(lnn2gnn[0]), &(owner_range[0]), comm);
342  else
343  mesh = new Mesh<real_t>(NNodes, NElements, &(ENList[0]), &(x[0]), &(y[0]), &(z[0]), &(lnn2gnn[0]), &(owner_range[0]), comm);
344  }
345 
346  if(nparts==1){ // If nparts!=1, then the mesh has been created already by the code a few lines above.
347  if(ndims==2)
348  mesh = new Mesh<real_t>(NNodes, NElements, &(ENList[0]), &(x[0]), &(y[0]));
349  else
350  mesh = new Mesh<real_t>(NNodes, NElements, &(ENList[0]), &(x[0]), &(y[0]), &(z[0]));
351  }
352 
353  return mesh;
354  }
355 
356  static void export_vtu(const char *basename, const Mesh<real_t> *mesh, const real_t *psi=NULL){
357  size_t NElements = mesh->get_number_elements();
358  size_t ndims = mesh->get_number_dimensions();
359 
360  // Set the orientation of elements.
361  ElementProperty<real_t> *property = NULL;
362  for(size_t i=0;i<NElements;i++){
363  const int *n=mesh->get_element(i);
364  assert(n[0]>=0);
365 
366  if(ndims==2)
367  property = new ElementProperty<real_t>(mesh->get_coords(n[0]),
368  mesh->get_coords(n[1]),
369  mesh->get_coords(n[2]));
370  else
371  property = new ElementProperty<real_t>(mesh->get_coords(n[0]),
372  mesh->get_coords(n[1]),
373  mesh->get_coords(n[2]),
374  mesh->get_coords(n[3]));
375  break;
376  }
377 
378  // Create VTU object to write out.
379  vtkSmartPointer<vtkUnstructuredGrid> ug = vtkSmartPointer<vtkUnstructuredGrid>::New();
380 
381  vtkSmartPointer<vtkPoints> vtk_points = vtkSmartPointer<vtkPoints>::New();
382  size_t NNodes = mesh->get_number_nodes();
383  vtk_points->SetNumberOfPoints(NNodes);
384 
385  vtkSmartPointer<vtkDoubleArray> vtk_psi = NULL;
386 
387  if(psi!=NULL){
388  vtk_psi = vtkSmartPointer<vtkDoubleArray>::New();
389  vtk_psi->SetNumberOfComponents(1);
390  vtk_psi->SetNumberOfTuples(NNodes);
391  vtk_psi->SetName("psi");
392  }
393 
394  vtkSmartPointer<vtkIntArray> vtk_node_numbering = vtkSmartPointer<vtkIntArray>::New();
395  vtk_node_numbering->SetNumberOfComponents(1);
396  vtk_node_numbering->SetNumberOfTuples(NNodes);
397  vtk_node_numbering->SetName("nid");
398 
399  vtkSmartPointer<vtkDoubleArray> vtk_metric = vtkSmartPointer<vtkDoubleArray>::New();
400  vtk_metric->SetNumberOfComponents(ndims*ndims);
401  vtk_metric->SetNumberOfTuples(NNodes);
402  vtk_metric->SetName("Metric");
403 
404  vtkSmartPointer<vtkDoubleArray> vtk_edge_length = vtkSmartPointer<vtkDoubleArray>::New();
405  vtk_edge_length->SetNumberOfComponents(1);
406  vtk_edge_length->SetNumberOfTuples(NNodes);
407  vtk_edge_length->SetName("mean_edge_length");
408 
409  vtkSmartPointer<vtkDoubleArray> vtk_max_desired_length = vtkSmartPointer<vtkDoubleArray>::New();
410  vtk_max_desired_length->SetNumberOfComponents(1);
411  vtk_max_desired_length->SetNumberOfTuples(NNodes);
412  vtk_max_desired_length->SetName("max_desired_edge_length");
413 
414  vtkSmartPointer<vtkDoubleArray> vtk_min_desired_length = vtkSmartPointer<vtkDoubleArray>::New();
415  vtk_min_desired_length->SetNumberOfComponents(1);
416  vtk_min_desired_length->SetNumberOfTuples(NNodes);
417  vtk_min_desired_length->SetName("min_desired_edge_length");
418 
419 #pragma omp parallel for
420  for(size_t i=0;i<NNodes;i++){
421  const real_t *r = mesh->get_coords(i);
422  const double *m = mesh->get_metric(i);
423 
424  if(vtk_psi!=NULL)
425  vtk_psi->SetTuple1(i, psi[i]);
426  vtk_node_numbering->SetTuple1(i, i);
427  if(ndims==2){
428  vtk_points->SetPoint(i, r[0], r[1], 0.0);
429  vtk_metric->SetTuple4(i,
430  m[0], m[1],
431  m[1], m[2]);
432  }else{
433  vtk_points->SetPoint(i, r[0], r[1], r[2]);
434  vtk_metric->SetTuple9(i,
435  m[0], m[1], m[2],
436  m[1], m[3], m[4],
437  m[2], m[4], m[5]);
438  }
439  int nedges=mesh->NNList[i].size();
440  real_t mean_edge_length=0;
441  real_t max_desired_edge_length=0;
442  real_t min_desired_edge_length=DBL_MAX;
443 
444  if(ndims==2)
445  for(typename std::vector<index_t>::const_iterator it=mesh->NNList[i].begin();it!=mesh->NNList[i].end();++it){
446  double length = mesh->calc_edge_length(i, *it);
447  mean_edge_length += length;
448 
449  MetricTensor<double,2> M(m);
450 
451  max_desired_edge_length = std::max(max_desired_edge_length, M.max_length());
452  min_desired_edge_length = std::min(min_desired_edge_length, M.min_length());
453  }
454  else if(ndims==3)
455  for(typename std::vector<index_t>::const_iterator it=mesh->NNList[i].begin();it!=mesh->NNList[i].end();++it){
456  double length = mesh->calc_edge_length(i, *it);
457  mean_edge_length += length;
458 
459  MetricTensor<double,3> M(m);
460 
461  max_desired_edge_length = std::max(max_desired_edge_length, M.max_length());
462  min_desired_edge_length = std::min(min_desired_edge_length, M.min_length());
463  }
464 
465 
466  mean_edge_length/=nedges;
467  vtk_edge_length->SetTuple1(i, mean_edge_length);
468  vtk_max_desired_length->SetTuple1(i, max_desired_edge_length);
469  vtk_min_desired_length->SetTuple1(i, min_desired_edge_length);
470  }
471 
472  ug->SetPoints(vtk_points);
473  if(vtk_psi!=NULL){
474  ug->GetPointData()->AddArray(vtk_psi);
475  }
476  ug->GetPointData()->AddArray(vtk_node_numbering);
477  ug->GetPointData()->AddArray(vtk_metric);
478  ug->GetPointData()->AddArray(vtk_edge_length);
479  ug->GetPointData()->AddArray(vtk_max_desired_length);
480  ug->GetPointData()->AddArray(vtk_min_desired_length);
481 
482  // Create a point data array to illustrate the boundary nodes.
483  vtkSmartPointer<vtkIntArray> vtk_boundary_nodes = vtkSmartPointer<vtkIntArray>::New();
484  vtk_boundary_nodes->SetNumberOfComponents(1);
485  vtk_boundary_nodes->SetNumberOfTuples(NNodes);
486  vtk_boundary_nodes->SetName("BoundaryNodes");
487  for(size_t i=0;i<NNodes;i++)
488  vtk_boundary_nodes->SetTuple1(i, -1);
489 
490  for(size_t i=0;i<NElements;i++){
491  const int *n=mesh->get_element(i);
492  if(n[0]==-1)
493  continue;
494 
495  if(ndims==2){
496  for(int j=0;j<3;j++){
497  if(mesh->boundary[i*3+j]>0){
498  vtk_boundary_nodes->SetTuple1(n[(j+1)%3], mesh->boundary[i*3+j]);
499  vtk_boundary_nodes->SetTuple1(n[(j+2)%3], mesh->boundary[i*3+j]);
500  }
501  }
502  }else{
503  for(int j=0;j<4;j++){
504  if(mesh->boundary[i*4+j]>0){
505  vtk_boundary_nodes->SetTuple1(n[(j+1)%4], mesh->boundary[i*4+j]);
506  vtk_boundary_nodes->SetTuple1(n[(j+2)%4], mesh->boundary[i*4+j]);
507  vtk_boundary_nodes->SetTuple1(n[(j+3)%4], mesh->boundary[i*4+j]);
508  }
509  }
510  }
511  }
512  ug->GetPointData()->AddArray(vtk_boundary_nodes);
513 
514  vtkSmartPointer<vtkIntArray> vtk_cell_numbering = vtkSmartPointer<vtkIntArray>::New();
515  vtk_cell_numbering->SetNumberOfComponents(1);
516  vtk_cell_numbering->SetNumberOfTuples(NElements);
517  vtk_cell_numbering->SetName("eid");
518 
519  vtkSmartPointer<vtkIntArray> vtk_boundary = vtkSmartPointer<vtkIntArray>::New();
520  if(ndims==2)
521  vtk_boundary->SetNumberOfComponents(3);
522  else
523  vtk_boundary->SetNumberOfComponents(4);
524  vtk_boundary->SetNumberOfTuples(NElements);
525  vtk_boundary->SetName("Boundary");
526 
527  vtkSmartPointer<vtkDoubleArray> vtk_quality = vtkSmartPointer<vtkDoubleArray>::New();
528  vtk_quality->SetNumberOfComponents(1);
529  vtk_quality->SetNumberOfTuples(NElements);
530  vtk_quality->SetName("quality");
531 
532  for(size_t i=0, k=0;i<NElements;i++){
533  const index_t *n = mesh->get_element(i);
534  assert(n[0]>=0);
535 
536  if(ndims==2){
537  vtkIdType pts[] = {n[0], n[1], n[2]};
538  ug->InsertNextCell(VTK_TRIANGLE, 3, pts);
539  vtk_boundary->SetTuple3(i, mesh->boundary[i*3], mesh->boundary[i*3+1], mesh->boundary[i*3+2]);
540 
541  vtk_quality->SetTuple1(k, property->lipnikov(mesh->get_coords(n[0]), mesh->get_coords(n[1]), mesh->get_coords(n[2]),
542  mesh->get_metric(n[0]), mesh->get_metric(n[1]), mesh->get_metric(n[2])));
543  }else{
544  vtkIdType pts[] = {n[0], n[1], n[2], n[3]};
545  ug->InsertNextCell(VTK_TETRA, 4, pts);
546  vtk_boundary->SetTuple4(i, mesh->boundary[i*4], mesh->boundary[i*4+1], mesh->boundary[i*4+2], mesh->boundary[i*4+3]);
547 
548  vtk_quality->SetTuple1(k, property->lipnikov(mesh->get_coords(n[0]), mesh->get_coords(n[1]), mesh->get_coords(n[2]), mesh->get_coords(n[3]),
549  mesh->get_metric(n[0]), mesh->get_metric(n[1]), mesh->get_metric(n[2]), mesh->get_metric(n[3])));
550  }
551 
552  vtk_cell_numbering->SetTuple1(k, i);
553 
554  k++;
555  }
556 
557  ug->GetCellData()->AddArray(vtk_cell_numbering);
558  ug->GetCellData()->AddArray(vtk_quality);
559  ug->GetCellData()->AddArray(vtk_boundary);
560 
561  int rank=0, nparts=1;
562 #ifdef HAVE_MPI
563  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
564  MPI_Comm_size(MPI_COMM_WORLD, &nparts);
565 #endif
566 
567  if(nparts==1){
568  vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
569  std::string filename = std::string(basename)+std::string(".vtu");
570  writer->SetFileName(filename.c_str());
571 #if VTK_MAJOR_VERSION < 6
572  writer->SetInput(ug);
573 #else
574  writer->SetInputData(ug);
575 #endif
576  writer->Write();
577  }
578 #ifdef HAVE_MPI
579  else{
580  // Set ghost levels
581  vtkSmartPointer<vtkUnsignedCharArray> vtk_ghost = vtkSmartPointer<vtkUnsignedCharArray>::New();
582  vtk_ghost->SetNumberOfComponents(1);
583  vtk_ghost->SetNumberOfTuples(NElements);
584  vtk_ghost->SetName("vtkGhostLevels");
585 
586  for(size_t i=0;i<NElements;i++){
587  const index_t *n = mesh->get_element(i);
588  int owner;
589  if(ndims==2)
590  owner=std::min(mesh->node_owner[n[0]], std::min(mesh->node_owner[n[1]], mesh->node_owner[n[2]]));
591  else
592  owner=std::min(std::min(mesh->node_owner[n[0]], mesh->node_owner[n[1]]), std::min(mesh->node_owner[n[2]], mesh->node_owner[n[3]]));
593 
594  if(owner==rank){
595  vtk_ghost->SetTuple1(i, 0);
596  }else{
597  vtk_ghost->SetTuple1(i, 1);
598  }
599  }
600  ug->GetCellData()->AddArray(vtk_ghost);
601 
602  // Set GlobalIds
603  vtkSmartPointer<vtkIdTypeArray> vtk_gnn = vtkSmartPointer<vtkIdTypeArray>::New();
604  vtk_gnn->SetNumberOfComponents(1);
605  vtk_gnn->SetNumberOfTuples(NNodes);
606  vtk_gnn->SetName("GlobalId");
607 
608  for(size_t i=0;i<NNodes;i++){
609  vtk_gnn->SetTuple1(i, mesh->lnn2gnn[i]);
610  }
611  // ug->GetPointData()->AddArray(vtk_gnn);
612  // ug->GetPointData()->SetActiveGlobalIds("GlobalId");
613  ug->GetPointData()->SetGlobalIds(vtk_gnn);
614 
615  vtkSmartPointer<vtkXMLPUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
616  std::string filename = std::string(basename)+std::string(".pvtu");
617  writer->SetFileName(filename.c_str());
618  writer->SetNumberOfPieces(nparts);
619  writer->SetGhostLevel(1);
620  writer->SetStartPiece(rank);
621  writer->SetEndPiece(rank);
622 #if VTK_MAJOR_VERSION < 6
623  writer->SetInput(ug);
624 #else
625  writer->SetInputData(ug);
626 #endif
627  writer->Write();
628  }
629 #endif
630 
631  delete property;
632 
633  return;
634  }
635 
636 };
637 #endif
VTKTools::import_vtu
static Mesh< real_t > * import_vtu(std::string filename)
Definition: VTKTools.h:97
MetricTensor
symmetric metric tensor class.
Definition: MetricTensor.h:48
MetricTensor::max_length
treal_t max_length() const
Definition: MetricTensor.h:392
Mesh::get_coords
const real_t * get_coords(index_t nid) const
Return positions vector.
Definition: Mesh.h:384
index_t
int index_t
Definition: PragmaticTypes.h:41
ElementProperty
Calculates a number of element properties.
Definition: ElementProperty.h:65
Mesh
Manages mesh data.
Definition: Mesh.h:70
generate_Steiner_ellipse_3d.M
tuple M
Definition: generate_Steiner_ellipse_3d.py:15
VTKTools::export_vtu
static void export_vtu(const char *basename, const Mesh< real_t > *mesh, const real_t *psi=NULL)
Definition: VTKTools.h:356
Mesh::calc_edge_length
real_t calc_edge_length(index_t nid0, index_t nid1) const
Calculates the edge lengths in metric space.
Definition: Mesh.h:911
Mesh::get_number_elements
size_t get_number_elements() const
Return the number of elements in the mesh.
Definition: Mesh.h:374
mpi_tools.h
ellipse_convergence.mesh
tuple mesh
Definition: ellipse_convergence.py:108
vtkFloatingPointType
#define vtkFloatingPointType
Definition: VTKTools.h:66
mpi_type_wrapper
Definition: mpi_tools.h:45
mpi_type_wrapper::mpi_type
const MPI_Datatype mpi_type
Definition: mpi_tools.h:46
MetricTensor::min_length
double min_length() const
Definition: MetricTensor.h:408
Mesh::get_metric
const double * get_metric(index_t nid) const
Return metric at that vertex.
Definition: Mesh.h:396
Mesh::get_number_dimensions
size_t get_number_dimensions() const
Return the number of spatial dimensions.
Definition: Mesh.h:379
VTKTools
Toolkit for importing and exporting VTK files. This is mostly used as part of the internal unit testi...
Definition: VTKTools.h:95
Mesh::get_number_nodes
size_t get_number_nodes() const
Return the number of nodes in the mesh.
Definition: Mesh.h:369
Mesh::get_element
const index_t * get_element(size_t eid) const
Return a pointer to the element-node list.
Definition: Mesh.h:358
Generated on Mon Feb 23 2015 22:50:44 for PRAgMaTIc by   doxygen 1.8.7