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 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
339 const int * newConstituentValue = _constituent->_nodal->getValue();
340 const int * newConstituentIndex = _constituent->_nodal->getIndex();
342 const int * newReverseDescendingIndex =
343 _reverseDescendingConnectivity->getIndex();
345 const int * newDescendingIndex = _descending->getIndex();
346 // const int * newDescendingValue = _descending->getValue();
348 // loop on all family,
349 // for all constituent in family, we get it's old connectivity
350 // (with oldCconstituentValue and oldConstituentIndex)
351 // and search the constituent in newConstituentValue with class
354 // when a new face is found, replace old constituent
355 // number in family with new one
356 // If face is not rigth oriented, we must change _descending attribute
357 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
359 // Voila a toi de jouer Nadir :-)
361 // First we get the renumbering from the oldCconstituentValue and
362 // oldConstituentIndex in the the new one, newConstituentValue and
363 // newConstituentIndex with a negative sign if the face is not
366 int * renumberingFromOldToNew = new int [oldNumberOfFace];
370 _constituent->calculateReverseNodalConnectivity();
372 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
376 renumberingFromOldToNew[iOldFace] = iOldFace+1;
377 // renumberingFromOldToNew[iOldFace] = 999999;
379 int face_it_beginOld = oldConstituentIndex[iOldFace];
380 int face_it_endOld = oldConstituentIndex[iOldFace+1];
381 int face_size_itOld = face_it_endOld - face_it_beginOld;
383 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
386 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
387 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
389 // set an array wich contains faces numbers arround first node
390 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
391 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
392 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
394 int * FacesList = new int[NumberOfFacesInList];
396 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
397 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
399 // foreach node in sub cell, we search elements which are in common
400 // at the end, we must have only one !
402 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
404 int NewNumberOfFacesInList = 0;
405 int * NewFacesList = new int[NumberOfFacesInList];
407 for (int l1=0; l1<NumberOfFacesInList; l1++) {
408 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
409 if (FacesList[l1]<reverseFaceNodal[l2-1])
410 // increasing order : FacesList[l1] are not in elements list of node l
412 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
414 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
415 NewNumberOfFacesInList++;
420 NumberOfFacesInList = NewNumberOfFacesInList;
422 FacesList = NewFacesList;
425 if (!NumberOfFacesInList==0)
427 if (NumberOfFacesInList>1)
428 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
430 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
432 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
433 int face_it_endNew = newConstituentIndex[FacesList[0]];
434 face_size_itNew = face_it_endNew - face_it_beginNew;
436 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
437 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
439 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
441 //SCRUTE(retCompareNewOld);
443 // Real new face found
445 if(retCompareNewOld == 1)
447 renumberingFromOldToNew[iOldFace] = FacesList[0];
452 // Reverse new face found
454 if(retCompareNewOld == -1)
456 renumberingFromOldToNew[iOldFace] = FacesList[0];
460 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
461 int face_it_end = newReverseDescendingIndex[FacesList[0]];
462 int face_size_it = face_it_end - face_it_begin;
464 if (face_size_it == 1)
465 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
467 if (face_size_it > 2)
468 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
470 // we have always 2 neighbourings
471 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
472 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
473 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
475 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
477 if (cell2 != 0) { // we are not on border !!!!
479 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
480 // Updating _constituent->_nodal because of reversity
481 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
482 for(int iarray=1;iarray<=face_size_itNew;iarray++){
483 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
486 // Updating _reverseDescendingConnectivity
489 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
491 // Updating _descending for cell1 and cell2
492 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
493 if (_descending->getIndexValue(iface)==FacesList[0])
494 _descending->setIndexValue(iface,-FacesList[0]);
495 else if (_descending->getIndexValue(iface)==-FacesList[0])
496 _descending->setIndexValue(iface,FacesList[0]);
497 // else nothing to do
499 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];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 } else {// else we are on border and we do nothing !!!!!!!!
507 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
508 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
509 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
515 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
516 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
522 MESSAGE(LOC<<"The Renumbering is finished and the status is");
524 // Updating the Family
526 for(int i=0; i<numberOfFamilies; i++) {
527 FAMILY * myFamily = myFamilies[i];
529 MEDSKYLINEARRAY * number = myFamily->getnumber();
530 int numberOfLines_skyline = number->getNumberOf();
531 const int * index_skyline = number->getIndex();
533 for (int i=0;i<numberOfLines_skyline;i++) {
534 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
535 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
538 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
541 delete oldConstituent ;
542 delete [] renumberingFromOldToNew;
550 // meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
551 // l'heritage car les pointeurs sont dans un conteneur.
552 void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
553 //------------------------------------------------------------//
555 const char * LOC = "CONNECTIVITY::updateFamily(vector<GROUP*>) ";
558 int numberOfFamilies = myFamilies.size();
559 if (numberOfFamilies == 0 ) {
560 MESSAGE(LOC<<"No family");
563 // does we do an update ?
564 if ((_constituent != NULL)&(_descending != NULL)) {
565 MESSAGE(LOC<<"Constituent is already defined");
569 if ((_constituent != NULL)&(_descending == NULL)) {
570 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
571 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
575 // well we could go !
576 CONNECTIVITY * oldConstituent = _constituent;
578 // for(int i=0; i<numberOfFamilies; i++) {
579 // FAMILY * myFamily = myFamilies[i];
580 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
583 _constituent = (CONNECTIVITY *)NULL;
584 // for instance we must have nodal connectivity in constituent :
585 if (oldConstituent->_nodal == NULL)
587 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
588 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
590 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
591 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
592 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
594 calculateDescendingConnectivity();
596 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
597 const int * newConstituentValue = _constituent->_nodal->getValue();
598 const int * newConstituentIndex = _constituent->_nodal->getIndex();
600 const int * newReverseDescendingIndex =
601 _reverseDescendingConnectivity->getIndex();
603 const int * newDescendingIndex = _descending->getIndex();
604 // const int * newDescendingValue = _descending->getValue();
606 // loop on all family,
607 // for all constituent in family, we get it's old connectivity
608 // (with oldCconstituentValue and oldConstituentIndex)
609 // and search the constituent in newConstituentValue with class
612 // when a new face is found, replace old constituent
613 // number in family with new one
614 // If face is not rigth oriented, we must change _descending attribute
615 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
617 // Voila a toi de jouer Nadir :-)
619 // First we get the renumbering from the oldCconstituentValue and
620 // oldConstituentIndex in the the new one, newConstituentValue and
621 // newConstituentIndex with a negative sign if the face is not
624 int * renumberingFromOldToNew = new int [oldNumberOfFace];
628 _constituent->calculateReverseNodalConnectivity();
630 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
634 renumberingFromOldToNew[iOldFace] = iOldFace+1;
635 // renumberingFromOldToNew[iOldFace] = 999999;
637 int face_it_beginOld = oldConstituentIndex[iOldFace];
638 int face_it_endOld = oldConstituentIndex[iOldFace+1];
639 int face_size_itOld = face_it_endOld - face_it_beginOld;
641 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
644 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
645 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
647 // set an array wich contains faces numbers arround first node
648 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
649 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
650 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
652 int * FacesList = new int[NumberOfFacesInList];
654 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
655 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
657 // foreach node in sub cell, we search elements which are in common
658 // at the end, we must have only one !
660 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
662 int NewNumberOfFacesInList = 0;
663 int * NewFacesList = new int[NumberOfFacesInList];
665 for (int l1=0; l1<NumberOfFacesInList; l1++) {
666 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
667 if (FacesList[l1]<reverseFaceNodal[l2-1])
668 // increasing order : FacesList[l1] are not in elements list of node l
670 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
672 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
673 NewNumberOfFacesInList++;
678 NumberOfFacesInList = NewNumberOfFacesInList;
680 FacesList = NewFacesList;
683 if (!NumberOfFacesInList==0)
685 if (NumberOfFacesInList>1)
686 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
688 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
690 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
691 int face_it_endNew = newConstituentIndex[FacesList[0]];
692 face_size_itNew = face_it_endNew - face_it_beginNew;
694 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
695 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
697 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
699 //SCRUTE(retCompareNewOld);
701 // Real new face found
703 if(retCompareNewOld == 1)
705 renumberingFromOldToNew[iOldFace] = FacesList[0];
710 // Reverse new face found
712 if(retCompareNewOld == -1)
714 renumberingFromOldToNew[iOldFace] = FacesList[0];
718 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
719 int face_it_end = newReverseDescendingIndex[FacesList[0]];
720 int face_size_it = face_it_end - face_it_begin;
722 if (face_size_it == 1)
723 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
725 if (face_size_it > 2)
726 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
728 // we have always 2 neighbourings
729 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
730 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
731 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
733 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
735 if (cell2 != 0) { // we are not on border !!!!
737 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
738 // Updating _constituent->_nodal because of reversity
739 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
740 for(int iarray=1;iarray<=face_size_itNew;iarray++){
741 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
744 // Updating _reverseDescendingConnectivity
747 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
749 // Updating _descending for cell1 and cell2
750 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
751 if (_descending->getIndexValue(iface)==FacesList[0])
752 _descending->setIndexValue(iface,-FacesList[0]);
753 else if (_descending->getIndexValue(iface)==-FacesList[0])
754 _descending->setIndexValue(iface,FacesList[0]);
755 // else nothing to do
757 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
758 if (_descending->getIndexValue(iface)==FacesList[0])
759 _descending->setIndexValue(iface,-FacesList[0]);
760 else if (_descending->getIndexValue(iface)==-FacesList[0])
761 _descending->setIndexValue(iface,FacesList[0]);
762 // else nothing to do
764 } else {// else we are on border and we do nothing !!!!!!!!
765 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
766 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
767 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
773 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
774 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
780 MESSAGE(LOC<<"The Renumbering is finished and the status is");
782 // Updating the Family
784 for(int i=0; i<numberOfFamilies; i++) {
785 GROUP * myFamily = myFamilies[i];
787 MEDSKYLINEARRAY * number = myFamily->getnumber();
788 int numberOfLines_skyline = number->getNumberOf();
789 const int * index_skyline = number->getIndex();
791 for (int i=0;i<numberOfLines_skyline;i++) {
792 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
793 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
796 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
799 delete oldConstituent ;
800 delete [] renumberingFromOldToNew;
808 //------------------------------------------------------------------------------------------------------------------//
809 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
810 //------------------------------------------------------------------------------------------------------------------//
812 const char * LOC = "CONNECTIVITY::getConnectivity";
815 MEDSKYLINEARRAY * Connectivity;
816 if (Entity==_entity) {
818 if (ConnectivityType==MED_NODAL)
820 calculateNodalConnectivity();
825 calculateDescendingConnectivity();
826 Connectivity=_descending;
829 if (Connectivity!=NULL)
830 if (Type==MED_ALL_ELEMENTS)
831 return Connectivity->getValue();
833 for (med_int i=0; i<_numberOfTypes; i++)
834 if (_geometricTypes[i]==Type)
835 //return Connectivity->getI(i+1);
836 return Connectivity->getI(_count[i]);
837 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
840 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
842 if (_constituent != NULL)
843 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
845 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
848 /*! Give morse index array to use with
849 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
851 Each value give start index for corresponding entity in connectivity array.
853 Example : i-th element, j-th node of it :
854 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
855 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
856 //-----------------------------------------------------------------------------------------------//
857 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
858 //----0000000--------------------------------------------------------------------------------------------//
860 const char * LOC = "CONNECTIVITY::getConnectivityIndex";
863 MEDSKYLINEARRAY * Connectivity;
864 if (Entity==_entity) {
866 if (ConnectivityType==MED_NODAL)
869 Connectivity=_descending;
871 if (Connectivity!=NULL)
872 return Connectivity->getIndex();
874 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
877 if (_constituent != NULL)
878 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
880 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
884 //--------------------------------------------------------------//
885 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
886 //--------------------------------------------------------------//
888 const char * LOC = "CONNECTIVITY::getType";
891 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
892 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
893 for (med_int i=0; i<_numberOfTypes; i++)
894 if (_geometricTypes[i]==Type)
896 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
899 /*! Returns the number of elements of type <med geometrie element>.
900 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
901 //------------------------------------------------------------------------//
902 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
903 //------------------------------------------------------------------------//
905 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
908 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
909 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
910 for (med_int i=0; i<_numberOfTypes; i++)
911 if (_geometricTypes[i]==Type)
912 return _type[i].getNumberOfNodes();
913 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
916 /*! Returns the number of geometric sub cells of <med geometrie element> type.
917 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
918 //------------------------------------------------------------------------//
919 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
920 //------------------------------------------------------------------------//
922 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
923 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
924 for (med_int i=0; i<_numberOfTypes; i++)
925 if (_geometricTypes[i]==Type)
926 return _type[i].getNumberOfConstituents(1);
927 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
930 /*! Returns the number of elements of type <med geometrie element>.
933 - Implemented for MED_ALL_ELEMENTS
934 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
935 - Not implemented for MED_NONE */
936 //-----------------------------------------------------------------------------------//
937 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
938 //-----------------------------------------------------------------------------------//
940 const char * LOC = "CONNECTIVITY::getNumberOf";
943 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
945 if (Entity==_entity) {
947 return 0; // not defined !
948 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
949 if (Type==MED_ALL_ELEMENTS)
950 return _count[_numberOfTypes]-1;
951 for (med_int i=0; i<_numberOfTypes; i++)
952 if (_geometricTypes[i]==Type)
953 return (_count[i+1] - _count[i]);
954 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
956 if (_constituent != NULL)
957 return _constituent->getNumberOf(Entity,Type);
959 return 0; // valid if they are nothing else !
960 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
964 //--------------------------------------------------------------//
965 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
966 medGeometryElement Type)
967 //--------------------------------------------------------------//
969 if (TypeConnectivity == MED_NODAL)
971 calculateNodalConnectivity();
972 if (Type==MED_ALL_ELEMENTS)
973 return _nodal->getValue();
974 for (med_int i=0; i<_numberOfTypes; i++)
975 if (_geometricTypes[i]==Type)
976 return _nodal->getI(_count[i]);
980 calculateDescendingConnectivity();
981 if (Type==MED_ALL_ELEMENTS)
982 return _descending->getValue();
983 for (med_int i=0; i<_numberOfTypes; i++)
984 if (_geometricTypes[i]==Type)
985 return _descending->getI(Type);
987 throw MEDEXCEPTION("Not found");
991 //---------------------------------------------------------------------//
992 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
993 //---------------------------------------------------------------------//
995 if (TypeConnectivity == MED_NODAL)
997 calculateNodalConnectivity();
998 return _nodal->getIndex();
1002 calculateDescendingConnectivity();
1003 return _descending->getIndex();
1007 /*! Not yet implemented */
1008 //----------------------------------------------//
1009 const med_int* CONNECTIVITY:: getNeighbourhood() const
1010 //----------------------------------------------//
1012 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1015 /*! Returns an array which contains, for each node, all cells
1017 //-------------------------------------------------//
1018 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
1019 //-------------------------------------------------//
1021 calculateReverseNodalConnectivity();
1022 return _reverseNodalConnectivity->getValue();
1025 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
1026 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
1027 //-------------------------------------------------------//
1028 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1029 //-------------------------------------------------------//
1031 calculateReverseNodalConnectivity();
1032 return _reverseNodalConnectivity->getIndex();
1035 /*! Returns an array which contains, for each face (or edge),
1036 the 2 cells of each side. First is cell which face normal is outgoing.
1038 //------------------------------------------------------//
1039 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
1040 //------------------------------------------------------//
1042 // it is in _constituent connectivity only if we are in MED_CELL
1043 // (we could not for instance calculate face-edge connectivity !)
1044 if (_entity!=MED_CELL)
1045 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1047 // we need descending connectivity
1048 calculateDescendingConnectivity();
1049 return _reverseDescendingConnectivity->getValue();
1052 /*! calculate the reverse descending Connectivity
1053 and returns the index ( A DOCUMENTER MIEUX)*/
1054 //-----------------------------------------------------------//
1055 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1056 //-----------------------------------------------------------//
1058 // it is in _constituent connectivity only if we are in MED_CELL
1059 if (_entity!=MED_CELL)
1060 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1062 // we need descending connectivity
1063 calculateDescendingConnectivity();
1064 return _reverseDescendingConnectivity->getIndex();
1067 /*! A DOCUMENTER (et a finir ???) */
1068 //--------------------------------------------//
1069 void CONNECTIVITY::calculateNodalConnectivity()
1070 //--------------------------------------------//
1074 if (_descending==NULL)
1075 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1076 // calculate _nodal from _descending
1080 /*! If not yet done, calculate the nodal Connectivity
1081 and the reverse nodal Connectivity*/
1082 //---------------------------------------------------//
1083 void CONNECTIVITY::calculateReverseNodalConnectivity()
1084 //---------------------------------------------------//
1086 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1090 calculateNodalConnectivity();
1092 if(_reverseNodalConnectivity==NULL) {
1094 med_int node_number = 0;
1095 vector <vector <med_int> > reverse_connectivity;
1096 reverse_connectivity.resize(_numberOfNodes+1);
1098 // Treat all cells types
1100 for (med_int j = 0; j < _numberOfTypes; j++)
1102 // node number of the cell type
1103 node_number = _type[j].getNumberOfNodes();
1104 // treat all cells of a particular type
1105 for (med_int k = _count[j]; k < _count[j+1]; k++)
1106 // treat all nodes of the cell type
1107 for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1109 med_int global_node = _nodal->getIJ(k,local_node_number);
1110 reverse_connectivity[global_node].push_back(k);
1114 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1116 //calculate size of reverse_nodal_connectivity array
1117 med_int size_reverse_nodal_connectivity = 0;
1118 for (med_int i = 1; i < _numberOfNodes+1; i++)
1119 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1121 //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
1122 med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
1123 med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
1124 //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
1125 //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
1127 reverse_nodal_connectivity_index[0] = 1;
1128 for (med_int i = 1; i < _numberOfNodes+1; i++)
1130 med_int size = reverse_connectivity[i].size();
1131 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1132 for (med_int j = 0; j < size; j++)
1133 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1136 //_reverseNodalConnectivity = ReverseConnectivity;
1137 _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1138 reverse_nodal_connectivity_index,
1139 reverse_nodal_connectivity);
1140 delete [] reverse_nodal_connectivity_index;
1141 delete [] reverse_nodal_connectivity;
1145 /*! If not yet done, calculate the Descending Connectivity */
1146 //-------------------------------------------------//
1147 void CONNECTIVITY::calculateDescendingConnectivity()
1148 //-------------------------------------------------//
1150 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1153 if (_descending==NULL)
1156 MESSAGE(LOC<<"No connectivity found !");
1157 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1159 // calcul _descending from _nodal
1160 // we need CONNECTIVITY for constituent
1162 if (_constituent != NULL)
1163 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1164 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1166 if (_entityDimension == 3)
1167 _constituent = new CONNECTIVITY(MED_FACE);
1168 else if (_entityDimension == 2)
1169 _constituent = new CONNECTIVITY(MED_EDGE);
1171 MESSAGE(LOC<<"We are in 1D");
1175 _constituent->_typeConnectivity = MED_NODAL;
1176 _constituent->_numberOfNodes = _numberOfNodes;
1177 // foreach cells, we built array of constituent
1178 int DescendingSize = 0;
1179 for(int i=0; i<_numberOfTypes; i++)
1180 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1181 //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1182 //const int * descend_connectivity = _descending->getValue();
1183 int * descend_connectivity = new int[DescendingSize];
1184 for (int i=0; i<DescendingSize; i++)
1185 descend_connectivity[i]=0;
1186 //const int * descend_connectivity_index = _descending->getIndex();
1187 int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1188 descend_connectivity_index[0]=1;
1189 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1190 ConstituentsTypes[0]=MED_NONE;
1191 ConstituentsTypes[1]=MED_NONE;
1192 int * NumberOfConstituentsForeachType = new int[2];
1193 NumberOfConstituentsForeachType[0]=0;
1194 NumberOfConstituentsForeachType[1]=0;
1195 for(int i=0; i<_numberOfTypes; i++) {
1196 // initialize descend_connectivity_index array :
1197 int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1198 for (int j=_count[i];j<_count[i+1];j++){
1199 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1200 // compute number of constituent of all cell for each type
1201 for(int k=1;k<NumberOfConstituents+1;k++) {
1202 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1203 if (ConstituentsTypes[0]==MED_NONE) {
1204 ConstituentsTypes[0]=MEDType;
1205 NumberOfConstituentsForeachType[0]++;
1206 } else if (ConstituentsTypes[0]==MEDType)
1207 NumberOfConstituentsForeachType[0]++;
1208 else if (ConstituentsTypes[1]==MED_NONE) {
1209 ConstituentsTypes[1]=MEDType;
1210 NumberOfConstituentsForeachType[1]++;
1211 } else if (ConstituentsTypes[1]==MEDType)
1212 NumberOfConstituentsForeachType[1]++;
1214 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1219 // we could built _constituent !
1220 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1221 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1223 //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1225 // we use _constituent->_nodal
1226 //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1227 //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1228 int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1229 int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1230 ConstituentNodalConnectivityIndex[0]=1;
1232 _constituent->_entityDimension=ConstituentsTypes[0]/100;
1233 if (ConstituentsTypes[1]==MED_NONE)
1234 _constituent->_numberOfTypes = 1;
1236 _constituent->_numberOfTypes = 2;
1237 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1238 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1239 _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1240 _constituent->_count[0]=1;
1241 int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1242 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1243 for (int i=0; i<_constituent->_numberOfTypes;i++) {
1244 _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1245 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1246 CELLMODEL Type(ConstituentsTypes[i]);
1247 _constituent->_type[i]=Type;
1248 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1249 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1250 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1252 delete [] ConstituentsTypes;
1253 delete [] NumberOfConstituentsForeachType;
1255 // we need reverse nodal connectivity
1256 if (! _reverseNodalConnectivity)
1257 calculateReverseNodalConnectivity();
1258 const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1259 const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1261 // array to keep reverse descending connectivity
1262 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1264 int TotalNumberOfSubCell = 0;
1265 for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
1266 CELLMODEL Type = _type[i];
1267 // int NumberOfNodesPerCell = Type.getNumberOfNodes();
1268 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1269 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1270 for (int k=1; k<=NumberOfConstituentPerCell; k++)
1271 { // we loop on all sub cell of it
1272 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1273 // it is a new sub cell !
1274 // TotalNumberOfSubCell++;
1276 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1277 tmp_NumberOfConstituentsForeachType[0]++;
1278 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1280 tmp_NumberOfConstituentsForeachType[1]++;
1281 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1283 //we have maximum two types
1285 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1286 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1287 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1289 int * NodesLists = new int[NumberOfNodesPerConstituent];
1290 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1291 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1292 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1294 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1296 // all elements which contains first node of sub cell :
1297 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1298 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1299 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1301 if (NumberOfCellsInList > 0) { // we could have no element !
1302 int * CellsList = new int[NumberOfCellsInList];
1303 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1304 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1306 // foreach node in sub cell, we search elements which are in common
1307 // at the end, we must have only one !
1309 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1310 int NewNumberOfCellsInList = 0;
1311 int * NewCellsList = new int[NumberOfCellsInList];
1312 for (int l1=0; l1<NumberOfCellsInList; l1++)
1313 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
1314 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1315 // increasing order : CellsList[l1] are not in elements list of node l
1317 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1318 // we have found one
1319 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1320 NewNumberOfCellsInList++;
1324 NumberOfCellsInList = NewNumberOfCellsInList;
1326 delete [] CellsList;
1327 CellsList = NewCellsList;
1330 if (NumberOfCellsInList > 0) { // We have found some elements !
1331 int CellNumber = CellsList[0];
1332 //delete [] CellsList;
1333 if (NumberOfCellsInList>1)
1334 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1336 if (NumberOfCellsInList==1) {
1337 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1339 // we search sub cell number in this cell to not calculate it another time
1342 for (int l=0; l<_numberOfTypes; l++)
1343 if (CellNumber < _count[l+1]) {
1347 //int sub_cell_count2 = Type2.get_entities_count(1);
1348 //int nodes_cell_count2 = Type2.get_nodes_count();
1350 for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1352 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1353 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
1354 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1357 if (counter==Type.getConstituentType(1,k)%100) {
1358 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1365 MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1368 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1370 delete [] CellsList;
1373 delete [] NodesLists;
1377 // we adjust _constituent
1378 int NumberOfConstituent=0;
1379 int SizeOfConstituentNodal=0;
1380 for (int i=0;i<_constituent->_numberOfTypes; i++) {
1381 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1382 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1384 // we built new _nodal attribute in _constituent
1385 //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1386 //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1387 //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1388 int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1389 int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1390 ConstituentNodalIndex[0]=1;
1391 // we build _reverseDescendingConnectivity
1392 //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1393 //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1394 //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1395 int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1396 int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1397 reverseDescendingConnectivityIndex[0]=1;
1399 // first constituent type
1400 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1401 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1402 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1403 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1405 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1406 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1407 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1410 // second type if any
1411 if (_constituent->_numberOfTypes==2) {
1412 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1413 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1414 int offset2=offset*2;
1415 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1416 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1417 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1418 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1419 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1421 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1422 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1423 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1426 _constituent->_count[2]=NumberOfConstituent+1;
1427 // we correct _descending to adjust face number
1428 for(int j=0;j<DescendingSize;j++)
1429 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1430 descend_connectivity[j]-=offset;
1433 delete [] ConstituentNodalConnectivityIndex;
1434 delete [] ConstituentNodalConnectivity;
1436 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1438 descend_connectivity_index,
1439 descend_connectivity);
1440 delete [] descend_connectivity_index;
1441 delete [] descend_connectivity;
1442 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1443 2*NumberOfConstituent,
1444 reverseDescendingConnectivityIndex,
1445 reverseDescendingConnectivityValue);
1446 delete [] reverseDescendingConnectivityIndex;
1447 delete [] reverseDescendingConnectivityValue;
1449 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1450 delete [] tmp_NumberOfConstituentsForeachType;
1452 //delete _constituent->_nodal;
1453 _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1454 SizeOfConstituentNodal,
1455 ConstituentNodalIndex,
1456 ConstituentNodalValue);
1458 delete [] ConstituentNodalIndex;
1459 delete [] ConstituentNodalValue;
1461 delete [] ReverseDescendingConnectivityValue;
1466 //-----------------------------------------------------------------------------------------//
1467 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1468 //-----------------------------------------------------------------------------------------//
1470 // if (_entity==MED_CELL)
1471 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1472 // Connectivity : we could not do that with MED_CELL entity !");
1473 // if (myConnectivity->getEntity()!=MED_CELL)
1474 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1475 // Connectivity : we need MED_CELL connectivity !");
1479 /*! Not implemented yet */
1480 //--------------------------------------------------------------------//
1481 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1482 //--------------------------------------------------------------------//
1484 const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1487 MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1496 Returns the geometry of an element given by its entity type & its global number.
1498 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1500 //--------------------------------------------------------------------//
1501 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1502 //--------------------------------------------------------------------//
1504 const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1507 // The globalNumber must verify : 1 <= globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1508 if ( (globalNumber < 1) || (globalNumber > _count[_numberOfTypes]-1 ) )
1509 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1510 << 1 <<"| and <= |" << _count[_numberOfTypes]-1 ));
1512 if (_entity==Entity) {
1513 for(int i=1; i<=_numberOfTypes;i++)
1514 if (globalNumber<_count[i])
1515 return _geometricTypes[i-1];
1517 else if (_constituent!=NULL)
1518 return _constituent->getElementType(Entity,globalNumber);
1520 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1521 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1527 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1529 os << endl << "------------- Entity = ";
1544 case MED_ALL_ENTITIES:
1545 os << "MED_ALL_ENTITIES";
1551 os << " -------------\n\nMedConnectivity : ";
1552 switch (co._typeConnectivity)
1555 os << "MED_NODAL\n";
1557 case MED_DESCENDING:
1558 os << "MED_DESCENDING\n";
1563 os << "Entity dimension : " << co._entityDimension << endl;
1564 os << "Number of nodes : " << co._numberOfNodes << endl;
1565 os << "Number of types : " << co._numberOfTypes << endl;
1566 for (int i=0; i!=co._numberOfTypes ; ++i)
1567 os << " -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1568 << co._count[i+1]-co._count[i] << " elements" << endl;
1570 if (co._typeConnectivity == MED_NODAL )
1572 for (int i=0; i<co._numberOfTypes; i++)
1574 os << endl << co._type[i].getName() << " : " << endl;
1575 int numberofelements = co._count[i+1]-co._count[i];
1576 const med_int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1577 int numberofnodespercell = co._geometricTypes[i]%100;
1578 for (int j=0;j<numberofelements;j++)
1580 os << setw(6) << j+1 << " : ";
1581 for (int k=0;k<numberofnodespercell;k++)
1582 os << connectivity[j*numberofnodespercell+k]<<" ";
1587 else if (co._typeConnectivity == MED_DESCENDING)
1589 int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1590 const med_int *connectivity = co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1591 const int *connectivity_index = co.getConnectivityIndex( co._typeConnectivity, co._entity );
1593 for ( int j=0; j!=numberofelements; j++ )
1595 os << "Element " << j+1 << " : ";
1596 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1597 os << connectivity[k-1] << " ";
1602 if (co._constituent)
1603 os << endl << *co._constituent << endl;