src/Plugins/MedReader/IO/vtkMedUnstructuredGrid.cxx

   1 // Copyright (C) 2010-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "vtkMedUnstructuredGrid.h"
  21
  22 #include "vtkObjectFactory.h"
  23 #include "vtkSmartPointer.h"
  24 #include "vtkDataArray.h"
  25 #include "vtkUnstructuredGrid.h"
  26 #include "vtkPointData.h"
  27 #include "vtkCellData.h"
  28 #include "vtkIdList.h"
  29 #include "vtkCellType.h"
  30 #include "vtkInformation.h"
  31
  32 #include "vtkMedUtilities.h"
  33 #include "vtkMedEntityArray.h"
  34 #include "vtkMedMesh.h"
  35 #include "vtkMedFile.h"
  36 #include "vtkMedDriver.h"
  37 #include "vtkMedFamilyOnEntityOnProfile.h"
  38 #include "vtkMedFamilyOnEntity.h"
  39 #include "vtkMedProfile.h"
  40 #include "vtkMedUtilities.h"
  41 #include "vtkMedStructElement.h"
  42 #include "vtkMedVariableAttribute.h"
  43
  44 #include "vtkMultiProcessController.h"
  45
  46 vtkCxxSetObjectMacro(vtkMedUnstructuredGrid,Coordinates,vtkDataArray);
  47
  48 // vtkCxxRevisionMacro(vtkMedUnstructuredGrid, "$Revision$")
  49 vtkStandardNewMacro(vtkMedUnstructuredGrid)
  50
  51 vtkMedUnstructuredGrid::vtkMedUnstructuredGrid()
  52 {
  53   this->Coordinates = NULL;
  54   this->NumberOfPoints = 0;
  55 }
  56
  57 vtkMedUnstructuredGrid::~vtkMedUnstructuredGrid()
  58 {
  59   this->SetCoordinates(NULL);
  60 }
  61
  62 int vtkMedUnstructuredGrid::IsCoordinatesLoaded()
  63 {
  64   return this->Coordinates != NULL && this->Coordinates->GetNumberOfTuples()
  65      == this->NumberOfPoints;
  66 }
  67
  68 void  vtkMedUnstructuredGrid::InitializeCellGlobalIds()
  69 {
  70   vtkIdType ncells = 0;
  71
  72   for(int id = 0; id < this->EntityArray->size(); id++)
  73     {
  74     vtkMedEntityArray* array = this->EntityArray->at(id);
  75
  76     if(array == NULL)
  77       continue;
  78
  79     if(array->GetEntity().EntityType == MED_NODE)
  80       continue;
  81
  82     array->SetInitialGlobalId(ncells + 1);
  83     ncells += array->GetNumberOfEntity();
  84     }
  85 }
  86
  87 void  vtkMedUnstructuredGrid::ClearMedSupports()
  88 {
  89   this->Superclass::ClearMedSupports();
  90   for(int id = 0; id < this->EntityArray->size(); id++)
  91     {
  92     vtkMedEntityArray* array = this->EntityArray->at(id);
  93     array->Initialize();
  94     }
  95 }
  96
  97 void  vtkMedUnstructuredGrid::LoadCoordinates()
  98 {
  99   this->GetParentMesh()->GetParentFile()->GetMedDriver()->LoadCoordinates(this);
 100 }
 101
 102 double* vtkMedUnstructuredGrid::GetCoordTuple(med_int index)
 103 {
 104   return this->Coordinates->GetTuple(index);
 105 }
 106
 107 vtkDataSet* vtkMedUnstructuredGrid::CreateVTKDataSet(
 108     vtkMedFamilyOnEntityOnProfile* foep)
 109 {
 110   vtkMultiProcessController* controller =
 111       vtkMultiProcessController::GetGlobalController();
 112
 113   vtkMedFamilyOnEntity* foe = foep->GetFamilyOnEntity();
 114   vtkMedEntityArray* array = foe->GetEntityArray();
 115   if(foep->GetValid() == 0)
 116     {
 117     return NULL;
 118     }
 119
 120   vtkUnstructuredGrid* vtkugrid = vtkUnstructuredGrid::New();
 121
 122   // now we copy all the flagged nodes in the grid, shallow copy if possible
 123   bool shallowCopyPoints=true;
 124
 125   if(this->GetParentMesh()->GetSpaceDimension()!=3)
 126     {
 127     shallowCopyPoints=false;
 128     }
 129
 130   shallowCopyPoints = shallowCopyPoints && foep->CanShallowCopyPointField(NULL);
 131
 132   foe->GetParentGrid()->LoadCoordinates();
 133
 134   vtkIdType numberOfPoints;
 135   vtkPoints* points=vtkPoints::New(this->GetCoordinates()->GetDataType());
 136   vtkugrid->SetPoints(points);
 137   points->Delete();
 138
 139   vtkIdTypeArray* pointGlobalIds=vtkIdTypeArray::New();
 140   pointGlobalIds->SetName("MED_POINT_ID");
 141   pointGlobalIds->SetNumberOfComponents(1);
 142   vtkugrid->GetPointData()->SetGlobalIds(pointGlobalIds);
 143   pointGlobalIds->Delete();
 144
 145   vtkIdTypeArray* cellGlobalIds=vtkIdTypeArray::New();
 146   cellGlobalIds->SetName("MED_CELL_ID");
 147   cellGlobalIds->SetNumberOfComponents(1);
 148   vtkugrid->GetCellData()->SetGlobalIds(cellGlobalIds);
 149   cellGlobalIds->Delete();
 150
 151   if (shallowCopyPoints)
 152     {
 153     vtkugrid->GetPoints()->SetDataType(this->GetCoordinates()->GetDataType());
 154     vtkugrid->GetPoints()->SetData(this->GetCoordinates());
 155     // add global ids
 156     numberOfPoints=this->GetNumberOfPoints();
 157     pointGlobalIds->SetNumberOfTuples(numberOfPoints);
 158     vtkIdType* ptr=pointGlobalIds->GetPointer(0);
 159     for(int pid=0; pid<numberOfPoints; pid++)
 160       ptr[pid]=pid+1;
 161     }
 162   else
 163     {
 164     vtkIdType currentIndex=0;
 165     for(vtkIdType index=0; index<this->GetNumberOfPoints(); index++)
 166       {
 167       if(!foep->KeepPoint(index))
 168         continue;
 169
 170       double coord[3]={0.0, 0.0, 0.0};
 171       double * tuple=this->GetCoordinates()->GetTuple(index);
 172       for(int dim=0; dim<this->GetParentMesh()->GetSpaceDimension()&&dim<3; dim++)
 173         {
 174         coord[dim]=tuple[dim];
 175         }
 176       vtkugrid->GetPoints()->InsertPoint(currentIndex, coord);
 177       pointGlobalIds->InsertNextValue(index+1);
 178       currentIndex++;
 179       }
 180     vtkugrid->GetPoints()->Squeeze();
 181     pointGlobalIds->Squeeze();
 182     numberOfPoints=currentIndex;
 183     }
 184
 185   vtkMedStructElement* structelem = NULL;
 186   int vtkType = VTK_EMPTY_CELL;
 187   int nsupportcell = 1;
 188   vtkSmartPointer<vtkIdTypeArray> supportIndex = vtkSmartPointer<vtkIdTypeArray>::New();
 189   supportIndex->SetName("STRUCT_ELEMENT_INDEX");
 190   supportIndex->SetNumberOfComponents(2);
 191   supportIndex->GetInformation()->Set(vtkMedUtilities::STRUCT_ELEMENT_INDEX(), 1);
 192   supportIndex->GetInformation()->Set(vtkAbstractArray::GUI_HIDE(), 1);
 193   supportIndex->SetComponentName(0, "GLOBAL_ID");
 194   supportIndex->SetComponentName(1, "LOCAL_ID");
 195
 196   if(array->GetEntity().EntityType == MED_STRUCT_ELEMENT)
 197     {
 198     structelem = array->GetStructElement();
 199     if(structelem == NULL)
 200       {
 201       foep->SetValid(0);
 202       }
 203     else
 204       {
 205       vtkType = vtkMedUtilities::GetVTKCellType(
 206           structelem->GetSupportGeometryType());
 207       nsupportcell = structelem->GetSupportNumberOfCell();
 208       if(structelem->GetSupportEntityType() != MED_CELL
 209          || strcmp(structelem->GetName(), MED_PARTICLE_NAME) == 0)
 210         {
 211         // Special case : the support connectivity is implicit
 212         // (for particles)
 213         // map this to points
 214         vtkType = VTK_VERTEX;
 215         nsupportcell = structelem->GetSupportNumberOfNode();
 216         }
 217
 218       vtkugrid->GetInformation()->Set(vtkMedUtilities::STRUCT_ELEMENT(), structelem);
 219       std::cout << "structelem->GetNumberOfVariableAttribute() = " << structelem->GetNumberOfVariableAttribute() << std::endl;
 220       for(int varattid = 0; varattid<structelem->GetNumberOfVariableAttribute(); varattid++)
 221         {
 222         vtkMedVariableAttribute* varatt = structelem->GetVariableAttribute(varattid);
 223         varatt->Load(array);
 224         vtkAbstractArray* values = array->GetVariableAttributeValue(varatt);
 225         vtkugrid->GetFieldData()->AddArray(values);
 226         }
 227       vtkugrid->GetCellData()->AddArray(supportIndex);
 228       }
 229     }
 230   else
 231     {
 232     vtkType = vtkMedUtilities::GetVTKCellType(
 233         array->GetEntity().GeometryType);
 234     }
 235
 236   vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
 237   vtkSmartPointer<vtkIdList> vtkpts = vtkSmartPointer<vtkIdList>::New();
 238   vtkIdType intialGlobalId = array->GetInitialGlobalId();
 239   vtkMedIntArray* pids = (foep->GetProfile()!=NULL?
 240                            foep->GetProfile()->GetIds() : NULL);
 241   vtkIdType maxId = (pids!=NULL?
 242                      pids->GetNumberOfTuples():
 243                      array->GetNumberOfEntity());
 244
 245   int valid = foep->GetValid();
 246   if (controller != NULL)
 247     if (controller->GetNumberOfProcesses() > 1)
 248     valid = 1;
 249
 250   for (vtkIdType pindex = 0; pindex<maxId && valid; pindex++)
 251     {
 252     vtkIdType realIndex = (pids!=NULL?
 253                            pids->GetValue(pindex)-1:
 254                            pindex);
 255
 256     if (!foep->KeepCell(realIndex))
 257       continue;
 258
 259     array->GetCellVertices(realIndex, pts);
 260
 261     for(int sid = 0; sid < nsupportcell; sid++)
 262       {
 263       cellGlobalIds->InsertNextValue(intialGlobalId+pindex);
 264       }
 265
 266     // The supportIndex array has 2 component, the first to give the index
 267     // of the med cell in the variable attributes array, the second to give
 268     // the index of the cell in the n cells composing the structural element
 269     if(array->GetEntity().EntityType == MED_STRUCT_ELEMENT)
 270       {
 271       for(int sid = 0; sid < nsupportcell; sid++)
 272         {
 273         supportIndex->InsertNextTuple2(pindex, sid);
 274         }
 275       }
 276
 277     if (array->GetEntity().GeometryType==MED_POLYHEDRON)
 278       {
 279       if(vtkMedUtilities::FormatPolyhedronForVTK(foep, pindex, pts))
 280         {
 281         vtkugrid->InsertNextCell(VTK_POLYHEDRON, pts);
 282         }
 283       else
 284         {
 285         foep->SetValid(0);
 286         vtkugrid->Delete();
 287         return NULL;
 288         }
 289       }
 290     else
 291       {
 292       vtkpts->Initialize();
 293       vtkpts->SetNumberOfIds(pts->GetNumberOfIds());
 294
 295       for(vtkIdType node=0; node<pts->GetNumberOfIds(); node++)
 296         {
 297         vtkIdType pid = pts->GetId(node);
 298         vtkIdType ptid=foep->GetVTKPointIndex(pid);
 299         if(ptid < 0 || ptid >= vtkugrid->GetNumberOfPoints())
 300           {
 301           vtkDebugMacro("Index error, this cell"  <<
 302                         " is not on this profile");
 303 #ifndef MedReader_HAVE_PARALLEL_INFRASTRUCTURE
 304           vtkugrid->Delete();
 305           foep->SetValid(0);
 306           return NULL;
 307 #else
 308           return vtkugrid;
 309 #endif
 310           }
 311         vtkpts->SetId(vtkMedUtilities::MedToVTKIndex(vtkType, node), ptid);
 312         }
 313
 314       // for strutural elements, insert nsupportcell instead of only one
 315       if(nsupportcell > 1)
 316         {
 317         vtkIdType npts = vtkpts->GetNumberOfIds() / nsupportcell;
 318         for(int sid=0; sid < nsupportcell; sid++)
 319           {
 320           vtkIdType* ptids = vtkpts->GetPointer(sid * npts);
 321           vtkugrid->InsertNextCell(vtkType, npts, ptids);
 322           }
 323         }
 324       else
 325         {
 326         vtkugrid->InsertNextCell(vtkType, vtkpts);
 327         }
 328       }
 329     }
 330
 331   if(vtkugrid->GetNumberOfCells() == vtkugrid->GetNumberOfPoints())
 332     {
 333     for(int fieldId = 0; fieldId < vtkugrid->GetFieldData()->GetNumberOfArrays(); fieldId++)
 334       {
 335       vtkDataArray* fieldData = vtkugrid->GetFieldData()->GetArray(fieldId);
 336
 337       if(fieldData->GetNumberOfTuples() == vtkugrid->GetNumberOfCells())
 338         {
 339         vtkugrid->GetPointData()->AddArray(fieldData);
 340         }
 341       else
 342         {
 343         vtkDataArray* real_fieldData = fieldData->NewInstance();
 344         real_fieldData->SetName(fieldData->GetName());
 345         real_fieldData->SetNumberOfComponents(fieldData->GetNumberOfComponents());
 346
 347         for(vtkIdType cellId = 0; cellId < vtkugrid->GetNumberOfPoints(); cellId++)
 348           {
 349           vtkIdType supportId = static_cast<int>(supportIndex->GetTuple2(cellId)[0]);
 350           real_fieldData->InsertNextTuple(fieldData->GetTuple(supportId));
 351           }
 352
 353         vtkugrid->GetPointData()->AddArray(real_fieldData);
 354         real_fieldData->Delete();
 355         }
 356       }
 357     }
 358
 359   return vtkugrid;
 360 }
 361
 362 void vtkMedUnstructuredGrid::PrintSelf(ostream& os, vtkIndent indent)
 363 {
 364   this->Superclass::PrintSelf(os, indent);
 365   PRINT_IVAR(os, indent, NumberOfPoints);
 366 }