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"
14 Default Constructor. /n
15 Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
16 //--------------------------------------------------------------//
17 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
18 //--------------------------------------------------------------//
20 _typeConnectivity(MED_NODAL),
22 _geometricTypes((medGeometryElement*)NULL),
23 _type((CELLMODEL*)NULL),
26 _nodal((MEDSKYLINEARRAY*)NULL),
27 _descending((MEDSKYLINEARRAY*)NULL),
28 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
29 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
30 _neighbourhood((MEDSKYLINEARRAY*)NULL),
31 _constituent((CONNECTIVITY*)NULL)
33 MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
38 Default for Entity is MED_CELL */
39 //------------------------------------------------------------------------------//
40 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
41 //------------------------------------------------------------------------------//
43 _typeConnectivity(MED_NODAL),
44 _numberOfTypes(numberOfTypes),
46 _nodal((MEDSKYLINEARRAY*)NULL),
47 _descending((MEDSKYLINEARRAY*)NULL),
48 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
49 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
50 _neighbourhood((MEDSKYLINEARRAY*)NULL),
51 _constituent((CONNECTIVITY*)NULL)
53 MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
54 _geometricTypes = new medGeometryElement[numberOfTypes];
55 _type = new CELLMODEL[numberOfTypes];
56 _count = new int[numberOfTypes+1];
62 //------------------------------------------------------//
63 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
64 //----------------------------------------------------//
66 _typeConnectivity (m._typeConnectivity),
67 _numberOfTypes (m._numberOfTypes),
68 _entityDimension (m._entityDimension),
69 _numberOfNodes (m._numberOfNodes)
71 if (m._geometricTypes != NULL)
73 _geometricTypes = new medGeometryElement[_numberOfTypes];
74 memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
77 _geometricTypes = (medGeometryElement *) NULL;
81 _type = new CELLMODEL[_numberOfTypes];
82 for (int i=0;i<_numberOfTypes;i++)
83 _type[i] = CELLMODEL(m._type[i]);
86 _type = (CELLMODEL *) NULL;
90 _count = new med_int[_numberOfTypes+1];
91 memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(med_int));
94 _count = (med_int *) NULL;
97 _nodal = new MEDSKYLINEARRAY(* m._nodal);
99 _nodal = (MEDSKYLINEARRAY *) NULL;
101 if (m._descending != NULL)
102 _descending = new MEDSKYLINEARRAY(* m._descending);
104 _descending = (MEDSKYLINEARRAY *) NULL;
106 if (m._reverseNodalConnectivity != NULL)
107 _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
109 _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
111 if (m._reverseDescendingConnectivity != NULL)
112 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
114 _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
116 if (m._neighbourhood != NULL)
117 _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
119 _neighbourhood = (MEDSKYLINEARRAY *) NULL;
121 if (m._constituent != NULL)
122 _constituent = new CONNECTIVITY(* m._constituent);
124 _constituent = (CONNECTIVITY *) NULL;
129 desallocates existing pointers */
130 //----------------------------//
131 CONNECTIVITY::~CONNECTIVITY()
132 //----------------------------//
134 MESSAGE("Destructeur de CONNECTIVITY()");
136 if (_geometricTypes != NULL)
137 delete [] _geometricTypes;
144 if (_descending != NULL)
146 if (_reverseNodalConnectivity != NULL)
147 delete _reverseNodalConnectivity;
148 if (_reverseDescendingConnectivity != NULL)
149 delete _reverseDescendingConnectivity;
150 if (_constituent != NULL)
155 set _constituent to Constituent
156 be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
157 throws an exception if Constituent = MED_CELL
160 //----------------------------------------------------------//
161 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
163 //----------------------------------------------------------//
165 medEntityMesh Entity = Constituent->getEntity();
166 if (Entity == MED_CELL)
167 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
169 if ((Entity == MED_EDGE)&(_entityDimension == 3))
171 if (_constituent == NULL)
172 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
173 _constituent->setConstituent(Constituent);
176 _constituent = Constituent;
179 /*! Duplicated Types array in CONNECTIVITY object. */
180 //--------------------------------------------------------------------//
181 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
182 const medEntityMesh Entity)
184 //--------------------------------------------------------------------//
186 if (Entity == _entity)
187 for (int i=0; i<_numberOfTypes; i++)
189 _geometricTypes[i] = Types[i];
190 _type[i] = CELLMODEL(_geometricTypes[i]);
194 if (_constituent == NULL)
195 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
196 _constituent->setGeometricTypes(Types,Entity);
201 //--------------------------------------------------------------------//
202 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
204 //--------------------------------------------------------------------//
206 if (Entity == _entity)
208 int * index = new int[Count[_numberOfTypes]];
211 for (int i=0; i<_numberOfTypes; i++) {
212 _count[i+1] = Count[i+1];
213 int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
214 for (int j=_count[i]; j<_count[i+1]; j++)
215 index[j] = index[j-1]+NumberOfNodesPerElement;
218 if (_nodal != NULL) delete _nodal;
219 _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
220 _nodal->setIndex(index);
225 if (_constituent == NULL)
226 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
227 _constituent->setCount(Count,Entity);
231 //--------------------------------------------------------------------//
232 void CONNECTIVITY::setNodal(const int * Connectivity,
233 const medEntityMesh Entity,
234 const medGeometryElement Type)
236 //--------------------------------------------------------------------//
238 if (_entity == Entity)
240 // find geometric type
242 for (int i=0; i<_numberOfTypes; i++)
243 if (_geometricTypes[i] == Type) {
245 int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
246 //_nodal->setI(i+1,Connectivity);
247 for( int j=_count[i];j<_count[i+1]; j++)
248 _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
251 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
254 if (_constituent == NULL)
255 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
256 _constituent->setNodal(Connectivity,Entity,Type);
261 //------------------------------------------------------------------------------------------//
262 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
263 //------------------------------------------------------------------------------------------//
265 MESSAGE("CONNECTIVITY::calculateConnectivity");
267 // a temporary limitation due to calculteDescendingConnectivity function !!!
268 if ((_entityDimension==3) & (Entity==MED_EDGE))
269 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
272 if (ConnectivityType==MED_NODAL)
273 calculateNodalConnectivity();
275 if (Entity==MED_CELL)
276 calculateDescendingConnectivity();
278 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
279 if (Entity!=_entity) {
280 calculateDescendingConnectivity();
281 if (_entityDimension == 2 || _entityDimension == 3)
282 _constituent->calculateConnectivity(ConnectivityType,Entity);
286 /*! Give, in full or no interlace mode (for nodal connectivity),
287 descending or nodal connectivity.
289 You must give a <medEntityMesh> (ie:MED_EDGE)
290 and a <medGeometryElement> (ie:MED_SEG3). */
292 //------------------------------------------------------------//
293 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
294 //------------------------------------------------------------//
296 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
299 int numberOfFamilies = myFamilies.size();
300 if (numberOfFamilies == 0 ) {
301 MESSAGE(LOC<<"No family");
304 // does we do an update ?
305 if ((_constituent != NULL)&(_descending != NULL)) {
306 MESSAGE(LOC<<"Constituent is already defined");
310 if ((_constituent != NULL)&(_descending == NULL)) {
311 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
312 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
316 // well we could go !
317 CONNECTIVITY * oldConstituent = _constituent;
319 // for(int i=0; i<numberOfFamilies; i++) {
320 // FAMILY * myFamily = myFamilies[i];
321 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
324 _constituent = (CONNECTIVITY *)NULL;
325 // for instance we must have nodal connectivity in constituent :
326 if (oldConstituent->_nodal == NULL)
328 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
329 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
331 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
332 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
333 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
335 calculateDescendingConnectivity();
337 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
338 const int * newConstituentValue = _constituent->_nodal->getValue();
339 const int * newConstituentIndex = _constituent->_nodal->getIndex();
341 const int * newReverseDescendingIndex =
342 _reverseDescendingConnectivity->getIndex();
344 const int * newDescendingIndex = _descending->getIndex();
345 // const int * newDescendingValue = _descending->getValue();
347 // loop on all family,
348 // for all constituent in family, we get it's old connectivity
349 // (with oldCconstituentValue and oldConstituentIndex)
350 // and search the constituent in newConstituentValue with class
353 // when a new face is found, replace old constituent
354 // number in family with new one
355 // If face is not rigth oriented, we must change _descending attribute
356 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
358 // Voila a toi de jouer Nadir :-)
360 // First we get the renumbering from the oldCconstituentValue and
361 // oldConstituentIndex in the the new one, newConstituentValue and
362 // newConstituentIndex with a negative sign if the face is not
365 int * renumberingFromOldToNew = new int [oldNumberOfFace];
369 _constituent->calculateReverseNodalConnectivity();
371 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
375 renumberingFromOldToNew[iOldFace] = iOldFace+1;
376 // renumberingFromOldToNew[iOldFace] = 999999;
378 int face_it_beginOld = oldConstituentIndex[iOldFace];
379 int face_it_endOld = oldConstituentIndex[iOldFace+1];
380 int face_size_itOld = face_it_endOld - face_it_beginOld;
382 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
385 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
386 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
388 // set an array wich contains faces numbers arround first node
389 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
390 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
391 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
393 int * FacesList = new int[NumberOfFacesInList];
395 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
396 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
398 // foreach node in sub cell, we search elements which are in common
399 // at the end, we must have only one !
401 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
403 int NewNumberOfFacesInList = 0;
404 int * NewFacesList = new int[NumberOfFacesInList];
406 for (int l1=0; l1<NumberOfFacesInList; l1++) {
407 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
408 if (FacesList[l1]<reverseFaceNodal[l2-1])
409 // increasing order : FacesList[l1] are not in elements list of node l
411 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
413 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
414 NewNumberOfFacesInList++;
419 NumberOfFacesInList = NewNumberOfFacesInList;
421 FacesList = NewFacesList;
424 if (!NumberOfFacesInList==0)
426 if (NumberOfFacesInList>1)
427 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
429 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
431 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
432 int face_it_endNew = newConstituentIndex[FacesList[0]];
433 face_size_itNew = face_it_endNew - face_it_beginNew;
435 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
436 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
438 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
440 //SCRUTE(retCompareNewOld);
442 // Real new face found
444 if(retCompareNewOld == 1)
446 renumberingFromOldToNew[iOldFace] = FacesList[0];
451 // Reverse new face found
453 if(retCompareNewOld == -1)
455 renumberingFromOldToNew[iOldFace] = FacesList[0];
459 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
460 int face_it_end = newReverseDescendingIndex[FacesList[0]];
461 int face_size_it = face_it_end - face_it_begin;
463 if (face_size_it == 1)
464 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
466 if (face_size_it > 2)
467 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
469 // we have always 2 neighbourings
470 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
471 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
472 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
474 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
476 if (cell2 != 0) { // we are not on border !!!!
478 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
479 // Updating _constituent->_nodal because of reversity
480 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
481 for(int iarray=1;iarray<=face_size_itNew;iarray++){
482 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
485 // Updating _reverseDescendingConnectivity
488 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
490 // Updating _descending for cell1 and cell2
491 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
492 if (_descending->getIndexValue(iface)==FacesList[0])
493 _descending->setIndexValue(iface,-FacesList[0]);
494 else if (_descending->getIndexValue(iface)==-FacesList[0])
495 _descending->setIndexValue(iface,FacesList[0]);
496 // else nothing to do
498 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
499 if (_descending->getIndexValue(iface)==FacesList[0])
500 _descending->setIndexValue(iface,-FacesList[0]);
501 else if (_descending->getIndexValue(iface)==-FacesList[0])
502 _descending->setIndexValue(iface,FacesList[0]);
503 // else nothing to do
505 } else {// else we are on border and we do nothing !!!!!!!!
506 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
507 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
508 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
514 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
515 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
521 MESSAGE(LOC<<"The Renumbering is finished and the status is");
523 // Updating the Family
525 for(int i=0; i<numberOfFamilies; i++) {
526 FAMILY * myFamily = myFamilies[i];
528 MEDSKYLINEARRAY * number = myFamily->getnumber();
529 int numberOfLines_skyline = number->getNumberOf();
530 const int * index_skyline = number->getIndex();
532 for (int i=0;i<numberOfLines_skyline;i++) {
533 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
534 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
537 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
540 delete oldConstituent ;
541 delete [] renumberingFromOldToNew;
549 // meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
550 // l'heritage car les pointeurs sont dans un conteneur.
551 void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
552 //------------------------------------------------------------//
554 const char * LOC = "CONNECTIVITY::updateFamily(vector<GROUP*>) ";
557 int numberOfFamilies = myFamilies.size();
558 if (numberOfFamilies == 0 ) {
559 MESSAGE(LOC<<"No family");
562 // does we do an update ?
563 if ((_constituent != NULL)&(_descending != NULL)) {
564 MESSAGE(LOC<<"Constituent is already defined");
568 if ((_constituent != NULL)&(_descending == NULL)) {
569 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
570 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
574 // well we could go !
575 CONNECTIVITY * oldConstituent = _constituent;
577 // for(int i=0; i<numberOfFamilies; i++) {
578 // FAMILY * myFamily = myFamilies[i];
579 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
582 _constituent = (CONNECTIVITY *)NULL;
583 // for instance we must have nodal connectivity in constituent :
584 if (oldConstituent->_nodal == NULL)
586 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
587 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
589 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
590 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
591 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
593 calculateDescendingConnectivity();
595 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
596 const int * newConstituentValue = _constituent->_nodal->getValue();
597 const int * newConstituentIndex = _constituent->_nodal->getIndex();
599 const int * newReverseDescendingIndex =
600 _reverseDescendingConnectivity->getIndex();
602 const int * newDescendingIndex = _descending->getIndex();
603 // const int * newDescendingValue = _descending->getValue();
605 // loop on all family,
606 // for all constituent in family, we get it's old connectivity
607 // (with oldCconstituentValue and oldConstituentIndex)
608 // and search the constituent in newConstituentValue with class
611 // when a new face is found, replace old constituent
612 // number in family with new one
613 // If face is not rigth oriented, we must change _descending attribute
614 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
616 // Voila a toi de jouer Nadir :-)
618 // First we get the renumbering from the oldCconstituentValue and
619 // oldConstituentIndex in the the new one, newConstituentValue and
620 // newConstituentIndex with a negative sign if the face is not
623 int * renumberingFromOldToNew = new int [oldNumberOfFace];
627 _constituent->calculateReverseNodalConnectivity();
629 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
633 renumberingFromOldToNew[iOldFace] = iOldFace+1;
634 // renumberingFromOldToNew[iOldFace] = 999999;
636 int face_it_beginOld = oldConstituentIndex[iOldFace];
637 int face_it_endOld = oldConstituentIndex[iOldFace+1];
638 int face_size_itOld = face_it_endOld - face_it_beginOld;
640 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
643 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
644 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
646 // set an array wich contains faces numbers arround first node
647 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
648 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
649 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
651 int * FacesList = new int[NumberOfFacesInList];
653 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
654 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
656 // foreach node in sub cell, we search elements which are in common
657 // at the end, we must have only one !
659 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
661 int NewNumberOfFacesInList = 0;
662 int * NewFacesList = new int[NumberOfFacesInList];
664 for (int l1=0; l1<NumberOfFacesInList; l1++) {
665 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
666 if (FacesList[l1]<reverseFaceNodal[l2-1])
667 // increasing order : FacesList[l1] are not in elements list of node l
669 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
671 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
672 NewNumberOfFacesInList++;
677 NumberOfFacesInList = NewNumberOfFacesInList;
679 FacesList = NewFacesList;
682 if (!NumberOfFacesInList==0)
684 if (NumberOfFacesInList>1)
685 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
687 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
689 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
690 int face_it_endNew = newConstituentIndex[FacesList[0]];
691 face_size_itNew = face_it_endNew - face_it_beginNew;
693 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
694 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
696 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
698 //SCRUTE(retCompareNewOld);
700 // Real new face found
702 if(retCompareNewOld == 1)
704 renumberingFromOldToNew[iOldFace] = FacesList[0];
709 // Reverse new face found
711 if(retCompareNewOld == -1)
713 renumberingFromOldToNew[iOldFace] = FacesList[0];
717 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
718 int face_it_end = newReverseDescendingIndex[FacesList[0]];
719 int face_size_it = face_it_end - face_it_begin;
721 if (face_size_it == 1)
722 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
724 if (face_size_it > 2)
725 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
727 // we have always 2 neighbourings
728 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
729 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
730 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
732 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
734 if (cell2 != 0) { // we are not on border !!!!
736 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
737 // Updating _constituent->_nodal because of reversity
738 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
739 for(int iarray=1;iarray<=face_size_itNew;iarray++){
740 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
743 // Updating _reverseDescendingConnectivity
746 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
748 // Updating _descending for cell1 and cell2
749 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
750 if (_descending->getIndexValue(iface)==FacesList[0])
751 _descending->setIndexValue(iface,-FacesList[0]);
752 else if (_descending->getIndexValue(iface)==-FacesList[0])
753 _descending->setIndexValue(iface,FacesList[0]);
754 // else nothing to do
756 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
757 if (_descending->getIndexValue(iface)==FacesList[0])
758 _descending->setIndexValue(iface,-FacesList[0]);
759 else if (_descending->getIndexValue(iface)==-FacesList[0])
760 _descending->setIndexValue(iface,FacesList[0]);
761 // else nothing to do
763 } else {// else we are on border and we do nothing !!!!!!!!
764 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
765 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
766 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
772 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
773 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
779 MESSAGE(LOC<<"The Renumbering is finished and the status is");
781 // Updating the Family
783 for(int i=0; i<numberOfFamilies; i++) {
784 GROUP * myFamily = myFamilies[i];
786 MEDSKYLINEARRAY * number = myFamily->getnumber();
787 int numberOfLines_skyline = number->getNumberOf();
788 const int * index_skyline = number->getIndex();
790 for (int i=0;i<numberOfLines_skyline;i++) {
791 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
792 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
795 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
798 delete oldConstituent ;
799 delete [] renumberingFromOldToNew;
807 //------------------------------------------------------------------------------------------------------------------//
808 const med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
809 //------------------------------------------------------------------------------------------------------------------//
811 const char * LOC = "CONNECTIVITY::getConnectivity";
814 MEDSKYLINEARRAY * Connectivity;
815 if (Entity==_entity) {
817 if (ConnectivityType==MED_NODAL)
819 calculateNodalConnectivity();
824 calculateDescendingConnectivity();
825 Connectivity=_descending;
828 if (Connectivity!=NULL)
829 if (Type==MED_ALL_ELEMENTS)
830 return Connectivity->getValue();
832 for (med_int i=0; i<_numberOfTypes; i++)
833 if (_geometricTypes[i]==Type)
834 //return Connectivity->getI(i+1);
835 return Connectivity->getI(_count[i]);
836 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
839 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
841 if (_constituent != NULL)
842 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
844 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
847 /*! Give morse index array to use with
848 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
850 Each value give start index for corresponding entity in connectivity array.
852 Example : i-th element, j-th node of it :
853 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
854 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
855 //-----------------------------------------------------------------------------------------------//
856 const med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
857 //----0000000--------------------------------------------------------------------------------------------//
859 const char * LOC = "CONNECTIVITY::getConnectivityIndex";
862 MEDSKYLINEARRAY * Connectivity;
863 if (Entity==_entity) {
865 if (ConnectivityType==MED_NODAL)
868 Connectivity=_descending;
870 if (Connectivity!=NULL)
871 return Connectivity->getIndex();
873 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
876 if (_constituent != NULL)
877 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
879 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
883 //--------------------------------------------------------------//
884 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
885 //--------------------------------------------------------------//
887 const char * LOC = "CONNECTIVITY::getType";
890 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
891 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
892 for (med_int i=0; i<_numberOfTypes; i++)
893 if (_geometricTypes[i]==Type)
895 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
898 /*! Returns the number of elements of type <med geometrie element>.
899 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
900 //------------------------------------------------------------------------//
901 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
902 //------------------------------------------------------------------------//
904 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
907 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
908 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
909 for (med_int i=0; i<_numberOfTypes; i++)
910 if (_geometricTypes[i]==Type)
911 return _type[i].getNumberOfNodes();
912 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
915 /*! Returns the number of geometric sub cells of <med geometrie element> type.
916 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
917 //------------------------------------------------------------------------//
918 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
919 //------------------------------------------------------------------------//
921 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
922 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
923 for (med_int i=0; i<_numberOfTypes; i++)
924 if (_geometricTypes[i]==Type)
925 return _type[i].getNumberOfConstituents(1);
926 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
929 /*! Returns the number of elements of type <med geometrie element>.
932 - Implemented for MED_ALL_ELEMENTS
933 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
934 - Not implemented for MED_NONE */
935 //-----------------------------------------------------------------------------------//
936 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
937 //-----------------------------------------------------------------------------------//
939 const char * LOC = "CONNECTIVITY::getNumberOf";
942 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
944 if (Entity==_entity) {
946 return 0; // not defined !
947 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
948 if (Type==MED_ALL_ELEMENTS)
949 return _count[_numberOfTypes]-1;
950 for (med_int i=0; i<_numberOfTypes; i++)
951 if (_geometricTypes[i]==Type)
952 return (_count[i+1] - _count[i]);
953 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
955 if (_constituent != NULL)
956 return _constituent->getNumberOf(Entity,Type);
958 return 0; // valid if they are nothing else !
959 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
963 //--------------------------------------------------------------//
964 const med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
965 medGeometryElement Type)
966 //--------------------------------------------------------------//
968 if (TypeConnectivity == MED_NODAL)
970 calculateNodalConnectivity();
971 if (Type==MED_ALL_ELEMENTS)
972 return _nodal->getValue();
973 for (med_int i=0; i<_numberOfTypes; i++)
974 if (_geometricTypes[i]==Type)
975 return _nodal->getI(_count[i]);
979 calculateDescendingConnectivity();
980 if (Type==MED_ALL_ELEMENTS)
981 return _descending->getValue();
982 for (med_int i=0; i<_numberOfTypes; i++)
983 if (_geometricTypes[i]==Type)
984 return _descending->getI(Type);
986 throw MEDEXCEPTION("Not found");
990 //---------------------------------------------------------------------//
991 const med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
992 //---------------------------------------------------------------------//
994 if (TypeConnectivity == MED_NODAL)
996 calculateNodalConnectivity();
997 return _nodal->getIndex();
1001 calculateDescendingConnectivity();
1002 return _descending->getIndex();
1006 /*! Not yet implemented */
1007 //----------------------------------------------//
1008 const med_int* CONNECTIVITY:: getNeighbourhood() const
1009 //----------------------------------------------//
1011 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1014 /*! Returns an array which contains, for each node, all cells
1016 //-------------------------------------------------//
1017 const med_int* CONNECTIVITY::getReverseNodalConnectivity()
1018 //-------------------------------------------------//
1020 calculateReverseNodalConnectivity();
1021 return _reverseNodalConnectivity->getValue();
1024 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
1025 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
1026 //-------------------------------------------------------//
1027 const med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1028 //-------------------------------------------------------//
1030 calculateReverseNodalConnectivity();
1031 return _reverseNodalConnectivity->getIndex();
1034 /*! Returns an array which contains, for each face (or edge),
1035 the 2 cells of each side. First is cell which face normal is outgoing.
1037 //------------------------------------------------------//
1038 const med_int* CONNECTIVITY::getReverseDescendingConnectivity()
1039 //------------------------------------------------------//
1041 // it is in _constituent connectivity only if we are in MED_CELL
1042 // (we could not for instance calculate face-edge connectivity !)
1043 if (_entity!=MED_CELL)
1044 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1046 // we need descending connectivity
1047 calculateDescendingConnectivity();
1048 return _reverseDescendingConnectivity->getValue();
1051 /*! calculate the reverse descending Connectivity
1052 and returns the index ( A DOCUMENTER MIEUX)*/
1053 //-----------------------------------------------------------//
1054 const med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1055 //-----------------------------------------------------------//
1057 // it is in _constituent connectivity only if we are in MED_CELL
1058 if (_entity!=MED_CELL)
1059 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1061 // we need descending connectivity
1062 calculateDescendingConnectivity();
1063 return _reverseDescendingConnectivity->getIndex();
1066 /*! A DOCUMENTER (et a finir ???) */
1067 //--------------------------------------------//
1068 void CONNECTIVITY::calculateNodalConnectivity()
1069 //--------------------------------------------//
1073 if (_descending==NULL)
1074 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1075 // calculate _nodal from _descending
1079 /*! If not yet done, calculate the nodal Connectivity
1080 and the reverse nodal Connectivity*/
1081 //---------------------------------------------------//
1082 void CONNECTIVITY::calculateReverseNodalConnectivity()
1083 //---------------------------------------------------//
1085 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1089 calculateNodalConnectivity();
1091 if(_reverseNodalConnectivity==NULL) {
1093 med_int node_number = 0;
1094 vector <vector <med_int> > reverse_connectivity;
1095 reverse_connectivity.resize(_numberOfNodes+1);
1097 // Treat all cells types
1099 for (med_int j = 0; j < _numberOfTypes; j++)
1101 // node number of the cell type
1102 node_number = _type[j].getNumberOfNodes();
1103 // treat all cells of a particular type
1104 for (med_int k = _count[j]; k < _count[j+1]; k++)
1105 // treat all nodes of the cell type
1106 for (med_int local_node_number = 1; local_node_number < node_number+1; local_node_number++)
1108 med_int global_node = _nodal->getIJ(k,local_node_number);
1109 reverse_connectivity[global_node].push_back(k);
1113 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1115 //calculate size of reverse_nodal_connectivity array
1116 med_int size_reverse_nodal_connectivity = 0;
1117 for (med_int i = 1; i < _numberOfNodes+1; i++)
1118 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1120 //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
1121 med_int * reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
1122 med_int * reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
1123 //const med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
1124 //const med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
1126 reverse_nodal_connectivity_index[0] = 1;
1127 for (med_int i = 1; i < _numberOfNodes+1; i++)
1129 med_int size = reverse_connectivity[i].size();
1130 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
1131 for (med_int j = 0; j < size; j++)
1132 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1135 //_reverseNodalConnectivity = ReverseConnectivity;
1136 _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1137 reverse_nodal_connectivity_index,
1138 reverse_nodal_connectivity);
1139 delete [] reverse_nodal_connectivity_index;
1140 delete [] reverse_nodal_connectivity;
1144 /*! If not yet done, calculate the Descending Connectivity */
1145 //-------------------------------------------------//
1146 void CONNECTIVITY::calculateDescendingConnectivity()
1147 //-------------------------------------------------//
1149 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1152 if (_descending==NULL)
1155 MESSAGE(LOC<<"No connectivity found !");
1156 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1158 // calcul _descending from _nodal
1159 // we need CONNECTIVITY for constituent
1161 if (_constituent != NULL)
1162 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1163 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1165 if (_entityDimension == 3)
1166 _constituent = new CONNECTIVITY(MED_FACE);
1167 else if (_entityDimension == 2)
1168 _constituent = new CONNECTIVITY(MED_EDGE);
1170 MESSAGE(LOC<<"We are in 1D");
1174 _constituent->_typeConnectivity = MED_NODAL;
1175 _constituent->_numberOfNodes = _numberOfNodes;
1176 // foreach cells, we built array of constituent
1177 int DescendingSize = 0;
1178 for(int i=0; i<_numberOfTypes; i++)
1179 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1180 //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1181 //const int * descend_connectivity = _descending->getValue();
1182 int * descend_connectivity = new int[DescendingSize];
1183 for (int i=0; i<DescendingSize; i++)
1184 descend_connectivity[i]=0;
1185 //const int * descend_connectivity_index = _descending->getIndex();
1186 int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1187 descend_connectivity_index[0]=1;
1188 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1189 ConstituentsTypes[0]=MED_NONE;
1190 ConstituentsTypes[1]=MED_NONE;
1191 int * NumberOfConstituentsForeachType = new int[2];
1192 NumberOfConstituentsForeachType[0]=0;
1193 NumberOfConstituentsForeachType[1]=0;
1194 for(int i=0; i<_numberOfTypes; i++) {
1195 // initialize descend_connectivity_index array :
1196 int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1197 for (int j=_count[i];j<_count[i+1];j++){
1198 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1199 // compute number of constituent of all cell for each type
1200 for(int k=1;k<NumberOfConstituents+1;k++) {
1201 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1202 if (ConstituentsTypes[0]==MED_NONE) {
1203 ConstituentsTypes[0]=MEDType;
1204 NumberOfConstituentsForeachType[0]++;
1205 } else if (ConstituentsTypes[0]==MEDType)
1206 NumberOfConstituentsForeachType[0]++;
1207 else if (ConstituentsTypes[1]==MED_NONE) {
1208 ConstituentsTypes[1]=MEDType;
1209 NumberOfConstituentsForeachType[1]++;
1210 } else if (ConstituentsTypes[1]==MEDType)
1211 NumberOfConstituentsForeachType[1]++;
1213 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1218 // we could built _constituent !
1219 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1220 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1222 //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1224 // we use _constituent->_nodal
1225 //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1226 //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1227 int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1228 int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1229 ConstituentNodalConnectivityIndex[0]=1;
1231 _constituent->_entityDimension=ConstituentsTypes[0]/100;
1232 if (ConstituentsTypes[1]==MED_NONE)
1233 _constituent->_numberOfTypes = 1;
1235 _constituent->_numberOfTypes = 2;
1236 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1237 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1238 _constituent->_count = new med_int[_constituent->_numberOfTypes+1];
1239 _constituent->_count[0]=1;
1240 int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1];
1241 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1242 for (int i=0; i<_constituent->_numberOfTypes;i++) {
1243 _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1244 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1245 CELLMODEL Type(ConstituentsTypes[i]);
1246 _constituent->_type[i]=Type;
1247 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1248 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1249 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1251 delete [] ConstituentsTypes;
1252 delete [] NumberOfConstituentsForeachType;
1254 // we need reverse nodal connectivity
1255 if (! _reverseNodalConnectivity)
1256 calculateReverseNodalConnectivity();
1257 const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1258 const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1260 // array to keep reverse descending connectivity
1261 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1263 int TotalNumberOfSubCell = 0;
1264 for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
1265 CELLMODEL Type = _type[i];
1266 // int NumberOfNodesPerCell = Type.getNumberOfNodes();
1267 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1268 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1269 for (int k=1; k<=NumberOfConstituentPerCell; k++)
1270 { // we loop on all sub cell of it
1271 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
1272 // it is a new sub cell !
1273 // TotalNumberOfSubCell++;
1275 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
1276 tmp_NumberOfConstituentsForeachType[0]++;
1277 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1279 tmp_NumberOfConstituentsForeachType[1]++;
1280 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1282 //we have maximum two types
1284 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1285 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1286 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1288 int * NodesLists = new int[NumberOfNodesPerConstituent];
1289 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
1290 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1291 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1293 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1295 // all elements which contains first node of sub cell :
1296 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1297 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]];
1298 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1300 if (NumberOfCellsInList > 0) { // we could have no element !
1301 int * CellsList = new int[NumberOfCellsInList];
1302 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1303 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1305 // foreach node in sub cell, we search elements which are in common
1306 // at the end, we must have only one !
1308 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
1309 int NewNumberOfCellsInList = 0;
1310 int * NewCellsList = new int[NumberOfCellsInList];
1311 for (int l1=0; l1<NumberOfCellsInList; l1++)
1312 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
1313 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1314 // increasing order : CellsList[l1] are not in elements list of node l
1316 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1317 // we have found one
1318 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1319 NewNumberOfCellsInList++;
1323 NumberOfCellsInList = NewNumberOfCellsInList;
1325 delete [] CellsList;
1326 CellsList = NewCellsList;
1329 if (NumberOfCellsInList > 0) { // We have found some elements !
1330 int CellNumber = CellsList[0];
1331 //delete [] CellsList;
1332 if (NumberOfCellsInList>1)
1333 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1335 if (NumberOfCellsInList==1) {
1336 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1338 // we search sub cell number in this cell to not calculate it another time
1341 for (int l=0; l<_numberOfTypes; l++)
1342 if (CellNumber < _count[l+1]) {
1346 //int sub_cell_count2 = Type2.get_entities_count(1);
1347 //int nodes_cell_count2 = Type2.get_nodes_count();
1349 for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) { // on all sub cell
1351 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1352 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
1353 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1356 if (counter==Type.getConstituentType(1,k)%100) {
1357 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1364 MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1367 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1369 delete [] CellsList;
1372 delete [] NodesLists;
1376 // we adjust _constituent
1377 int NumberOfConstituent=0;
1378 int SizeOfConstituentNodal=0;
1379 for (int i=0;i<_constituent->_numberOfTypes; i++) {
1380 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1381 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1383 // we built new _nodal attribute in _constituent
1384 //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1385 //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1386 //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1387 int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1388 int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1389 ConstituentNodalIndex[0]=1;
1390 // we build _reverseDescendingConnectivity
1391 //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1392 //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1393 //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1394 int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1395 int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1396 reverseDescendingConnectivityIndex[0]=1;
1398 // first constituent type
1399 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
1400 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1401 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1402 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1404 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1405 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1406 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1409 // second type if any
1410 if (_constituent->_numberOfTypes==2) {
1411 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1412 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1413 int offset2=offset*2;
1414 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1415 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1416 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1417 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1418 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1420 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1421 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1422 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1425 _constituent->_count[2]=NumberOfConstituent+1;
1426 // we correct _descending to adjust face number
1427 for(int j=0;j<DescendingSize;j++)
1428 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1429 descend_connectivity[j]-=offset;
1432 delete [] ConstituentNodalConnectivityIndex;
1433 delete [] ConstituentNodalConnectivity;
1435 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1437 descend_connectivity_index,
1438 descend_connectivity);
1439 delete [] descend_connectivity_index;
1440 delete [] descend_connectivity;
1441 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1442 2*NumberOfConstituent,
1443 reverseDescendingConnectivityIndex,
1444 reverseDescendingConnectivityValue);
1445 delete [] reverseDescendingConnectivityIndex;
1446 delete [] reverseDescendingConnectivityValue;
1448 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1449 delete [] tmp_NumberOfConstituentsForeachType;
1451 //delete _constituent->_nodal;
1452 _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1453 SizeOfConstituentNodal,
1454 ConstituentNodalIndex,
1455 ConstituentNodalValue);
1457 delete [] ConstituentNodalIndex;
1458 delete [] ConstituentNodalValue;
1460 delete [] ReverseDescendingConnectivityValue;
1465 //-----------------------------------------------------------------------------------------//
1466 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1467 //-----------------------------------------------------------------------------------------//
1469 // if (_entity==MED_CELL)
1470 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1471 // Connectivity : we could not do that with MED_CELL entity !");
1472 // if (myConnectivity->getEntity()!=MED_CELL)
1473 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1474 // Connectivity : we need MED_CELL connectivity !");
1478 /*! Not implemented yet */
1479 //--------------------------------------------------------------------//
1480 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1481 //--------------------------------------------------------------------//
1483 const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1486 MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1495 Returns the geometry of an element given by its entity type & its global number.
1497 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1499 //--------------------------------------------------------------------//
1500 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1501 //--------------------------------------------------------------------//
1503 const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1506 // The globalNumber must verify : 1 <= globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1507 if ( (globalNumber < 1) || (globalNumber > _count[_numberOfTypes]-1 ) )
1508 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1509 << 1 <<"| and <= |" << _count[_numberOfTypes]-1 ));
1511 if (_entity==Entity) {
1512 for(int i=1; i<=_numberOfTypes;i++)
1513 if (globalNumber<_count[i])
1514 return _geometricTypes[i-1];
1516 else if (_constituent!=NULL)
1517 return _constituent->getElementType(Entity,globalNumber);
1519 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1520 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1526 ostream & operator<<(ostream &os, CONNECTIVITY &co)
1528 os << endl << "------------- Entity = ";
1543 case MED_ALL_ENTITIES:
1544 os << "MED_ALL_ENTITIES";
1550 os << " -------------\n\nMedConnectivity : ";
1551 switch (co._typeConnectivity)
1554 os << "MED_NODAL\n";
1556 case MED_DESCENDING:
1557 os << "MED_DESCENDING\n";
1562 os << "Entity dimension : " << co._entityDimension << endl;
1563 os << "Number of nodes : " << co._numberOfNodes << endl;
1564 os << "Number of types : " << co._numberOfTypes << endl;
1565 for (int i=0; i!=co._numberOfTypes ; ++i)
1566 os << " -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1567 << co._count[i+1]-co._count[i] << " elements" << endl;
1569 if (co._typeConnectivity == MED_NODAL )
1571 for (int i=0; i<co._numberOfTypes; i++)
1573 os << endl << co._type[i].getName() << " : " << endl;
1574 int numberofelements = co._count[i+1]-co._count[i];
1575 const med_int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1576 int numberofnodespercell = co._geometricTypes[i]%100;
1577 for (int j=0;j<numberofelements;j++)
1579 os << setw(6) << j+1 << " : ";
1580 for (int k=0;k<numberofnodespercell;k++)
1581 os << connectivity[j*numberofnodespercell+k]<<" ";
1586 else if (co._typeConnectivity == MED_DESCENDING)
1588 int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1589 const med_int *connectivity = co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1590 const int *connectivity_index = co.getConnectivityIndex( co._typeConnectivity, co._entity );
1592 for ( int j=0; j!=numberofelements; j++ )
1594 os << "Element " << j+1 << " : ";
1595 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1596 os << connectivity[k-1] << " ";
1601 if (co._constituent)
1602 os << endl << *co._constituent << endl;