src/Plugins/MedReader/IO/vtkMedFamilyOnEntityOnProfile.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 "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 grid is structured, the entity is on cell
 273     // and there is at most 1 family on his entity, then all points are used
 274     if(vtkMedUnstructuredGrid::SafeDownCast(grid) == NULL &&
 275        this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell &&
 276        this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() <= 1)
 277       {
 278       this->UseAllPoints = true;
 279       return;
 280       }
 281     }
 282
 283   vtkSmartPointer<vtkBitArray> flag = vtkSmartPointer<vtkBitArray>::New();
 284   flag->SetNumberOfTuples(grid->GetNumberOfPoints());
 285
 286   // initialize the array to false
 287   for(vtkIdType pid = 0; pid < flag->GetNumberOfTuples(); pid++)
 288     flag->SetValue(pid, false);
 289
 290   // for each cell, flag the used points
 291   if(this->Profile)
 292     this->Profile->Load();
 293
 294   vtkMedIntArray* pids=(this->Profile!=NULL?this->Profile->GetIds():NULL);
 295
 296   med_int famId = this->FamilyOnEntity->GetFamily()->GetId();
 297   vtkMedEntityArray* array = this->FamilyOnEntity->GetEntityArray();
 298   vtkSmartPointer<vtkIdList> ids = vtkSmartPointer<vtkIdList>::New();
 299
 300   array->LoadConnectivity();
 301
 302   vtkIdType pflsize = (pids != NULL ? pids->GetNumberOfTuples():array->GetNumberOfEntity());
 303   for(vtkIdType pindex=0; pindex<pflsize; pindex++)
 304     {
 305     med_int pid = (pids != NULL ? pids->GetValue(pindex)-1 : pindex);
 306     med_int fid = array->GetFamilyId(pid);
 307     if(famId==fid)
 308       {
 309       // this cell is of the family and on the profile.
 310       // --> flag all vertices of this cell
 311       array->GetCellVertices(pid, ids);
 312       for(int id = 0; id<ids->GetNumberOfIds(); id++)
 313         {
 314         vtkIdType subid = ids->GetId(id);
 315         if(subid < 0 || subid >= flag->GetNumberOfTuples())
 316           {
 317           vtkDebugMacro("invalid sub id : " << subid);
 318           this->SetValid(0);
 319           break;
 320           }
 321         flag->SetValue(subid, 1);
 322         }
 323       }
 324     }
 325
 326   // now, the flag array contains all vertices used by this support
 327   this->UseAllPoints = true;
 328   for(vtkIdType pid = 0; pid<flag->GetNumberOfTuples(); pid++)
 329     {
 330     if(flag->GetValue(pid) == false)
 331       {
 332       this->UseAllPoints = false;
 333       break;
 334       }
 335     }
 336
 337   if(!this->UseAllPoints)
 338     {
 339     // If all points are not used, I compute the index mapping
 340     vtkIdType vtk_index = 0;
 341
 342     for(vtkIdType pid=0; pid < flag->GetNumberOfTuples(); pid++)
 343       {
 344       if(flag->GetValue(pid) == true)
 345         {
 346         this->MedToVTKPointIndexMap[pid] = vtk_index;
 347         vtk_index++;
 348         }
 349       }
 350     }
 351 }
 352
 353 void vtkMedFamilyOnEntityOnProfile::ComputeCellFamilyVsCellProfileMatch()
 354 {
 355   if(this->MatchComputed)
 356     return;
 357
 358   // this computes the UseAllPoints flag.
 359   this->ComputeUsedPoints();
 360
 361   if(this->Profile == NULL)
 362     {
 363     // If there is no profile, then the match is exact if and only
 364     // if there is 1 cell family on this entity
 365     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() == 1)
 366       {
 367       this->IntersectionStatus =
 368           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 369       }
 370     else
 371       {
 372       this->IntersectionStatus =
 373           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 374       }
 375     return;
 376     }
 377
 378   this->Profile->Load();
 379   vtkMedIntArray* pids=this->Profile->GetIds();
 380
 381   if(pids==NULL)
 382     {
 383     vtkErrorMacro("Could not load profile indices!");
 384     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
 385     this->UseAllPoints = false;
 386     return;
 387     }
 388
 389   med_int famId = this->GetFamilyOnEntity()->GetFamily()->GetId();
 390   vtkIdType pindex=0;
 391   vtkMedEntityArray* array=this->FamilyOnEntity->GetEntityArray();
 392   bool profile_included = true;
 393   bool profile_intersect = false;
 394   for(int pindex=0; pindex<pids->GetNumberOfTuples(); pindex++)
 395     {
 396     med_int pid=pids->GetValue(pindex)-1;
 397     med_int fid=array->GetFamilyId(pid);
 398     if(famId==fid)
 399       {// the current cell is on the familyand on the profile
 400       // --> there is an overlap
 401       profile_intersect = true;
 402       }
 403     else
 404       {
 405       // the cell is on the profile but not on the family --> no inclusion
 406       profile_included=false;
 407       }
 408     }
 409
 410   if(profile_included && profile_intersect)
 411     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 412   else if(profile_intersect)
 413     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 414   else
 415     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
 416 }
 417
 418 void vtkMedFamilyOnEntityOnProfile::
 419     ComputeCellFamilyVsPointProfileMatch(vtkMedProfile* profile)
 420 {
 421   // first test if the cache is already set
 422   if(this->PointProfileMatch.find(profile) != this->PointProfileMatch.end())
 423     return;
 424
 425   // this will compute the cell match cache, as well as the UseAllPoints flag.
 426   this->ComputeCellFamilyVsCellProfileMatch();
 427
 428   if(profile == NULL)
 429     {
 430     // If there is no profile, then the match is at least partial.
 431     // It is exact if and only
 432     // if the cell family uses all points
 433     int match =  (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
 434     this->PointProfileMatch[profile] =
 435         (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
 436     return;
 437     }
 438
 439   // if profile is not NULL and I use all points --> BadOrNoIntersection
 440   if(this->UseAllPoints)
 441     {
 442     this->PointProfileMatch[profile] = BadOrNoIntersection;
 443     }
 444
 445   profile->Load();
 446
 447   vtkMedIntArray* pids=profile->GetIds();
 448
 449   if(pids == NULL)
 450     {
 451     vtkErrorMacro("profile indices could not be loaded!");
 452     this->PointProfileMatch[profile] = BadOrNoIntersection;
 453     return;
 454     }
 455
 456   med_int pindex = 0;
 457   bool exact_match = true;
 458   vtkIdType numberOfUsedPoints = pids->GetNumberOfTuples();
 459   for(med_int pindex=0; pindex < pids->GetNumberOfTuples(); pindex++)
 460     {
 461     med_int id = pids->GetValue(pindex);
 462     if(this->MedToVTKPointIndexMap.find(id-1) == this->MedToVTKPointIndexMap.end())
 463       {
 464       // The given point profile index is not used by this support.
 465       // the superposition is at most partial.
 466       exact_match = false;
 467       numberOfUsedPoints--;
 468       }
 469     }
 470
 471   // if this profile is smaller than the number of points, I can't match
 472   // the profile to this support
 473   if(numberOfUsedPoints < this->MedToVTKPointIndexMap.size())
 474     {
 475     this->PointProfileMatch[profile] = BadOrNoIntersection;
 476     }
 477   else if(exact_match)
 478     {
 479     this->PointProfileMatch[profile] = ProfileEqualsSupport;
 480     }
 481   else
 482     {
 483     this->PointProfileMatch[profile] = ProfileLargerThanSupport;
 484     }
 485 }
 486
 487 void  vtkMedFamilyOnEntityOnProfile::ComputePointFamilyVsPointProfileMatch()
 488 {
 489   if(this->MatchComputed)
 490     return;
 491
 492   this->ComputeUsedPoints();
 493
 494   if(this->Profile == NULL)
 495     {
 496     // If there is no profile, then the match is exact if there is at most
 497     // 1 point family on the grid
 498     if(this->FamilyOnEntity->GetParentGrid()->GetParentMesh()
 499       ->GetNumberOfPointFamily() <= 1)
 500       {
 501       this->IntersectionStatus = ProfileIncludedInFamily;
 502       }
 503     }
 504
 505   // there is a profile, we have to compute the match between the family and
 506   // the profile
 507   vtkMedFamilyOnEntity* foe = this->GetFamilyOnEntity();
 508   vtkMedEntityArray* pea = foe->GetEntityArray();
 509   vtkMedIntArray* pIds = NULL;
 510
 511   if(this->Profile)
 512     {
 513     this->Profile->Load();
 514     pIds=this->Profile->GetIds();
 515     }
 516
 517   if(pIds == NULL)
 518     {
 519     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() > 1)
 520       {
 521       this->IntersectionStatus =
 522           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 523       }
 524     else
 525       {
 526       this->IntersectionStatus =
 527           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 528       }
 529     return;
 530     }
 531
 532   bool profile_intersects=false;
 533   bool profile_included=true;
 534   med_int famId=this->FamilyOnEntity->GetFamily()->GetId();
 535   for(vtkIdType pindex=0; pindex<pIds->GetNumberOfTuples(); pindex++)
 536     {
 537     med_int pid=pIds->GetValue(pindex)-1;
 538     med_int fid = pea->GetFamilyId(pid);
 539    // med_int fid=grid->GetPointFamilyId(pid);
 540     if(fid==famId)
 541       {// the family of the current point is the good one
 542       profile_intersects=true;
 543       }
 544     else
 545       {
 546       // we are on the profile and not on the family -->
 547       // no exact match, but the the profile might be larger than the family.
 548       profile_included=false;
 549       }
 550     }
 551
 552   if(!profile_intersects)
 553     {
 554     this->IntersectionStatus =
 555         vtkMedFamilyOnEntityOnProfile::NoIntersection;
 556     }
 557   else if(profile_included)
 558     {
 559     this->IntersectionStatus =
 560         vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
 561     }
 562   else
 563     {
 564     this->IntersectionStatus =
 565         vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
 566     }
 567 }
 568
 569 vtkIdType vtkMedFamilyOnEntityOnProfile::GetVTKPointIndex(vtkIdType medCIndex)
 570 {
 571   if(this->IntersectionStatus == NotComputed)
 572     this->ComputeIntersection(NULL);
 573
 574   if(this->UseAllPoints)
 575     return medCIndex;
 576
 577   if(this->MedToVTKPointIndexMap.find(medCIndex)
 578     == this->MedToVTKPointIndexMap.end())
 579     {
 580     vtkDebugMacro("GetVTKPointIndex asked for "
 581                   << medCIndex << " which has not been mapped");
 582     return -1;
 583     }
 584
 585   return this->MedToVTKPointIndexMap[medCIndex];
 586 }
 587
 588 void vtkMedFamilyOnEntityOnProfile::PrintSelf(ostream& os, vtkIndent indent)
 589 {
 590   this->Superclass::PrintSelf(os, indent);
 591 }