1 #include "MEDMEM_Connectivity.hxx"
2 #include "MEDMEM_Family.hxx"
3 #include "MEDMEM_Group.hxx"
4 #include "MEDMEM_CellModel.hxx"
6 #include "MEDMEM_SkyLineArray.hxx"
7 #include "MEDMEM_ModulusArray.hxx"
9 #include "MEDMEM_STRING.hxx"
13 using namespace MEDMEM;
14 using namespace MED_EN;
16 Default Constructor. \n
17 Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
18 //--------------------------------------------------------------//
19 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
20 //--------------------------------------------------------------//
22 _typeConnectivity(MED_NODAL),
24 _geometricTypes((medGeometryElement*)NULL),
25 _type((CELLMODEL*)NULL),
28 _nodal((MEDSKYLINEARRAY*)NULL),
29 _descending((MEDSKYLINEARRAY*)NULL),
30 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
31 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
32 _neighbourhood((MEDSKYLINEARRAY*)NULL),
33 _constituent((CONNECTIVITY*)NULL)
35 MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
40 Default for Entity is MED_CELL */
41 //------------------------------------------------------------------------------//
42 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
43 //------------------------------------------------------------------------------//
45 _typeConnectivity(MED_NODAL),
46 _numberOfTypes(numberOfTypes),
48 _nodal((MEDSKYLINEARRAY*)NULL),
49 _descending((MEDSKYLINEARRAY*)NULL),
50 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
51 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
52 _neighbourhood((MEDSKYLINEARRAY*)NULL),
53 _constituent((CONNECTIVITY*)NULL)
55 MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
56 _geometricTypes = new medGeometryElement[numberOfTypes];
57 _type = new CELLMODEL[numberOfTypes];
58 _count = new int[numberOfTypes+1];
64 //------------------------------------------------------//
65 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
66 //----------------------------------------------------//
68 _typeConnectivity (m._typeConnectivity),
69 _numberOfTypes (m._numberOfTypes),
70 _entityDimension (m._entityDimension),
71 _numberOfNodes (m._numberOfNodes)
73 if (m._geometricTypes != NULL)
75 _geometricTypes = new medGeometryElement[_numberOfTypes];
76 memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
79 _geometricTypes = (medGeometryElement *) NULL;
83 _type = new CELLMODEL[_numberOfTypes];
84 for (int i=0;i<_numberOfTypes;i++)
85 _type[i] = CELLMODEL(m._type[i]);
88 _type = (CELLMODEL *) NULL;
92 _count = new int[_numberOfTypes+1];
93 memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
96 _count = (int *) NULL;
99 _nodal = new MEDSKYLINEARRAY(* m._nodal);
101 _nodal = (MEDSKYLINEARRAY *) NULL;
103 if (m._descending != NULL)
104 _descending = new MEDSKYLINEARRAY(* m._descending);
106 _descending = (MEDSKYLINEARRAY *) NULL;
108 if (m._reverseNodalConnectivity != NULL)
109 _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
111 _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
113 if (m._reverseDescendingConnectivity != NULL)
114 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
116 _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
118 if (m._neighbourhood != NULL)
119 _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
121 _neighbourhood = (MEDSKYLINEARRAY *) NULL;
123 if (m._constituent != NULL)
124 _constituent = new CONNECTIVITY(* m._constituent);
126 _constituent = (CONNECTIVITY *) NULL;
131 desallocates existing pointers */
132 //----------------------------//
133 CONNECTIVITY::~CONNECTIVITY()
134 //----------------------------//
136 MESSAGE("Destructeur de CONNECTIVITY()");
138 if (_geometricTypes != NULL)
139 delete [] _geometricTypes;
146 if (_descending != NULL)
148 if (_reverseNodalConnectivity != NULL)
149 delete _reverseNodalConnectivity;
150 if (_reverseDescendingConnectivity != NULL)
151 delete _reverseDescendingConnectivity;
152 if (_constituent != NULL)
157 set _constituent to Constituent
158 be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
159 throws an exception if Constituent = MED_CELL
162 //----------------------------------------------------------//
163 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
165 //----------------------------------------------------------//
167 medEntityMesh Entity = Constituent->getEntity();
168 if (Entity == MED_CELL)
169 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
171 if ((Entity == MED_EDGE)&(_entityDimension == 3))
173 if (_constituent == NULL)
174 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
175 _constituent->setConstituent(Constituent);
178 _constituent = Constituent;
181 /*! Duplicated Types array in CONNECTIVITY object. */
182 //--------------------------------------------------------------------//
183 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
184 const medEntityMesh Entity)
186 //--------------------------------------------------------------------//
188 if (Entity == _entity)
189 for (int i=0; i<_numberOfTypes; i++)
191 _geometricTypes[i] = Types[i];
192 _type[i] = CELLMODEL(_geometricTypes[i]);
196 if (_constituent == NULL)
197 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
198 _constituent->setGeometricTypes(Types,Entity);
203 //--------------------------------------------------------------------//
204 void CONNECTIVITY::setCount(const int * Count, const medEntityMesh Entity)
206 //--------------------------------------------------------------------//
208 if (Entity == _entity)
210 int * index = new int[Count[_numberOfTypes]];
213 for (int i=0; i<_numberOfTypes; i++) {
214 _count[i+1] = Count[i+1];
215 int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
216 for (int j=_count[i]; j<_count[i+1]; j++)
217 index[j] = index[j-1]+NumberOfNodesPerElement;
220 if (_nodal != NULL) delete _nodal;
221 _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
222 _nodal->setIndex(index);
227 if (_constituent == NULL)
228 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
229 _constituent->setCount(Count,Entity);
233 //--------------------------------------------------------------------//
234 void CONNECTIVITY::setNodal(const int * Connectivity,
235 const medEntityMesh Entity,
236 const medGeometryElement Type)
238 //--------------------------------------------------------------------//
240 if (_entity == Entity)
242 // find geometric type
244 for (int i=0; i<_numberOfTypes; i++)
245 if (_geometricTypes[i] == Type) {
247 int NumberOfNodePerElements = _type[i].getNumberOfNodes() ;
248 //_nodal->setI(i+1,Connectivity);
249 for( int j=_count[i];j<_count[i+1]; j++)
250 _nodal->setI(j,Connectivity+(j-_count[i])*NumberOfNodePerElements);
253 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
256 if (_constituent == NULL)
257 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
258 _constituent->setNodal(Connectivity,Entity,Type);
263 //------------------------------------------------------------------------------------------//
264 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
265 //------------------------------------------------------------------------------------------//
267 MESSAGE("CONNECTIVITY::calculateConnectivity");
269 // a temporary limitation due to calculteDescendingConnectivity function !!!
270 if ((_entityDimension==3) & (Entity==MED_EDGE))
271 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
274 if (ConnectivityType==MED_NODAL)
275 calculateNodalConnectivity();
277 if (Entity==MED_CELL)
278 calculateDescendingConnectivity();
280 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
281 if (Entity!=_entity) {
282 calculateDescendingConnectivity();
283 if (_entityDimension == 2 || _entityDimension == 3)
284 _constituent->calculateConnectivity(ConnectivityType,Entity);
288 /*! Give, in full or no interlace mode (for nodal connectivity),
289 descending or nodal connectivity.
291 You must give a %medEntityMesh (ie:MED_EDGE)
292 and a %medGeometryElement (ie:MED_SEG3). */
294 //------------------------------------------------------------//
295 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
296 //------------------------------------------------------------//
298 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
301 int numberOfFamilies = myFamilies.size();
302 if (numberOfFamilies == 0 ) {
303 MESSAGE(LOC<<"No family");
306 // does we do an update ?
307 if ((_constituent != NULL)&(_descending != NULL)) {
308 MESSAGE(LOC<<"Constituent is already defined");
312 if ((_constituent != NULL)&(_descending == NULL)) {
313 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
314 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
318 // well we could go !
319 CONNECTIVITY * oldConstituent = _constituent;
321 // for(int i=0; i<numberOfFamilies; i++) {
322 // FAMILY * myFamily = myFamilies[i];
323 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
326 _constituent = (CONNECTIVITY *)NULL;
327 // for instance we must have nodal connectivity in constituent :
328 if (oldConstituent->_nodal == NULL)
330 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
331 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
333 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
334 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
335 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
337 calculateDescendingConnectivity();
339 MESSAGE(LOC << " Right after the call to calculateDescendingConnectivity");
340 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
341 const int * newConstituentValue = _constituent->_nodal->getValue();
342 const int * newConstituentIndex = _constituent->_nodal->getIndex();
344 const int * newReverseDescendingIndex =
345 _reverseDescendingConnectivity->getIndex();
347 const int * newDescendingIndex = _descending->getIndex();
348 // const int * newDescendingValue = _descending->getValue();
350 // loop on all family,
351 // for all constituent in family, we get it's old connectivity
352 // (with oldCconstituentValue and oldConstituentIndex)
353 // and search the constituent in newConstituentValue with class
356 // when a new face is found, replace old constituent
357 // number in family with new one
358 // If face is not rigth oriented, we must change _descending attribute
359 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
361 // Voila a toi de jouer Nadir :-)
363 // First we get the renumbering from the oldCconstituentValue and
364 // oldConstituentIndex in the the new one, newConstituentValue and
365 // newConstituentIndex with a negative sign if the face is not
368 int * renumberingFromOldToNew = new int [oldNumberOfFace];
372 MESSAGE(LOC << " Right before the call to _constituent->calculateReverseNodalConnectivity");
373 _constituent->calculateReverseNodalConnectivity();
374 MESSAGE(LOC << " Right after the call to _constituent->calculateReverseNodalConnectivity");
377 SCRUTE(oldNumberOfFace);
380 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
384 renumberingFromOldToNew[iOldFace] = iOldFace+1;
385 // renumberingFromOldToNew[iOldFace] = 999999;
387 int face_it_beginOld = oldConstituentIndex[iOldFace];
388 int face_it_endOld = oldConstituentIndex[iOldFace+1];
389 int face_size_itOld = face_it_endOld - face_it_beginOld;
391 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
394 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
395 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
397 // set an array wich contains faces numbers arround first node
398 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
399 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
400 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
402 int * FacesList = new int[NumberOfFacesInList];
404 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
405 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
407 // foreach node in sub cell, we search elements which are in common
408 // at the end, we must have only one !
410 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
412 int NewNumberOfFacesInList = 0;
413 int * NewFacesList = new int[NumberOfFacesInList];
415 for (int l1=0; l1<NumberOfFacesInList; l1++) {
416 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
417 if (FacesList[l1]<reverseFaceNodal[l2-1])
418 // increasing order : FacesList[l1] are not in elements list of node l
420 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
422 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
423 NewNumberOfFacesInList++;
428 NumberOfFacesInList = NewNumberOfFacesInList;
430 FacesList = NewFacesList;
433 if (!NumberOfFacesInList==0)
435 if (NumberOfFacesInList>1)
436 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
438 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
440 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
441 int face_it_endNew = newConstituentIndex[FacesList[0]];
442 face_size_itNew = face_it_endNew - face_it_beginNew;
444 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
445 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
447 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
449 //SCRUTE(retCompareNewOld);
451 // Real new face found
453 if(retCompareNewOld == 1)
455 renumberingFromOldToNew[iOldFace] = FacesList[0];
460 // Reverse new face found
462 if(retCompareNewOld == -1)
464 renumberingFromOldToNew[iOldFace] = FacesList[0];
468 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
469 int face_it_end = newReverseDescendingIndex[FacesList[0]];
470 int face_size_it = face_it_end - face_it_begin;
472 if (face_size_it == 1)
473 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
475 if (face_size_it > 2)
476 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
478 // we have always 2 neighbourings
479 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
480 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
481 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
483 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
485 if (cell2 != 0) { // we are not on border !!!!
487 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
488 // Updating _constituent->_nodal because of reversity
489 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
490 for(int iarray=1;iarray<=face_size_itNew;iarray++){
491 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
494 // Updating _reverseDescendingConnectivity
497 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
499 // Updating _descending for cell1 and cell2
500 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
501 if (_descending->getIndexValue(iface)==FacesList[0])
502 _descending->setIndexValue(iface,-FacesList[0]);
503 else if (_descending->getIndexValue(iface)==-FacesList[0])
504 _descending->setIndexValue(iface,FacesList[0]);
505 // else nothing to do
507 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
508 if (_descending->getIndexValue(iface)==FacesList[0])
509 _descending->setIndexValue(iface,-FacesList[0]);
510 else if (_descending->getIndexValue(iface)==-FacesList[0])
511 _descending->setIndexValue(iface,FacesList[0]);
512 // else nothing to do
514 } else {// else we are on border and we do nothing !!!!!!!!
515 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
516 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
517 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
523 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
524 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
530 MESSAGE(LOC<<"The Renumbering is finished and the status is");
532 // Updating the Family
534 SCRUTE(numberOfFamilies);
536 for(int i=0; i<numberOfFamilies; i++) {
537 FAMILY * myFamily = myFamilies[i];
540 // MESSAGE(LOC<<(*myFamily));
542 if (myFamily->isOnAllElements()) {
544 // we must have more constituent ?
545 if (oldNumberOfFace==_constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS))
546 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a family which is already in all constituent !"));
547 myFamily->setAll(false);
549 int * values = new int[oldNumberOfFace] ;
550 for (int ind=0;ind<oldNumberOfFace;ind++)
553 int NumberOfTypes = myFamily->getNumberOfTypes();
554 const int * count = oldConstituent->getGlobalNumberingIndex(_constituent->getEntity());
555 int * index = new int[NumberOfTypes+1] ;
556 memcpy(index,count,(NumberOfTypes+1)*sizeof(int));
557 // build new number attribut
558 myFamily->setNumber(index,values);
561 // MESSAGE(LOC<<(*myFamily));
562 MEDSKYLINEARRAY * number = myFamily->getnumber();
564 int numberOfLines_skyline = number->getNumberOf();
567 SCRUTE(numberOfLines_skyline);
569 const int * index_skyline = number->getIndex();
571 SCRUTE(index_skyline);
573 for (int i=0;i<numberOfLines_skyline;i++) {
574 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
575 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
578 // MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
581 delete oldConstituent ;
582 delete [] renumberingFromOldToNew;
590 // meme methode que updateFamily, mais avec des groupes. Il n'est pas possible d'utiliser
591 // l'heritage car les pointeurs sont dans un conteneur.
592 void CONNECTIVITY::updateGroup(vector<GROUP*> myFamilies)
593 //------------------------------------------------------------//
595 const char * LOC = "CONNECTIVITY::updateGroup(vector<GROUP*>) ";
598 int numberOfFamilies = myFamilies.size();
599 if (numberOfFamilies == 0 ) {
600 MESSAGE(LOC<<"No family");
603 // does we do an update ?
604 if ((_constituent != NULL)&(_descending != NULL)) {
605 MESSAGE(LOC<<"Constituent is already defined");
609 if ((_constituent != NULL)&(_descending == NULL)) {
610 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
611 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing");
615 // well we could go !
616 CONNECTIVITY * oldConstituent = _constituent;
618 // for(int i=0; i<numberOfFamilies; i++) {
619 // FAMILY * myFamily = myFamilies[i];
620 // MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
623 _constituent = (CONNECTIVITY *)NULL;
624 // for instance we must have nodal connectivity in constituent :
625 if (oldConstituent->_nodal == NULL)
627 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
628 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
630 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
631 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
632 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
634 calculateDescendingConnectivity();
636 // int newNumberOfFace = _constituent->_nodal->getNumberOf();
637 const int * newConstituentValue = _constituent->_nodal->getValue();
638 const int * newConstituentIndex = _constituent->_nodal->getIndex();
640 const int * newReverseDescendingIndex =
641 _reverseDescendingConnectivity->getIndex();
643 const int * newDescendingIndex = _descending->getIndex();
644 // const int * newDescendingValue = _descending->getValue();
646 // loop on all family,
647 // for all constituent in family, we get it's old connectivity
648 // (with oldCconstituentValue and oldConstituentIndex)
649 // and search the constituent in newConstituentValue with class
652 // when a new face is found, replace old constituent
653 // number in family with new one
654 // If face is not rigth oriented, we must change _descending attribute
655 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
657 // Voila a toi de jouer Nadir :-)
659 // First we get the renumbering from the oldCconstituentValue and
660 // oldConstituentIndex in the the new one, newConstituentValue and
661 // newConstituentIndex with a negative sign if the face is not
664 int * renumberingFromOldToNew = new int [oldNumberOfFace];
668 _constituent->calculateReverseNodalConnectivity();
670 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
674 renumberingFromOldToNew[iOldFace] = iOldFace+1;
675 // renumberingFromOldToNew[iOldFace] = 999999;
677 int face_it_beginOld = oldConstituentIndex[iOldFace];
678 int face_it_endOld = oldConstituentIndex[iOldFace+1];
679 int face_size_itOld = face_it_endOld - face_it_beginOld;
681 const int* NodesLists = oldConstituentValue + (face_it_beginOld-1);
684 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
685 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
687 // set an array wich contains faces numbers arround first node
688 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1];
689 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]];
690 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
692 int * FacesList = new int[NumberOfFacesInList];
694 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++){
695 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1];
697 // foreach node in sub cell, we search elements which are in common
698 // at the end, we must have only one !
700 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
702 int NewNumberOfFacesInList = 0;
703 int * NewFacesList = new int[NumberOfFacesInList];
705 for (int l1=0; l1<NumberOfFacesInList; l1++) {
706 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
707 if (FacesList[l1]<reverseFaceNodal[l2-1])
708 // increasing order : FacesList[l1] are not in elements list of node l
710 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
712 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
713 NewNumberOfFacesInList++;
718 NumberOfFacesInList = NewNumberOfFacesInList;
720 FacesList = NewFacesList;
723 if (!NumberOfFacesInList==0)
725 if (NumberOfFacesInList>1)
726 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
728 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
730 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
731 int face_it_endNew = newConstituentIndex[FacesList[0]];
732 face_size_itNew = face_it_endNew - face_it_beginNew;
734 const int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1;
735 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newNodesLists);
737 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
739 //SCRUTE(retCompareNewOld);
741 // Real new face found
743 if(retCompareNewOld == 1)
745 renumberingFromOldToNew[iOldFace] = FacesList[0];
750 // Reverse new face found
752 if(retCompareNewOld == -1)
754 renumberingFromOldToNew[iOldFace] = FacesList[0];
758 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
759 int face_it_end = newReverseDescendingIndex[FacesList[0]];
760 int face_size_it = face_it_end - face_it_begin;
762 if (face_size_it == 1)
763 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
765 if (face_size_it > 2)
766 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
768 // we have always 2 neighbourings
769 int cell1 = _reverseDescendingConnectivity->getIJ(FacesList[0],1);
770 int cell2 = _reverseDescendingConnectivity->getIJ(FacesList[0],2);
771 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
773 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
775 if (cell2 != 0) { // we are not on border !!!!
777 _reverseDescendingConnectivity->setIJ(FacesList[0],1,cell2);
778 // Updating _constituent->_nodal because of reversity
779 const int * oldArray = oldConstituentValue+face_it_beginOld-1;
780 for(int iarray=1;iarray<=face_size_itNew;iarray++){
781 _constituent->_nodal->setIJ(FacesList[0],iarray,oldArray[iarray-1]);
784 // Updating _reverseDescendingConnectivity
787 _reverseDescendingConnectivity->setIJ(FacesList[0],2,cell1);
789 // Updating _descending for cell1 and cell2
790 for(int iface=newDescendingIndex[cell1-1];iface<=newDescendingIndex[cell1];iface++)
791 if (_descending->getIndexValue(iface)==FacesList[0])
792 _descending->setIndexValue(iface,-FacesList[0]);
793 else if (_descending->getIndexValue(iface)==-FacesList[0])
794 _descending->setIndexValue(iface,FacesList[0]);
795 // else nothing to do
797 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
798 if (_descending->getIndexValue(iface)==FacesList[0])
799 _descending->setIndexValue(iface,-FacesList[0]);
800 else if (_descending->getIndexValue(iface)==-FacesList[0])
801 _descending->setIndexValue(iface,FacesList[0]);
802 // else nothing to do
804 } else {// else we are on border and we do nothing !!!!!!!!
805 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
806 MESSAGE(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !");
807 MESSAGE("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
813 MESSAGE(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
814 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
820 MESSAGE(LOC<<"The Renumbering is finished and the status is");
822 // Updating the Family
824 for(int i=0; i<numberOfFamilies; i++) {
825 GROUP * myFamily = myFamilies[i];
827 MEDSKYLINEARRAY * number = myFamily->getnumber();
828 int numberOfLines_skyline = number->getNumberOf();
829 const int * index_skyline = number->getIndex();
831 for (int i=0;i<numberOfLines_skyline;i++) {
832 for (int j=index_skyline[i]; j<index_skyline[i+1];j++) {
833 number->setIndexValue(j,renumberingFromOldToNew[number->getIndexValue(j)-1]);
836 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
839 delete oldConstituent ;
840 delete [] renumberingFromOldToNew;
848 //------------------------------------------------------------------------------------------------------------------//
849 const int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
850 //------------------------------------------------------------------------------------------------------------------//
852 const char * LOC = "CONNECTIVITY::getConnectivity";
855 MEDSKYLINEARRAY * Connectivity;
856 if (Entity==_entity) {
858 if (ConnectivityType==MED_NODAL)
860 calculateNodalConnectivity();
865 calculateDescendingConnectivity();
866 Connectivity=_descending;
869 if (Connectivity!=NULL)
870 if (Type==MED_ALL_ELEMENTS)
871 return Connectivity->getValue();
873 for (int i=0; i<_numberOfTypes; i++)
874 if (_geometricTypes[i]==Type)
875 //return Connectivity->getI(i+1);
876 return Connectivity->getI(_count[i]);
877 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
880 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
882 if (_constituent != NULL)
883 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
885 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
888 /*! Give morse index array to use with
889 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
891 Each value give start index for corresponding entity in connectivity array.
893 Example : i-th element, j-th node of it :
894 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
895 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
896 //-----------------------------------------------------------------------------------------------//
897 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
898 //----0000000--------------------------------------------------------------------------------------------//
900 const char * LOC = "CONNECTIVITY::getConnectivityIndex";
903 MEDSKYLINEARRAY * Connectivity;
904 if (Entity==_entity) {
906 if (ConnectivityType==MED_NODAL)
909 Connectivity=_descending;
911 if (Connectivity!=NULL)
912 return Connectivity->getIndex();
914 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
917 if (_constituent != NULL)
918 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
920 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
924 //--------------------------------------------------------------//
925 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
926 //--------------------------------------------------------------//
928 const char * LOC = "CONNECTIVITY::getType";
931 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
932 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !"));
933 for (int i=0; i<_numberOfTypes; i++)
934 if (_geometricTypes[i]==Type)
936 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
939 /*! Returns the number of elements of type %medGeometryElement.
940 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
941 //------------------------------------------------------------------------//
942 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
943 //------------------------------------------------------------------------//
945 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
948 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
949 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
950 for (int i=0; i<_numberOfTypes; i++)
951 if (_geometricTypes[i]==Type)
952 return _type[i].getNumberOfNodes();
953 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
956 /*! Returns the number of geometric sub cells of %medGeometryElement type.
957 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
958 //------------------------------------------------------------------------//
959 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
960 //------------------------------------------------------------------------//
962 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
963 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
964 for (int i=0; i<_numberOfTypes; i++)
965 if (_geometricTypes[i]==Type)
966 return _type[i].getNumberOfConstituents(1);
967 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
970 /*! Returns the number of elements of type %medGeometryElement.
973 - Implemented for MED_ALL_ELEMENTS
974 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
975 - Not implemented for MED_NONE */
976 //-----------------------------------------------------------------------------------//
977 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
978 //-----------------------------------------------------------------------------------//
980 const char * LOC = "CONNECTIVITY::getNumberOf";
983 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
985 if (Entity==_entity) {
987 return 0; // not defined !
988 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
989 if (Type==MED_ALL_ELEMENTS)
990 return _count[_numberOfTypes]-1;
991 for (int i=0; i<_numberOfTypes; i++)
992 if (_geometricTypes[i]==Type)
993 return (_count[i+1] - _count[i]);
994 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
996 if (_constituent != NULL)
997 return _constituent->getNumberOf(Entity,Type);
999 return 0; // valid if they are nothing else !
1000 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
1004 //--------------------------------------------------------------//
1005 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
1006 medGeometryElement Type)
1007 //--------------------------------------------------------------//
1009 if (TypeConnectivity == MED_NODAL)
1011 calculateNodalConnectivity();
1012 if (Type==MED_ALL_ELEMENTS)
1013 return _nodal->getValue();
1014 for (int i=0; i<_numberOfTypes; i++)
1015 if (_geometricTypes[i]==Type)
1016 return _nodal->getI(_count[i]);
1020 calculateDescendingConnectivity();
1021 if (Type==MED_ALL_ELEMENTS)
1022 return _descending->getValue();
1023 for (int i=0; i<_numberOfTypes; i++)
1024 if (_geometricTypes[i]==Type)
1025 return _descending->getI(Type);
1027 throw MEDEXCEPTION("Not found");
1031 //---------------------------------------------------------------------//
1032 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
1033 //---------------------------------------------------------------------//
1035 if (TypeConnectivity == MED_NODAL)
1037 calculateNodalConnectivity();
1038 return _nodal->getIndex();
1042 calculateDescendingConnectivity();
1043 return _descending->getIndex();
1047 /*! Not yet implemented */
1048 //----------------------------------------------//
1049 const int* CONNECTIVITY:: getNeighbourhood() const
1050 //----------------------------------------------//
1052 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
1055 /*! Returns an array which contains, for each node, all cells
1057 //-------------------------------------------------//
1058 const int* CONNECTIVITY::getReverseNodalConnectivity()
1059 //-------------------------------------------------//
1061 calculateReverseNodalConnectivity();
1062 return _reverseNodalConnectivity->getValue();
1065 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
1066 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
1067 //-------------------------------------------------------//
1068 const int* CONNECTIVITY::getReverseNodalConnectivityIndex()
1069 //-------------------------------------------------------//
1071 calculateReverseNodalConnectivity();
1072 return _reverseNodalConnectivity->getIndex();
1075 /*! Returns an array which contains, for each face (or edge),
1076 the 2 cells of each side. First is cell which face normal is outgoing.
1078 //------------------------------------------------------//
1079 const int* CONNECTIVITY::getReverseDescendingConnectivity()
1080 //------------------------------------------------------//
1082 // it is in _constituent connectivity only if we are in MED_CELL
1083 // (we could not for instance calculate face-edge connectivity !)
1084 if (_entity!=MED_CELL)
1085 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
1087 // we need descending connectivity
1088 calculateDescendingConnectivity();
1089 return _reverseDescendingConnectivity->getValue();
1092 /*! calculate the reverse descending Connectivity
1093 and returns the index ( A DOCUMENTER MIEUX)*/
1094 //-----------------------------------------------------------//
1095 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
1096 //-----------------------------------------------------------//
1098 // it is in _constituent connectivity only if we are in MED_CELL
1099 if (_entity!=MED_CELL)
1100 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
1102 // we need descending connectivity
1103 calculateDescendingConnectivity();
1104 return _reverseDescendingConnectivity->getIndex();
1107 /*! A DOCUMENTER (et a finir ???) */
1108 //--------------------------------------------//
1109 void CONNECTIVITY::calculateNodalConnectivity()
1110 //--------------------------------------------//
1114 if (_descending==NULL)
1115 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
1116 // calculate _nodal from _descending
1120 /*! If not yet done, calculate the nodal Connectivity
1121 and the reverse nodal Connectivity*/
1122 //---------------------------------------------------//
1123 void CONNECTIVITY::calculateReverseNodalConnectivity()
1124 //---------------------------------------------------//
1126 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
1130 SCRUTE(_reverseNodalConnectivity);
1134 calculateNodalConnectivity();
1136 if(_reverseNodalConnectivity==NULL)
1138 int node_number = 0;
1139 vector <vector <int> > reverse_connectivity;
1140 reverse_connectivity.resize(_numberOfNodes+1);
1142 // Treat all cells types
1144 for (int j = 0; j < _numberOfTypes; j++)
1146 // node number of the cell type
1147 node_number = _type[j].getNumberOfNodes();
1148 // treat all cells of a particular type
1149 for (int k = _count[j]; k < _count[j+1]; k++)
1150 // treat all nodes of the cell type
1151 for (int local_node_number = 1;
1152 local_node_number < node_number+1;
1153 local_node_number++)
1155 int global_node = _nodal->getIJ(k,local_node_number);
1156 reverse_connectivity[global_node].push_back(k);
1160 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
1162 //calculate size of reverse_nodal_connectivity array
1163 int size_reverse_nodal_connectivity = 0;
1164 for (int i = 1; i < _numberOfNodes+1; i++)
1165 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
1167 //MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity);
1168 int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
1169 int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
1170 //const int * reverse_nodal_connectivity = ReverseConnectivity->getValue();
1171 //const int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex();
1173 reverse_nodal_connectivity_index[0] = 1;
1174 for (int i = 1; i < _numberOfNodes+1; i++)
1176 int size = reverse_connectivity[i].size();
1177 reverse_nodal_connectivity_index[i] =
1178 reverse_nodal_connectivity_index[i-1] + size;
1179 for (int j = 0; j < size; j++)
1180 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
1183 //_reverseNodalConnectivity = ReverseConnectivity;
1184 _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
1185 reverse_nodal_connectivity_index,
1186 reverse_nodal_connectivity);
1187 delete [] reverse_nodal_connectivity_index;
1188 delete [] reverse_nodal_connectivity;
1193 /*! If not yet done, calculate the Descending Connectivity */
1194 //-------------------------------------------------//
1195 void CONNECTIVITY::calculateDescendingConnectivity()
1196 //-------------------------------------------------//
1198 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
1201 if (_descending==NULL)
1205 MESSAGE(LOC<<"No connectivity found !");
1206 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1208 // calcul _descending from _nodal
1209 // we need CONNECTIVITY for constituent
1211 if (_constituent != NULL)
1212 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
1213 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
1215 if (_entityDimension == 3)
1216 _constituent = new CONNECTIVITY(MED_FACE);
1217 else if (_entityDimension == 2)
1218 _constituent = new CONNECTIVITY(MED_EDGE);
1221 MESSAGE(LOC<<"We are in 1D");
1225 _constituent->_typeConnectivity = MED_NODAL;
1226 _constituent->_numberOfNodes = _numberOfNodes;
1227 // foreach cells, we built array of constituent
1228 int DescendingSize = 0;
1229 for(int i=0; i<_numberOfTypes; i++)
1230 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1231 //_descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize);
1232 //const int * descend_connectivity = _descending->getValue();
1233 int * descend_connectivity = new int[DescendingSize];
1234 for (int i=0; i<DescendingSize; i++)
1235 descend_connectivity[i]=0;
1236 //const int * descend_connectivity_index = _descending->getIndex();
1237 int * descend_connectivity_index = new int[_count[_numberOfTypes]];
1238 descend_connectivity_index[0]=1;
1239 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
1240 ConstituentsTypes[0]=MED_NONE;
1241 ConstituentsTypes[1]=MED_NONE;
1242 int * NumberOfConstituentsForeachType = new int[2];
1243 NumberOfConstituentsForeachType[0]=0;
1244 NumberOfConstituentsForeachType[1]=0;
1245 for(int i=0; i<_numberOfTypes; i++)
1247 // initialize descend_connectivity_index array :
1248 int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1249 for (int j=_count[i];j<_count[i+1];j++)
1251 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1252 // compute number of constituent of all cell for each type
1253 for(int k=1;k<NumberOfConstituents+1;k++)
1255 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1256 if (ConstituentsTypes[0]==MED_NONE)
1258 ConstituentsTypes[0]=MEDType;
1259 NumberOfConstituentsForeachType[0]++;
1260 } else if (ConstituentsTypes[0]==MEDType)
1261 NumberOfConstituentsForeachType[0]++;
1262 else if (ConstituentsTypes[1]==MED_NONE)
1264 ConstituentsTypes[1]=MEDType;
1265 NumberOfConstituentsForeachType[1]++;
1266 } else if (ConstituentsTypes[1]==MEDType)
1267 NumberOfConstituentsForeachType[1]++;
1269 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
1274 // we could built _constituent !
1275 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1276 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1278 //_constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes);
1280 // we use _constituent->_nodal
1281 //const int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
1282 //const int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
1283 int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1284 int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1285 ConstituentNodalConnectivityIndex[0]=1;
1287 _constituent->_entityDimension=ConstituentsTypes[0]/100;
1288 if (ConstituentsTypes[1]==MED_NONE)
1289 _constituent->_numberOfTypes = 1;
1291 _constituent->_numberOfTypes = 2;
1292 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1293 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1294 _constituent->_count = new int[_constituent->_numberOfTypes+1];
1295 _constituent->_count[0]=1;
1296 int* tmp_NumberOfConstituentsForeachType = new int[_constituent->_numberOfTypes+1];
1297 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1298 for (int i=0; i<_constituent->_numberOfTypes;i++) {
1299 _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1300 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1301 CELLMODEL Type(ConstituentsTypes[i]);
1302 _constituent->_type[i]=Type;
1303 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1304 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1305 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1307 delete [] ConstituentsTypes;
1308 delete [] NumberOfConstituentsForeachType;
1310 // we need reverse nodal connectivity
1311 if (! _reverseNodalConnectivity)
1312 calculateReverseNodalConnectivity();
1313 const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1314 const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1316 // array to keep reverse descending connectivity
1317 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1319 int TotalNumberOfSubCell = 0;
1320 for (int i=0; i<_numberOfTypes; i++)
1321 { // loop on all cell type
1322 CELLMODEL Type = _type[i];
1323 // int NumberOfNodesPerCell = Type.getNumberOfNodes();
1324 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1325 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1326 for (int k=1; k<=NumberOfConstituentPerCell; k++)
1327 { // we loop on all sub cell of it
1328 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0)
1330 // it is a new sub cell !
1331 // TotalNumberOfSubCell++;
1333 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0])
1335 tmp_NumberOfConstituentsForeachType[0]++;
1336 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1340 tmp_NumberOfConstituentsForeachType[1]++;
1341 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1343 //we have maximum two types
1345 descend_connectivity[descend_connectivity_index[j-1]+k-2]=
1346 TotalNumberOfSubCell;
1347 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1349 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1351 int * NodesLists = new int[NumberOfNodesPerConstituent];
1352 for (int l=0; l<NumberOfNodesPerConstituent; l++)
1354 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1355 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l] =
1358 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1360 // all elements which contains first node of sub cell :
1361 int ReverseNodalConnectivityIndex_0 =
1362 ReverseNodalConnectivityIndex[NodesLists[0]-1];
1363 int ReverseNodalConnectivityIndex_1 =
1364 ReverseNodalConnectivityIndex[NodesLists[0]];
1365 int NumberOfCellsInList =
1366 ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1368 if (NumberOfCellsInList > 0)
1369 { // we could have no element !
1370 int * CellsList = new int[NumberOfCellsInList];
1371 for (int l=ReverseNodalConnectivityIndex_0;
1372 l<ReverseNodalConnectivityIndex_1; l++)
1373 CellsList[l-ReverseNodalConnectivityIndex_0]=
1374 ReverseNodalConnectivityValue[l-1];
1376 // foreach node in sub cell, we search elements which are in common
1377 // at the end, we must have only one !
1379 for (int l=1; l<NumberOfNodesPerConstituent; l++)
1381 int NewNumberOfCellsInList = 0;
1382 int * NewCellsList = new int[NumberOfCellsInList];
1383 for (int l1=0; l1<NumberOfCellsInList; l1++)
1384 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1];
1385 l2<ReverseNodalConnectivityIndex[NodesLists[l]];
1388 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1389 // increasing order : CellsList[l1] are not in elements list of node l
1391 if ((CellsList[l1]==
1392 ReverseNodalConnectivityValue[l2-1])&
1395 // we have found one
1396 NewCellsList[NewNumberOfCellsInList] =
1398 NewNumberOfCellsInList++;
1402 NumberOfCellsInList = NewNumberOfCellsInList;
1404 delete [] CellsList;
1405 CellsList = NewCellsList;
1408 if (NumberOfCellsInList > 0)
1409 { // We have found some elements !
1410 if (NumberOfCellsInList>1)
1411 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1413 int CellNumber = CellsList[0];
1415 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1] =
1418 // we search sub cell number in this cell to not calculate it another time
1421 for (int l=0; l<_numberOfTypes; l++)
1422 if (CellNumber < _count[l+1]) {
1426 //int sub_cell_count2 = Type2.get_entities_count(1);
1427 //int nodes_cell_count2 = Type2.get_nodes_count();
1429 for (int l=1; l<=Type2.getNumberOfConstituents(1);l++)
1430 { // on all sub cell
1432 for (int m=1; m<=Type2.getConstituentType(1,l)%100;
1434 for (int n=1; n<=Type.getConstituentType(1,k)%100;
1437 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) ==
1441 if (counter==Type.getConstituentType(1,k)%100)
1443 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=
1444 -1*TotalNumberOfSubCell; // because, see it in other side !
1452 MESSAGE(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1457 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1459 delete [] CellsList;
1462 delete [] NodesLists;
1466 // we adjust _constituent
1467 int NumberOfConstituent=0;
1468 int SizeOfConstituentNodal=0;
1469 for (int i=0;i<_constituent->_numberOfTypes; i++)
1471 NumberOfConstituent +=
1472 tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1473 SizeOfConstituentNodal +=
1474 (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1476 // we built new _nodal attribute in _constituent
1477 //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1478 //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1479 //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1480 int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1481 int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1482 ConstituentNodalIndex[0]=1;
1483 // we build _reverseDescendingConnectivity
1484 //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1485 //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1486 //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1487 int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1488 int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1489 reverseDescendingConnectivityIndex[0]=1;
1491 // first constituent type
1492 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++)
1494 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1495 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
1497 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1499 reverseDescendingConnectivityIndex[j+1] =
1500 reverseDescendingConnectivityIndex[j]+2;
1501 for(int k=reverseDescendingConnectivityIndex[j];
1502 k<reverseDescendingConnectivityIndex[j+1]; k++)
1504 reverseDescendingConnectivityValue[k-1] =
1505 ReverseDescendingConnectivityValue[k-1];
1508 // second type if any
1509 if (_constituent->_numberOfTypes==2)
1511 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1512 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1513 int offset2=offset*2;
1514 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1515 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++)
1517 ConstituentNodalIndex[j+1]=
1518 ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1519 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
1521 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1523 reverseDescendingConnectivityIndex[j+1] =
1524 reverseDescendingConnectivityIndex[j]+2;
1525 for(int k=reverseDescendingConnectivityIndex[j];
1526 k<reverseDescendingConnectivityIndex[j+1]; k++)
1528 reverseDescendingConnectivityValue[k-1] =
1529 ReverseDescendingConnectivityValue[offset2+k-1];
1532 _constituent->_count[2]=NumberOfConstituent+1;
1533 // we correct _descending to adjust face number
1534 for(int j=0;j<DescendingSize;j++)
1536 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1537 descend_connectivity[j]-=offset;
1538 else if (descend_connectivity[j]<-tmp_NumberOfConstituentsForeachType[0])
1539 descend_connectivity[j]+=offset;
1543 delete [] ConstituentNodalConnectivityIndex;
1544 delete [] ConstituentNodalConnectivity;
1546 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1548 descend_connectivity_index,
1549 descend_connectivity);
1550 delete [] descend_connectivity_index;
1551 delete [] descend_connectivity;
1552 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,
1553 2*NumberOfConstituent,
1554 reverseDescendingConnectivityIndex,
1555 reverseDescendingConnectivityValue);
1556 delete [] reverseDescendingConnectivityIndex;
1557 delete [] reverseDescendingConnectivityValue;
1559 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1560 delete [] tmp_NumberOfConstituentsForeachType;
1562 //delete _constituent->_nodal;
1563 _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1564 SizeOfConstituentNodal,
1565 ConstituentNodalIndex,
1566 ConstituentNodalValue);
1568 delete [] ConstituentNodalIndex;
1569 delete [] ConstituentNodalValue;
1571 delete [] ReverseDescendingConnectivityValue;
1576 //-----------------------------------------------------------------------------------------//
1577 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1578 //-----------------------------------------------------------------------------------------//
1580 // if (_entity==MED_CELL)
1581 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1582 // Connectivity : we could not do that with MED_CELL entity !");
1583 // if (myConnectivity->getEntity()!=MED_CELL)
1584 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1585 // Connectivity : we need MED_CELL connectivity !");
1589 /*! Not implemented yet */
1590 //--------------------------------------------------------------------//
1591 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1592 //--------------------------------------------------------------------//
1594 const char * LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1597 MESSAGE(LOC<<"method not yet implemented " << myConnectivity._entity);
1606 Returns the geometry of an element given by its entity type & its global number.
1608 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1610 //--------------------------------------------------------------------//
1611 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1612 //--------------------------------------------------------------------//
1614 const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1616 int globalNumberMin = 1;
1617 int globalNumberMax ;
1619 if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
1620 else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
1622 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1624 // The globalNumber must verify : 1 <= globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1626 if ( (globalNumber < globalNumberMin) || (globalNumber > globalNumberMax-1 ) )
1627 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1628 << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
1630 if (_entity==Entity) {
1631 for(int i=1; i<=_numberOfTypes;i++)
1632 if (globalNumber<_count[i])
1633 return _geometricTypes[i-1];
1635 else if (_constituent!=NULL)
1636 return _constituent->getElementType(Entity,globalNumber);
1638 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1639 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1645 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1647 os << endl << "------------- Entity = ";
1662 case MED_ALL_ENTITIES:
1663 os << "MED_ALL_ENTITIES";
1669 os << " -------------\n\nMedConnectivity : ";
1670 switch (co._typeConnectivity)
1673 os << "MED_NODAL\n";
1675 case MED_DESCENDING:
1676 os << "MED_DESCENDING\n";
1681 os << "Entity dimension : " << co._entityDimension << endl;
1682 os << "Number of nodes : " << co._numberOfNodes << endl;
1683 os << "Number of types : " << co._numberOfTypes << endl;
1684 for (int i=0; i!=co._numberOfTypes ; ++i)
1685 os << " -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1686 << co._count[i+1]-co._count[i] << " elements" << endl;
1688 if (co._typeConnectivity == MED_NODAL )
1690 for (int i=0; i<co._numberOfTypes; i++)
1692 os << endl << co._type[i].getName() << " : " << endl;
1693 int numberofelements = co._count[i+1]-co._count[i];
1694 const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, co._geometricTypes[i]);
1695 int numberofnodespercell = co._geometricTypes[i]%100;
1696 for (int j=0;j<numberofelements;j++)
1698 os << setw(6) << j+1 << " : ";
1699 for (int k=0;k<numberofnodespercell;k++)
1700 os << connectivity[j*numberofnodespercell+k]<<" ";
1705 else if (co._typeConnectivity == MED_DESCENDING)
1707 int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1708 const int *connectivity = co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1709 const int *connectivity_index = co.getConnectivityIndex( co._typeConnectivity, co._entity );
1711 for ( int j=0; j!=numberofelements; j++ )
1713 os << "Element " << j+1 << " : ";
1714 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1715 os << connectivity[k-1] << " ";
1720 if (co._constituent)
1721 os << endl << *co._constituent << endl;