2 #include "MEDMEM_Connectivity.hxx"
3 #include "MEDMEM_Family.hxx"
4 #include "MEDMEM_Group.hxx"
5 #include "MEDMEM_CellModel.hxx"
7 #include "MEDMEM_SkyLineArray.hxx"
8 #include "MEDMEM_ModulusArray.hxx"
10 #include "MEDMEM_STRING.hxx"
13 using namespace MEDMEM;
15 Default Constructor. /n
16 Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
17 //--------------------------------------------------------------//
18 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
19 //--------------------------------------------------------------//
21 _typeConnectivity(MED_NODAL),
23 _geometricTypes((medGeometryElement*)NULL),
24 _type((CELLMODEL*)NULL),
27 _nodal((MEDSKYLINEARRAY*)NULL),
28 _descending((MEDSKYLINEARRAY*)NULL),
29 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
30 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
31 _neighbourhood((MEDSKYLINEARRAY*)NULL),
32 _constituent((CONNECTIVITY*)NULL)
34 MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
39 Default for Entity is MED_CELL */
40 //------------------------------------------------------------------------------//
41 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
42 //------------------------------------------------------------------------------//
44 _typeConnectivity(MED_NODAL),
45 _numberOfTypes(numberOfTypes),
47 _nodal((MEDSKYLINEARRAY*)NULL),
48 _descending((MEDSKYLINEARRAY*)NULL),
49 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
50 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
51 _neighbourhood((MEDSKYLINEARRAY*)NULL),
52 _constituent((CONNECTIVITY*)NULL)
54 MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
55 _geometricTypes = new medGeometryElement[numberOfTypes];
56 _type = new CELLMODEL[numberOfTypes];
57 _count = new int[numberOfTypes+1];
63 //------------------------------------------------------//
64 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
65 //----------------------------------------------------//
67 _typeConnectivity (m._typeConnectivity),
68 _numberOfTypes (m._numberOfTypes),
69 _entityDimension (m._entityDimension),
70 _numberOfNodes (m._numberOfNodes)
72 if (m._geometricTypes != NULL)
74 _geometricTypes = new medGeometryElement[_numberOfTypes];
75 memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
78 _geometricTypes = (medGeometryElement *) NULL;
82 _type = new CELLMODEL[_numberOfTypes];
83 for (int i=0;i<_numberOfTypes;i++)
84 _type[i] = CELLMODEL(m._type[i]);
87 _type = (CELLMODEL *) NULL;
91 _count = new med_int[_numberOfTypes+1];
92 memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(med_int));
95 _count = (med_int *) NULL;
98 _nodal = new MEDSKYLINEARRAY(* m._nodal);
100 _nodal = (MEDSKYLINEARRAY *) NULL;
102 if (m._descending != NULL)
103 _descending = new MEDSKYLINEARRAY(* m._descending);
105 _descending = (MEDSKYLINEARRAY *) NULL;
107 if (m._reverseNodalConnectivity != NULL)
108 _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
110 _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
112 if (m._reverseDescendingConnectivity != NULL)
113 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
115 _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
117 if (m._neighbourhood != NULL)
118 _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
120 _neighbourhood = (MEDSKYLINEARRAY *) NULL;
122 if (m._constituent != NULL)
123 _constituent = new CONNECTIVITY(* m._constituent);
125 _constituent = (CONNECTIVITY *) NULL;
130 desallocates existing pointers */
131 //----------------------------//
132 CONNECTIVITY::~CONNECTIVITY()
133 //----------------------------//
135 MESSAGE("Destructeur de CONNECTIVITY()");
137 if (_geometricTypes != NULL)
138 delete [] _geometricTypes;
145 if (_descending != NULL)
147 if (_reverseNodalConnectivity != NULL)
148 delete _reverseNodalConnectivity;
149 if (_reverseDescendingConnectivity != NULL)
150 delete _reverseDescendingConnectivity;
151 if (_constituent != NULL)
156 set _constituent to Constituent
157 be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
158 throws an exception if Constituent = MED_CELL
161 //----------------------------------------------------------//
162 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
164 //----------------------------------------------------------//
166 medEntityMesh Entity = Constituent->getEntity();
167 if (Entity == MED_CELL)
168 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
170 if ((Entity == MED_EDGE)&(_entityDimension == 3))
172 if (_constituent == NULL)
173 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
174 _constituent->setConstituent(Constituent);
177 _constituent = Constituent;
180 /*! Duplicated Types array in CONNECTIVITY object. */
181 //--------------------------------------------------------------------//
182 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
183 const medEntityMesh Entity)
185 //--------------------------------------------------------------------//
187 if (Entity == _entity)
188 for (int i=0; i<_numberOfTypes; i++)
190 _geometricTypes[i] = Types[i];
191 _type[i] = CELLMODEL(_geometricTypes[i]);
195 if (_constituent == NULL)
196 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
197 _constituent->setGeometricTypes(Types,Entity);
202 //--------------------------------------------------------------------//
203 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
205 //--------------------------------------------------------------------//
207 if (Entity == _entity)
209 int * index = new int[Count[_numberOfTypes]];
212 for (int i=0; i<_numberOfTypes; i++) {
213 _count[i+1] = Count[i+1];
214 int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
215 for (int j=_count[i]; j<_count[i+1]; j++)
216 index[j] = index[j-1]+NumberOfNodesPerElement;
219 if (_nodal != NULL) delete _nodal;
220 _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
221 _nodal->setIndex(index);
226 if (_constituent == NULL)
227 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
228 _constituent->setCount(Count,Entity);
232 //--------------------------------------------------------------------//
233 void CONNECTIVITY::setNodal(const int * Connectivity,
234 const medEntityMesh Entity,
235 const medGeometryElement Type)
237 //--------------------------------------------------------------------//
239 if (_entity == Entity)
241 // find geometric type
243 for (int i=0; i<_numberOfTypes; i++)
244 if (_geometricTypes[i] == Type) {
246 int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
247 //_nodal->setI(i+1,Connectivity);
248 for( int j=_count[i];j<_count[i+1]; j++)
249 _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
252 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
255 if (_constituent == NULL)
256 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
257 _constituent->setNodal(Connectivity,Entity,Type);
262 //------------------------------------------------------------------------------------------//
263 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
264 //------------------------------------------------------------------------------------------//
266 MESSAGE("CONNECTIVITY::calculateConnectivity");
268 // a temporary limitation due to calculteDescendingConnectivity function !!!
269 if ((_entityDimension==3) & (Entity==MED_EDGE))
270 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
273 if (ConnectivityType==MED_NODAL)
274 calculateNodalConnectivity();
276 if (Entity==MED_CELL)
277 calculateDescendingConnectivity();
279 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
280 if (Entity!=_entity) {
281 calculateDescendingConnectivity();
282 if (_entityDimension == 2 || _entityDimension == 3)
283 _constituent->calculateConnectivity(ConnectivityType,Entity);
287 /*! Give, in full or no interlace mode (for nodal connectivity),
288 descending or nodal connectivity.
290 You must give a %medEntityMesh (ie:MED_EDGE)
291 and a %medGeometryElement (ie:MED_SEG3). */
293 //------------------------------------------------------------//
294 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
295 //------------------------------------------------------------//
297 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
300 int numberOfFamilies = myFamilies.size();
301 if (numberOfFamilies == 0 ) {
302 MESSAGE(LOC<<"No family");
305 // does we do an update ?
306 if ((_constituent != NULL)&(_descending != NULL)) {
307 MESSAGE(LOC<<"Constituent is already defined");
311 if ((_constituent != NULL)&(_descending == NULL)) {
312 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
313 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
317 // well we could go !
318 CONNECTIVITY * oldConstituent = _constituent;
320 // for(int i=0; i<numberOfFamilies; i++) {
321 // FAMILY * myFamily = myFamilies[i];
322 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
325 _constituent = (CONNECTIVITY *)NULL;
326 // for instance we must have nodal connectivity in constituent :
327 if (oldConstituent->_nodal == NULL)
329 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
330 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
332 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
333 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
334 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
336 calculateDescendingConnectivity();
338 MESSAGE(LOC << " Right after the call to calculateDescendingConnectivity");
339 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
340 const int * newConstituentValue = _constituent->_nodal->getValue();
341 const int * newConstituentIndex = _constituent->_nodal->getIndex();
343 const int * newReverseDescendingIndex =
344 _reverseDescendingConnectivity->getIndex();
346 const int * newDescendingIndex = _descending->getIndex();
347 // const int * newDescendingValue = _descending->getValue();
349 // loop on all family,
350 // for all constituent in family, we get it's old connectivity
351 // (with oldCconstituentValue and oldConstituentIndex)
352 // and search the constituent in newConstituentValue with class
355 // when a new face is found, replace old constituent
356 // number in family with new one
357 // If face is not rigth oriented, we must change _descending attribute
358 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
360 // Voila a toi de jouer Nadir :-)
362 // First we get the renumbering from the oldCconstituentValue and
363 // oldConstituentIndex in the the new one, newConstituentValue and
364 // newConstituentIndex with a negative sign if the face is not
367 int * renumberingFromOldToNew = new int [oldNumberOfFace];
371 MESSAGE(LOC << " Right before the call to _constituent->calculateReverseNodalConnectivity");
372 _constituent->calculateReverseNodalConnectivity();
373 MESSAGE(LOC << " Right after the call to _constituent->calculateReverseNodalConnectivity");
376 SCRUTE(oldNumberOfFace);
379 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
383 renumberingFromOldToNew[iOldFace] = iOldFace+1;
384 // renumberingFromOldToNew[iOldFace] = 999999;
386 int face_it_beginOld = oldConstituentIndex[iOldFace];
387 int face_it_endOld = oldConstituentIndex[iOldFace+1];
388 int face_size_itOld = face_it_endOld - face_it_beginOld;
390 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
393 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
394 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
396 // set an array wich contains faces numbers arround first node
397 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
398 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
399 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
401 int * FacesList = new int[NumberOfFacesInList];
403 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
404 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
406 // foreach node in sub cell, we search elements which are in common
407 // at the end, we must have only one !
409 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
411 int NewNumberOfFacesInList = 0;
412 int * NewFacesList = new int[NumberOfFacesInList];
414 for (int l1=0; l1<NumberOfFacesInList; l1++) {
415 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
416 if (FacesList[l1]<reverseFaceNodal[l2-1])
417 // increasing order : FacesList[l1] are not in elements list of node l
419 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
421 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
422 NewNumberOfFacesInList++;
427 NumberOfFacesInList = NewNumberOfFacesInList;
429 FacesList = NewFacesList;
432 if (!NumberOfFacesInList==0)
434 if (NumberOfFacesInList>1)
435 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
437 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
439 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
440 int face_it_endNew = newConstituentIndex[FacesList[0]];
441 face_size_itNew = face_it_endNew - face_it_beginNew;
443 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
444 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
446 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
448 //SCRUTE(retCompareNewOld);
450 // Real new face found
452 if(retCompareNewOld == 1)
454 renumberingFromOldToNew[iOldFace] = FacesList[0];
459 // Reverse new face found
461 if(retCompareNewOld == -1)
463 renumberingFromOldToNew[iOldFace] = FacesList[0];
467 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
468 int face_it_end = newReverseDescendingIndex[FacesList[0]];
469 int face_size_it = face_it_end - face_it_begin;
471 if (face_size_it == 1)
472 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
474 if (face_size_it > 2)
475 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
477 // we have always 2 neighbourings
478 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
479 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
480 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
482 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
484 if (cell2 != 0) { // we are not on border !!!!
486 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
487 // Updating _constituent->_nodal because of reversity
488 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
489 for(int iarray=1;iarray<=face_size_itNew;iarray++){
490 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
493 // Updating _reverseDescendingConnectivity
496 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
498 // Updating _descending for cell1 and cell2
499 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
500 if (_descending->getIndexValue(iface)==FacesList[0])
501 _descending->setIndexValue(iface,-FacesList[0]);
502 else if (_descending->getIndexValue(iface)==-FacesList[0])
503 _descending->setIndexValue(iface,FacesList[0]);
504 // else nothing to do
506 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
507 if (_descending->getIndexValue(iface)==FacesList[0])
508 _descending->setIndexValue(iface,-FacesList[0]);
509 else if (_descending->getIndexValue(iface)==-FacesList[0])
510 _descending->setIndexValue(iface,FacesList[0]);
511 // else nothing to do
513 } else {// else we are on border and we do nothing !!!!!!!!
514 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
515 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
516 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
522 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
523 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
529 MESSAGE(LOC<<"The Renumbering is finished and the status is");
531 // Updating the Family
533 SCRUTE(numberOfFamilies);
535 for(int i=0; i<numberOfFamilies; i++) {
536 FAMILY * myFamily = myFamilies[i];
539 // MESSAGE(LOC<<(*myFamily));
541 if (myFamily->isOnAllElements()) {
543 // we must have more constituent ?
544 if (oldNumberOfFace==_constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS))
545 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a family which is already in all constituent !"));
546 myFamily->setAll(false);
548 int * values = new int[oldNumberOfFace] ;
549 for (int ind=0;ind<oldNumberOfFace;ind++)
552 int NumberOfTypes = myFamily->getNumberOfTypes();
553 const int * count = oldConstituent->getGlobalNumberingIndex(_constituent->getEntity());
554 int * index = new int[NumberOfTypes+1] ;
555 memcpy(index,count,(NumberOfTypes+1)*sizeof(int));
556 // build new number attribut
557 myFamily->setNumber(index,values);
560 // MESSAGE(LOC<<(*myFamily));
561 MEDSKYLINEARRAY * number = myFamily->getnumber();
563 int numberOfLines_skyline = number->getNumberOf();
566 SCRUTE(numberOfLines_skyline);
568 const int * index_skyline = number->getIndex();
570 SCRUTE(index_skyline);
572 for (int i=0;i<numberOfLines_skyline;i++) {
573 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
574 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
577 // MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
580 delete oldConstituent ;
581 delete [] renumberingFromOldToNew;
589 // meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
590 // l'heritage car les pointeurs sont dans un conteneur.
591 void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
592 //------------------------------------------------------------//
594 const char * LOC = "CONNECTIVITY::updateGroup(vector<GROUP*>) ";
597 int numberOfFamilies = myFamilies.size();
598 if (numberOfFamilies == 0 ) {
599 MESSAGE(LOC<<"No family");
602 // does we do an update ?
603 if ((_constituent != NULL)&(_descending != NULL)) {
604 MESSAGE(LOC<<"Constituent is already defined");
608 if ((_constituent != NULL)&(_descending == NULL)) {
609 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
610 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
614 // well we could go !
615 CONNECTIVITY * oldConstituent = _constituent;
617 // for(int i=0; i<numberOfFamilies; i++) {
618 // FAMILY * myFamily = myFamilies[i];
619 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
622 _constituent = (CONNECTIVITY *)NULL;
623 // for instance we must have nodal connectivity in constituent :
624 if (oldConstituent->_nodal == NULL)
626 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
627 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
629 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
630 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
631 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
633 calculateDescendingConnectivity();
635 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
636 const int * newConstituentValue = _constituent->_nodal->getValue();
637 const int * newConstituentIndex = _constituent->_nodal->getIndex();
639 const int * newReverseDescendingIndex =
640 _reverseDescendingConnectivity->getIndex();
642 const int * newDescendingIndex = _descending->getIndex();
643 // const int * newDescendingValue = _descending->getValue();
645 // loop on all family,
646 // for all constituent in family, we get it's old connectivity
647 // (with oldCconstituentValue and oldConstituentIndex)
648 // and search the constituent in newConstituentValue with class
651 // when a new face is found, replace old constituent
652 // number in family with new one
653 // If face is not rigth oriented, we must change _descending attribute
654 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
656 // Voila a toi de jouer Nadir :-)
658 // First we get the renumbering from the oldCconstituentValue and
659 // oldConstituentIndex in the the new one, newConstituentValue and
660 // newConstituentIndex with a negative sign if the face is not
663 int * renumberingFromOldToNew = new int [oldNumberOfFace];
667 _constituent->calculateReverseNodalConnectivity();
669 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
673 renumberingFromOldToNew[iOldFace] = iOldFace+1;
674 // renumberingFromOldToNew[iOldFace] = 999999;
676 int face_it_beginOld = oldConstituentIndex[iOldFace];
677 int face_it_endOld = oldConstituentIndex[iOldFace+1];
678 int face_size_itOld = face_it_endOld - face_it_beginOld;
680 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
683 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
684 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
686 // set an array wich contains faces numbers arround first node
687 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
688 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
689 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
691 int * FacesList = new int[NumberOfFacesInList];
693 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
694 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
696 // foreach node in sub cell, we search elements which are in common
697 // at the end, we must have only one !
699 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
701 int NewNumberOfFacesInList = 0;
702 int * NewFacesList = new int[NumberOfFacesInList];
704 for (int l1=0; l1<NumberOfFacesInList; l1++) {
705 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
706 if (FacesList[l1]<reverseFaceNodal[l2-1])
707 // increasing order : FacesList[l1] are not in elements list of node l
709 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
711 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
712 NewNumberOfFacesInList++;
717 NumberOfFacesInList = NewNumberOfFacesInList;
719 FacesList = NewFacesList;
722 if (!NumberOfFacesInList==0)
724 if (NumberOfFacesInList>1)
725 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
727 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
729 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
730 int face_it_endNew = newConstituentIndex[FacesList[0]];
731 face_size_itNew = face_it_endNew - face_it_beginNew;
733 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
734 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
736 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
738 //SCRUTE(retCompareNewOld);
740 // Real new face found
742 if(retCompareNewOld == 1)
744 renumberingFromOldToNew[iOldFace] = FacesList[0];
749 // Reverse new face found
751 if(retCompareNewOld == -1)
753 renumberingFromOldToNew[iOldFace] = FacesList[0];
757 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
758 int face_it_end = newReverseDescendingIndex[FacesList[0]];
759 int face_size_it = face_it_end - face_it_begin;
761 if (face_size_it == 1)
762 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
764 if (face_size_it > 2)
765 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
767 // we have always 2 neighbourings
768 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
769 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
770 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
772 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
774 if (cell2 != 0) { // we are not on border !!!!
776 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
777 // Updating _constituent->_nodal because of reversity
778 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
779 for(int iarray=1;iarray<=face_size_itNew;iarray++){
780 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
783 // Updating _reverseDescendingConnectivity
786 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
788 // Updating _descending for cell1 and cell2
789 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
790 if (_descending->getIndexValue(iface)==FacesList[0])
791 _descending->setIndexValue(iface,-FacesList[0]);
792 else if (_descending->getIndexValue(iface)==-FacesList[0])
793 _descending->setIndexValue(iface,FacesList[0]);
794 // else nothing to do
796 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
797 if (_descending->getIndexValue(iface)==FacesList[0])
798 _descending->setIndexValue(iface,-FacesList[0]);
799 else if (_descending->getIndexValue(iface)==-FacesList[0])
800 _descending->setIndexValue(iface,FacesList[0]);
801 // else nothing to do
803 } else {// else we are on border and we do nothing !!!!!!!!
804 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
805 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
806 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
812 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
813 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
819 MESSAGE(LOC<<"The Renumbering is finished and the status is");
821 // Updating the Family
823 for(int i=0; i<numberOfFamilies; i++) {
824 GROUP * myFamily = myFamilies[i];
826 MEDSKYLINEARRAY * number = myFamily->getnumber();
827 int numberOfLines_skyline = number->getNumberOf();
828 const int * index_skyline = number->getIndex();
830 for (int i=0;i<numberOfLines_skyline;i++) {
831 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
832 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
835 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
838 delete oldConstituent ;
839 delete [] renumberingFromOldToNew;
847 //------------------------------------------------------------------------------------------------------------------//
848 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
849 //------------------------------------------------------------------------------------------------------------------//
851 const char * LOC = "CONNECTIVITY::getConnectivity";
854 MEDSKYLINEARRAY * Connectivity;
855 if (Entity==_entity) {
857 if (ConnectivityType==MED_NODAL)
859 calculateNodalConnectivity();
864 calculateDescendingConnectivity();
865 Connectivity=_descending;
868 if (Connectivity!=NULL)
869 if (Type==MED_ALL_ELEMENTS)
870 return Connectivity->getValue();
872 for (med_int i=0; i<_numberOfTypes; i++)
873 if (_geometricTypes[i]==Type)
874 //return Connectivity->getI(i+1);
875 return Connectivity->getI(_count[i]);
876 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
879 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
881 if (_constituent != NULL)
882 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
884 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
887 /*! Give morse index array to use with
888 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
890 Each value give start index for corresponding entity in connectivity array.
892 Example : i-th element, j-th node of it :
893 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
894 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
895 //-----------------------------------------------------------------------------------------------//
896 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
897 //----0000000--------------------------------------------------------------------------------------------//
899 const char * LOC = "CONNECTIVITY::getConnectivityIndex";
902 MEDSKYLINEARRAY * Connectivity;
903 if (Entity==_entity) {
905 if (ConnectivityType==MED_NODAL)
908 Connectivity=_descending;
910 if (Connectivity!=NULL)
911 return Connectivity->getIndex();
913 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
916 if (_constituent != NULL)
917 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
919 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
923 //--------------------------------------------------------------//
924 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
925 //--------------------------------------------------------------//
927 const char * LOC = "CONNECTIVITY::getType";
930 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
931 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
932 for (med_int i=0; i<_numberOfTypes; i++)
933 if (_geometricTypes[i]==Type)
935 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
938 /*! Returns the number of elements of type %medGeometryElement.
939 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
940 //------------------------------------------------------------------------//
941 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
942 //------------------------------------------------------------------------//
944 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
947 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
948 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
949 for (med_int i=0; i<_numberOfTypes; i++)
950 if (_geometricTypes[i]==Type)
951 return _type[i].getNumberOfNodes();
952 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
955 /*! Returns the number of geometric sub cells of %medGeometryElement type.
956 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
957 //------------------------------------------------------------------------//
958 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
959 //------------------------------------------------------------------------//
961 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
962 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
963 for (med_int i=0; i<_numberOfTypes; i++)
964 if (_geometricTypes[i]==Type)
965 return _type[i].getNumberOfConstituents(1);
966 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
969 /*! Returns the number of elements of type %medGeometryElement.
972 - Implemented for MED_ALL_ELEMENTS
973 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
974 - Not implemented for MED_NONE */
975 //-----------------------------------------------------------------------------------//
976 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
977 //-----------------------------------------------------------------------------------//
979 const char * LOC = "CONNECTIVITY::getNumberOf";
982 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
984 if (Entity==_entity) {
986 return 0; // not defined !
987 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
988 if (Type==MED_ALL_ELEMENTS)
989 return _count[_numberOfTypes]-1;
990 for (med_int i=0; i<_numberOfTypes; i++)
991 if (_geometricTypes[i]==Type)
992 return (_count[i+1] - _count[i]);
993 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
995 if (_constituent != NULL)
996 return _constituent->getNumberOf(Entity,Type);
998 return 0; // valid if they are nothing else !
999 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
1003 //--------------------------------------------------------------//
1004 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
1005 medGeometryElement Type)
1006 //--------------------------------------------------------------//
1008 if (TypeConnectivity == MED_NODAL)
1010 calculateNodalConnectivity();
1011 if (Type==MED_ALL_ELEMENTS)
1012 return _nodal->getValue();
1013 for (med_int i=0; i<_numberOfTypes; i++)
1014 if (_geometricTypes[i]==Type)
1015 return _nodal->getI(_count[i]);
1019 calculateDescendingConnectivity();
1020 if (Type==MED_ALL_ELEMENTS)
1021 return _descending->getValue();
1022 for (med_int i=0; i<_numberOfTypes; i++)
1023 if (_geometricTypes[i]==Type)
1024 return _descending->getI(Type);
1026 throw MEDEXCEPTION("Not found");
1030 //---------------------------------------------------------------------//
1031 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
1032 //---------------------------------------------------------------------//
1034 if (TypeConnectivity == MED_NODAL)
1036 calculateNodalConnectivity();
1037 return _nodal->getIndex();
1041 calculateDescendingConnectivity();
1042 return _descending->getIndex();
1046 /*! Not yet implemented */
1047 //----------------------------------------------//
1048 const med_int* CONNECTIVITY:: getNeighbourhood() const
1049 //----------------------------------------------//
1051 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1054 /*! Returns an array which contains, for each node, all cells
1056 //-------------------------------------------------//
1057 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
1058 //-------------------------------------------------//
1060 calculateReverseNodalConnectivity();
1061 return _reverseNodalConnectivity->getValue();
1064 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
1065 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
1066 //-------------------------------------------------------//
1067 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1068 //-------------------------------------------------------//
1070 calculateReverseNodalConnectivity();
1071 return _reverseNodalConnectivity->getIndex();
1074 /*! Returns an array which contains, for each face (or edge),
1075 the 2 cells of each side. First is cell which face normal is outgoing.
1077 //------------------------------------------------------//
1078 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
1079 //------------------------------------------------------//
1081 // it is in _constituent connectivity only if we are in MED_CELL
1082 // (we could not for instance calculate face-edge connectivity !)
1083 if (_entity!=MED_CELL)
1084 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1086 // we need descending connectivity
1087 calculateDescendingConnectivity();
1088 return _reverseDescendingConnectivity->getValue();
1091 /*! calculate the reverse descending Connectivity
1092 and returns the index ( A DOCUMENTER MIEUX)*/
1093 //-----------------------------------------------------------//
1094 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1095 //-----------------------------------------------------------//
1097 // it is in _constituent connectivity only if we are in MED_CELL
1098 if (_entity!=MED_CELL)
1099 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1101 // we need descending connectivity
1102 calculateDescendingConnectivity();
1103 return _reverseDescendingConnectivity->getIndex();
1106 /*! A DOCUMENTER (et a finir ???) */
1107 //--------------------------------------------//
1108 void CONNECTIVITY::calculateNodalConnectivity()
1109 //--------------------------------------------//
1113 if (_descending==NULL)
1114 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1115 // calculate _nodal from _descending
1119 /*! If not yet done, calculate the nodal Connectivity
1120 and the reverse nodal Connectivity*/
1121 //---------------------------------------------------//
1122 void CONNECTIVITY::calculateReverseNodalConnectivity()
1123 //---------------------------------------------------//
1125 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1129 SCRUTE(_reverseNodalConnectivity);
1133 calculateNodalConnectivity();
1135 if(_reverseNodalConnectivity==NULL) {
1137 med_int node_number = 0;
1138 vector <vector <med_int> > reverse_connectivity;
1139 reverse_connectivity.resize(_numberOfNodes+1);
1141 // Treat all cells types
1143 for (med_int j = 0; j < _numberOfTypes; j++)
1145 // node number of the cell type
1146 node_number = _type[j].getNumberOfNodes();
1147 // treat all cells of a particular type
1148 for (med_int k = _count[j]; k < _count[j+1]; k++)
1149 // treat all nodes of the cell type
1150 for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1152 med_int global_node = _nodal->getIJ(k,local_node_number);
1153 reverse_connectivity[global_node].push_back(k);
1157 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1159 //calculate size of reverse_nodal_connectivity array
1160 med_int size_reverse_nodal_connectivity = 0;
1161 for (med_int i = 1; i < _numberOfNodes+1; i++)
1162 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1164 //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
1165 med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
1166 med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
1167 //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
1168 //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
1170 reverse_nodal_connectivity_index[0] = 1;
1171 for (med_int i = 1; i < _numberOfNodes+1; i++)
1173 med_int size = reverse_connectivity[i].size();
1174 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1175 for (med_int j = 0; j < size; j++)
1176 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1179 //_reverseNodalConnectivity = ReverseConnectivity;
1180 _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1181 reverse_nodal_connectivity_index,
1182 reverse_nodal_connectivity);
1183 delete [] reverse_nodal_connectivity_index;
1184 delete [] reverse_nodal_connectivity;
1189 /*! If not yet done, calculate the Descending Connectivity */
1190 //-------------------------------------------------//
1191 void CONNECTIVITY::calculateDescendingConnectivity()
1192 //-------------------------------------------------//
1194 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1197 if (_descending==NULL)
1201 MESSAGE(LOC<<"No connectivity found !");
1202 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1204 // calcul _descending from _nodal
1205 // we need CONNECTIVITY for constituent
1207 if (_constituent != NULL)
1208 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1209 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1211 if (_entityDimension == 3)
1212 _constituent = new CONNECTIVITY(MED_FACE);
1213 else if (_entityDimension == 2)
1214 _constituent = new CONNECTIVITY(MED_EDGE);
1217 MESSAGE(LOC<<"We are in 1D");
1221 _constituent->_typeConnectivity = MED_NODAL;
1222 _constituent->_numberOfNodes = _numberOfNodes;
1223 // foreach cells, we built array of constituent
1224 int DescendingSize = 0;
1225 for(int i=0; i<_numberOfTypes; i++)
1226 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1227 //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1228 //const int * descend_connectivity = _descending->getValue();
1229 int * descend_connectivity = new int[DescendingSize];
1230 for (int i=0; i<DescendingSize; i++)
1231 descend_connectivity[i]=0;
1232 //const int * descend_connectivity_index = _descending->getIndex();
1233 int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1234 descend_connectivity_index[0]=1;
1235 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1236 ConstituentsTypes[0]=MED_NONE;
1237 ConstituentsTypes[1]=MED_NONE;
1238 int * NumberOfConstituentsForeachType = new int[2];
1239 NumberOfConstituentsForeachType[0]=0;
1240 NumberOfConstituentsForeachType[1]=0;
1241 for(int i=0; i<_numberOfTypes; i++)
1243 // initialize descend_connectivity_index array :
1244 int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1245 for (int j=_count[i];j<_count[i+1];j++)
1247 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1248 // compute number of constituent of all cell for each type
1249 for(int k=1;k<NumberOfConstituents+1;k++)
1251 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1252 if (ConstituentsTypes[0]==MED_NONE)
1254 ConstituentsTypes[0]=MEDType;
1255 NumberOfConstituentsForeachType[0]++;
1256 } else if (ConstituentsTypes[0]==MEDType)
1257 NumberOfConstituentsForeachType[0]++;
1258 else if (ConstituentsTypes[1]==MED_NONE)
1260 ConstituentsTypes[1]=MEDType;
1261 NumberOfConstituentsForeachType[1]++;
1262 } else if (ConstituentsTypes[1]==MEDType)
1263 NumberOfConstituentsForeachType[1]++;
1265 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1270 // we could built _constituent !
1271 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1272 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1274 //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1276 // we use _constituent->_nodal
1277 //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1278 //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1279 int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1280 int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1281 ConstituentNodalConnectivityIndex[0]=1;
1283 _constituent->_entityDimension=ConstituentsTypes[0]/100;
1284 if (ConstituentsTypes[1]==MED_NONE)
1285 _constituent->_numberOfTypes = 1;
1287 _constituent->_numberOfTypes = 2;
1288 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1289 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1290 _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1291 _constituent->_count[0]=1;
1292 int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1293 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1294 for (int i=0; i<_constituent->_numberOfTypes;i++) {
1295 _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1296 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1297 CELLMODEL Type(ConstituentsTypes[i]);
1298 _constituent->_type[i]=Type;
1299 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1300 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1301 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1303 delete [] ConstituentsTypes;
1304 delete [] NumberOfConstituentsForeachType;
1306 // we need reverse nodal connectivity
1307 if (! _reverseNodalConnectivity)
1308 calculateReverseNodalConnectivity();
1309 const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1310 const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1312 // array to keep reverse descending connectivity
1313 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1315 int TotalNumberOfSubCell = 0;
1316 for (int i=0; i<_numberOfTypes; i++)
1317 { // loop on all cell type
1318 CELLMODEL Type = _type[i];
1319 // int NumberOfNodesPerCell = Type.getNumberOfNodes();
1320 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1321 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1322 for (int k=1; k<=NumberOfConstituentPerCell; k++)
1323 { // we loop on all sub cell of it
1324 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1325 // it is a new sub cell !
1326 // TotalNumberOfSubCell++;
1328 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1329 tmp_NumberOfConstituentsForeachType[0]++;
1330 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1332 tmp_NumberOfConstituentsForeachType[1]++;
1333 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1335 //we have maximum two types
1337 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1338 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1339 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1341 int * NodesLists = new int[NumberOfNodesPerConstituent];
1342 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1343 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1344 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1346 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1348 // all elements which contains first node of sub cell :
1349 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1350 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1351 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1353 if (NumberOfCellsInList > 0)
1354 { // we could have no element !
1355 int * CellsList = new int[NumberOfCellsInList];
1356 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1357 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1359 // foreach node in sub cell, we search elements which are in common
1360 // at the end, we must have only one !
1362 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1363 int NewNumberOfCellsInList = 0;
1364 int * NewCellsList = new int[NumberOfCellsInList];
1365 for (int l1=0; l1<NumberOfCellsInList; l1++)
1366 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
1368 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1369 // increasing order : CellsList[l1] are not in elements list of node l
1371 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1372 // we have found one
1373 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1374 NewNumberOfCellsInList++;
1378 NumberOfCellsInList = NewNumberOfCellsInList;
1380 delete [] CellsList;
1381 CellsList = NewCellsList;
1384 if (NumberOfCellsInList > 0) { // We have found some elements !
1385 int CellNumber = CellsList[0];
1386 //delete [] CellsList;
1387 if (NumberOfCellsInList>1)
1388 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1390 if (NumberOfCellsInList==1)
1392 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1394 // we search sub cell number in this cell to not calculate it another time
1397 for (int l=0; l<_numberOfTypes; l++)
1398 if (CellNumber < _count[l+1]) {
1402 //int sub_cell_count2 = Type2.get_entities_count(1);
1403 //int nodes_cell_count2 = Type2.get_nodes_count();
1405 for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1407 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1408 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
1410 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1413 if (counter==Type.getConstituentType(1,k)%100)
1415 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1422 MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1425 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1427 delete [] CellsList;
1430 delete [] NodesLists;
1434 // we adjust _constituent
1435 int NumberOfConstituent=0;
1436 int SizeOfConstituentNodal=0;
1437 for (int i=0;i<_constituent->_numberOfTypes; i++) {
1438 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1439 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1441 // we built new _nodal attribute in _constituent
1442 //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1443 //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1444 //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1445 int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1446 int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1447 ConstituentNodalIndex[0]=1;
1448 // we build _reverseDescendingConnectivity
1449 //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1450 //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1451 //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1452 int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1453 int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1454 reverseDescendingConnectivityIndex[0]=1;
1456 // first constituent type
1457 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1458 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1459 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1460 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1462 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1463 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1464 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1467 // second type if any
1468 if (_constituent->_numberOfTypes==2) {
1469 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1470 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1471 int offset2=offset*2;
1472 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1473 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1474 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1475 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1476 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1478 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1479 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1480 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1483 _constituent->_count[2]=NumberOfConstituent+1;
1484 // we correct _descending to adjust face number
1485 for(int j=0;j<DescendingSize;j++)
1486 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1487 descend_connectivity[j]-=offset;
1490 delete [] ConstituentNodalConnectivityIndex;
1491 delete [] ConstituentNodalConnectivity;
1493 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1495 descend_connectivity_index,
1496 descend_connectivity);
1497 delete [] descend_connectivity_index;
1498 delete [] descend_connectivity;
1499 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1500 2*NumberOfConstituent,
1501 reverseDescendingConnectivityIndex,
1502 reverseDescendingConnectivityValue);
1503 delete [] reverseDescendingConnectivityIndex;
1504 delete [] reverseDescendingConnectivityValue;
1506 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1507 delete [] tmp_NumberOfConstituentsForeachType;
1509 //delete _constituent->_nodal;
1510 _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1511 SizeOfConstituentNodal,
1512 ConstituentNodalIndex,
1513 ConstituentNodalValue);
1515 delete [] ConstituentNodalIndex;
1516 delete [] ConstituentNodalValue;
1518 delete [] ReverseDescendingConnectivityValue;
1523 //-----------------------------------------------------------------------------------------//
1524 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1525 //-----------------------------------------------------------------------------------------//
1527 // if (_entity==MED_CELL)
1528 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1529 // Connectivity : we could not do that with MED_CELL entity !");
1530 // if (myConnectivity->getEntity()!=MED_CELL)
1531 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1532 // Connectivity : we need MED_CELL connectivity !");
1536 /*! Not implemented yet */
1537 //--------------------------------------------------------------------//
1538 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1539 //--------------------------------------------------------------------//
1541 const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1544 MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1553 Returns the geometry of an element given by its entity type & its global number.
1555 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1557 //--------------------------------------------------------------------//
1558 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1559 //--------------------------------------------------------------------//
1561 const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1563 int globalNumberMin = 1;
1564 int globalNumberMax ;
1566 if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
1567 else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
1569 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1571 // The globalNumber must verify : 1 <= globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1573 if ( (globalNumber < globalNumberMin) || (globalNumber > globalNumberMax-1 ) )
1574 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1575 << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
1577 if (_entity==Entity) {
1578 for(int i=1; i<=_numberOfTypes;i++)
1579 if (globalNumber<_count[i])
1580 return _geometricTypes[i-1];
1582 else if (_constituent!=NULL)
1583 return _constituent->getElementType(Entity,globalNumber);
1585 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1586 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1592 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1594 os << endl << "------------- Entity = ";
1609 case MED_ALL_ENTITIES:
1610 os << "MED_ALL_ENTITIES";
1616 os << " -------------\n\nMedConnectivity : ";
1617 switch (co._typeConnectivity)
1620 os << "MED_NODAL\n";
1622 case MED_DESCENDING:
1623 os << "MED_DESCENDING\n";
1628 os << "Entity dimension : " << co._entityDimension << endl;
1629 os << "Number of nodes : " << co._numberOfNodes << endl;
1630 os << "Number of types : " << co._numberOfTypes << endl;
1631 for (int i=0; i!=co._numberOfTypes ; ++i)
1632 os << " -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1633 << co._count[i+1]-co._count[i] << " elements" << endl;
1635 if (co._typeConnectivity == MED_NODAL )
1637 for (int i=0; i<co._numberOfTypes; i++)
1639 os << endl << co._type[i].getName() << " : " << endl;
1640 int numberofelements = co._count[i+1]-co._count[i];
1641 const med_int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1642 int numberofnodespercell = co._geometricTypes[i]%100;
1643 for (int j=0;j<numberofelements;j++)
1645 os << setw(6) << j+1 << " : ";
1646 for (int k=0;k<numberofnodespercell;k++)
1647 os << connectivity[j*numberofnodespercell+k]<<" ";
1652 else if (co._typeConnectivity == MED_DESCENDING)
1654 int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1655 const med_int *connectivity = co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1656 const int *connectivity_index = co.getConnectivityIndex( co._typeConnectivity, co._entity );
1658 for ( int j=0; j!=numberofelements; j++ )
1660 os << "Element " << j+1 << " : ";
1661 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1662 os << connectivity[k-1] << " ";
1667 if (co._constituent)
1668 os << endl << *co._constituent << endl;