Back to index

salome-paravis  6.5.0
vtkMedFamilyOnEntityOnProfile.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2010-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "vtkMedFamilyOnEntityOnProfile.h"
00021 
00022 #include "vtkObjectFactory.h"
00023 #include "vtkMedUtilities.h"
00024 #include "vtkMedFamilyOnEntity.h"
00025 #include "vtkMedProfile.h"
00026 #include "vtkMedEntityArray.h"
00027 #include "vtkMedFamily.h"
00028 #include "vtkMedField.h"
00029 #include "vtkMedGrid.h"
00030 #include "vtkMedMesh.h"
00031 #include "vtkMedFieldOnProfile.h"
00032 #include "vtkMedFieldOverEntity.h"
00033 #include "vtkMedFieldStep.h"
00034 #include "vtkMedFile.h"
00035 #include "vtkMedGrid.h"
00036 #include "vtkMedDriver.h"
00037 #include "vtkMedUnstructuredGrid.h"
00038 #include "vtkMedIntArray.h"
00039 
00040 #include "vtkBitArray.h"
00041 #include "vtkIdList.h"
00042 
00043 #include "vtkMultiProcessController.h"
00044 
00045 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile,FamilyOnEntity, vtkMedFamilyOnEntity);
00046 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile, Profile, vtkMedProfile);
00047 
00048 vtkCxxRevisionMacro(vtkMedFamilyOnEntityOnProfile, "$Revision: 1.1.4.7.2.1 $");
00049 vtkStandardNewMacro(vtkMedFamilyOnEntityOnProfile)
00050 
00051 vtkMedFamilyOnEntityOnProfile::vtkMedFamilyOnEntityOnProfile()
00052 {
00053   this->FamilyOnEntity = NULL;
00054   this->Profile = NULL;
00055   this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
00056   this->UseAllPoints = false;
00057   this->MatchComputed = false;
00058   this->Valid = true;
00059 }
00060 
00061 vtkMedFamilyOnEntityOnProfile::~vtkMedFamilyOnEntityOnProfile()
00062 {
00063   this->SetFamilyOnEntity(NULL);
00064   this->SetProfile(NULL);
00065 }
00066 
00067 bool vtkMedFamilyOnEntityOnProfile::KeepPoint(med_int index)
00068 {
00069   if(this->IntersectionStatus == NotComputed)
00070     this->ComputeIntersection(NULL);
00071 
00072   if(this->UseAllPoints)
00073     return true;
00074 
00075   if(this->MedToVTKPointIndexMap.find(index)
00076     == this->MedToVTKPointIndexMap.end())
00077     return false;
00078 
00079   return true;
00080 }
00081 
00082 bool vtkMedFamilyOnEntityOnProfile::KeepCell(med_int index)
00083 {
00084   if(this->FamilyOnEntity->GetEntityArray()->GetFamilyId(index)
00085     != this->FamilyOnEntity->GetFamily()->GetId())
00086     return false;
00087   return true;
00088 }
00089 
00090 int vtkMedFamilyOnEntityOnProfile::CanMapField(vtkMedFieldOnProfile* fop)
00091 {
00092   // only point fields can be mapped on point supports.
00093   if(this->GetFamilyOnEntity()->GetEntity().EntityType == MED_NODE &&
00094      fop->GetParentFieldOverEntity()->GetEntity().EntityType != MED_NODE)
00095     return false;
00096 
00097   // if it is a cell-centered field, the geometry need to be the same
00098   if(fop->GetParentFieldOverEntity()->GetEntity().EntityType != MED_NODE
00099      && fop->GetParentFieldOverEntity()->GetEntity().GeometryType !=
00100      this->GetFamilyOnEntity()->GetEntity().GeometryType)
00101     return false;
00102 
00103   int numProc = 1;
00104   vtkMultiProcessController* controller =
00105                 vtkMultiProcessController::GetGlobalController();
00106   if (controller != NULL)
00107     {
00108     numProc = controller->GetNumberOfProcesses();
00109     }
00110 
00111   if ((this->GetValid() == 0) && numProc == 1)
00112     return false;
00113 
00114   this->ComputeIntersection(fop);
00115 
00116   if(this->IntersectionStatus == vtkMedFamilyOnEntityOnProfile::NoIntersection)
00117     return false;
00118 
00119   if(fop != NULL &&
00120      this->GetFamilyOnEntity()->GetEntity().EntityType != MED_NODE &&
00121      fop->GetParentFieldOverEntity()->GetEntity().EntityType == MED_NODE &&
00122      this->PointProfileMatch[fop->GetProfile()] == BadOrNoIntersection)
00123     return false;
00124 
00125   return true;
00126 }
00127 
00128 int vtkMedFamilyOnEntityOnProfile::CanShallowCopy(vtkMedFieldOnProfile *fop)
00129 {
00130   if(fop == NULL)
00131     {
00132     bool shallow_on_points = this->CanShallowCopyPointField(NULL);
00133     bool shallow_on_cells = this->CanShallowCopyCellField(NULL);
00134     if(shallow_on_points && shallow_on_cells)
00135       return true;
00136     if(!shallow_on_points && !shallow_on_cells)
00137       return false;
00138     vtkErrorMacro("CanShallowCopy cannot answer : is it a point or a cell field?");
00139     return false;
00140     }
00141 
00142   if(fop->GetParentFieldOverEntity()->GetParentStep()->GetParentField()
00143     ->GetFieldType() == vtkMedField::PointField)
00144     return this->CanShallowCopyPointField(fop);
00145   else
00146     return this->CanShallowCopyCellField(fop);
00147 }
00148 
00149 void vtkMedFamilyOnEntityOnProfile::ComputeIntersection(vtkMedFieldOnProfile* fop)
00150 {
00151   int nodeOrCellSupport=this->GetFamilyOnEntity()->GetPointOrCell();
00152   int fieldType;
00153   if(fop)
00154     {
00155     fieldType = fop->GetParentFieldOverEntity()->GetParentStep()->
00156                   GetParentField()->GetFieldType();
00157     }
00158   else
00159     {
00160     fieldType = (nodeOrCellSupport==vtkMedUtilities::OnPoint?vtkMedField::PointField:vtkMedField::CellField);
00161     }
00162   // Cell fields cannot match point supports
00163   if(fieldType != vtkMedField::PointField
00164      && nodeOrCellSupport == vtkMedUtilities::OnPoint)
00165     {
00166     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
00167     this->UseAllPoints = false;
00168     }
00169   else if(fieldType != vtkMedField::PointField
00170      && nodeOrCellSupport ==vtkMedUtilities::OnCell)
00171     {
00172     this->ComputeCellFamilyVsCellProfileMatch();
00173     }
00174   else if(fieldType == vtkMedField::PointField
00175      && nodeOrCellSupport ==vtkMedUtilities::OnPoint)
00176     {
00177     vtkMedProfile* profile = NULL;
00178     if(fop != NULL)
00179       {
00180       profile = fop->GetProfile();
00181       }
00182     // point fields must share the same profile as the point support.
00183     this->ComputePointFamilyVsPointProfileMatch();
00184 
00185     }
00186   else if(fieldType == vtkMedField::PointField
00187      && nodeOrCellSupport == vtkMedUtilities::OnCell)
00188     {
00189     vtkMedProfile* profile = NULL;
00190     if(fop != NULL)
00191       {
00192       profile = fop->GetProfile();
00193       }
00194     this->ComputeCellFamilyVsPointProfileMatch(profile);
00195     }
00196 
00197   this->MatchComputed = true;
00198 }
00199 
00200 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyPointField(vtkMedFieldOnProfile* fop)
00201 {
00202   vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
00203   if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell)
00204     {
00205     if(this->PointProfileMatch.find(profile) == this->PointProfileMatch.end())
00206       {
00207       this->ComputeCellFamilyVsPointProfileMatch(profile);
00208       }
00209     int match = this->PointProfileMatch[profile];
00210     return match
00211         == vtkMedFamilyOnEntityOnProfile::ProfileEqualsSupport;
00212     }
00213   else
00214     {
00215     // this is a point support.
00216     // The only case when there is shallow copy is if there is at most 1 family
00217     // and the profile is shared.
00218     if(this->Profile == profile &&
00219        this->GetFamilyOnEntity()->GetEntityArray()
00220         ->GetNumberOfFamilyOnEntity() <= 1)
00221       {
00222       return true;
00223       }
00224     else
00225       {
00226       return false;
00227       }
00228     }
00229 }
00230 
00231 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyCellField(vtkMedFieldOnProfile* fop)
00232 {
00233   vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
00234   // cell fields cannot be mapped to cell supports
00235   if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnPoint)
00236     {
00237     return false;
00238     }
00239 
00240   // this is a cell support.
00241   if(this->Profile != profile)
00242     {
00243     return false;
00244     }
00245 
00246   // the only case I can shallow copy is if there is only one family
00247   // defined on those cells.
00248   if(this->Profile == NULL &&
00249      this->GetFamilyOnEntity()->GetEntityArray()
00250     ->GetNumberOfFamilyOnEntity() <= 1)
00251     return true;
00252 
00253   return false;
00254 }
00255 
00256 void  vtkMedFamilyOnEntityOnProfile::ComputeUsedPoints()
00257 {
00258   this->MedToVTKPointIndexMap.clear();
00259 
00260   //first test a few special cases where no heavy computing is necessary
00261   vtkMedGrid* grid = this->FamilyOnEntity->GetParentGrid();
00262   if(this->Profile == NULL)
00263     {
00264     // If there is no profile, the entity is on points and there
00265     // at most 1 point family, then all points are used.
00266     if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnPoint &&
00267        this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() <= 1)
00268       {
00269       this->UseAllPoints = true;
00270       return;
00271       }
00272     // if there is no profile, the grid is structured, the entity is on cell
00273     // and there is at most 1 family on his entity, then all points are used
00274     if(vtkMedUnstructuredGrid::SafeDownCast(grid) == NULL &&
00275        this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell &&
00276        this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() <= 1)
00277       {
00278       this->UseAllPoints = true;
00279       return;
00280       }
00281     }
00282 
00283   vtkSmartPointer<vtkBitArray> flag = vtkSmartPointer<vtkBitArray>::New();
00284   flag->SetNumberOfTuples(grid->GetNumberOfPoints());
00285 
00286   // initialize the array to false
00287   for(vtkIdType pid = 0; pid < flag->GetNumberOfTuples(); pid++)
00288     flag->SetValue(pid, false);
00289 
00290   // for each cell, flag the used points
00291   if(this->Profile)
00292     this->Profile->Load();
00293 
00294   vtkMedIntArray* pids=(this->Profile!=NULL?this->Profile->GetIds():NULL);
00295 
00296   med_int famId = this->FamilyOnEntity->GetFamily()->GetId();
00297   vtkMedEntityArray* array = this->FamilyOnEntity->GetEntityArray();
00298   vtkSmartPointer<vtkIdList> ids = vtkSmartPointer<vtkIdList>::New();
00299 
00300   array->LoadConnectivity();
00301 
00302   vtkIdType pflsize = (pids != NULL ? pids->GetNumberOfTuples():array->GetNumberOfEntity());
00303   for(vtkIdType pindex=0; pindex<pflsize; pindex++)
00304     {
00305     med_int pid = (pids != NULL ? pids->GetValue(pindex)-1 : pindex);
00306     med_int fid = array->GetFamilyId(pid);
00307     if(famId==fid)
00308       {
00309       // this cell is of the family and on the profile.
00310       // --> flag all vertices of this cell
00311       array->GetCellVertices(pid, ids);
00312       for(int id = 0; id<ids->GetNumberOfIds(); id++)
00313         {
00314         vtkIdType subid = ids->GetId(id);
00315         if(subid < 0 || subid >= flag->GetNumberOfTuples())
00316           {
00317           vtkDebugMacro("invalid sub id : " << subid);
00318           this->SetValid(0);
00319           break;
00320           }
00321         flag->SetValue(subid, 1);
00322         }
00323       }
00324     }
00325 
00326   // now, the flag array contains all vertices used by this support
00327   this->UseAllPoints = true;
00328   for(vtkIdType pid = 0; pid<flag->GetNumberOfTuples(); pid++)
00329     {
00330     if(flag->GetValue(pid) == false)
00331       {
00332       this->UseAllPoints = false;
00333       break;
00334       }
00335     }
00336 
00337   if(!this->UseAllPoints)
00338     {
00339     // If all points are not used, I compute the index mapping
00340     vtkIdType vtk_index = 0;
00341 
00342     for(vtkIdType pid=0; pid < flag->GetNumberOfTuples(); pid++)
00343       {
00344       if(flag->GetValue(pid) == true)
00345         {
00346         this->MedToVTKPointIndexMap[pid] = vtk_index;
00347         vtk_index++;
00348         }
00349       }
00350     }
00351 }
00352 
00353 void vtkMedFamilyOnEntityOnProfile::ComputeCellFamilyVsCellProfileMatch()
00354 {
00355   if(this->MatchComputed)
00356     return;
00357 
00358   // this computes the UseAllPoints flag.
00359   this->ComputeUsedPoints();
00360 
00361   if(this->Profile == NULL)
00362     {
00363     // If there is no profile, then the match is exact if and only
00364     // if there is 1 cell family on this entity
00365     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() == 1)
00366       {
00367       this->IntersectionStatus =
00368           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
00369       }
00370     else
00371       {
00372       this->IntersectionStatus =
00373           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
00374       }
00375     return;
00376     }
00377 
00378   this->Profile->Load();
00379   vtkMedIntArray* pids=this->Profile->GetIds();
00380 
00381   if(pids==NULL)
00382     {
00383     vtkErrorMacro("Could not load profile indices!");
00384     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
00385     this->UseAllPoints = false;
00386     return;
00387     }
00388 
00389   med_int famId = this->GetFamilyOnEntity()->GetFamily()->GetId();
00390   vtkIdType pindex=0;
00391   vtkMedEntityArray* array=this->FamilyOnEntity->GetEntityArray();
00392   bool profile_included = true;
00393   bool profile_intersect = false;
00394   for(int pindex=0; pindex<pids->GetNumberOfTuples(); pindex++)
00395     {
00396     med_int pid=pids->GetValue(pindex)-1;
00397     med_int fid=array->GetFamilyId(pid);
00398     if(famId==fid)
00399       {// the current cell is on the familyand on the profile
00400       // --> there is an overlap
00401       profile_intersect = true;
00402       }
00403     else
00404       {
00405       // the cell is on the profile but not on the family --> no inclusion
00406       profile_included=false;
00407       }
00408     }
00409 
00410   if(profile_included && profile_intersect)
00411     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
00412   else if(profile_intersect)
00413     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
00414   else
00415     this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
00416 }
00417 
00418 void vtkMedFamilyOnEntityOnProfile::
00419     ComputeCellFamilyVsPointProfileMatch(vtkMedProfile* profile)
00420 {
00421   // first test if the cache is already set
00422   if(this->PointProfileMatch.find(profile) != this->PointProfileMatch.end())
00423     return;
00424 
00425   // this will compute the cell match cache, as well as the UseAllPoints flag.
00426   this->ComputeCellFamilyVsCellProfileMatch();
00427 
00428   if(profile == NULL)
00429     {
00430     // If there is no profile, then the match is at least partial.
00431     // It is exact if and only
00432     // if the cell family uses all points
00433     int match =  (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
00434     this->PointProfileMatch[profile] =
00435         (this->UseAllPoints? ProfileEqualsSupport : ProfileLargerThanSupport);
00436     return;
00437     }
00438 
00439   // if profile is not NULL and I use all points --> BadOrNoIntersection
00440   if(this->UseAllPoints)
00441     {
00442     this->PointProfileMatch[profile] = BadOrNoIntersection;
00443     }
00444 
00445   profile->Load();
00446 
00447   vtkMedIntArray* pids=profile->GetIds();
00448 
00449   if(pids == NULL)
00450     {
00451     vtkErrorMacro("profile indices could not be loaded!");
00452     this->PointProfileMatch[profile] = BadOrNoIntersection;
00453     return;
00454     }
00455 
00456   med_int pindex = 0;
00457   bool exact_match = true;
00458   vtkIdType numberOfUsedPoints = pids->GetNumberOfTuples();
00459   for(med_int pindex=0; pindex < pids->GetNumberOfTuples(); pindex++)
00460     {
00461     med_int id = pids->GetValue(pindex);
00462     if(this->MedToVTKPointIndexMap.find(id-1) == this->MedToVTKPointIndexMap.end())
00463       {
00464       // The given point profile index is not used by this support.
00465       // the superposition is at most partial.
00466       exact_match = false;
00467       numberOfUsedPoints--;
00468       }
00469     }
00470 
00471   // if this profile is smaller than the number of points, I can't match
00472   // the profile to this support
00473   if(numberOfUsedPoints < this->MedToVTKPointIndexMap.size())
00474     {
00475     this->PointProfileMatch[profile] = BadOrNoIntersection;
00476     }
00477   else if(exact_match)
00478     {
00479     this->PointProfileMatch[profile] = ProfileEqualsSupport;
00480     }
00481   else
00482     {
00483     this->PointProfileMatch[profile] = ProfileLargerThanSupport;
00484     }
00485 }
00486 
00487 void  vtkMedFamilyOnEntityOnProfile::ComputePointFamilyVsPointProfileMatch()
00488 {
00489   if(this->MatchComputed)
00490     return;
00491 
00492   this->ComputeUsedPoints();
00493 
00494   if(this->Profile == NULL)
00495     {
00496     // If there is no profile, then the match is exact if there is at most
00497     // 1 point family on the grid
00498     if(this->FamilyOnEntity->GetParentGrid()->GetParentMesh()
00499       ->GetNumberOfPointFamily() <= 1)
00500       {
00501       this->IntersectionStatus = ProfileIncludedInFamily;
00502       }
00503     }
00504 
00505   // there is a profile, we have to compute the match between the family and
00506   // the profile
00507   vtkMedFamilyOnEntity* foe = this->GetFamilyOnEntity();
00508   vtkMedEntityArray* pea = foe->GetEntityArray();
00509   vtkMedIntArray* pIds = NULL;
00510 
00511   if(this->Profile)
00512     {
00513     this->Profile->Load();
00514     pIds=this->Profile->GetIds();
00515     }
00516 
00517   if(pIds == NULL)
00518     {
00519     if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() > 1)
00520       {
00521       this->IntersectionStatus =
00522           vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
00523       }
00524     else
00525       {
00526       this->IntersectionStatus =
00527           vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
00528       }
00529     return;
00530     }
00531 
00532   bool profile_intersects=false;
00533   bool profile_included=true;
00534   med_int famId=this->FamilyOnEntity->GetFamily()->GetId();
00535   for(vtkIdType pindex=0; pindex<pIds->GetNumberOfTuples(); pindex++)
00536     {
00537     med_int pid=pIds->GetValue(pindex)-1;
00538     med_int fid = pea->GetFamilyId(pid);
00539    // med_int fid=grid->GetPointFamilyId(pid);
00540     if(fid==famId)
00541       {// the family of the current point is the good one
00542       profile_intersects=true;
00543       }
00544     else
00545       {
00546       // we are on the profile and not on the family -->
00547       // no exact match, but the the profile might be larger than the family.
00548       profile_included=false;
00549       }
00550     }
00551 
00552   if(!profile_intersects)
00553     {
00554     this->IntersectionStatus =
00555         vtkMedFamilyOnEntityOnProfile::NoIntersection;
00556     }
00557   else if(profile_included)
00558     {
00559     this->IntersectionStatus =
00560         vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
00561     }
00562   else
00563     {
00564     this->IntersectionStatus =
00565         vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
00566     }
00567 }
00568 
00569 vtkIdType vtkMedFamilyOnEntityOnProfile::GetVTKPointIndex(vtkIdType medCIndex)
00570 {
00571   if(this->IntersectionStatus == NotComputed)
00572     this->ComputeIntersection(NULL);
00573 
00574   if(this->UseAllPoints)
00575     return medCIndex;
00576 
00577   if(this->MedToVTKPointIndexMap.find(medCIndex)
00578     == this->MedToVTKPointIndexMap.end())
00579     {
00580     vtkDebugMacro("GetVTKPointIndex asked for "
00581                   << medCIndex << " which has not been mapped");
00582     return -1;
00583     }
00584 
00585   return this->MedToVTKPointIndexMap[medCIndex];
00586 }
00587 
00588 void vtkMedFamilyOnEntityOnProfile::PrintSelf(ostream& os, vtkIndent indent)
00589 {
00590   this->Superclass::PrintSelf(os, indent);
00591 }