src/Plugins/MedReader/IO/vtkMedEntityArray.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 "vtkMedEntityArray.h"
  21
  22 #include "vtkObjectFactory.h"
  23 #include "vtkMedIntArray.h"
  24 #include "vtkMedUtilities.h"
  25 #include "vtkMedFamily.h"
  26 #include "vtkMedFamilyOnEntity.h"
  27 #include "vtkMedMesh.h"
  28 #include "vtkMedGrid.h"
  29 #include "vtkMedUnstructuredGrid.h"
  30 #include "vtkMedRegularGrid.h"
  31 #include "vtkMedCurvilinearGrid.h"
  32 #include "vtkMedFile.h"
  33 #include "vtkMedDriver.h"
  34 #include "vtkMedStructElement.h"
  35
  36 #include "vtkIdList.h"
  37
  38 #include <set>
  39 using std::set;
  40
  41 vtkCxxSetObjectVectorMacro(vtkMedEntityArray, FamilyOnEntity, vtkMedFamilyOnEntity);
  42 vtkCxxGetObjectVectorMacro(vtkMedEntityArray, FamilyOnEntity, vtkMedFamilyOnEntity);
  43 vtkCxxSetObjectMacro(vtkMedEntityArray,FamilyIds,vtkMedIntArray);
  44 vtkCxxSetObjectMacro(vtkMedEntityArray,GlobalIds,vtkMedIntArray);
  45 vtkCxxSetObjectMacro(vtkMedEntityArray,ConnectivityArray,vtkMedIntArray);
  46 vtkCxxSetObjectMacro(vtkMedEntityArray,FaceIndex,vtkMedIntArray);
  47 vtkCxxSetObjectMacro(vtkMedEntityArray,NodeIndex,vtkMedIntArray);
  48
  49 vtkCxxSetObjectMacro(vtkMedEntityArray,ParentGrid,vtkMedGrid);
  50 vtkCxxSetObjectMacro(vtkMedEntityArray,StructElement,vtkMedStructElement);
  51
  52 // vtkCxxRevisionMacro(vtkMedEntityArray, "$Revision$");
  53 vtkStandardNewMacro(vtkMedEntityArray);
  54
  55 vtkMedEntityArray::vtkMedEntityArray()
  56 {
  57   this->NumberOfEntity = 0;
  58   this->Connectivity = MED_NODAL;
  59   this->FamilyIds = NULL;
  60   this->GlobalIds = NULL;
  61   this->ConnectivityArray = NULL;
  62   this->FaceIndex = NULL;
  63   this->NodeIndex = NULL;
  64   this->InitialGlobalId = 0;
  65   this->FamilyOnEntity = new vtkObjectVector<vtkMedFamilyOnEntity> ();
  66   this->FamilyIdStatus = vtkMedEntityArray::FAMILY_ID_NOT_LOADED;
  67   this->ParentGrid = NULL;
  68   this->StructElement = NULL;
  69   //this->Filter = NULL;
  70 }
  71
  72 vtkMedEntityArray::~vtkMedEntityArray()
  73 {
  74   this->SetFamilyIds(NULL);
  75   this->SetGlobalIds(NULL);
  76   this->SetConnectivityArray(NULL);
  77   this->SetFaceIndex(NULL);
  78   this->SetNodeIndex(NULL);
  79   delete this->FamilyOnEntity;
  80   this->SetParentGrid(NULL);
  81   this->SetStructElement(NULL);
  82   //this->SetFilter(NULL);
  83 }
  84
  85 void vtkMedEntityArray::Initialize()
  86 {
  87   this->SetFamilyIds(NULL);
  88   this->SetGlobalIds(NULL);
  89   this->SetConnectivityArray(NULL);
  90   this->SetFaceIndex(NULL);
  91   this->SetNodeIndex(NULL);
  92   this->FamilyOnEntity->clear();
  93   this->FamilyIdStatus = FAMILY_ID_NOT_LOADED;
  94 }
  95
  96 void vtkMedEntityArray::ComputeFamilies()
  97 {
  98   this->FamilyOnEntity->clear();
  99   vtkMedMesh* mesh = this->ParentGrid->GetParentMesh();
 100
 101   if(this->FamilyIds == NULL)
 102     {
 103     vtkMedFamilyOnEntity* foe = vtkMedFamilyOnEntity::New();
 104     foe->SetParentGrid(this->ParentGrid);
 105     this->AppendFamilyOnEntity(foe);
 106     foe->Delete();
 107     if(this->GetEntity().EntityType != MED_NODE)
 108       {
 109       foe->SetFamily(mesh->GetOrCreateCellFamilyById(0));
 110       }
 111     else
 112       {
 113       foe->SetFamily(mesh->GetOrCreatePointFamilyById(0));
 114       }
 115     foe->SetEntityArray(this);
 116     this->FamilyIdStatus = vtkMedEntityArray::FAMILY_ID_IMPLICIT;
 117     return;
 118     }
 119
 120   this->FamilyIdStatus = vtkMedEntityArray::FAMILY_ID_EXPLICIT;
 121
 122   set<med_int> idset;
 123   for (vtkIdType index = 0; index < this->FamilyIds->GetNumberOfTuples(); index++)
 124     {
 125     med_int id = this->FamilyIds->GetValue(index);
 126     idset.insert(id);
 127     }
 128
 129   for (set<med_int>::iterator it = idset.begin(); it != idset.end(); it++)
 130     {
 131     vtkMedFamilyOnEntity* foe = vtkMedFamilyOnEntity::New();
 132     foe->SetParentGrid(this->ParentGrid);
 133     this->AppendFamilyOnEntity(foe);
 134     foe->Delete();
 135     if(this->GetEntity().EntityType != MED_NODE)
 136       {
 137       foe->SetFamily(mesh->GetOrCreateCellFamilyById(*it));
 138       }
 139     else
 140       {
 141       foe->SetFamily(mesh->GetOrCreatePointFamilyById(*it));
 142       }
 143     foe->SetEntityArray(this);
 144     }
 145 }
 146
 147 med_int vtkMedEntityArray::GetFamilyId(med_int id)
 148 {
 149   if(this->FamilyIdStatus == FAMILY_ID_IMPLICIT)
 150     return 0;
 151   if(this->FamilyIdStatus == FAMILY_ID_NOT_LOADED)
 152     {
 153     vtkErrorMacro("You have to load family ids before asking for it!");
 154     }
 155   return this->FamilyIds->GetValue(id);
 156 }
 157
 158 int vtkMedEntityArray::HasFamily(vtkMedFamily* family)
 159 {
 160   for (int i = 0; i < this->FamilyOnEntity->size(); i++)
 161     {
 162     vtkMedFamilyOnEntity* foe = this->FamilyOnEntity->at(i);
 163     if(foe->GetFamily() == family)
 164       return 1;
 165     }
 166   return 0;
 167 }
 168
 169 int vtkMedEntityArray::IsConnectivityLoaded()
 170 {
 171   // Entity Arrays representing something else than cells
 172   // have no connectivity
 173
 174   if(vtkMedUnstructuredGrid::SafeDownCast(this->GetParentGrid()) == NULL)
 175     return 1;
 176
 177   if( this->Entity.EntityType != MED_CELL &&
 178       this->Entity.EntityType != MED_DESCENDING_FACE &&
 179       this->Entity.EntityType != MED_DESCENDING_EDGE &&
 180       this->Entity.EntityType != MED_STRUCT_ELEMENT)
 181     return 1;
 182
 183   if(this->ConnectivityArray == NULL)
 184     return 0;
 185
 186   if(this->Connectivity == MED_NODAL && this->Entity.EntityType != MED_STRUCT_ELEMENT)
 187     {
 188     vtkIdType connSize = this->NumberOfEntity
 189         * vtkMedUtilities::GetNumberOfPoint(this->Entity.GeometryType);
 190
 191     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 192     }
 193   else if (this->Connectivity == MED_NODAL && this->Entity.EntityType == MED_STRUCT_ELEMENT)
 194     {
 195     if(this->StructElement == NULL)
 196       return 1;
 197
 198     vtkIdType connSize = this->NumberOfEntity
 199                          * this->StructElement->GetConnectivitySize();
 200
 201     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 202     }
 203   else
 204     {
 205     vtkIdType connSize = this->NumberOfEntity
 206         * vtkMedUtilities::GetNumberOfSubEntity(this->Entity.GeometryType);
 207
 208     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 209     }
 210 }
 211
 212 int vtkMedEntityArray::IsFamilyIdsLoaded()
 213 {
 214   return this->FamilyIdStatus != vtkMedEntityArray::FAMILY_ID_NOT_LOADED;;
 215 }
 216
 217 int vtkMedEntityArray::IsGlobalIdsLoaded()
 218 {
 219   return this->GlobalIds != NULL && this->GlobalIds->GetNumberOfTuples()
 220       == this->NumberOfEntity;
 221 }
 222
 223 void vtkMedEntityArray::GetCellVertices(vtkIdType index, vtkIdList* ids)
 224 {
 225   ids->Initialize();
 226
 227   if(this->Entity.EntityType == MED_NODE)
 228     {
 229     ids->InsertNextId(index);
 230     return;
 231     }
 232
 233   if( this->Entity.EntityType != MED_CELL &&
 234       this->Entity.EntityType != MED_DESCENDING_FACE &&
 235       this->Entity.EntityType != MED_DESCENDING_EDGE &&
 236       this->Entity.EntityType != MED_STRUCT_ELEMENT)
 237     {
 238     vtkErrorMacro("This reader is not compatible with those entities (yet)...");
 239     return;
 240     }
 241
 242   if(vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid) == NULL)
 243     {
 244     // this is a structured grid, connectivity is implicit...
 245
 246     if(this->Entity.GeometryType == MED_POINT1)
 247       {
 248       // degenerate case if there is only one point
 249       ids->InsertNextId(0);
 250       return;
 251       }
 252     if(this->Entity.GeometryType == MED_SEG2)
 253       {
 254       // line
 255       ids->InsertNextId(index);
 256       ids->InsertNextId(index+1);
 257       return;
 258       }
 259     vtkMedRegularGrid* vtkrgrid = vtkMedRegularGrid::SafeDownCast(
 260         this->GetParentGrid());
 261     vtkMedCurvilinearGrid* vtkcgrid = vtkMedCurvilinearGrid::SafeDownCast(
 262         this->GetParentGrid());
 263     vtkIdType xncell = 0;
 264     vtkIdType yncell = 0;
 265     vtkIdType zncell = 0;
 266     vtkIdType xnpts = 1;
 267     vtkIdType ynpts = 1;
 268     vtkIdType znpts = 1;
 269     if(vtkrgrid!=NULL)
 270       {
 271       xncell = vtkrgrid->GetAxisSize(0)-1;
 272       yncell = vtkrgrid->GetAxisSize(1)-1;
 273       zncell = vtkrgrid->GetAxisSize(2)-1;
 274       xnpts = vtkrgrid->GetAxisSize(0);
 275       ynpts = vtkrgrid->GetAxisSize(1);
 276       znpts = vtkrgrid->GetAxisSize(2);
 277       }
 278     if(vtkcgrid != NULL)
 279       {
 280       xncell = vtkcgrid->GetAxisSize(0)-1;
 281       yncell = vtkcgrid->GetAxisSize(1)-1;
 282       zncell = vtkcgrid->GetAxisSize(2)-1;
 283       xnpts = vtkcgrid->GetAxisSize(0);
 284       ynpts = vtkcgrid->GetAxisSize(1);
 285       znpts = vtkcgrid->GetAxisSize(2);
 286       }
 287     vtkIdType xindex = index % xncell;
 288     if(xncell <= 0)
 289       return;
 290
 291     vtkIdType yindex = index / xncell;
 292
 293     if(this->Entity.GeometryType == MED_QUAD4)
 294       {
 295       // plane
 296
 297       ids->InsertNextId(xindex + yindex*xnpts);
 298       ids->InsertNextId(xindex + 1 + yindex*xnpts);
 299       ids->InsertNextId(xindex + yindex*xnpts);
 300       ids->InsertNextId(xindex + 1 + (yindex + 1)*xnpts);
 301       return;
 302       }
 303
 304     if(yncell <= 0)
 305       return;
 306
 307     vtkIdType zindex = index / (xncell*yncell);
 308
 309     if(this->Entity.GeometryType == MED_HEXA8)
 310       {
 311       // volume
 312       ids->InsertNextId(xindex   + (yindex  )*xnpts + (zindex  )*xnpts*ynpts);
 313       ids->InsertNextId(xindex+1 + (yindex  )*xnpts + (zindex  )*xnpts*ynpts);
 314       ids->InsertNextId(xindex   + (yindex+1)*xnpts + (zindex  )*xnpts*ynpts);
 315       ids->InsertNextId(xindex+1 + (yindex+1)*xnpts + (zindex  )*xnpts*ynpts);
 316       ids->InsertNextId(xindex   + (yindex  )*xnpts + (zindex+1)*xnpts*ynpts);
 317       ids->InsertNextId(xindex+1 + (yindex  )*xnpts + (zindex+1)*xnpts*ynpts);
 318       ids->InsertNextId(xindex   + (yindex+1)*xnpts + (zindex+1)*xnpts*ynpts);
 319       ids->InsertNextId(xindex+1 + (yindex+1)*xnpts + (zindex+1)*xnpts*ynpts);
 320       return;
 321       }
 322     return;
 323     }
 324
 325   this->LoadConnectivity();
 326
 327   if (this->GetEntity().GeometryType==MED_POLYHEDRON)
 328     {
 329     vtkMedIntArray* conn = this->GetConnectivityArray();
 330     vtkMedIntArray* faceIndex = this->GetFaceIndex();
 331     vtkMedIntArray* nodeIndex = this->GetNodeIndex();
 332     med_int start = faceIndex->GetValue(index)-1;
 333     med_int end = faceIndex->GetValue(index+1)-1;
 334     // the use of a set loses the order, but VTK do not support this order anyway.
 335     if (this->GetConnectivity()==MED_NODAL)
 336       {
 337       for (int ff = start; ff<end; ff++)
 338         {
 339         med_int fstart = nodeIndex->GetValue(ff)-1;
 340         med_int fend = nodeIndex->GetValue(ff+1)-1;
 341         for (int pt = fstart; pt<fend; pt++)
 342           {
 343           med_int ptid = conn->GetValue(pt)-1;
 344           ids->InsertNextId(ptid);
 345           }
 346         }
 347       }
 348     else // MED_DESCENDING
 349       {
 350       vtkMedUnstructuredGrid* ugrid =
 351           vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 352       if (!ugrid)
 353         {
 354         vtkErrorMacro(
 355         "MED_DESCENDING connectivity is only supported on unstructured grids");
 356         return;
 357         }
 358       set<med_int> pts;
 359       vtkIdList* subIds = vtkIdList::New();
 360       for (int ff = start; ff<end; ff++)
 361         {
 362         med_int fid = conn->GetValue(ff)-1;
 363         vtkMedEntity entity;
 364         entity.GeometryType = (med_geometry_type) NodeIndex->GetValue(ff);
 365         entity.EntityType = MED_DESCENDING_FACE;
 366         vtkMedEntityArray* subarray = ugrid->GetEntityArray(entity);
 367         subarray->GetCellVertices(fid, subIds);
 368         for (int id = 0; id<subIds->GetNumberOfIds(); id++)
 369           {
 370           med_int ptid = subIds->GetId(id);
 371           if(pts.find(ptid) == pts.end())
 372             {
 373             ids->InsertNextId(ptid);
 374             pts.insert(ptid);
 375             }
 376           }
 377         }
 378       subIds->Delete();
 379       }
 380     }//end polyhedron
 381   else if (this->GetEntity().GeometryType==MED_POLYGON)
 382     {
 383     vtkMedIntArray* conn = this->GetConnectivityArray();
 384     vtkMedIntArray* nids = this->GetFaceIndex();
 385     med_int start = nids->GetValue(index)-1;
 386     med_int end = nids->GetValue(index+1)-1;
 387     if (this->GetConnectivity()==MED_NODAL)
 388       {
 389       for (int pt = start; pt<end; pt++)
 390         {
 391         ids->InsertNextId(conn->GetValue(pt)-1);
 392         }
 393       }
 394     else // MED_DESCENDING
 395       {
 396       vtkIdList* subpts=vtkIdList::New();
 397       vtkMedUnstructuredGrid* ugrid =
 398           vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 399       if (!ugrid)
 400         {
 401         vtkErrorMacro("MED_DESCENDING connectivity is only "
 402                       << "supported on unstructured grids");
 403         return;
 404         }
 405       set<med_int> pts;
 406       for (int sub = start; sub<end; sub++)
 407         {
 408         med_int subid = conn->GetValue(sub)-1;
 409         vtkMedEntity subentity;
 410         subentity.GeometryType = MED_SEG2;
 411         subentity.EntityType = MED_DESCENDING_EDGE;
 412         vtkMedEntityArray* subarray = ugrid->GetEntityArray(subentity);
 413         subarray->GetCellVertices(subid, subpts);
 414         for(int id=0; id<subpts->GetNumberOfIds(); id++)
 415           {
 416           med_int ptid = subpts->GetId(id);
 417           if(pts.find(ptid) != pts.end())
 418             {
 419             pts.insert(ptid);
 420             ids->InsertNextId(ptid);
 421             }
 422           }
 423         }
 424       subpts->Delete();
 425       }
 426     }//end poygon
 427   else if (this->GetConnectivity()==MED_NODAL ||
 428            vtkMedUtilities::GetDimension(this->GetEntity().GeometryType)<1)
 429     {
 430     int npts = 0;
 431     if(this->GetEntity().EntityType == MED_STRUCT_ELEMENT)
 432       {
 433       if(this->StructElement != NULL)
 434         {
 435         npts = this->StructElement->GetConnectivitySize();
 436         }
 437       }
 438     else
 439       {
 440       npts = vtkMedUtilities::GetNumberOfPoint(this->GetEntity().GeometryType);
 441       }
 442     vtkMedIntArray* conn = this->GetConnectivityArray();
 443     for (int i = 0; i<npts; i++)
 444       {
 445       vtkIdType ptid = conn->GetValue(npts*index+i)-1;
 446       ids->InsertNextId(ptid);
 447       }
 448     }//end nodal case
 449   else
 450     {
 451     vtkIdList* subpts=vtkIdList::New();
 452     int nsub=vtkMedUtilities::GetNumberOfSubEntity(
 453         this->GetEntity().GeometryType);
 454     vtkMedUnstructuredGrid* ugrid =
 455         vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 456     if (!ugrid)
 457       {
 458       vtkErrorMacro(
 459         "MED_DESCENDING connectivity is only supported on unstructured grids");
 460       return;
 461       }
 462     vtkMedIntArray* conn=this->GetConnectivityArray();
 463     ids->SetNumberOfIds(vtkMedUtilities::GetNumberOfPoint(
 464         this->GetEntity().GeometryType));
 465     for (int sub = 0; sub<nsub; sub++)
 466       {
 467       med_int subid = conn->GetValue(nsub*index+sub);
 468       bool invert = false;
 469       if(subid < 0)
 470         {
 471         subid = -subid;
 472         invert = true;
 473         }
 474       subid = subid-1;
 475
 476       vtkMedEntity subentity;
 477       subentity.GeometryType = vtkMedUtilities::GetSubGeometry(
 478           this->GetEntity().GeometryType, sub);
 479       subentity.EntityType = vtkMedUtilities::GetSubType(
 480           this->GetEntity().EntityType);
 481       vtkMedEntityArray* subarray = ugrid->GetEntityArray(subentity);
 482       if(subarray == NULL)
 483         {
 484         subentity.EntityType = MED_CELL;
 485         subarray = ugrid->GetEntityArray(subentity);
 486         }
 487       if(subarray == NULL)
 488         {
 489         vtkDebugMacro( << "Missing sub entity array " << subentity.GeometryType);
 490         this->Valid = false;
 491         break;
 492         }
 493       subarray->GetCellVertices(subid, subpts);
 494       vtkMedUtilities::ProjectConnectivity(this->GetEntity().GeometryType, ids, subpts,
 495           sub, invert);
 496       }
 497     subpts->Delete();
 498     }
 499 }
 500
 501 void  vtkMedEntityArray::LoadConnectivity()
 502 {
 503   if(this->IsConnectivityLoaded())
 504     return;
 505
 506   this->GetParentGrid()->GetParentMesh()->GetParentFile()->GetMedDriver()
 507       ->LoadConnectivity(this);
 508 }
 509
 510 void  vtkMedEntityArray::SetVariableAttributeValues(
 511     vtkMedVariableAttribute* varatt, vtkAbstractArray* value)
 512 {
 513   this->VariableAttributeValue[varatt] = value;
 514 }
 515
 516 vtkAbstractArray* vtkMedEntityArray::GetVariableAttributeValue(
 517     vtkMedVariableAttribute* varatt)
 518 {
 519   if(this->VariableAttributeValue.find(varatt)
 520     == this->VariableAttributeValue.end())
 521     return NULL;
 522
 523   return this->VariableAttributeValue[varatt];
 524 }
 525
 526 void vtkMedEntityArray::PrintSelf(ostream& os, vtkIndent indent)
 527 {
 528   this->Superclass::PrintSelf(os, indent);
 529   PRINT_IVAR(os, indent, NumberOfEntity)
 530   PRINT_IVAR(os, indent, Connectivity)
 531   PRINT_IVAR(os, indent, InitialGlobalId)
 532 }