src/Plugins/MedReader/IO/vtkMedFamilyOnEntityOnProfile.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 "vtkMedFamilyOnEntityOnProfile.h"
  21
  22 #include "vtkObjectFactory.h"
  23 #include "vtkMedUtilities.h"
  24 #include "vtkMedFamilyOnEntity.h"
  25 #include "vtkMedProfile.h"
  26 #include "vtkMedEntityArray.h"
  27 #include "vtkMedFamily.h"
  28 #include "vtkMedField.h"
  29 #include "vtkMedGrid.h"
  30 #include "vtkMedMesh.h"
  31 #include "vtkMedFieldOnProfile.h"
  32 #include "vtkMedFieldOverEntity.h"
  33 #include "vtkMedFieldStep.h"
  34 #include "vtkMedFile.h"
  35 #include "vtkMedGrid.h"
  36 #include "vtkMedDriver.h"
  37 #include "vtkMedUnstructuredGrid.h"
  38 #include "vtkMedIntArray.h"
  39
  40 #include "vtkBitArray.h"
  41 #include "vtkIdList.h"
  42
  43 #include "vtkMultiProcessController.h"
  44
  45 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile,FamilyOnEntity, vtkMedFamilyOnEntity);
  46 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile, Profile, vtkMedProfile);
  47
  48 vtkCxxRevisionMacro(vtkMedFamilyOnEntityOnProfile, "$Revision$");
  49 vtkStandardNewMacro(vtkMedFamilyOnEntityOnProfile)
  50
  51 vtkMedFamilyOnEntityOnProfile::vtkMedFamilyOnEntityOnProfile()
  52 {
  53   this->FamilyOnEntity = NULL;
  54   this->Profile = NULL;
  55   this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
  56   this->UseAllPoints = false;
  57   this->MatchComputed = false;
  58   this->Valid = true;
  59 }
  60
  61 vtkMedFamilyOnEntityOnProfile::~vtkMedFamilyOnEntityOnProfile()
  62 {
  63   this->SetFamilyOnEntity(NULL);
  64   this->SetProfile(NULL);
  65 }
  66
  67 bool vtkMedFamilyOnEntityOnProfile::KeepPoint(med_int index)
  68 {
  69   if(this->IntersectionStatus == NotComputed)
  70     this->ComputeIntersection(NULL);
  71
  72   if(this->UseAllPoints)
  73     return true;
  74
  75   if(this->MedToVTKPointIndexMap.find(index)
  76     == this->MedToVTKPointIndexMap.end())
  77     return false;
  78
  79   return true;
  80 }
  81
  82 bool vtkMedFamilyOnEntityOnProfile::KeepCell(med_int index)
  83 {
  84   if(this->FamilyOnEntity->GetEntityArray()->GetFamilyId(index)
  85     != this->FamilyOnEntity->GetFamily()->GetId())
  86     return false;
  87   return true;
  88 }
  89
  90 int vtkMedFamilyOnEntityOnProfile::CanMapField(vtkMedFieldOnProfile* fop)
  91 {
  92   // only point fields can be mapped on point supports.
  93   if(this->GetFamilyOnEntity()->GetEntity().EntityType == MED_NODE &&
  94      fop->GetParentFieldOverEntity()->GetEntity().EntityType != MED_NODE)
  95     return false;
  96
  97   // if it is a cell-centered field, the geometry need to be the same
  98   if(fop->GetParentFieldOverEntity()->GetEntity().EntityType != MED_NODE
  99      && fop->GetParentFieldOverEntity()->GetEntity().GeometryType !=
 100      this->GetFamilyOnEntity()->GetEntity().GeometryType)
 101     return false;
 102
 103   int numProc = 1;
 104   vtkMultiProcessController* controller =
 105                   vtkMultiProcessController::GetGlobalController();
 106   if (controller != NULL)
 107     {
 108     numProc = controller->GetNumberOfProcesses();
 109     }
 110
 111   if ((this->GetValid() == 0) && numProc == 1)
 112     return false;
 113
 114   this->ComputeIntersection(fop);
 115
 116   if(this->IntersectionStatus == vtkMedFamilyOnEntityOnProfile::NoIntersection)
 117     return false;
 118
 119   if(fop != NULL &&
 120      this->GetFamilyOnEntity()->GetEntity().EntityType != MED_NODE &&
 121      fop->GetParentFieldOverEntity()->GetEntity().EntityType == MED_NODE &&
 122      this->PointProfileMatch[fop->GetProfile()] == BadOrNoIntersection)
 123     return false;
 124
 125   return true;
 126 }
 127
 128 int vtkMedFamilyOnEntityOnProfile::CanShallowCopy(vtkMedFieldOnProfile *fop)
 129 {
 130   if(fop == NULL)
 131     {
 132     bool shallow_on_points = this->CanShallowCopyPointField(NULL);
 133     bool shallow_on_cells = this->CanShallowCopyCellField(NULL);
 134     if(shallow_on_points && shallow_on_cells)
 135       return true;
 136     if(!shallow_on_points && !shallow_on_cells)
 137       return false;
 138     vtkErrorMacro("CanShallowCopy cannot answer : is it a point or a cell field?");
 139     return false;
 140     }
 141
 142   if(fop->GetParentFieldOverEntity()->GetParentStep()->GetParentField()
 143     ->GetFieldType() == vtkMedField::PointField)
 144     return this->CanShallowCopyPointField(fop);
 145   else
 146     return this->CanShallowCopyCellField(fop);
 147 }
 148
 149 void vtkMedFamilyOnEntityOnProfile::ComputeIntersection(vtkMedFieldOnProfile* fop)
 150 {
 151   int nodeOrCellSupport=this->GetFamilyOnEntity()->GetPointOrCell();
 152   int fieldType;
 153   if(fop)
 154     {
 155     fieldType = fop->GetParentFieldOverEntity()->GetParentStep()->
 156                   GetParentField()->GetFieldType();
 157     }
 158   else
 159     {
 160     fieldType = (nodeOrCellSupport==vtkMedUtilities::OnPoint?vtkMedField::PointField:vtkMedField::CellField);
 161     }
 162   // Cell fields cannot match point supports
 163   if(fieldType != vtkMedField::PointField
 164      && nodeOrCellSupport == vtkMedUtilities::OnPoint)
 165     {
 166     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
 167     this->UseAllPoints = false;
 168     }
 169   else if(fieldType != vtkMedField::PointField
 170      && nodeOrCellSupport ==vtkMedUtilities::OnCell)
 171     {
 172     this->ComputeCellFamilyVsCellProfileMatch();
 173     }
 174   else if(fieldType == vtkMedField::PointField
 175      && nodeOrCellSupport ==vtkMedUtilities::OnPoint)
 176     {
 177     vtkMedProfile* profile = NULL;
 178     if(fop != NULL)
 179       {
 180       profile = fop->GetProfile();
 181       }
 182     // point fields must share the same profile as the point support.
 183     this->ComputePointFamilyVsPointProfileMatch();
 184
 185     }
 186   else if(fieldType == vtkMedField::PointField
 187      && nodeOrCellSupport == vtkMedUtilities::OnCell)
 188     {
 189     vtkMedProfile* profile = NULL;
 190     if(fop != NULL)
 191       {
 192       profile = fop->GetProfile();
 193       }
 194     this->ComputeCellFamilyVsPointProfileMatch(profile);
 195     }
 196
 197   this->MatchComputed = true;
 198 }
 199
 200 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyPointField(vtkMedFieldOnProfile* fop)
 201 {
 202   vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
 203   if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell)
 204     {
 205     if(this->PointProfileMatch.find(profile) == this->PointProfileMatch.end())
 206       {
 207       this->ComputeCellFamilyVsPointProfileMatch(profile);
 208       }
 209     int match = this->PointProfileMatch[profile];
 210     return match
 211         == vtkMedFamilyOnEntityOnProfile::ProfileEqualsSupport;
 212     }
 213   else
 214     {
 215     // this is a point support.
 216     // The only case when there is shallow copy is if there is at most 1 family
 217     // and the profile is shared.
 218     if(this->Profile == profile &&
 219        this->GetFamilyOnEntity()->GetEntityArray()
 220         ->GetNumberOfFamilyOnEntity() <= 1)
 221       {
 222       return true;
 223       }
 224     else
 225       {
 226       return false;
 227       }
 228     }
 229 }
 230
 231 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyCellField(vtkMedFieldOnProfile* fop)
 232 {
 233   vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
 234   // cell fields cannot be mapped to cell supports
 235   if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnPoint)
 236     {
 237     return false;
 238     }
 239
 240   // this is a cell support.
 241   if(this->Profile != profile)
 242     {
 243     return false;
 244     }
 245
 246   // the only case I can shallow copy is if there is only one family
 247   // defined on those cells.
 248   if(this->Profile == NULL &&
 249      this->GetFamilyOnEntity()->GetEntityArray()
 250     ->GetNumberOfFamilyOnEntity() <= 1)
 251     return true;
 252
 253   return false;
 254 }
 255
 256 void  vtkMedFamilyOnEntityOnProfile::ComputeUsedPoints()
 257 {
 258   this->MedToVTKPointIndexMap.clear();
 259
 260   //first test a few special cases where no heavy computing is necessary
 261   vtkMedGrid* grid = this->FamilyOnEntity->GetParentGrid();
 262   if(this->Profile == NULL)
 263     {
 264     // If there is no profile, the entity is on points and there
 265     // at most 1 point family, then all points are used.
 266     if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnPoint &&
 267        this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() <= 1)
 268       {
 269       this->UseAllPoints = true;
 270       return;
 271       }
 272     // if there is no profile, the entity is on cell
 273     // and there is at most 1 family on his entity, then all points are used
 274     if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell &&
 275        this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() <= 1)
 276       {
 277       this->UseAllPoints = true;
 278       return;
 279       }
 280     }
 281
 282   vtkSmartPointer<vtkBitArray> flag = vtkSmartPointer<vtkBitArray>::New();
 283   flag->SetNumberOfTuples(grid->GetNumberOfPoints());
 284
 285   // initialize the array to false
 286   for(vtkIdType pid = 0; pid < flag->GetNumberOfTuples(); pid++)
 287     flag->SetValue(pid, false);
 288
 289   // for each cell, flag the used points
 290   if(this->Profile)
 291     this->Profile->Load();
 292
 293   vtkMedIntArray* pids=(this->Profile!=NULL?this->Profile->GetIds():NULL);
 294
 295   med_int famId = this->FamilyOnEntity->GetFamily()->GetId();
 296   vtkMedEntityArray* array = this->FamilyOnEntity->GetEntityArray();
 297   vtkSmartPointer<vtkIdList> ids = vtkSmartPointer<vtkIdList>::New();
 298
 299   array->LoadConnectivity();
 300
 301   vtkIdType pflsize = (pids != NULL ? pids->GetNumberOfTuples():array->GetNumberOfEntity());
 302   for(vtkIdType pindex=0; pindex<pflsize; pindex++)
 303     {
 304     med_int pid = (pids != NULL ? pids->GetValue(pindex)-1 : pindex);
 305     med_int fid = array->GetFamilyId(pid);
 306     if(famId==fid)
 307       {
 308       // this cell is of the family and on the profile.
 309       // --> flag all vertices of this cell
 310       array->GetCellVertices(pid, ids);
 311       for(int id = 0; id<ids->GetNumberOfIds(); id++)
 312         {
 313         vtkIdType subid = ids->GetId(id);
 314         if(subid < 0 || subid >= flag->GetNumberOfTuples())
 315           {
 316           vtkDebugMacro("invalid sub id : " << subid);
 317           this->SetValid(0);
 318           break;
 319           }
 320         flag->SetValue(subid, 1);
 321         }
 322       }
 323     }
 324
 325   // now, the flag array contains all vertices used by this support
 326   this->UseAllPoints = true;
 327   for(vtkIdType pid = 0; pid<flag->GetNumberOfTuples(); pid++)
 328     {
 329     if(flag->GetValue(pid) == false)
 330       {
 331       this->UseAllPoints = false;
 332       break;
 333       }
 334     }
 335
 336   if(!this->UseAllPoints)
 337     {
 338     // If all points are not used, I compute the index mapping
 339     vtkIdType vtk_index = 0;
 340
 341     for(vtkIdType pid=0; pid < flag->GetNumberOfTuples(); pid++)
 342       {
 343       if(flag->GetValue(pid) == true)
 344         {
 345         this->MedToVTKPointIndexMap[pid] = vtk_index;
 346         vtk_index++;
 347         }
 348       }
 349     }
 350 }
 351
 352 void vtkMedFamilyOnEntityOnProfile::ComputeCellFamilyVsCellProfileMatch()
 353 {
 354   if(this->MatchComputed)
 355     return;
 356
 357   // this computes the UseAllPoints flag.
 358   this->ComputeUsedPoints();
 359
 360   if(this->Profile == NULL)
 361     {
 362     // If there is no profile, then the match is exact if and only
 363     // if there is 1 cell family on this entity
 364     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() == 1)
 365       {
 366       this->IntersectionStatus =
 367           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 368       }
 369     else
 370       {
 371       this->IntersectionStatus =
 372           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 373       }
 374     return;
 375     }
 376
 377   this->Profile->Load();
 378   vtkMedIntArray* pids=this->Profile->GetIds();
 379
 380   if(pids==NULL)
 381     {
 382     vtkErrorMacro("Could not load profile indices!");
 383     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
 384     this->UseAllPoints = false;
 385     return;
 386     }
 387
 388   med_int famId = this->GetFamilyOnEntity()->GetFamily()->GetId();
 389   vtkIdType pindex=0;
 390   vtkMedEntityArray* array=this->FamilyOnEntity->GetEntityArray();
 391   bool profile_included = true;
 392   bool profile_intersect = false;
 393   for(int pindex=0; pindex<pids->GetNumberOfTuples(); pindex++)
 394     {
 395     med_int pid=pids->GetValue(pindex)-1;
 396     med_int fid=array->GetFamilyId(pid);
 397     if(famId==fid)
 398       {// the current cell is on the familyand on the profile
 399       // --> there is an overlap
 400       profile_intersect = true;
 401       }
 402     else
 403       {
 404       // the cell is on the profile but not on the family --> no inclusion
 405       profile_included=false;
 406       }
 407     }
 408
 409   if(profile_included && profile_intersect)
 410     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 411   else if(profile_intersect)
 412     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 413   else
 414     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
 415 }
 416
 417 void vtkMedFamilyOnEntityOnProfile::
 418     ComputeCellFamilyVsPointProfileMatch(vtkMedProfile* profile)
 419 {
 420   // first test if the cache is already set
 421   if(this->PointProfileMatch.find(profile) != this->PointProfileMatch.end())
 422     return;
 423
 424   // this will compute the cell match cache, as well as the UseAllPoints flag.
 425   this->ComputeCellFamilyVsCellProfileMatch();
 426
 427   if(profile == NULL)
 428     {
 429     // If there is no profile, then the match is at least partial.
 430     // It is exact if and only
 431     // if the cell family uses all points
 432     int match =  (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
 433     this->PointProfileMatch[profile] =
 434         (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
 435     return;
 436     }
 437
 438   // if profile is not NULL and I use all points --> BadOrNoIntersection
 439   if(this->UseAllPoints)
 440     {
 441     this->PointProfileMatch[profile] = BadOrNoIntersection;
 442     }
 443
 444   profile->Load();
 445
 446   vtkMedIntArray* pids=profile->GetIds();
 447
 448   if(pids == NULL)
 449     {
 450     vtkErrorMacro("profile indices could not be loaded!");
 451     this->PointProfileMatch[profile] = BadOrNoIntersection;
 452     return;
 453     }
 454
 455   med_int pindex = 0;
 456   bool exact_match = true;
 457   vtkIdType numberOfUsedPoints = pids->GetNumberOfTuples();
 458   for(med_int pindex=0; pindex < pids->GetNumberOfTuples(); pindex++)
 459     {
 460     med_int id = pids->GetValue(pindex);
 461     if(this->MedToVTKPointIndexMap.find(id-1) == this->MedToVTKPointIndexMap.end())
 462       {
 463       // The given point profile index is not used by this support.
 464       // the superposition is at most partial.
 465       exact_match = false;
 466       numberOfUsedPoints--;
 467       }
 468     }
 469
 470   // if this profile is smaller than the number of points, I can't match
 471   // the profile to this support
 472   if(numberOfUsedPoints < this->MedToVTKPointIndexMap.size())
 473     {
 474     this->PointProfileMatch[profile] = BadOrNoIntersection;
 475     }
 476   else if(exact_match)
 477     {
 478     this->PointProfileMatch[profile] = ProfileEqualsSupport;
 479     }
 480   else
 481     {
 482     this->PointProfileMatch[profile] = ProfileLargerThanSupport;
 483     }
 484 }
 485
 486 void  vtkMedFamilyOnEntityOnProfile::ComputePointFamilyVsPointProfileMatch()
 487 {
 488   if(this->MatchComputed)
 489     return;
 490
 491   this->ComputeUsedPoints();
 492
 493   if(this->Profile == NULL)
 494     {
 495     // If there is no profile, then the match is exact if there is at most
 496     // 1 point family on the grid
 497     if(this->FamilyOnEntity->GetParentGrid()->GetParentMesh()
 498       ->GetNumberOfPointFamily() <= 1)
 499       {
 500       this->IntersectionStatus = ProfileIncludedInFamily;
 501       }
 502     }
 503
 504   // there is a profile, we have to compute the match between the family and
 505   // the profile
 506   vtkMedFamilyOnEntity* foe = this->GetFamilyOnEntity();
 507   vtkMedEntityArray* pea = foe->GetEntityArray();
 508   vtkMedIntArray* pIds = NULL;
 509
 510   if(this->Profile)
 511     {
 512     this->Profile->Load();
 513     pIds=this->Profile->GetIds();
 514     }
 515
 516   if(pIds == NULL)
 517     {
 518     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() > 1)
 519       {
 520       this->IntersectionStatus =
 521           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 522       }
 523     else
 524       {
 525       this->IntersectionStatus =
 526           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 527       }
 528     return;
 529     }
 530
 531   bool profile_intersects=false;
 532   bool profile_included=true;
 533   med_int famId=this->FamilyOnEntity->GetFamily()->GetId();
 534   for(vtkIdType pindex=0; pindex<pIds->GetNumberOfTuples(); pindex++)
 535     {
 536     med_int pid=pIds->GetValue(pindex)-1;
 537     med_int fid = pea->GetFamilyId(pid);
 538    // med_int fid=grid->GetPointFamilyId(pid);
 539     if(fid==famId)
 540       {// the family of the current point is the good one
 541       profile_intersects=true;
 542       }
 543     else
 544       {
 545       // we are on the profile and not on the family -->
 546       // no exact match, but the the profile might be larger than the family.
 547       profile_included=false;
 548       }
 549     }
 550
 551   if(!profile_intersects)
 552     {
 553     this->IntersectionStatus =
 554         vtkMedFamilyOnEntityOnProfile::NoIntersection;
 555     }
 556   else if(profile_included)
 557     {
 558     this->IntersectionStatus =
 559         vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 560     }
 561   else
 562     {
 563     this->IntersectionStatus =
 564         vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 565     }
 566 }
 567
 568 vtkIdType vtkMedFamilyOnEntityOnProfile::GetVTKPointIndex(vtkIdType medCIndex)
 569 {
 570   if(this->IntersectionStatus == NotComputed)
 571     this->ComputeIntersection(NULL);
 572
 573   if(this->UseAllPoints)
 574     return medCIndex;
 575
 576   if(this->MedToVTKPointIndexMap.find(medCIndex)
 577     == this->MedToVTKPointIndexMap.end())
 578     {
 579     vtkDebugMacro("GetVTKPointIndex asked for "
 580                   << medCIndex << " which has not been mapped");
 581     return -1;
 582     }
 583
 584   return this->MedToVTKPointIndexMap[medCIndex];
 585 }
 586
 587 void vtkMedFamilyOnEntityOnProfile::PrintSelf(ostream& os, vtkIndent indent)
 588 {
 589   this->Superclass::PrintSelf(os, indent);
 590 }