1 // Copyright (C) 2010-2011 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "vtkMedFamilyOnEntityOnProfile.h"
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"
40 #include "vtkBitArray.h"
41 #include "vtkIdList.h"
43 #include "vtkMultiProcessController.h"
45 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile,FamilyOnEntity, vtkMedFamilyOnEntity);
46 vtkCxxSetObjectMacro(vtkMedFamilyOnEntityOnProfile, Profile, vtkMedProfile);
48 // vtkCxxRevisionMacro(vtkMedFamilyOnEntityOnProfile, "$Revision$");
49 vtkStandardNewMacro(vtkMedFamilyOnEntityOnProfile)
51 vtkMedFamilyOnEntityOnProfile::vtkMedFamilyOnEntityOnProfile()
53 this->FamilyOnEntity = NULL;
55 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
56 this->UseAllPoints = false;
57 this->MatchComputed = false;
61 vtkMedFamilyOnEntityOnProfile::~vtkMedFamilyOnEntityOnProfile()
63 this->SetFamilyOnEntity(NULL);
64 this->SetProfile(NULL);
67 bool vtkMedFamilyOnEntityOnProfile::KeepPoint(med_int index)
69 if(this->IntersectionStatus == NotComputed)
70 this->ComputeIntersection(NULL);
72 if(this->UseAllPoints)
75 if(this->MedToVTKPointIndexMap.find(index)
76 == this->MedToVTKPointIndexMap.end())
82 bool vtkMedFamilyOnEntityOnProfile::KeepCell(med_int index)
84 if(this->FamilyOnEntity->GetEntityArray()->GetFamilyId(index)
85 != this->FamilyOnEntity->GetFamily()->GetId())
90 int vtkMedFamilyOnEntityOnProfile::CanMapField(vtkMedFieldOnProfile* fop)
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)
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)
104 vtkMultiProcessController* controller =
105 vtkMultiProcessController::GetGlobalController();
106 if (controller != NULL)
108 numProc = controller->GetNumberOfProcesses();
111 if ((this->GetValid() == 0) && numProc == 1)
114 this->ComputeIntersection(fop);
116 if(this->IntersectionStatus == vtkMedFamilyOnEntityOnProfile::NoIntersection)
120 this->GetFamilyOnEntity()->GetEntity().EntityType != MED_NODE &&
121 fop->GetParentFieldOverEntity()->GetEntity().EntityType == MED_NODE &&
122 this->PointProfileMatch[fop->GetProfile()] == BadOrNoIntersection)
128 int vtkMedFamilyOnEntityOnProfile::CanShallowCopy(vtkMedFieldOnProfile *fop)
132 bool shallow_on_points = this->CanShallowCopyPointField(NULL);
133 bool shallow_on_cells = this->CanShallowCopyCellField(NULL);
134 if(shallow_on_points && shallow_on_cells)
136 if(!shallow_on_points && !shallow_on_cells)
138 vtkErrorMacro("CanShallowCopy cannot answer : is it a point or a cell field?");
142 if(fop->GetParentFieldOverEntity()->GetParentStep()->GetParentField()
143 ->GetFieldType() == vtkMedField::PointField)
144 return this->CanShallowCopyPointField(fop);
146 return this->CanShallowCopyCellField(fop);
149 void vtkMedFamilyOnEntityOnProfile::ComputeIntersection(vtkMedFieldOnProfile* fop)
151 int nodeOrCellSupport=this->GetFamilyOnEntity()->GetPointOrCell();
155 fieldType = fop->GetParentFieldOverEntity()->GetParentStep()->
156 GetParentField()->GetFieldType();
160 fieldType = (nodeOrCellSupport==vtkMedUtilities::OnPoint?vtkMedField::PointField:vtkMedField::CellField);
162 // Cell fields cannot match point supports
163 if(fieldType != vtkMedField::PointField
164 && nodeOrCellSupport == vtkMedUtilities::OnPoint)
166 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
167 this->UseAllPoints = false;
169 else if(fieldType != vtkMedField::PointField
170 && nodeOrCellSupport ==vtkMedUtilities::OnCell)
172 this->ComputeCellFamilyVsCellProfileMatch();
174 else if(fieldType == vtkMedField::PointField
175 && nodeOrCellSupport ==vtkMedUtilities::OnPoint)
177 vtkMedProfile* profile = NULL;
180 profile = fop->GetProfile();
182 // point fields must share the same profile as the point support.
183 this->ComputePointFamilyVsPointProfileMatch();
186 else if(fieldType == vtkMedField::PointField
187 && nodeOrCellSupport == vtkMedUtilities::OnCell)
189 vtkMedProfile* profile = NULL;
192 profile = fop->GetProfile();
194 this->ComputeCellFamilyVsPointProfileMatch(profile);
197 this->MatchComputed = true;
200 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyPointField(vtkMedFieldOnProfile* fop)
202 vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
203 if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnCell)
205 if(this->PointProfileMatch.find(profile) == this->PointProfileMatch.end())
207 this->ComputeCellFamilyVsPointProfileMatch(profile);
209 int match = this->PointProfileMatch[profile];
211 == vtkMedFamilyOnEntityOnProfile::ProfileEqualsSupport;
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)
231 int vtkMedFamilyOnEntityOnProfile::CanShallowCopyCellField(vtkMedFieldOnProfile* fop)
233 vtkMedProfile* profile = (fop != NULL?fop->GetProfile(): NULL);
234 // cell fields cannot be mapped to cell supports
235 if(this->FamilyOnEntity->GetPointOrCell() == vtkMedUtilities::OnPoint)
240 // this is a cell support.
241 if(this->Profile != profile)
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)
256 void vtkMedFamilyOnEntityOnProfile::ComputeUsedPoints()
258 this->MedToVTKPointIndexMap.clear();
260 //first test a few special cases where no heavy computing is necessary
261 vtkMedGrid* grid = this->FamilyOnEntity->GetParentGrid();
262 if(this->Profile == NULL)
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)
269 this->UseAllPoints = true;
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)
278 this->UseAllPoints = true;
283 vtkSmartPointer<vtkBitArray> flag = vtkSmartPointer<vtkBitArray>::New();
284 flag->SetNumberOfTuples(grid->GetNumberOfPoints());
286 // initialize the array to false
287 for(vtkIdType pid = 0; pid < flag->GetNumberOfTuples(); pid++)
288 flag->SetValue(pid, false);
290 // for each cell, flag the used points
292 this->Profile->Load();
294 vtkMedIntArray* pids=(this->Profile!=NULL?this->Profile->GetIds():NULL);
296 med_int famId = this->FamilyOnEntity->GetFamily()->GetId();
297 vtkMedEntityArray* array = this->FamilyOnEntity->GetEntityArray();
298 vtkSmartPointer<vtkIdList> ids = vtkSmartPointer<vtkIdList>::New();
300 array->LoadConnectivity();
302 vtkIdType pflsize = (pids != NULL ? pids->GetNumberOfTuples():array->GetNumberOfEntity());
303 for(vtkIdType pindex=0; pindex<pflsize; pindex++)
305 med_int pid = (pids != NULL ? pids->GetValue(pindex)-1 : pindex);
306 med_int fid = array->GetFamilyId(pid);
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++)
314 vtkIdType subid = ids->GetId(id);
315 if(subid < 0 || subid >= flag->GetNumberOfTuples())
317 vtkDebugMacro("invalid sub id : " << subid);
321 flag->SetValue(subid, 1);
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++)
330 if(flag->GetValue(pid) == false)
332 this->UseAllPoints = false;
337 if(!this->UseAllPoints)
339 // If all points are not used, I compute the index mapping
340 vtkIdType vtk_index = 0;
342 for(vtkIdType pid=0; pid < flag->GetNumberOfTuples(); pid++)
344 if(flag->GetValue(pid) == true)
346 this->MedToVTKPointIndexMap[pid] = vtk_index;
353 void vtkMedFamilyOnEntityOnProfile::ComputeCellFamilyVsCellProfileMatch()
355 if(this->MatchComputed)
358 // this computes the UseAllPoints flag.
359 this->ComputeUsedPoints();
361 if(this->Profile == NULL)
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)
367 this->IntersectionStatus =
368 vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
372 this->IntersectionStatus =
373 vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
378 this->Profile->Load();
379 vtkMedIntArray* pids=this->Profile->GetIds();
383 vtkErrorMacro("Could not load profile indices!");
384 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NotComputed;
385 this->UseAllPoints = false;
389 med_int famId = this->GetFamilyOnEntity()->GetFamily()->GetId();
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++)
396 med_int pid=pids->GetValue(pindex)-1;
397 med_int fid=array->GetFamilyId(pid);
399 {// the current cell is on the familyand on the profile
400 // --> there is an overlap
401 profile_intersect = true;
405 // the cell is on the profile but not on the family --> no inclusion
406 profile_included=false;
410 if(profile_included && profile_intersect)
411 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
412 else if(profile_intersect)
413 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
415 this->IntersectionStatus = vtkMedFamilyOnEntityOnProfile::NoIntersection;
418 void vtkMedFamilyOnEntityOnProfile::
419 ComputeCellFamilyVsPointProfileMatch(vtkMedProfile* profile)
421 // first test if the cache is already set
422 if(this->PointProfileMatch.find(profile) != this->PointProfileMatch.end())
425 // this will compute the cell match cache, as well as the UseAllPoints flag.
426 this->ComputeCellFamilyVsCellProfileMatch();
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);
439 // if profile is not NULL and I use all points --> BadOrNoIntersection
440 if(this->UseAllPoints)
442 this->PointProfileMatch[profile] = BadOrNoIntersection;
447 vtkMedIntArray* pids=profile->GetIds();
451 vtkErrorMacro("profile indices could not be loaded!");
452 this->PointProfileMatch[profile] = BadOrNoIntersection;
457 bool exact_match = true;
458 vtkIdType numberOfUsedPoints = pids->GetNumberOfTuples();
459 for(med_int pindex=0; pindex < pids->GetNumberOfTuples(); pindex++)
461 med_int id = pids->GetValue(pindex);
462 if(this->MedToVTKPointIndexMap.find(id-1) == this->MedToVTKPointIndexMap.end())
464 // The given point profile index is not used by this support.
465 // the superposition is at most partial.
467 numberOfUsedPoints--;
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())
475 this->PointProfileMatch[profile] = BadOrNoIntersection;
479 this->PointProfileMatch[profile] = ProfileEqualsSupport;
483 this->PointProfileMatch[profile] = ProfileLargerThanSupport;
487 void vtkMedFamilyOnEntityOnProfile::ComputePointFamilyVsPointProfileMatch()
489 if(this->MatchComputed)
492 this->ComputeUsedPoints();
494 if(this->Profile == NULL)
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)
501 this->IntersectionStatus = ProfileIncludedInFamily;
505 // there is a profile, we have to compute the match between the family and
507 vtkMedFamilyOnEntity* foe = this->GetFamilyOnEntity();
508 vtkMedEntityArray* pea = foe->GetEntityArray();
509 vtkMedIntArray* pIds = NULL;
513 this->Profile->Load();
514 pIds=this->Profile->GetIds();
519 if(this->FamilyOnEntity->GetEntityArray()->GetNumberOfFamilyOnEntity() > 1)
521 this->IntersectionStatus =
522 vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
526 this->IntersectionStatus =
527 vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
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++)
537 med_int pid=pIds->GetValue(pindex)-1;
538 med_int fid = pea->GetFamilyId(pid);
539 // med_int fid=grid->GetPointFamilyId(pid);
541 {// the family of the current point is the good one
542 profile_intersects=true;
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;
552 if(!profile_intersects)
554 this->IntersectionStatus =
555 vtkMedFamilyOnEntityOnProfile::NoIntersection;
557 else if(profile_included)
559 this->IntersectionStatus =
560 vtkMedFamilyOnEntityOnProfile::ProfileIncludedInFamily;
564 this->IntersectionStatus =
565 vtkMedFamilyOnEntityOnProfile::ProfileIntersectsFamily;
569 vtkIdType vtkMedFamilyOnEntityOnProfile::GetVTKPointIndex(vtkIdType medCIndex)
571 if(this->IntersectionStatus == NotComputed)
572 this->ComputeIntersection(NULL);
574 if(this->UseAllPoints)
577 if(this->MedToVTKPointIndexMap.find(medCIndex)
578 == this->MedToVTKPointIndexMap.end())
580 vtkDebugMacro("GetVTKPointIndex asked for "
581 << medCIndex << " which has not been mapped");
585 return this->MedToVTKPointIndexMap[medCIndex];
588 void vtkMedFamilyOnEntityOnProfile::PrintSelf(ostream& os, vtkIndent indent)
590 this->Superclass::PrintSelf(os, indent);