src/Plugins/MedReader/IO/vtkMedUnstructuredGrid.cxx

   1 // Copyright (C) 2010-2012  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       for(int varattid = 0; varattid<structelem->GetNumberOfVariableAttribute(); varattid++)
 220         {
 221         vtkMedVariableAttribute* varatt = structelem->GetVariableAttribute(varattid);
 222         varatt->Load(array);
 223         vtkAbstractArray* values = array->GetVariableAttributeValue(varatt);
 224         vtkugrid->GetFieldData()->AddArray(values);
 225         }
 226       vtkugrid->GetCellData()->AddArray(supportIndex);
 227       }
 228     }
 229   else
 230     {
 231     vtkType = vtkMedUtilities::GetVTKCellType(
 232         array->GetEntity().GeometryType);
 233     }
 234
 235   vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
 236   vtkSmartPointer<vtkIdList> vtkpts = vtkSmartPointer<vtkIdList>::New();
 237   vtkIdType intialGlobalId = array->GetInitialGlobalId();
 238   vtkMedIntArray* pids = (foep->GetProfile()!=NULL?
 239                            foep->GetProfile()->GetIds() : NULL);
 240   vtkIdType maxId = (pids!=NULL?
 241                      pids->GetNumberOfTuples():
 242                      array->GetNumberOfEntity());
 243
 244   int valid = foep->GetValid();
 245   if (controller != NULL)
 246     if (controller->GetNumberOfProcesses() > 1)
 247     valid = 1;
 248
 249   for (vtkIdType pindex = 0; pindex<maxId && valid; pindex++)
 250     {
 251     vtkIdType realIndex = (pids!=NULL?
 252                            pids->GetValue(pindex)-1:
 253                            pindex);
 254
 255     if (!foep->KeepCell(realIndex))
 256       continue;
 257
 258     array->GetCellVertices(realIndex, pts);
 259
 260     for(int sid = 0; sid < nsupportcell; sid++)
 261       {
 262       cellGlobalIds->InsertNextValue(intialGlobalId+pindex);
 263       }
 264
 265     // The supportIndex array has 2 component, the first to give the index
 266     // of the med cell in the variable attributes array, the second to give
 267     // the index of the cell in the n cells composing the structural element
 268     if(array->GetEntity().EntityType == MED_STRUCT_ELEMENT)
 269       {
 270       for(int sid = 0; sid < nsupportcell; sid++)
 271         {
 272         supportIndex->InsertNextTuple2(pindex, sid);
 273         }
 274       }
 275
 276     if (array->GetEntity().GeometryType==MED_POLYHEDRON)
 277       {
 278       if(vtkMedUtilities::FormatPolyhedronForVTK(foep, pindex, pts))
 279         {
 280         vtkugrid->InsertNextCell(VTK_POLYHEDRON, pts);
 281         }
 282       else
 283         {
 284         foep->SetValid(0);
 285         vtkugrid->Delete();
 286         return NULL;
 287         }
 288       }
 289     else
 290       {
 291       vtkpts->Initialize();
 292       vtkpts->SetNumberOfIds(pts->GetNumberOfIds());
 293
 294       for(vtkIdType node=0; node<pts->GetNumberOfIds(); node++)
 295         {
 296         vtkIdType pid = pts->GetId(node);
 297         vtkIdType ptid=foep->GetVTKPointIndex(pid);
 298         if(ptid < 0 || ptid >= vtkugrid->GetNumberOfPoints())
 299           {
 300           vtkDebugMacro("Index error, this cell"  <<
 301                         " is not on this profile");
 302 #ifndef MedReader_HAVE_PARALLEL_INFRASTRUCTURE
 303           vtkugrid->Delete();
 304           foep->SetValid(0);
 305           return NULL;
 306 #else
 307           return vtkugrid;
 308 #endif
 309           }
 310         vtkpts->SetId(vtkMedUtilities::MedToVTKIndex(vtkType, node), ptid);
 311         }
 312
 313       // for strutural elements, insert nsupportcell instead of only one
 314       if(nsupportcell > 1)
 315         {
 316         vtkIdType npts = vtkpts->GetNumberOfIds() / nsupportcell;
 317         for(int sid=0; sid < nsupportcell; sid++)
 318           {
 319           vtkIdType* ptids = vtkpts->GetPointer(sid * npts);
 320           vtkugrid->InsertNextCell(vtkType, npts, ptids);
 321           }
 322         }
 323       else
 324         {
 325         vtkugrid->InsertNextCell(vtkType, vtkpts);
 326         }
 327       }
 328     }
 329
 330   if(vtkugrid->GetNumberOfCells() == vtkugrid->GetNumberOfPoints())
 331     {
 332     for(int fieldId = 0; fieldId < vtkugrid->GetFieldData()->GetNumberOfArrays(); fieldId++)
 333       {
 334       vtkDataArray* fieldData = vtkugrid->GetFieldData()->GetArray(fieldId);
 335
 336       if(fieldData->GetNumberOfTuples() == vtkugrid->GetNumberOfCells())
 337         {
 338         vtkugrid->GetPointData()->AddArray(fieldData);
 339         }
 340       else
 341         {
 342         vtkDataArray* real_fieldData = fieldData->NewInstance();
 343         real_fieldData->SetName(fieldData->GetName());
 344         real_fieldData->SetNumberOfComponents(fieldData->GetNumberOfComponents());
 345
 346         for(vtkIdType cellId = 0; cellId < vtkugrid->GetNumberOfPoints(); cellId++)
 347           {
 348           vtkIdType supportId = static_cast<int>(supportIndex->GetTuple2(cellId)[0]);
 349           real_fieldData->InsertNextTuple(fieldData->GetTuple(supportId));
 350           }
 351
 352         vtkugrid->GetPointData()->AddArray(real_fieldData);
 353         real_fieldData->Delete();
 354         }
 355       }
 356     }
 357
 358   return vtkugrid;
 359 }
 360
 361 void vtkMedUnstructuredGrid::PrintSelf(ostream& os, vtkIndent indent)
 362 {
 363   this->Superclass::PrintSelf(os, indent);
 364   PRINT_IVAR(os, indent, NumberOfPoints);
 365 }