XMLIOSERVER 0.4
Serveur d'Entrées/Sorties parallèles
|
00001 00002 #ifndef LSCE_EXPORTS 00003 # include "lscereader.hpp" 00004 #else 00005 # include "vtkLSCEReader.h" 00006 #endif //LSCE_EXPORTS 00007 00008 00009 // Bibliothèque VTK 00010 #include "vtkObjectFactory.h" 00011 #include "vtkInformation.h" 00012 #include "vtkInformationVector.h" 00013 #include "vtkStreamingDemandDrivenPipeline.h" 00014 00015 #include "vtkPointLocator.h" 00016 00017 #include <vtkProperty.h> 00018 #include <vtkDataSetMapper.h> 00019 #include <vtkActor.h> 00020 #include <vtkRenderWindow.h> 00021 #include <vtkRenderer.h> 00022 #include <vtkRenderWindowInteractor.h> 00023 #include <vtkLegendScaleActor.h> 00024 00025 #include "vtkMath.h" 00026 00027 #include "inetcdf4_impl.hpp" 00028 #include "inetcdf4_adv_impl.hpp" 00029 00030 // Bibliothèque C++ 00031 #include <algorithm> 00032 00033 // Type Fortran 00034 #define ARRAY(valuetype, numdims) boost::shared_ptr<boost::multi_array<valuetype, numdims> > 00035 00036 #define ARRAY_ASSIGN(value, valuetype, numdims, extent)\ 00037 value.reset(new boost::multi_array<valuetype, numdims>(boost::extents extent)) 00038 00039 #define ARRAY_CREATE(value, valuetype, numdims, extent)\ 00040 ARRAY(valuetype, numdims) value = \ 00041 ARRAY(valuetype, numdims)(new boost::multi_array<valuetype, numdims>(boost::extents extent)) 00042 00043 // Type C 00044 #define ARRAY_C_ASSIGN(value, valuetype, numdims, extent)\ 00045 value = ARRAY(valuetype, numdims) \ 00046 (new boost::multi_array<valuetype, numdims>(boost::extents extent, c_storage_order())) 00047 00048 #define ARRAY_C_CREATE(value, valuetype, numdims, extent)\ 00049 ARRAY_C_ASSIGN(ARRAY(valuetype, numdims) value, valuetype, numdims, extent) 00050 00051 00052 #ifndef LSCE_EXPORTS 00053 namespace xmlioserver 00054 { 00055 namespace vtk 00056 { 00057 #endif //LSCE_EXPORTS 00058 00060 vtkStandardNewMacro(vtkLSCEReader); 00061 00062 vtkLSCEReader::vtkLSCEReader(void) 00063 : FileName(), CurGridType(RECTILINEAR), VarNames() 00064 , A3D(true), ATemporal(true), ACell(true), AGridDef(false) 00065 , Reader() 00066 { 00067 this->SetNumberOfInputPorts (0); 00068 this->SetNumberOfOutputPorts(1); 00069 } 00070 00071 vtkLSCEReader::~vtkLSCEReader(void) 00072 { 00073 00074 } 00075 00077 00078 void vtkLSCEReader::SetFileName(const vtkStdString & fileName) 00079 { 00080 this->FileName = fileName; 00081 this->Reader = boost::shared_ptr<xmlioserver::io::CINetCDF4Adv> 00082 (new xmlioserver::io::CINetCDF4Adv(fileName)); 00083 this->Modified(); 00084 } 00085 00086 void vtkLSCEReader::SetVariable(const vtkStdString & variable) 00087 { 00088 this->RemoveAllSelectedVariables(); 00089 if (Reader->hasVariable(variable)) 00090 this->AddVariableToSelection(variable); 00091 this->Modified(); 00092 } 00093 00094 const char * vtkLSCEReader::GetFileName(void) const 00095 { 00096 return (this->FileName.c_str()); 00097 } 00098 00099 void vtkLSCEReader::PrintSelf(ostream& os, vtkIndent indent) 00100 { 00101 this->Superclass::PrintSelf(os,indent); 00102 } 00103 00104 //---------------------------------------------------------------- 00105 void vtkLSCEReader::AddVariableToSelection(const vtkStdString & varName) 00106 { 00107 if (this->Reader->isRectilinear(varName)) 00108 this->SetGridType (RECTILINEAR); 00109 else if (this->Reader->isCurvilinear(varName)) 00110 this->SetGridType (CURVILINEAR); 00111 else 00112 this->SetGridType (UNSTRUCTURED); 00113 00114 this->VarNames.insert(varName); 00115 } 00116 00117 void vtkLSCEReader::RemoveSelectedVariable(const vtkStdString & varName) 00118 { 00119 std::set<vtkStdString>::iterator it = this->VarNames.find(varName); 00120 if (it != this->VarNames.end()) this->VarNames.erase(it); 00121 } 00122 00123 void vtkLSCEReader::RemoveAllSelectedVariables(void) 00124 { 00125 this->VarNames.clear(); 00126 } 00127 00128 const std::set<vtkStdString> & 00129 vtkLSCEReader::GetSelectedVariables(void) const 00130 { 00131 return (this->VarNames); 00132 } 00133 00134 bool vtkLSCEReader::HasSelectedVariable(void) const 00135 { 00136 return (this->VarNames.size() != 0); 00137 } 00138 00139 //---------------------------------------------------------------- 00140 00141 void vtkLSCEReader::SetGridType (GridType type) 00142 { 00143 if ((this->IsUnstructured() && (type == CURVILINEAR)) || 00144 (this->IsUnstructured() && (type == RECTILINEAR)) || 00145 (this->IsCurvilinear() && (type == RECTILINEAR))) 00146 this->RemoveAllSelectedVariables(); 00147 this->CurGridType = type; 00148 } 00149 00150 bool vtkLSCEReader::IsUnstructured(void) const 00151 { 00152 return (this->CurGridType == UNSTRUCTURED); 00153 } 00154 00155 bool vtkLSCEReader::IsCurvilinear(void) const 00156 { 00157 return (this->CurGridType == CURVILINEAR); 00158 } 00159 00160 bool vtkLSCEReader::IsRectilinear(void) const 00161 { 00162 return (this->CurGridType == RECTILINEAR); 00163 } 00164 00165 //---------------------------------------------------------------- 00166 00167 void vtkLSCEReader::AcceptTemporalOnly(bool value) 00168 { 00169 this->ATemporal = value; 00170 } 00171 00172 void vtkLSCEReader::Accept3DOnly(bool value) 00173 { 00174 this->A3D = value; 00175 } 00176 00177 void vtkLSCEReader::AcceptCellOnly(bool value) 00178 { 00179 this->ACell = value; 00180 } 00181 00182 //---------------------------------------------------------------- 00183 00184 int vtkLSCEReader::ProcessRequest 00185 (vtkInformation * request, 00186 vtkInformationVector ** vtkNotUsed(inputVector), 00187 vtkInformationVector * outputVector) 00188 { 00189 try 00190 { 00191 00192 if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) 00193 { 00194 return (this->RequestData(NULL, NULL, outputVector)); 00195 } 00196 00197 if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA_OBJECT())) 00198 { 00199 // return (this->RequestDataObject(NULL, NULL, outputVector)); 00200 } 00201 00202 if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION())) 00203 { 00204 return (this->RequestInformation(NULL, NULL, outputVector)); 00205 } 00206 00207 if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT())) 00208 { 00209 return (this->RequestUpdateExtent(NULL, NULL, outputVector)); 00210 } 00211 } 00212 catch (xmlioserver::CException & _exception) 00213 { 00214 std::cerr << _exception.getMessage() << std::endl; 00215 return (-1); 00216 } 00217 return (1); 00218 } 00219 00220 //---------------------------------------------------------------- 00221 00222 void vtkLSCEReader::ShowVariable(const vtkStdString & filename, 00223 const vtkStdString & varname) 00224 { 00225 vtkSmartPointer<vtkLSCEReader> lscereader = 00226 vtkSmartPointer<vtkLSCEReader>::New(); 00227 00228 lscereader->SetFileName(filename); 00229 lscereader->AddVariableToSelection(varname); 00230 lscereader->Update(); 00231 00232 // Création d'un mappeur et d'un acteur. 00233 vtkSmartPointer<vtkDataSetMapper> mapper = vtkSmartPointer<vtkDataSetMapper>::New(); 00234 mapper->SetInputConnection(lscereader->GetOutput()->GetProducerPort()); 00235 mapper->SetScalarRange (1000, 10000); 00236 00237 vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New(); 00238 actor->SetMapper(mapper); 00239 actor->GetProperty()->SetRepresentationToWireframe(); 00240 00241 // Création d'une fenêtre de rendu. 00242 vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New(); 00243 vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New(); 00244 renderWindow->AddRenderer(renderer); 00245 vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = 00246 vtkSmartPointer<vtkRenderWindowInteractor>::New(); 00247 renderWindowInteractor->SetRenderWindow(renderWindow); 00248 00249 00250 vtkSmartPointer<vtkLegendScaleActor> legendScaleActor = 00251 vtkSmartPointer<vtkLegendScaleActor>::New(); 00252 00253 // Ajout de l'acteur à la scène. 00254 renderer->AddActor(actor); 00255 renderer->AddActor(legendScaleActor); 00256 00257 renderer->GradientBackgroundOn(); 00258 renderer->SetBackground(1,1,1); 00259 renderer->SetBackground2(0,0,0); 00260 00261 // Rendu et intéraction 00262 renderWindow->Render(); 00263 renderWindowInteractor->Start(); 00264 } 00265 00266 //---------------------------------------------------------------- 00267 00268 void vtkLSCEReader::CreateRectilinearGrid(vtkRectilinearGrid * grid, 00269 vtkInformation * outInfo, 00270 vtkFloatArray * xspacing, 00271 vtkFloatArray * yspacing, 00272 vtkFloatArray * zspacing, 00273 vtkIntArray * dimensions) 00274 { 00275 int extent[6] = { 0, dimensions->GetValue(0) - 1 00276 , 0, dimensions->GetValue(1) - 1 00277 , 0, dimensions->GetValue(2) - 1}; 00278 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6); 00279 grid->SetExtent (extent); 00280 grid->SetDimensions (dimensions->GetValue(0), 00281 dimensions->GetValue(1), 00282 dimensions->GetValue(2)); 00283 grid->SetXCoordinates (xspacing); 00284 grid->SetYCoordinates (yspacing); 00285 grid->SetZCoordinates (zspacing); 00286 } 00287 00288 void vtkLSCEReader::CreateStructuredGrid(vtkStructuredGrid * grid, 00289 vtkInformation * outInfo, 00290 vtkPoints * points, 00291 vtkIntArray * dimensions) 00292 { 00293 int extent[6] = { 0, dimensions->GetValue(0) - 1 00294 , 0, dimensions->GetValue(1) - 1 00295 , 0, dimensions->GetValue(2) - 1}; 00296 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6); 00297 grid->SetExtent (extent); 00298 grid->SetDimensions (dimensions->GetValue(0), 00299 dimensions->GetValue(1), 00300 dimensions->GetValue(2)); 00301 grid->SetPoints(points); 00302 } 00303 00304 void vtkLSCEReader::CreateUnstructuredGrid(vtkUnstructuredGrid * grid, 00305 vtkInformation * outInfo, 00306 vtkPoints * points, 00307 vtkCellArray * cells, 00308 int celltype) 00309 { 00310 grid->SetPoints(points); 00311 grid->SetCells(celltype, cells); 00312 } 00313 00314 //---------------------------------------------------------------- 00315 00316 void vtkLSCEReader::GetRectilinearConnectivity(int dimZ, int dimY, int dimX, vtkCellArray * cells) 00317 { 00318 if (dimZ == 1 && dimY > 1 && dimX > 1) // 2D 00319 { 00320 vtkIdType pts[4] ; 00321 00322 for (int i = 0, k = 0; i < (dimX - 1); i++) 00323 for (int j = 0; j < (dimY - 1); j++, k++) 00324 { 00325 pts[0] = i * dimY + j; 00326 pts[1] = i * dimY + j + 1; 00327 pts[2] = (i+1) * dimY + j + 1; 00328 pts[3] = (i+1) * dimY + j; 00329 cells->InsertNextCell(4, pts); 00330 } 00331 } 00332 else if (dimZ > 1 && dimY > 1 && dimX > 1) // 3D 00333 { 00334 vtkIdType pts[8]; 00335 for (int i = 0, l = 0; i < (dimX - 1); i++) 00336 for (int j = 0; j < (dimY - 1); j++) 00337 for (int k = 0; k < (dimZ - 1); k++, l++) 00338 { 00339 pts[0] = i * dimY * dimZ + j * dimZ + k; 00340 pts[1] = i * dimY * dimZ + j * dimZ + k + 1; 00341 pts[2] = (i+1) * dimY * dimZ + j * dimZ + k + 1; 00342 pts[3] = (i+1) * dimY * dimZ + j * dimZ + k; 00343 00344 pts[4] = i * dimY * dimZ + (j + 1) * dimZ + k ; 00345 pts[5] = i * dimY * dimZ + (j + 1) * dimZ + k +1; 00346 pts[6] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k + 1; 00347 pts[7] = (i+1) * dimY * dimZ + (j + 1) * dimZ + k ; 00348 cells->InsertNextCell(8, pts); 00349 } 00350 } 00351 } 00352 00353 //---------------------------------------------------------------- 00354 00355 void vtkLSCEReader::CreateSimpleGrid (int xi, int xf, 00356 int yi, int yf, 00357 int zi, int zf, 00358 vtkFloatArray * xspacing, 00359 vtkFloatArray * yspacing, 00360 vtkFloatArray * zspacing, 00361 vtkIntArray * dimensions) 00362 { 00363 dimensions->InsertNextValue (xf - xi + 1); // X 00364 dimensions->InsertNextValue (yf - yi + 1); // Y 00365 dimensions->InsertNextValue (zf - zi + 1); // Z 00366 for (int i = xi; i <= xf; i++) 00367 xspacing->InsertNextValue (i); 00368 for (int j = yi; j <= yf; j++) 00369 yspacing->InsertNextValue (j); 00370 for (int k = zi; k <= zf; k++) 00371 zspacing->InsertNextValue (k); 00372 } 00373 00374 void vtkLSCEReader::CreateSimpleGrid (int xi, int xf, 00375 int yi, int yf, 00376 int zi, int zf, 00377 vtkPoints * points, 00378 vtkIntArray * dimensions) 00379 { 00380 dimensions->InsertNextValue (xf - xi + 1); // X 00381 dimensions->InsertNextValue (yf - yi + 1); // Y 00382 dimensions->InsertNextValue (zf - zi + 1); // Z 00383 for (int i = xi; i <= xf; i++) 00384 for (int j = yi; j <= yf; j++) 00385 for (int k = zi; k <= zf; k++) 00386 points->InsertNextPoint(i, j, k); 00387 } 00388 00389 void vtkLSCEReader::CreateSimpleGrid (int xi, int xf, 00390 int yi, int yf, 00391 int zi, int zf, 00392 vtkPoints * points, 00393 vtkCellArray * cells, 00394 vtkIntArray * dimensions) 00395 { 00396 this->CreateSimpleGrid (xi, xf, yi, yf, zi, zf, points, dimensions); 00397 int dimX = dimensions->GetValue(0), 00398 dimY = dimensions->GetValue(1), 00399 dimZ = dimensions->GetValue(2); 00400 this->GetRectilinearConnectivity(dimZ, dimY, dimX, cells); 00401 } 00402 00403 //---------------------------------------------------------------- 00404 void vtkLSCEReader::GetSpacings(const vtkStdString & coordinate, 00405 bool bounds, vtkFloatArray * spacing) 00406 { 00407 std::size_t size = (Reader->getDimensions(&coordinate).begin())->second; 00408 ARRAY_CREATE(data_arr, float, 1, [1]); 00409 if (Reader->hasVariable(coordinate)) 00410 { 00411 if (bounds) 00412 { 00413 if (Reader->hasBounds(coordinate)) 00414 { // CellData, variable de coordonnée avec attributes bounds. 00415 vtkStdString boundid = Reader->getBoundsId(coordinate); 00416 float f; 00417 Reader->readData(*data_arr, boundid); 00418 spacing->InsertNextValue (f = (*data_arr)[0]); 00419 for (std::size_t j = 1; j < (2 * size); j++) // << Super laid 00420 if (f != (*data_arr)[j]) 00421 spacing->InsertNextValue (f = (*data_arr)[j]); 00422 } 00423 else 00424 { // CellData, variable de coordonnée sans attributes bounds. 00425 Reader->readData(*data_arr, coordinate); 00426 spacing->InsertNextValue ((*data_arr)[0] - ((*data_arr)[1] - (*data_arr)[0])/2); 00427 for (std::size_t j = 0; j < (size-1); j++) 00428 spacing->InsertNextValue ((*data_arr)[j] + ((*data_arr)[j+1]-(*data_arr)[j])/2); 00429 spacing->InsertNextValue ((*data_arr)[size-2] + ((*data_arr)[size-1] - (*data_arr)[size-2])/2); 00430 } 00431 } 00432 else 00433 { // PointData, variable de coordonnée. 00434 Reader->readData(*data_arr, coordinate); 00435 for (std::size_t j = 0; j < (size-1); j++) 00436 spacing->InsertNextValue ((*data_arr)[j]); 00437 } 00438 } 00439 else 00440 { 00441 if (bounds) // CellData, dimension de coordonnée. 00442 for (float i = -0.5; i < ((float)size + 1.); i = i + 1.) 00443 spacing->InsertNextValue (i); 00444 else // PointData, dimension de coordonnée. 00445 for (float j = 0.; j < (float)size; j = j + 1.) 00446 spacing->InsertNextValue (j); 00447 } 00448 } 00449 //---------------------------------------------------------------- 00450 00451 void vtkLSCEReader::AddPoint(vtkPoints * points, float * value, bool proj) 00452 { 00453 if (proj) 00454 { 00455 double h = 1.0; 00456 double cartesianCoord[3]; 00457 value[0] = 0.017453292519943295 * value[0]; 00458 value[1] = 0.017453292519943295 * value[1]; 00459 00460 cartesianCoord[0] = h * cos(value[0]) * cos(value[1]); 00461 cartesianCoord[1] = h * sin(value[0]) * cos(value[1]); 00462 cartesianCoord[2] = h * sin(value[1]); 00463 points->InsertNextPoint (cartesianCoord); 00464 } 00465 else 00466 { 00467 points->InsertNextPoint (value); 00468 } 00469 } 00470 00471 //---------------------------------------------------------------- 00472 00473 void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate, 00474 const vtkStdString & ycoordinate, 00475 const vtkStdString & zcoordinate, 00476 vtkIntArray * dimensions, bool bounds) 00477 { 00478 int a = (bounds) ? 1 : 0; 00479 dimensions->InsertNextValue 00480 ((Reader->getDimensions(&xcoordinate).begin())->second + a); 00481 dimensions->InsertNextValue 00482 ((Reader->getDimensions(&ycoordinate).begin())->second + a); 00483 if (zcoordinate.size() > 0) 00484 { 00485 dimensions->InsertNextValue 00486 ((Reader->getDimensions(&zcoordinate).begin())->second + a); 00487 } 00488 else 00489 { 00490 dimensions->InsertNextValue(1); 00491 } 00492 } 00493 00494 //---------------------------------------------------------------- 00495 00496 void vtkLSCEReader::GetDimensions(const vtkStdString & xcoordinate, 00497 const vtkStdString & zcoordinate, 00498 vtkIntArray * dimensions, bool bounds) 00499 { 00500 int a = (bounds) ? 1 : 0; 00501 dimensions->InsertNextValue 00502 (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)] + a); 00503 dimensions->InsertNextValue 00504 (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)] + a); 00505 if (zcoordinate.size() > 0) 00506 { 00507 dimensions->InsertNextValue 00508 ((Reader->getDimensions(&zcoordinate).begin())->second + a); 00509 } 00510 else 00511 { 00512 dimensions->InsertNextValue(1); 00513 } 00514 } 00515 00516 //---------------------------------------------------------------- 00517 void vtkLSCEReader::GetPoints(const vtkStdString & xcoordinate, 00518 const vtkStdString & ycoordinate, 00519 const vtkStdString & zcoordinate, 00520 bool bounds, bool proj, 00521 vtkPoints * points, vtkIntArray * dimensions) 00522 { 00523 float value[3]; 00524 bool is2D = true; 00525 ARRAY_CREATE(xvalue, float, 1, [1]); 00526 ARRAY_CREATE(yvalue, float, 1, [1]); 00527 vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New(); 00528 if (zcoordinate.size() > 0) 00529 { 00530 GetSpacings(zcoordinate, bounds, zspacing); 00531 is2D = false; 00532 } 00533 else 00534 zspacing->InsertNextValue(0); 00535 00536 dimensions->InsertNextValue 00537 (Reader->getDimensions(&xcoordinate)[Reader->getLonCoordName(xcoordinate)]); 00538 dimensions->InsertNextValue 00539 (Reader->getDimensions(&xcoordinate)[Reader->getLatCoordName(xcoordinate)]); 00540 dimensions->InsertNextValue(zspacing->GetNumberOfTuples()); 00541 00542 if (!bounds) 00543 { 00544 Reader->readData(*xvalue, xcoordinate); 00545 Reader->readData(*yvalue, ycoordinate); 00546 00547 for (int i = 0; i < dimensions->GetValue(0); i++) 00548 for (int j = 0; j < dimensions->GetValue(1); j++) 00549 for (int k = 0; k < dimensions->GetValue(2); k++) 00550 { 00551 value[0] = (*xvalue)[i*dimensions->GetValue(1) + j]; 00552 value[1] = (*yvalue)[i*dimensions->GetValue(1) + j]; 00553 value[2] = zspacing->GetValue(k); 00554 AddPoint(points, value, proj); 00555 } 00556 return; 00557 } 00558 00559 dimensions->SetValue(0, dimensions->GetValue(0) + 1); 00560 dimensions->SetValue(1, dimensions->GetValue(1) + 1); 00561 00562 if (Reader->hasBounds(xcoordinate)) 00563 { 00564 vtkStdString boundid = Reader->getBoundsId(xcoordinate); 00565 Reader->readData(*xvalue, boundid); 00566 } 00567 if (Reader->hasBounds(ycoordinate)) 00568 { 00569 vtkStdString boundid = Reader->getBoundsId(ycoordinate); 00570 Reader->readData(*yvalue, boundid); 00571 } 00572 int ydim = dimensions->GetValue(0) - 1 , 00573 xdim = dimensions->GetValue(1) - 1; 00574 00575 for (int i = 0; i <= xdim; i++) 00576 for (int j = 0; j <= ydim; j++) 00577 for (int k = 0; k < dimensions->GetValue(2); k++) 00578 { 00579 if ((j != ydim) && (i != xdim)) 00580 { 00581 value[0] = (*xvalue)[4 * (i * ydim + j) + 0]; 00582 value[1] = (*yvalue)[4 * (i * ydim + j) + 0]; 00583 value[2] = zspacing->GetValue(k); 00584 AddPoint(points, value, proj); 00585 00586 continue; 00587 } 00588 else if ((j == ydim) && (i != xdim)) 00589 { 00590 value[0] = (*xvalue)[4 * (i * ydim + j - 1) + 1]; 00591 value[1] = (*yvalue)[4 * (i * ydim + j - 1) + 1]; 00592 value[2] = zspacing->GetValue(k); 00593 AddPoint(points, value, proj); 00594 00595 continue; 00596 } 00597 else if ((j != ydim) && (i == xdim)) 00598 { 00599 value[0] = (*xvalue)[4 * ((i-1) * ydim + j) + 3]; 00600 value[1] = (*yvalue)[4 * ((i-1) * ydim + j) + 3]; 00601 value[2] = zspacing->GetValue(k); 00602 AddPoint(points, value, proj); 00603 00604 continue; 00605 } 00606 else if ((j == ydim) && (i == xdim)) 00607 { 00608 value[0] = (*xvalue)[4 * ((i-1) * ydim + j - 1) + 2]; 00609 value[1] = (*yvalue)[4 * ((i-1) * ydim + j - 1) + 2]; 00610 value[2] = zspacing->GetValue(k); 00611 AddPoint(points, value, proj); 00612 00613 continue; 00614 } 00615 } 00616 } 00617 00618 //---------------------------------------------------------------- 00619 00620 void vtkLSCEReader::GetCellsAndPoints(const vtkStdString & xcoordinate, 00621 const vtkStdString & ycoordinate, 00622 const vtkStdString & zcoordinate, 00623 bool bounds, bool proj, bool clean, std::size_t nbvertex, 00624 vtkCellArray * cells, vtkPoints * points, 00625 vtkIntArray * dimensions) 00626 { 00627 float value[3]; 00628 double h = 1.0; 00629 double cartesianCoord[3], vbounds[6], range[2]; 00630 vtkIdType pts[nbvertex] ; 00631 bool is2D = true; 00632 ARRAY_CREATE(xvalue, float, 1, [1]); 00633 ARRAY_CREATE(yvalue, float, 1, [1]); 00634 vtkSmartPointer<vtkPointLocator> locator = vtkSmartPointer<vtkPointLocator>::New(); 00635 vtkSmartPointer<vtkFloatArray> zspacing = vtkSmartPointer<vtkFloatArray>::New(); 00636 00637 if (zcoordinate.size() > 0) 00638 { 00639 GetSpacings(zcoordinate, bounds, zspacing); 00640 is2D = false; 00641 } 00642 else 00643 zspacing->InsertNextValue(0); 00644 00645 zspacing->GetRange (range); 00646 00647 std::size_t size = (*Reader->getDimensions(&xcoordinate).begin()).second; 00648 if (Reader->hasBounds(xcoordinate)) 00649 { 00650 vtkStdString boundid = Reader->getBoundsId(xcoordinate); 00651 Reader->readData(*xvalue, boundid); 00652 if (!proj) 00653 { 00654 vbounds[0] = *std::min_element(xvalue->begin(), xvalue->end()); 00655 vbounds[1] = *std::max_element(xvalue->begin(), xvalue->end()); 00656 } 00657 else 00658 { 00659 vbounds[0] = -h; 00660 vbounds[1] = h; 00661 } 00662 } 00663 if (Reader->hasBounds(ycoordinate)) 00664 { 00665 vtkStdString boundid = Reader->getBoundsId(ycoordinate); 00666 Reader->readData(*yvalue, boundid); 00667 if (!proj) 00668 { 00669 vbounds[2] = *std::min_element(yvalue->begin(), yvalue->end()); 00670 vbounds[3] = *std::max_element(yvalue->begin(), yvalue->end()); 00671 } 00672 else 00673 { 00674 vbounds[2] = -h; 00675 vbounds[3] = h; 00676 } 00677 } 00678 if (!proj) 00679 { 00680 vbounds[4] = range[0]; 00681 vbounds[5] = range[1]; 00682 } 00683 else 00684 { 00685 vbounds[4] = -h; 00686 vbounds[5] = h; 00687 } 00688 locator->InitPointInsertion(points, vbounds, size * zspacing->GetNumberOfTuples() * nbvertex); 00689 00690 for (std::size_t u = 0; u < size; u++) 00691 { 00692 for (int v = 0; v < zspacing->GetNumberOfTuples(); v++) 00693 { 00694 for (std::size_t w = 0; w < nbvertex; w++) 00695 { 00696 value[0] = (*xvalue)[nbvertex * u + w]; 00697 value[1] = (*yvalue)[nbvertex * u + w]; 00698 value[2] = zspacing->GetValue(v); 00699 pts[w] = nbvertex * u + w; 00700 00701 if (proj) 00702 { 00703 value[0] = 0.017453292519943295 * value[0]; 00704 value[1] = 0.017453292519943295 * value[1]; 00705 00706 cartesianCoord[0] = h * cos(value[0]) * cos(value[1]); 00707 cartesianCoord[1] = h * sin(value[0]) * cos(value[1]); 00708 cartesianCoord[2] = h * sin(value[1]); 00709 if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]); 00710 else points->InsertPoint (pts[w], cartesianCoord); 00711 } 00712 else 00713 { 00714 cartesianCoord[0] = value[0]; 00715 cartesianCoord[1] = value[1]; 00716 cartesianCoord[2] = value[2]; 00717 if (clean) locator->InsertUniquePoint (cartesianCoord, pts[w]); 00718 else points->InsertPoint (pts[w], cartesianCoord); 00719 } 00720 } 00721 cells->InsertNextCell(nbvertex, pts); 00722 } 00723 } 00724 } 00725 00726 //---------------------------------------------------------------- 00727 00728 void vtkLSCEReader::AddScalarData(vtkDataSet * output, 00729 const vtkStdString & varname, 00730 std::size_t record, bool bounds) 00731 { 00732 vtkSmartPointer<vtkDoubleArray> data = vtkDoubleArray::New(); 00733 data->SetName(varname.c_str()); 00734 ARRAY_CREATE(data_arr, double, 1, [1]); 00735 Reader->readData(*data_arr, varname, record); 00736 Reader->replaceMissingValue(varname, *data_arr, vtkMath::Nan()); 00737 00738 00739 data->SetNumberOfValues (data_arr->size()); 00740 for (std::size_t i = 0; i < data_arr->size(); i++) 00741 data->SetValue (i, (*data_arr)[i]); 00742 00743 if (bounds) 00744 { 00745 output->GetCellData()->AddArray(data); 00746 output->GetCellData()->SetActiveScalars (varname.c_str()); 00747 } 00748 else 00749 { 00750 output->GetPointData()->AddArray(data); 00751 output->GetPointData()->SetActiveScalars (varname.c_str()); 00752 } 00753 } 00754 00755 //---------------------------------------------------------------- 00756 00757 int vtkLSCEReader::RequestData 00758 (vtkInformation * vtkNotUsed(request), 00759 vtkInformationVector ** vtkNotUsed(inputVector), 00760 vtkInformationVector * outputVector) 00761 { 00762 vtkInformation * outInfo = outputVector->GetInformationObject(0); 00763 vtkDataSet * output = vtkDataSet::GetData(outInfo); 00764 00765 if (!output) 00766 { 00767 switch (this->CurGridType) 00768 { 00769 case(RECTILINEAR) : 00770 output = vtkRectilinearGrid::New(); 00771 break; 00772 case(CURVILINEAR) : 00773 output = vtkStructuredGrid::New(); 00774 break; 00775 default : 00776 output = vtkUnstructuredGrid::New(); 00777 break; 00778 } 00779 output->SetPipelineInformation(outInfo); 00780 output->Delete(); 00781 } 00782 00783 vtkUnstructuredGrid * uoutput = vtkUnstructuredGrid::SafeDownCast(output); 00784 vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output); 00785 vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output); 00786 00787 double time = 0.0; 00788 if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())) 00789 time = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0]; 00790 00791 if (routput && !this->VarNames.empty()) 00792 { // Grille rectiliniéaire 00793 vtkSmartPointer<vtkFloatArray> 00794 xspacing = vtkSmartPointer<vtkFloatArray>::New(), 00795 yspacing = vtkSmartPointer<vtkFloatArray>::New(), 00796 zspacing = vtkSmartPointer<vtkFloatArray>::New(); 00797 vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New(); 00798 00799 //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, xspacing, yspacing, zspacing, dims); 00800 GetSpacings(Reader->getLonCoordName(*VarNames.begin()), true, xspacing); 00801 GetSpacings(Reader->getLatCoordName(*VarNames.begin()), true, yspacing); 00802 GetSpacings(Reader->getVertCoordName(*VarNames.begin()), true, zspacing); 00803 /*std::cout << zspacing->GetNumberOfTuples() << std::endl;*/ 00804 00805 dims->InsertNextValue(xspacing->GetNumberOfTuples()); 00806 dims->InsertNextValue(yspacing->GetNumberOfTuples()); 00807 dims->InsertNextValue(zspacing->GetNumberOfTuples()); 00808 00809 CreateRectilinearGrid(routput, outInfo, xspacing, yspacing, zspacing, dims); 00810 AddScalarData(routput, *VarNames.begin(), 0, true); 00811 00812 this->AGridDef = true; 00813 } 00814 00815 if (soutput && !this->VarNames.empty()) 00816 { // Grille structurée 00817 00818 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New(); 00819 vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New(); 00820 00821 //CreateSimpleGrid (-1, 10, -1, 10, -1, 10, pts, dims); 00822 GetPoints(Reader->getLonCoordName(*VarNames.begin()), 00823 Reader->getLatCoordName(*VarNames.begin()), 00824 Reader->getVertCoordName(*VarNames.begin()), true, true, pts, dims); 00825 00826 CreateStructuredGrid(soutput, outInfo, pts, dims); 00827 AddScalarData(soutput, *VarNames.begin(), 0, true); 00828 00829 this->AGridDef = true; 00830 } 00831 00832 if (uoutput && !this->VarNames.empty()) 00833 { // Grille non-structurée 00834 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New(); 00835 vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New(); 00836 vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New(); 00837 // VTK_QUAD, VTK_HEXAHEDRON, VTK_POLYGON 00838 //CreateSimpleGrid (0, 12, 0, 12, 0, 12, pts, cells, dims); 00839 //std::cout << Reader->getVertCoordName(*VarNames.begin()) << std::endl; 00840 00841 GetCellsAndPoints(Reader->getLonCoordName(*VarNames.begin()), 00842 Reader->getLatCoordName(*VarNames.begin()), 00843 Reader->getVertCoordName(*VarNames.begin()), true, true, true, 00844 Reader->getNbVertex(*VarNames.begin()), cells, pts, dims); 00845 00846 CreateUnstructuredGrid(uoutput, outInfo, pts, cells, VTK_POLYGON); 00847 00848 AddScalarData(uoutput, *VarNames.begin(), this->GetRecord(*VarNames.begin(), time), true); 00849 00850 this->AGridDef = true; 00851 } 00852 return (1); 00853 } 00854 00855 //---------------------------------------------------------------- 00856 00857 void vtkLSCEReader::GetTimeInformation 00858 (const vtkStdString & variable, vtkDoubleArray * values, double * timeRange) 00859 { 00860 std::string timename = this->Reader->getTimeCoordName(variable); 00861 if (this->Reader->hasVariable(timename)) 00862 { 00863 ARRAY_CREATE(data_arr, double, 1, [1]); 00864 Reader->readData(*data_arr, timename); 00865 00866 values->SetNumberOfValues (data_arr->size()); 00867 for (std::size_t i = 0; i < data_arr->size(); i++) 00868 values->SetValue (i, (*data_arr)[i]); 00869 } 00870 else 00871 { 00872 std::size_t nbtimestep = this->Reader->getNbOfTimestep(); 00873 values->SetNumberOfValues (nbtimestep); 00874 for (std::size_t i = 0; i < nbtimestep; i++) 00875 values->SetValue (i, i); 00876 } 00877 timeRange[0] = values->GetValue(0); 00878 timeRange[1] = values->GetValue(values->GetNumberOfTuples()-1); 00879 } 00880 00881 std::size_t vtkLSCEReader::GetRecord(const vtkStdString & variable, const double & value) 00882 { 00883 std::string timename = this->Reader->getTimeCoordName(variable); 00884 if (this->Reader->hasVariable(timename)) 00885 { 00886 ARRAY_CREATE(data_arr, double, 1, [1]); 00887 Reader->readData(*data_arr, timename); 00888 return (std::find(data_arr->begin(), data_arr->end(), value) - data_arr->begin()); 00889 } 00890 else 00891 { 00892 return (value); 00893 } 00894 } 00895 00896 //---------------------------------------------------------------- 00897 00898 int vtkLSCEReader::RequestInformation 00899 (vtkInformation * vtkNotUsed(request), 00900 vtkInformationVector ** vtkNotUsed(inputVector), 00901 vtkInformationVector * outputVector) 00902 { 00903 vtkInformation * outInfo = outputVector->GetInformationObject(0); 00904 vtkDataSet * output = vtkDataSet::GetData(outInfo); 00905 00906 this->AGridDef = false; 00907 00908 if (!this->VarNames.empty()) 00909 { 00910 switch (this->CurGridType) 00911 { 00912 case(RECTILINEAR) : 00913 output = vtkRectilinearGrid::New(); 00914 break; 00915 case(CURVILINEAR) : 00916 output = vtkStructuredGrid::New(); 00917 break; 00918 default : 00919 output = vtkUnstructuredGrid::New(); 00920 break; 00921 } 00922 output->SetPipelineInformation(outInfo); 00923 output->Delete(); 00924 } 00925 00926 vtkRectilinearGrid * routput = vtkRectilinearGrid::SafeDownCast(output); 00927 vtkStructuredGrid * soutput = vtkStructuredGrid::SafeDownCast(output); 00928 00929 if (!this->VarNames.empty()) 00930 { // Informations temporelles 00931 if (this->Reader->hasTemporalDim() && this->Reader->isTemporal(*VarNames.begin())) 00932 { 00933 vtkSmartPointer<vtkDoubleArray> timeStepValues = vtkSmartPointer<vtkDoubleArray>::New(); 00934 double timeRange[2]; 00935 this->GetTimeInformation(*VarNames.begin(), timeStepValues, timeRange); 00936 00937 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), 00938 timeStepValues->GetPointer(0), 00939 timeStepValues->GetNumberOfTuples()); 00940 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2); 00941 } 00942 else 00943 { 00944 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); 00945 outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE()); 00946 } 00947 00948 } 00949 00950 if (routput && !this->VarNames.empty()) 00951 { // Grille rectiliniéaire 00952 vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New(); 00953 this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()), 00954 Reader->getLatCoordName(*VarNames.begin()), 00955 Reader->getVertCoordName(*VarNames.begin()), 00956 dims, true); 00957 00958 int extent[6] = { 0, dims->GetValue(0) - 1 00959 , 0, dims->GetValue(1) - 1 00960 , 0, dims->GetValue(2) - 1}; 00961 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6); 00962 } 00963 00964 if (soutput && !this->VarNames.empty()) 00965 { // Grille structurée 00966 00967 vtkSmartPointer<vtkIntArray> dims = vtkSmartPointer<vtkIntArray>::New(); 00968 this->GetDimensions(Reader->getLonCoordName(*VarNames.begin()), 00969 Reader->getVertCoordName(*VarNames.begin()), 00970 dims, true); 00971 00972 int extent[6] = { 0, dims->GetValue(0) - 1 00973 , 0, dims->GetValue(1) - 1 00974 , 0, dims->GetValue(2) - 1}; 00975 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),extent,6); 00976 } 00977 00978 return (1); 00979 } 00980 00982 00983 #ifndef LSCE_EXPORTS 00984 } // namespace vtk 00985 } // namespace xmlioserver 00986 #endif //LSCE_EXPORTS