src/Plugins/MedReader/IO/vtkMedEntityArray.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 "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 = 0;
 189     if(this->Entity.GeometryType == MED_POLYHEDRON || this->Entity.GeometryType == MED_POLYGON)
 190       {
 191       if(!this->NodeIndex || !this->FaceIndex)
 192         return 0;
 193       for(med_int i  = 0; i < this->NumberOfEntity; i++ )
 194         {
 195         med_int start = this->FaceIndex->GetValue(i)-1;
 196         med_int end = this->FaceIndex->GetValue(i+1)-1;
 197         for(med_int fi = start ; fi < end; fi++ )
 198           {
 199           med_int fstart = this->NodeIndex->GetValue(fi)-1;
 200           med_int fend = this->NodeIndex->GetValue(fi+1)-1;
 201           connSize += (fend-fstart);
 202           }
 203         }
 204       } else
 205         {
 206         connSize = this->NumberOfEntity *
 207         vtkMedUtilities::GetNumberOfPoint(this->Entity.GeometryType);
 208         }
 209     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 210     }
 211   else if (this->Connectivity == MED_NODAL && this->Entity.EntityType == MED_STRUCT_ELEMENT)
 212     {
 213     if(this->StructElement == NULL)
 214       return 1;
 215
 216     vtkIdType connSize = this->NumberOfEntity
 217                          * this->StructElement->GetConnectivitySize();
 218
 219     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 220     }
 221   else
 222     {
 223     vtkIdType connSize = this->NumberOfEntity
 224         * vtkMedUtilities::GetNumberOfSubEntity(this->Entity.GeometryType);
 225
 226     return connSize == this->ConnectivityArray->GetNumberOfTuples();
 227     }
 228 }
 229
 230 int vtkMedEntityArray::IsFamilyIdsLoaded()
 231 {
 232   return this->FamilyIdStatus != vtkMedEntityArray::FAMILY_ID_NOT_LOADED;;
 233 }
 234
 235 int vtkMedEntityArray::IsGlobalIdsLoaded()
 236 {
 237   return this->GlobalIds != NULL && this->GlobalIds->GetNumberOfTuples()
 238       == this->NumberOfEntity;
 239 }
 240
 241 void vtkMedEntityArray::GetCellVertices(vtkIdType index, vtkIdList* ids)
 242 {
 243   ids->Initialize();
 244
 245   if(this->Entity.EntityType == MED_NODE)
 246     {
 247     ids->InsertNextId(index);
 248     return;
 249     }
 250
 251   if( this->Entity.EntityType != MED_CELL &&
 252       this->Entity.EntityType != MED_DESCENDING_FACE &&
 253       this->Entity.EntityType != MED_DESCENDING_EDGE &&
 254       this->Entity.EntityType != MED_STRUCT_ELEMENT)
 255     {
 256     vtkErrorMacro("This reader is not compatible with those entities (yet)...");
 257     return;
 258     }
 259
 260   if(vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid) == NULL)
 261     {
 262     // this is a structured grid, connectivity is implicit...
 263
 264     if(this->Entity.GeometryType == MED_POINT1)
 265       {
 266       // degenerate case if there is only one point
 267       ids->InsertNextId(0);
 268       return;
 269       }
 270     if(this->Entity.GeometryType == MED_SEG2)
 271       {
 272       // line
 273       ids->InsertNextId(index);
 274       ids->InsertNextId(index+1);
 275       return;
 276       }
 277     vtkMedRegularGrid* vtkrgrid = vtkMedRegularGrid::SafeDownCast(
 278         this->GetParentGrid());
 279     vtkMedCurvilinearGrid* vtkcgrid = vtkMedCurvilinearGrid::SafeDownCast(
 280         this->GetParentGrid());
 281     vtkIdType xncell = 0;
 282     vtkIdType yncell = 0;
 283     vtkIdType zncell = 0;
 284     vtkIdType xnpts = 1;
 285     vtkIdType ynpts = 1;
 286     vtkIdType znpts = 1;
 287     if(vtkrgrid!=NULL)
 288       {
 289       xncell = vtkrgrid->GetAxisSize(0)-1;
 290       yncell = vtkrgrid->GetAxisSize(1)-1;
 291       zncell = vtkrgrid->GetAxisSize(2)-1;
 292       xnpts = vtkrgrid->GetAxisSize(0);
 293       ynpts = vtkrgrid->GetAxisSize(1);
 294       znpts = vtkrgrid->GetAxisSize(2);
 295       }
 296     if(vtkcgrid != NULL)
 297       {
 298       xncell = vtkcgrid->GetAxisSize(0)-1;
 299       yncell = vtkcgrid->GetAxisSize(1)-1;
 300       zncell = vtkcgrid->GetAxisSize(2)-1;
 301       xnpts = vtkcgrid->GetAxisSize(0);
 302       ynpts = vtkcgrid->GetAxisSize(1);
 303       znpts = vtkcgrid->GetAxisSize(2);
 304       }
 305     vtkIdType xindex = index % xncell;
 306     if(xncell <= 0)
 307       return;
 308
 309     vtkIdType yindex = index / xncell;
 310
 311     if(this->Entity.GeometryType == MED_QUAD4)
 312       {
 313       // plane
 314
 315       ids->InsertNextId(xindex + yindex*xnpts);
 316       ids->InsertNextId(xindex + 1 + yindex*xnpts);
 317       ids->InsertNextId(xindex + yindex*xnpts);
 318       ids->InsertNextId(xindex + 1 + (yindex + 1)*xnpts);
 319       return;
 320       }
 321
 322     if(yncell <= 0)
 323       return;
 324
 325     vtkIdType zindex = index / (xncell*yncell);
 326
 327     if(this->Entity.GeometryType == MED_HEXA8)
 328       {
 329       // volume
 330       ids->InsertNextId(xindex   + (yindex  )*xnpts + (zindex  )*xnpts*ynpts);
 331       ids->InsertNextId(xindex+1 + (yindex  )*xnpts + (zindex  )*xnpts*ynpts);
 332       ids->InsertNextId(xindex   + (yindex+1)*xnpts + (zindex  )*xnpts*ynpts);
 333       ids->InsertNextId(xindex+1 + (yindex+1)*xnpts + (zindex  )*xnpts*ynpts);
 334       ids->InsertNextId(xindex   + (yindex  )*xnpts + (zindex+1)*xnpts*ynpts);
 335       ids->InsertNextId(xindex+1 + (yindex  )*xnpts + (zindex+1)*xnpts*ynpts);
 336       ids->InsertNextId(xindex   + (yindex+1)*xnpts + (zindex+1)*xnpts*ynpts);
 337       ids->InsertNextId(xindex+1 + (yindex+1)*xnpts + (zindex+1)*xnpts*ynpts);
 338       return;
 339       }
 340     return;
 341     }
 342
 343   this->LoadConnectivity();
 344
 345   if (this->GetEntity().GeometryType==MED_POLYHEDRON)
 346     {
 347     vtkMedIntArray* conn = this->GetConnectivityArray();
 348     vtkMedIntArray* faceIndex = this->GetFaceIndex();
 349     vtkMedIntArray* nodeIndex = this->GetNodeIndex();
 350     med_int start = faceIndex->GetValue(index)-1;
 351     med_int end = faceIndex->GetValue(index+1)-1;
 352     // the use of a set loses the order, but VTK do not support this order anyway.
 353     if (this->GetConnectivity()==MED_NODAL)
 354       {
 355       for (int ff = start; ff<end; ff++)
 356         {
 357         med_int fstart = nodeIndex->GetValue(ff)-1;
 358         med_int fend = nodeIndex->GetValue(ff+1)-1;
 359         for (int pt = fstart; pt<fend; pt++)
 360           {
 361           med_int ptid = conn->GetValue(pt)-1;
 362           ids->InsertNextId(ptid);
 363           }
 364         }
 365       }
 366     else // MED_DESCENDING
 367       {
 368       vtkMedUnstructuredGrid* ugrid =
 369           vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 370       if (!ugrid)
 371         {
 372         vtkErrorMacro(
 373         "MED_DESCENDING connectivity is only supported on unstructured grids");
 374         return;
 375         }
 376       set<med_int> pts;
 377       vtkIdList* subIds = vtkIdList::New();
 378       for (int ff = start; ff<end; ff++)
 379         {
 380         med_int fid = conn->GetValue(ff)-1;
 381         vtkMedEntity entity;
 382         entity.GeometryType = (med_geometry_type) NodeIndex->GetValue(ff);
 383         entity.EntityType = MED_DESCENDING_FACE;
 384         vtkMedEntityArray* subarray = ugrid->GetEntityArray(entity);
 385         subarray->GetCellVertices(fid, subIds);
 386         for (int id = 0; id<subIds->GetNumberOfIds(); id++)
 387           {
 388           med_int ptid = subIds->GetId(id);
 389           if(pts.find(ptid) == pts.end())
 390             {
 391             ids->InsertNextId(ptid);
 392             pts.insert(ptid);
 393             }
 394           }
 395         }
 396       subIds->Delete();
 397       }
 398     }//end polyhedron
 399   else if (this->GetEntity().GeometryType==MED_POLYGON)
 400     {
 401     vtkMedIntArray* conn = this->GetConnectivityArray();
 402     vtkMedIntArray* nids = this->GetFaceIndex();
 403     med_int start = nids->GetValue(index)-1;
 404     med_int end = nids->GetValue(index+1)-1;
 405     if (this->GetConnectivity()==MED_NODAL)
 406       {
 407       for (int pt = start; pt<end; pt++)
 408         {
 409         ids->InsertNextId(conn->GetValue(pt)-1);
 410         }
 411       }
 412     else // MED_DESCENDING
 413       {
 414       vtkIdList* subpts=vtkIdList::New();
 415       vtkMedUnstructuredGrid* ugrid =
 416           vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 417       if (!ugrid)
 418         {
 419         vtkErrorMacro("MED_DESCENDING connectivity is only "
 420                       << "supported on unstructured grids");
 421         return;
 422         }
 423       set<med_int> pts;
 424       for (int sub = start; sub<end; sub++)
 425         {
 426         med_int subid = conn->GetValue(sub)-1;
 427         vtkMedEntity subentity;
 428         subentity.GeometryType = MED_SEG2;
 429         subentity.EntityType = MED_DESCENDING_EDGE;
 430         vtkMedEntityArray* subarray = ugrid->GetEntityArray(subentity);
 431         subarray->GetCellVertices(subid, subpts);
 432         for(int id=0; id<subpts->GetNumberOfIds(); id++)
 433           {
 434           med_int ptid = subpts->GetId(id);
 435           if(pts.find(ptid) != pts.end())
 436             {
 437             pts.insert(ptid);
 438             ids->InsertNextId(ptid);
 439             }
 440           }
 441         }
 442       subpts->Delete();
 443       }
 444     }//end poygon
 445   else if (this->GetConnectivity()==MED_NODAL ||
 446            vtkMedUtilities::GetDimension(this->GetEntity().GeometryType)<1)
 447     {
 448     int npts = 0;
 449     if(this->GetEntity().EntityType == MED_STRUCT_ELEMENT)
 450       {
 451       if(this->StructElement != NULL)
 452         {
 453         npts = this->StructElement->GetConnectivitySize();
 454         }
 455       }
 456     else
 457       {
 458       npts = vtkMedUtilities::GetNumberOfPoint(this->GetEntity().GeometryType);
 459       }
 460     vtkMedIntArray* conn = this->GetConnectivityArray();
 461     for (int i = 0; i<npts; i++)
 462       {
 463       vtkIdType ptid = conn->GetValue(npts*index+i)-1;
 464       ids->InsertNextId(ptid);
 465       }
 466     }//end nodal case
 467   else
 468     {
 469     vtkIdList* subpts=vtkIdList::New();
 470     int nsub=vtkMedUtilities::GetNumberOfSubEntity(
 471         this->GetEntity().GeometryType);
 472     vtkMedUnstructuredGrid* ugrid =
 473         vtkMedUnstructuredGrid::SafeDownCast(this->ParentGrid);
 474     if (!ugrid)
 475       {
 476       vtkErrorMacro(
 477         "MED_DESCENDING connectivity is only supported on unstructured grids");
 478       return;
 479       }
 480     vtkMedIntArray* conn=this->GetConnectivityArray();
 481     ids->SetNumberOfIds(vtkMedUtilities::GetNumberOfPoint(
 482         this->GetEntity().GeometryType));
 483     for (int sub = 0; sub<nsub; sub++)
 484       {
 485       med_int subid = conn->GetValue(nsub*index+sub)-1;
 486       bool invert = false;
 487       if(subid < 0)
 488         {
 489         subid = -subid;
 490         invert = true;
 491         }
 492
 493       vtkMedEntity subentity;
 494       subentity.GeometryType = vtkMedUtilities::GetSubGeometry(
 495           this->GetEntity().GeometryType, sub);
 496       subentity.EntityType = vtkMedUtilities::GetSubType(
 497           this->GetEntity().EntityType);
 498       vtkMedEntityArray* subarray = ugrid->GetEntityArray(subentity);
 499       if(subarray == NULL)
 500         {
 501         subentity.EntityType = MED_CELL;
 502         subarray = ugrid->GetEntityArray(subentity);
 503         }
 504       if(subarray == NULL)
 505         {
 506         vtkDebugMacro( << "Missing sub entity array " << subentity.GeometryType);
 507         this->Valid = false;
 508         break;
 509         }
 510       subarray->GetCellVertices(subid, subpts);
 511       vtkMedUtilities::ProjectConnectivity(this->GetEntity().GeometryType, ids, subpts,
 512           sub, invert);
 513       }
 514     subpts->Delete();
 515     }
 516 }
 517
 518 void  vtkMedEntityArray::LoadConnectivity()
 519 {
 520   if(this->IsConnectivityLoaded())
 521     return;
 522
 523   this->GetParentGrid()->GetParentMesh()->GetParentFile()->GetMedDriver()
 524       ->LoadConnectivity(this);
 525 }
 526
 527 void  vtkMedEntityArray::SetVariableAttributeValues(
 528     vtkMedVariableAttribute* varatt, vtkAbstractArray* value)
 529 {
 530   this->VariableAttributeValue[varatt] = value;
 531 }
 532
 533 vtkAbstractArray* vtkMedEntityArray::GetVariableAttributeValue(
 534     vtkMedVariableAttribute* varatt)
 535 {
 536   if(this->VariableAttributeValue.find(varatt)
 537     == this->VariableAttributeValue.end())
 538     return NULL;
 539
 540   return this->VariableAttributeValue[varatt];
 541 }
 542
 543 void vtkMedEntityArray::PrintSelf(ostream& os, vtkIndent indent)
 544 {
 545   this->Superclass::PrintSelf(os, indent);
 546   PRINT_IVAR(os, indent, NumberOfEntity)
 547   PRINT_IVAR(os, indent, Connectivity)
 548   PRINT_IVAR(os, indent, InitialGlobalId)
 549 }