source: XMLIO_V2/dev/dev_rv/src4/xmlio/vtk/lscereader.cpp @ 241

Last change on this file since 241 was 241, checked in by hozdoba, 11 years ago
File size: 38.0 KB
Line 
1
2#ifndef LSCE_EXPORTS
3#       include "lscereader.hpp"
4#else
5#       include "vtkLSCEReader.h"
6#endif //LSCE_EXPORTS
7
8
9// BibliothÚque VTK
10#include "vtkObjectFactory.h"
11#include "vtkInformation.h"
12#include "vtkInformationVector.h"
13#include "vtkStreamingDemandDrivenPipeline.h"
14
15#include "vtkPointLocator.h"
16
17#include <vtkProperty.h>
18#include <vtkDataSetMapper.h>
19#include <vtkActor.h>
20#include <vtkRenderWindow.h>
21#include <vtkRenderer.h>
22#include <vtkRenderWindowInteractor.h>
23#include <vtkLegendScaleActor.h>
24
25#include "vtkMath.h"
26
27#include "inetcdf4_impl.hpp"
28#include "inetcdf4_adv_impl.hpp"
29
30// BibliothÚque C++
31#include <algorithm>
32
33// Type Fortran
34#define ARRAY(valuetype, numdims) boost::shared_ptr<boost::multi_array<valuetype, numdims> >
35
36#define ARRAY_ASSIGN(value, valuetype, numdims, extent)\
37   value.reset(new boost::multi_array<valuetype, numdims>(boost::extents extent))
38
39#define ARRAY_CREATE(value, valuetype, numdims, extent)\
40   ARRAY(valuetype, numdims) value =                   \
41   ARRAY(valuetype, numdims)(new boost::multi_array<valuetype, numdims>(boost::extents extent))
42
43// Type C
44#define ARRAY_C_ASSIGN(value, valuetype, numdims, extent)\
45   value = ARRAY(valuetype, numdims)                     \
46   (new boost::multi_array<valuetype, numdims>(boost::extents extent, c_storage_order()))
47
48#define ARRAY_C_CREATE(value, valuetype, numdims, extent)\
49   ARRAY_C_ASSIGN(ARRAY(valuetype, numdims) value, valuetype, numdims, extent)
50
51
52#ifndef LSCE_EXPORTS
53namespace xmlioserver
54{
55   namespace vtk
56   {
57#endif //LSCE_EXPORTS
58       
59      /// ////////////////////// Définitions ////////////////////// ///
60      vtkStandardNewMacro(vtkLSCEReader);
61
62      vtkLSCEReader::vtkLSCEReader(void)
63         : FileName(), CurGridType(RECTILINEAR), VarNames()
64         , A3D(true), ATemporal(true), ACell(true), AGridDef(false)
65         , Reader()
66      {
67         this->SetNumberOfInputPorts (0);
68         this->SetNumberOfOutputPorts(1);
69      }
70
71      vtkLSCEReader::~vtkLSCEReader(void)
72      {
73
74      }
75
76      ///---------------------------------------------------------------
77
78      void vtkLSCEReader::SetFileName(const vtkStdString & fileName)
79      {
80         this->FileName = fileName;
81         this->Reader   = boost::shared_ptr<xmlioserver::io::CINetCDF4Adv>
82                                       (new xmlioserver::io::CINetCDF4Adv(fileName));
83         this->Modified();
84      }
85     
86      void vtkLSCEReader::SetVariable(const vtkStdString & variable)
87      {
88         this->RemoveAllSelectedVariables();
89         if (Reader->hasVariable(variable))
90            this->AddVariableToSelection(variable);
91         this->Modified();
92      }
93
94      const char * vtkLSCEReader::GetFileName(void) const
95      {
96         return (this->FileName.c_str());
97      }
98
99      void vtkLSCEReader::PrintSelf(ostream& os, vtkIndent indent)
100      {
101        this->Superclass::PrintSelf(os,indent);
102      }
103
104      //----------------------------------------------------------------
105      void vtkLSCEReader::AddVariableToSelection(const vtkStdString & varName)
106      {
107         if (this->Reader->isRectilinear(varName))
108            this->SetGridType (RECTILINEAR);
109         else if (this->Reader->isCurvilinear(varName))
110            this->SetGridType (CURVILINEAR);
111         else
112            this->SetGridType (UNSTRUCTURED);
113
114         this->VarNames.insert(varName);
115      }
116
117      void vtkLSCEReader::RemoveSelectedVariable(const vtkStdString & varName)
118      {
119         std::set<vtkStdString>::iterator it = this->VarNames.find(varName);
120         if (it != this->VarNames.end()) this->VarNames.erase(it);
121      }
122
123      void vtkLSCEReader::RemoveAllSelectedVariables(void)
124      {
125         this->VarNames.clear();
126      }
127
128      const std::set<vtkStdString> &
129            vtkLSCEReader::GetSelectedVariables(void) const
130      {
131         return (this->VarNames);
132      }
133
134      bool vtkLSCEReader::HasSelectedVariable(void) const
135      {
136         return (this->VarNames.size() != 0);
137      }
138
139      //----------------------------------------------------------------
140
141      void vtkLSCEReader::SetGridType (GridType type)
142      {
143         if ((this->IsUnstructured() && (type == CURVILINEAR)) ||
144             (this->IsUnstructured() && (type == RECTILINEAR)) ||
145             (this->IsCurvilinear()  && (type == RECTILINEAR)))
146            this->RemoveAllSelectedVariables();
147         this->CurGridType = type;
148      }
149
150      bool vtkLSCEReader::IsUnstructured(void) const
151      {
152         return (this->CurGridType == UNSTRUCTURED);
153      }
154
155      bool vtkLSCEReader::IsCurvilinear(void) const
156      {
157         return (this->CurGridType == CURVILINEAR);
158      }
159
160      bool vtkLSCEReader::IsRectilinear(void) const
161      {
162         return (this->CurGridType == RECTILINEAR);
163      }
164
165      //----------------------------------------------------------------
166
167      void vtkLSCEReader::AcceptTemporalOnly(bool value)
168      {
169         this->ATemporal = value;
170      }
171
172      void vtkLSCEReader::Accept3DOnly(bool value)
173      {
174         this->A3D = value;
175      }
176
177      void vtkLSCEReader::AcceptCellOnly(bool value)
178      {
179         this->ACell = value;
180      }
181
182      //----------------------------------------------------------------
183
184      int vtkLSCEReader::ProcessRequest
185                           (vtkInformation       *  request,
186                            vtkInformationVector ** vtkNotUsed(inputVector),
187                            vtkInformationVector *  outputVector)
188      {
189         try
190         {
191         
192            if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
193            {
194               return (this->RequestData(NULL, NULL, outputVector));
195            }
196
197            if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT()))
198            {
199            //   return (this->RequestDataObject(NULL, NULL, outputVector));
200            }
201
202            if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
203            {
204                 return (this->RequestInformation(NULL, NULL, outputVector));
205            }
206
207            if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
208            {
209               return (this->RequestUpdateExtent(NULL, NULL, outputVector));
210            }
211         }
212         catch (xmlioserver::CException & _exception)
213         {
214            std::cerr << _exception.getMessage() << std::endl;
215            return (-1);
216         }
217         return (1);
218      }
219
220      //----------------------------------------------------------------
221
222      void vtkLSCEReader::ShowVariable(const vtkStdString & filename,
223                                       const vtkStdString & varname)
224      {
225         vtkSmartPointer<vtkLSCEReader> lscereader =
226               vtkSmartPointer<vtkLSCEReader>::New();
227
228         lscereader->SetFileName(filename);
229         lscereader->AddVariableToSelection(varname);
230         lscereader->Update();
231
232         // Création d'un mappeur et d'un acteur.
233         vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New();
234         mapper->SetInputConnection(lscereader->GetOutput()->GetProducerPort());
235         mapper->SetScalarRange  (1000, 10000);
236
237         vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
238         actor->SetMapper(mapper);
239         actor->GetProperty()->SetRepresentationToWireframe();
240
241         // Création d'une fenêtre de rendu.
242         vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
243         vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
244         renderWindow->AddRenderer(renderer);
245         vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
246            vtkSmartPointer<vtkRenderWindowInteractor>::New();
247         renderWindowInteractor->SetRenderWindow(renderWindow);
248
249
250         vtkSmartPointer<vtkLegendScaleActor> legendScaleActor =
251            vtkSmartPointer<vtkLegendScaleActor>::New();
252
253         // Ajout de l'acteur à la scÚne.
254         renderer->AddActor(actor);
255         renderer->AddActor(legendScaleActor);
256
257         renderer->GradientBackgroundOn();
258         renderer->SetBackground(1,1,1);
259         renderer->SetBackground2(0,0,0);
260
261         // Rendu et intéraction
262         renderWindow->Render();
263         renderWindowInteractor->Start();
264      }
265
266      //----------------------------------------------------------------
267
268      void vtkLSCEReader::CreateRectilinearGrid(vtkRectilinearGrid * grid,
269                                                vtkInformation     * outInfo,
270                                                vtkFloatArray      * xspacing,
271                                                vtkFloatArray      * yspacing,
272                                                vtkFloatArray      * zspacing,
273                                                vtkIntArray        * dimensions)
274      {
275            int extent[6] = { 0, dimensions->GetValue(0) - 1
276                            , 0, dimensions->GetValue(1) - 1
277                            , 0, dimensions->GetValue(2) - 1};
278            outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
279            grid->SetExtent (extent);
280            grid->SetDimensions (dimensions->GetValue(0),
281                                 dimensions->GetValue(1),
282                                 dimensions->GetValue(2));
283            grid->SetXCoordinates (xspacing);
284            grid->SetYCoordinates (yspacing);
285            grid->SetZCoordinates (zspacing);
286      }
287
288      void vtkLSCEReader::CreateStructuredGrid(vtkStructuredGrid * grid,
289                                               vtkInformation    * outInfo,
290                                               vtkPoints         * points,
291                                               vtkIntArray       * dimensions)
292      {
293            int extent[6] = { 0, dimensions->GetValue(0) - 1
294                            , 0, dimensions->GetValue(1) - 1
295                            , 0, dimensions->GetValue(2) - 1};
296            outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
297            grid->SetExtent (extent);
298            grid->SetDimensions (dimensions->GetValue(0),
299                                 dimensions->GetValue(1),
300                                 dimensions->GetValue(2));
301            grid->SetPoints(points);
302      }
303
304      void vtkLSCEReader::CreateUnstructuredGrid(vtkUnstructuredGrid * grid,
305                                                 vtkInformation      * outInfo,
306                                                 vtkPoints           * points,
307                                                 vtkCellArray        * cells,
308                                                 int                   celltype)
309      {
310         grid->SetPoints(points);
311         grid->SetCells(celltype, cells);
312      }
313
314      //----------------------------------------------------------------
315
316      void vtkLSCEReader::GetRectilinearConnectivity(int dimZ, int dimY, int dimX, vtkCellArray * cells)
317      {
318         if (dimZ == 1 && dimY >  1 && dimX >  1) // 2D
319         {
320            vtkIdType pts[4] ;
321
322            for (int i = 0, k = 0; i < (dimX - 1); i++)
323               for (int j = 0; j < (dimY - 1); j++, k++)
324            {
325               pts[0] =  i * dimY + j;
326               pts[1] =  i * dimY + j + 1;
327               pts[2] =  (i+1) * dimY + j + 1;
328               pts[3] =  (i+1) * dimY + j;
329               cells->InsertNextCell(4,  pts);
330            }
331         }
332         else if (dimZ > 1 && dimY >  1 && dimX >  1) // 3D
333         {
334            vtkIdType pts[8];
335            for (int i = 0, l = 0; i < (dimX - 1); i++)
336               for (int j = 0; j < (dimY - 1); j++)
337                  for (int k = 0; k < (dimZ - 1); k++, l++)
338            {
339               pts[0] =  i * dimY * dimZ + j * dimZ + k;
340               pts[1] =  i * dimY * dimZ + j * dimZ + k + 1;
341               pts[2] =  (i+1) * dimY * dimZ + j * dimZ + k + 1;
342               pts[3] =  (i+1) * dimY * dimZ + j * dimZ + k;
343
344               pts[4] =  i * dimY * dimZ  + (j + 1) * dimZ + k ;
345               pts[5] =  i * dimY * dimZ  + (j + 1) * dimZ + k +1;
346               pts[6] =  (i+1) * dimY * dimZ + (j + 1) * dimZ + k + 1;
347               pts[7] =  (i+1) * dimY * dimZ + (j + 1) * dimZ + k ;
348               cells->InsertNextCell(8,  pts);
349            }
350         }
351      }
352
353      //----------------------------------------------------------------
354
355      void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
356                                            int yi, int yf,
357                                            int zi, int zf,
358                                            vtkFloatArray * xspacing,
359                                            vtkFloatArray * yspacing,
360                                            vtkFloatArray * zspacing,
361                                            vtkIntArray   * dimensions)
362      {
363         dimensions->InsertNextValue (xf - xi + 1); // X
364         dimensions->InsertNextValue (yf - yi + 1); // Y
365         dimensions->InsertNextValue (zf - zi + 1); // Z
366         for (int i = xi; i <= xf; i++)
367            xspacing->InsertNextValue (i);
368         for (int j = yi; j <= yf; j++)
369            yspacing->InsertNextValue (j);
370         for (int k = zi; k <= zf; k++)
371            zspacing->InsertNextValue (k);
372      }
373
374      void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
375                                            int yi, int yf,
376                                            int zi, int zf,
377                                            vtkPoints   * points,
378                                            vtkIntArray * dimensions)
379      {
380         dimensions->InsertNextValue (xf - xi + 1); // X
381         dimensions->InsertNextValue (yf - yi + 1); // Y
382         dimensions->InsertNextValue (zf - zi + 1); // Z
383         for (int i = xi; i <= xf; i++)
384            for (int j = yi; j <= yf; j++)
385               for (int k = zi; k <= zf; k++)
386                  points->InsertNextPoint(i, j, k);
387      }
388
389      void vtkLSCEReader::CreateSimpleGrid (int xi, int xf,
390                                            int yi, int yf,
391                                            int zi, int zf,
392                                            vtkPoints    * points,
393                                            vtkCellArray * cells,
394                                            vtkIntArray  * dimensions)
395      {
396         this->CreateSimpleGrid (xi, xf, yi, yf, zi, zf, points, dimensions);
397         int dimX = dimensions->GetValue(0),
398             dimY = dimensions->GetValue(1),
399             dimZ = dimensions->GetValue(2);
400         this->GetRectilinearConnectivity(dimZ, dimY, dimX, cells);
401      }
402
403      //----------------------------------------------------------------
404      void vtkLSCEReader::GetSpacings(const vtkStdString & coordinate,
405                                      bool bounds, vtkFloatArray * spacing)
406      {
407         std::size_t size = (Reader->getDimensions(&coordinate).begin())->second;
408         ARRAY_CREATE(data_arr, float, 1, [1]);
409         if (Reader->hasVariable(coordinate))
410         {
411            if (bounds)
412            {
413               if (Reader->hasBounds(coordinate))
414               { // CellData, variable de coordonnée avec attributes bounds.
415                  vtkStdString boundid = Reader->getBoundsId(coordinate);
416                  float f;
417                  Reader->readData(*data_arr, boundid);
418                  spacing->InsertNextValue (f = (*data_arr)[0]);
419                  for (std::size_t j = 1; j < (2 * size); j++) // << Super laid
420                     if (f != (*data_arr)[j])
421                        spacing->InsertNextValue (f = (*data_arr)[j]);
422               }
423               else
424               { // CellData, variable de coordonnée sans attributes bounds.
425                  Reader->readData(*data_arr, coordinate);
426                  spacing->InsertNextValue ((*data_arr)[0] - ((*data_arr)[1] - (*data_arr)[0])/2);
427                  for (std::size_t j = 0; j < (size-1); j++)
428                     spacing->InsertNextValue ((*data_arr)[j] + ((*data_arr)[j+1]-(*data_arr)[j])/2);
429                  spacing->InsertNextValue ((*data_arr)[size-2] + ((*data_arr)[size-1] - (*data_arr)[size-2])/2);
430               }
431            }
432            else
433            { // PointData, variable de coordonnée.
434               Reader->readData(*data_arr, coordinate);
435               for (std::size_t j = 0; j < (size-1); j++)
436                  spacing->InsertNextValue ((*data_arr)[j]);
437            }
438         }
439         else
440         {
441            if (bounds) // CellData, dimension de coordonnée.
442               for (float i = -0.5; i < ((float)size + 1.); i = i + 1.)
443                  spacing->InsertNextValue (i);
444            else // PointData, dimension de coordonnée.
445               for (float j = 0.; j < (float)size; j = j + 1.)
446                  spacing->InsertNextValue (j);
447         }
448      }
449      //----------------------------------------------------------------
450
451      void vtkLSCEReader::AddPoint(vtkPoints * points, float *  value, bool proj)
452      {
453         if (proj)
454         {
455            double h = 1.0;
456            double cartesianCoord[3];
457            value[0] =  0.017453292519943295 * value[0];
458            value[1] =  0.017453292519943295 * value[1];
459
460            cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
461            cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
462            cartesianCoord[2] = h * sin(value[1]);
463            points->InsertNextPoint (cartesianCoord);
464         }
465         else
466         {
467            points->InsertNextPoint (value);
468         }
469      }
470
471      //----------------------------------------------------------------
472     
473      void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
474                                        const vtkStdString & ycoordinate,
475                                        const vtkStdString & zcoordinate,
476                                        vtkIntArray * dimensions, bool bounds)
477      {
478         int a = (bounds) ? 1 : 0;
479         dimensions->InsertNextValue
480            ((Reader->getDimensions(&xcoordinate).begin())->second + a);
481         dimensions->InsertNextValue
482            ((Reader->getDimensions(&ycoordinate).begin())->second + a);
483         if (zcoordinate.size() > 0)
484         {
485            dimensions->InsertNextValue
486               ((Reader->getDimensions(&zcoordinate).begin())->second + a);
487         }
488         else
489         {
490            dimensions->InsertNextValue(1);
491         }
492      }
493
494      //----------------------------------------------------------------
495     
496      void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate,
497                                        const vtkStdString & zcoordinate,
498                                        vtkIntArray * dimensions, bool bounds)
499      {
500         int a = (bounds) ? 1 : 0;
501         dimensions->InsertNextValue
502            (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)] + a);
503         dimensions->InsertNextValue
504            (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)] + a);
505         if (zcoordinate.size() > 0)
506         {
507            dimensions->InsertNextValue
508               ((Reader->getDimensions(&zcoordinate).begin())->second + a);
509         }
510         else
511         {
512            dimensions->InsertNextValue(1);
513         }
514      }
515
516      //----------------------------------------------------------------
517      void vtkLSCEReader::GetPoints(const vtkStdString & xcoordinate,
518                                    const vtkStdString & ycoordinate,
519                                    const vtkStdString & zcoordinate,
520                                    bool bounds, bool proj,
521                                    vtkPoints * points, vtkIntArray * dimensions)
522      {
523         float value[3];
524         bool is2D = true;
525         ARRAY_CREATE(xvalue, float, 1, [1]);
526         ARRAY_CREATE(yvalue, float, 1, [1]);
527         vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New();
528         if (zcoordinate.size() > 0)
529         {
530            GetSpacings(zcoordinate, bounds, zspacing);
531            is2D = false;
532         }
533         else
534            zspacing->InsertNextValue(0);
535
536         dimensions->InsertNextValue
537            (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)]);
538         dimensions->InsertNextValue
539            (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)]);
540         dimensions->InsertNextValue(zspacing->GetNumberOfTuples());
541
542         if (!bounds)
543         {
544            Reader->readData(*xvalue, xcoordinate);
545            Reader->readData(*yvalue, ycoordinate);
546
547            for (int i = 0; i < dimensions->GetValue(0); i++)
548               for (int j = 0; j < dimensions->GetValue(1); j++)
549                  for (int k = 0; k < dimensions->GetValue(2); k++)
550                  {
551                     value[0] = (*xvalue)[i*dimensions->GetValue(1) + j];
552                     value[1] = (*yvalue)[i*dimensions->GetValue(1) + j];
553                     value[2] = zspacing->GetValue(k);
554                     AddPoint(points, value, proj);
555                  }
556            return;
557         }
558
559         dimensions->SetValue(0, dimensions->GetValue(0) + 1);
560         dimensions->SetValue(1, dimensions->GetValue(1) + 1);
561
562         if (Reader->hasBounds(xcoordinate))
563         {
564            vtkStdString boundid = Reader->getBoundsId(xcoordinate);
565            Reader->readData(*xvalue, boundid);
566         }
567         if (Reader->hasBounds(ycoordinate))
568         {
569            vtkStdString boundid = Reader->getBoundsId(ycoordinate);
570            Reader->readData(*yvalue, boundid);
571         }
572         int ydim = dimensions->GetValue(0) - 1 ,
573             xdim = dimensions->GetValue(1) - 1;
574
575         for (int i = 0; i <= xdim; i++)
576            for (int j = 0; j <= ydim; j++)
577               for (int k = 0; k < dimensions->GetValue(2); k++)
578               {
579                  if ((j != ydim) && (i != xdim))
580                  {
581                     value[0] = (*xvalue)[4 * (i * ydim + j) + 0];
582                     value[1] = (*yvalue)[4 * (i * ydim + j) + 0];
583                     value[2] = zspacing->GetValue(k);
584                     AddPoint(points, value, proj);
585
586                     continue;
587                  }
588                  else if ((j == ydim) && (i != xdim))
589                  {
590                     value[0] = (*xvalue)[4 * (i * ydim + j - 1) + 1];
591                     value[1] = (*yvalue)[4 * (i * ydim + j - 1) + 1];
592                     value[2] = zspacing->GetValue(k);
593                     AddPoint(points, value, proj);
594
595                     continue;
596                  }
597                  else if ((j != ydim) && (i == xdim))
598                  {
599                     value[0] = (*xvalue)[4 * ((i-1) * ydim + j) + 3];
600                     value[1] = (*yvalue)[4 * ((i-1) * ydim + j) + 3];
601                     value[2] = zspacing->GetValue(k);
602                     AddPoint(points, value, proj);
603
604                     continue;
605                  }
606                  else if ((j == ydim) && (i == xdim))
607                  {
608                     value[0] = (*xvalue)[4 * ((i-1) * ydim + j - 1) + 2];
609                     value[1] = (*yvalue)[4 * ((i-1) * ydim + j - 1) + 2];
610                     value[2] = zspacing->GetValue(k);
611                     AddPoint(points, value, proj);
612
613                     continue;
614                  }
615               }
616      }
617
618      //----------------------------------------------------------------
619
620      void vtkLSCEReader::GetCellsAndPoints(const vtkStdString & xcoordinate,
621                                            const vtkStdString & ycoordinate,
622                                            const vtkStdString & zcoordinate,
623                                            bool bounds, bool proj, bool clean, std::size_t nbvertex,
624                                            vtkCellArray * cells, vtkPoints * points,
625                                            vtkIntArray * dimensions)
626      {
627         float  value[3];
628         double h = 1.0;
629         double cartesianCoord[3], vbounds[6], range[2];
630         vtkIdType pts[nbvertex] ;
631         bool is2D = true;
632         ARRAY_CREATE(xvalue, float, 1, [1]);
633         ARRAY_CREATE(yvalue, float, 1, [1]);
634         vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New();
635         vtkSmartPointer<vtkFloatArray> zspacing  = vtkSmartPointer<vtkFloatArray>::New();
636         
637         if (zcoordinate.size() > 0)
638         {
639            GetSpacings(zcoordinate, bounds, zspacing);
640            is2D = false;
641         }
642         else
643            zspacing->InsertNextValue(0);
644         
645         zspacing->GetRange (range);
646
647         std::size_t size = (*Reader->getDimensions(&xcoordinate).begin()).second;
648         if (Reader->hasBounds(xcoordinate))
649         {
650            vtkStdString boundid = Reader->getBoundsId(xcoordinate);
651            Reader->readData(*xvalue, boundid);
652            if (!proj)
653            {
654               vbounds[0] = *std::min_element(xvalue->begin(), xvalue->end());
655               vbounds[1] = *std::max_element(xvalue->begin(), xvalue->end());
656            }
657            else
658            {
659               vbounds[0] = -h;
660               vbounds[1] = h;
661            }
662         }
663         if (Reader->hasBounds(ycoordinate))
664         {
665            vtkStdString boundid = Reader->getBoundsId(ycoordinate);
666            Reader->readData(*yvalue, boundid);
667            if (!proj)
668            {
669               vbounds[2] = *std::min_element(yvalue->begin(), yvalue->end());
670               vbounds[3] = *std::max_element(yvalue->begin(), yvalue->end());
671            }
672            else
673            {
674               vbounds[2] = -h;
675               vbounds[3] = h;
676            }
677         }
678         if (!proj)
679         {
680            vbounds[4] = range[0];
681            vbounds[5] = range[1];
682         }
683         else
684         {
685            vbounds[4] = -h;
686            vbounds[5] = h;
687         }
688         locator->InitPointInsertion(points, vbounds, size * zspacing->GetNumberOfTuples() * nbvertex); 
689         
690         for (std::size_t u = 0; u < size; u++)
691         {
692            for (int v = 0; v < zspacing->GetNumberOfTuples(); v++)
693            {
694               for (std::size_t w = 0; w < nbvertex; w++)
695               {
696                  value[0] = (*xvalue)[nbvertex * u + w];
697                  value[1] = (*yvalue)[nbvertex * u + w];
698                  value[2] = zspacing->GetValue(v);
699                  pts[w] = nbvertex * u + w;
700
701                  if (proj)
702                  {
703                     value[0] = 0.017453292519943295 * value[0];
704                     value[1] = 0.017453292519943295 * value[1];
705
706                     cartesianCoord[0] = h * cos(value[0]) * cos(value[1]);
707                     cartesianCoord[1] = h * sin(value[0]) * cos(value[1]);
708                     cartesianCoord[2] = h * sin(value[1]);
709                     if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
710                     else points->InsertPoint (pts[w], cartesianCoord);
711                  }
712                  else
713                  {
714                     cartesianCoord[0] = value[0];
715                     cartesianCoord[1] = value[1];
716                     cartesianCoord[2] = value[2];
717                     if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]);
718                     else points->InsertPoint (pts[w], cartesianCoord);
719                  }
720               }
721               cells->InsertNextCell(nbvertex,  pts);
722            }
723         }
724      }
725
726      //----------------------------------------------------------------
727
728      void vtkLSCEReader::AddScalarData(vtkDataSet * output,
729                                        const vtkStdString & varname,
730                                        std::size_t record, bool bounds)
731      {
732         vtkSmartPointer<vtkDoubleArray> data = vtkDoubleArray::New();
733         data->SetName(varname.c_str());
734         ARRAY_CREATE(data_arr, double, 1, [1]);
735         Reader->readData(*data_arr, varname, record);
736         Reader->replaceMissingValue(varname, *data_arr, vtkMath::Nan());
737         
738         
739         data->SetNumberOfValues (data_arr->size());
740         for (std::size_t i = 0; i < data_arr->size(); i++)
741            data->SetValue (i, (*data_arr)[i]);
742
743         if (bounds)
744         {
745            output->GetCellData()->AddArray(data);
746            output->GetCellData()->SetActiveScalars (varname.c_str());
747         }
748         else
749         {
750            output->GetPointData()->AddArray(data);
751            output->GetPointData()->SetActiveScalars (varname.c_str());
752         }
753      }
754
755      //----------------------------------------------------------------
756     
757      int vtkLSCEReader::RequestData
758                           (vtkInformation       *  vtkNotUsed(request),
759                            vtkInformationVector ** vtkNotUsed(inputVector),
760                            vtkInformationVector *  outputVector)
761      {
762         vtkInformation *  outInfo = outputVector->GetInformationObject(0);
763         vtkDataSet     *  output  = vtkDataSet::GetData(outInfo);
764         
765         if (!output)
766         {
767            switch (this->CurGridType)
768            {
769               case(RECTILINEAR)  :
770                  output = vtkRectilinearGrid::New();
771                  break;
772               case(CURVILINEAR)  :
773                  output = vtkStructuredGrid::New();
774                  break;
775               default :
776                  output = vtkUnstructuredGrid::New();
777                  break;
778            }
779            output->SetPipelineInformation(outInfo);
780            output->Delete();
781         }
782
783         vtkUnstructuredGrid * uoutput  = vtkUnstructuredGrid::SafeDownCast(output);
784         vtkRectilinearGrid  * routput  = vtkRectilinearGrid::SafeDownCast(output);
785         vtkStructuredGrid   * soutput  = vtkStructuredGrid::SafeDownCast(output);
786         
787        double time = 0.0;
788        if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
789           time = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0];
790       
791         if (routput && !this->VarNames.empty())
792         { // Grille rectiliniéaire
793            vtkSmartPointer<vtkFloatArray>
794                           xspacing = vtkSmartPointer<vtkFloatArray>::New(),
795                           yspacing = vtkSmartPointer<vtkFloatArray>::New(),
796                           zspacing = vtkSmartPointer<vtkFloatArray>::New();
797            vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
798
799            //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, xspacing, yspacing, zspacing, dims);
800            GetSpacings(Reader->getLonCoordName(*VarNames.begin()), true, xspacing);
801            GetSpacings(Reader->getLatCoordName(*VarNames.begin()), true, yspacing);
802            GetSpacings(Reader->getVertCoordName(*VarNames.begin()), true, zspacing);
803            /*std::cout << zspacing->GetNumberOfTuples() << std::endl;*/
804
805            dims->InsertNextValue(xspacing->GetNumberOfTuples());
806            dims->InsertNextValue(yspacing->GetNumberOfTuples());
807            dims->InsertNextValue(zspacing->GetNumberOfTuples());
808           
809            CreateRectilinearGrid(routput, outInfo, xspacing, yspacing, zspacing, dims);
810            AddScalarData(routput, *VarNames.begin(), 0, true);
811
812            this->AGridDef = true;
813         }
814
815         if (soutput && !this->VarNames.empty())
816         { // Grille structurée
817             
818            vtkSmartPointer<vtkPoints>   pts  = vtkSmartPointer<vtkPoints>::New();
819            vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
820           
821            //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, pts, dims);
822            GetPoints(Reader->getLonCoordName(*VarNames.begin()),
823                      Reader->getLatCoordName(*VarNames.begin()),
824                      Reader->getVertCoordName(*VarNames.begin()), true, true, pts, dims);
825
826            CreateStructuredGrid(soutput, outInfo, pts, dims);
827            AddScalarData(soutput, *VarNames.begin(), 0, true);
828
829            this->AGridDef = true;
830         }
831
832         if (uoutput  && !this->VarNames.empty())
833         { // Grille non-structurée
834            vtkSmartPointer<vtkPoints>    pts   = vtkSmartPointer<vtkPoints>::New();
835            vtkSmartPointer<vtkIntArray>  dims  = vtkSmartPointer<vtkIntArray>::New();
836            vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
837            // VTK_QUAD, VTK_HEXAHEDRON, VTK_POLYGON
838            //CreateSimpleGrid (0, 12, 0, 12, 0, 12, pts, cells, dims);
839            //std::cout << Reader->getVertCoordName(*VarNames.begin()) << std::endl;
840
841            GetCellsAndPoints(Reader->getLonCoordName(*VarNames.begin()),
842                              Reader->getLatCoordName(*VarNames.begin()),
843                              Reader->getVertCoordName(*VarNames.begin()), true, true, true,
844                              Reader->getNbVertex(*VarNames.begin()), cells, pts, dims);
845
846            CreateUnstructuredGrid(uoutput, outInfo, pts, cells, VTK_POLYGON);
847         
848            AddScalarData(uoutput, *VarNames.begin(), this->GetRecord(*VarNames.begin(), time), true);
849           
850            this->AGridDef = true;
851         }
852         return (1);
853      }
854
855      //----------------------------------------------------------------
856     
857      void vtkLSCEReader::GetTimeInformation
858               (const vtkStdString & variable, vtkDoubleArray * values, double * timeRange)
859      {
860         std::string timename = this->Reader->getTimeCoordName(variable);
861         if (this->Reader->hasVariable(timename))
862         {
863            ARRAY_CREATE(data_arr, double, 1, [1]);
864            Reader->readData(*data_arr, timename);
865           
866            values->SetNumberOfValues (data_arr->size());
867            for (std::size_t i = 0; i < data_arr->size(); i++)
868               values->SetValue (i, (*data_arr)[i]);
869         }
870         else
871         {
872            std::size_t nbtimestep = this->Reader->getNbOfTimestep();
873            values->SetNumberOfValues (nbtimestep);
874            for (std::size_t i = 0; i < nbtimestep; i++)
875               values->SetValue (i, i);
876         }
877         timeRange[0] = values->GetValue(0);
878         timeRange[1] = values->GetValue(values->GetNumberOfTuples()-1);
879      }
880     
881      std::size_t vtkLSCEReader::GetRecord(const vtkStdString & variable, const double & value)
882      {
883         std::string timename = this->Reader->getTimeCoordName(variable);
884         if (this->Reader->hasVariable(timename))
885         {
886            ARRAY_CREATE(data_arr, double, 1, [1]);
887            Reader->readData(*data_arr, timename);
888            return (std::find(data_arr->begin(), data_arr->end(), value) - data_arr->begin());
889         }
890         else
891         {
892            return (value);
893         }
894      }
895
896      //----------------------------------------------------------------
897     
898      int vtkLSCEReader::RequestInformation
899                           (vtkInformation       *  vtkNotUsed(request),
900                            vtkInformationVector ** vtkNotUsed(inputVector),
901                            vtkInformationVector *  outputVector)
902      {
903         vtkInformation *  outInfo = outputVector->GetInformationObject(0);
904         vtkDataSet     *  output  = vtkDataSet::GetData(outInfo);
905
906         this->AGridDef = false;
907         
908         if (!this->VarNames.empty())
909         {
910            switch (this->CurGridType)
911            {
912               case(RECTILINEAR)  :
913                  output = vtkRectilinearGrid::New();
914                  break;
915               case(CURVILINEAR)  :
916                  output = vtkStructuredGrid::New();
917                  break;
918               default :
919                  output = vtkUnstructuredGrid::New();
920                  break;
921            }
922            output->SetPipelineInformation(outInfo);
923            output->Delete();
924         }
925
926         vtkRectilinearGrid  * routput  = vtkRectilinearGrid::SafeDownCast(output);
927         vtkStructuredGrid   * soutput  = vtkStructuredGrid::SafeDownCast(output);
928         
929         if (!this->VarNames.empty())
930         { // Informations temporelles
931            if (this->Reader->hasTemporalDim() && this->Reader->isTemporal(*VarNames.begin()))
932            {
933               vtkSmartPointer<vtkDoubleArray> timeStepValues = vtkSmartPointer<vtkDoubleArray>::New();
934               double timeRange[2];
935               this->GetTimeInformation(*VarNames.begin(), timeStepValues, timeRange);
936
937               outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
938                            timeStepValues->GetPointer(0),
939                            timeStepValues->GetNumberOfTuples());
940               outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
941            }
942            else
943            {
944               outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
945               outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
946            }
947           
948         }
949         
950         if (routput && !this->VarNames.empty())
951         { // Grille rectiliniéaire
952            vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
953            this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
954                                Reader->getLatCoordName(*VarNames.begin()),
955                                Reader->getVertCoordName(*VarNames.begin()),
956                                dims, true);
957           
958            int extent[6] = { 0, dims->GetValue(0) - 1
959                            , 0, dims->GetValue(1) - 1
960                            , 0, dims->GetValue(2) - 1};
961            outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
962         }
963         
964         if (soutput && !this->VarNames.empty())
965         { // Grille structurée
966           
967            vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New();
968            this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()),
969                                Reader->getVertCoordName(*VarNames.begin()),
970                                dims, true);
971           
972            int extent[6] = { 0, dims->GetValue(0) - 1
973                            , 0, dims->GetValue(1) - 1
974                            , 0, dims->GetValue(2) - 1};
975            outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6);
976         }
977           
978         return (1);
979      }
980
981      ///---------------------------------------------------------------
982
983#ifndef LSCE_EXPORTS
984   } // namespace vtk
985} // namespace xmlioserver
986#endif //LSCE_EXPORTS
Note: See TracBrowser for help on using the repository browser.