2 #include "MEDMEM_Connectivity.hxx"
3 #include "MEDMEM_Family.hxx"
4 #include "MEDMEM_CellModel.hxx"
6 #include "MEDMEM_SkyLineArray.hxx"
7 #include "MEDMEM_ModulusArray.hxx"
9 #include "MEDMEM_STRING.hxx"
11 //------------------------------------------------------//
12 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity/*=MED_CELL*/):
13 //------------------------------------------------------//
15 _typeConnectivity(MED_NODAL),
17 _geometricTypes((medGeometryElement*)NULL),
18 _type((CELLMODEL*)NULL),
21 _nodal((MEDSKYLINEARRAY*)NULL),
22 _descending((MEDSKYLINEARRAY*)NULL),
23 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
24 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
25 _neighbourhood((MEDSKYLINEARRAY*)NULL),
26 _constituent((CONNECTIVITY*)NULL)
28 MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)") ;
31 //-------------------------------------------------------------------------//
32 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity/*=MED_CELL*/):
33 //-------------------------------------------------------------------------//
35 _typeConnectivity(MED_NODAL),
36 _numberOfTypes(numberOfTypes),
38 _nodal((MEDSKYLINEARRAY*)NULL),
39 _descending((MEDSKYLINEARRAY*)NULL),
40 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
41 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
42 _neighbourhood((MEDSKYLINEARRAY*)NULL),
43 _constituent((CONNECTIVITY*)NULL)
45 MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)") ;
46 _geometricTypes = new medGeometryElement[numberOfTypes];
47 _type = new CELLMODEL[numberOfTypes];
48 _count = new int[numberOfTypes];
51 //-------------------------------------------------------------------------//
52 CONNECTIVITY::CONNECTIVITY(CONNECTIVITY & m)
53 //-------------------------------------------------------------------------//
56 _typeConnectivity = m._typeConnectivity;
57 _numberOfTypes = m._numberOfTypes;
58 if (m._geometricTypes != NULL)
60 _geometricTypes = new medGeometryElement[m._numberOfTypes];
61 memcpy(_geometricTypes,m._geometricTypes,m._numberOfTypes*sizeof(medGeometryElement));
64 _geometricTypes = (medGeometryElement *) NULL;
66 _type = new CELLMODEL(* m._type);
68 _type = (CELLMODEL *) NULL;
69 _entityDimension = m._entityDimension;
70 _numberOfNodes = m._numberOfNodes;
73 _count = new med_int[m._numberOfTypes+1];
74 memcpy(_count,m._count,(m._numberOfTypes+1)*sizeof(med_int));
77 _count = (med_int *) NULL;
79 _nodal = new MEDSKYLINEARRAY(* m._nodal);
81 _nodal = (MEDSKYLINEARRAY *) NULL;
82 if (m._descending != NULL)
83 _descending = new MEDSKYLINEARRAY(* m._descending);
85 _descending = (MEDSKYLINEARRAY *) NULL;
86 if (m._reverseNodalConnectivity != NULL)
87 _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
89 _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
90 if (m._reverseDescendingConnectivity != NULL)
91 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
93 _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
94 if (m._neighbourhood != NULL)
95 _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
97 _neighbourhood = (MEDSKYLINEARRAY *) NULL;
98 if (m._constituent != NULL)
99 _constituent = new CONNECTIVITY(* m._constituent);
101 _constituent = (CONNECTIVITY *) NULL;
104 //----------------------------//
105 CONNECTIVITY::~CONNECTIVITY()
106 //----------------------------//
108 MESSAGE("Destructeur de CONNECTIVITY()") ;
109 if ( _geometricTypes != NULL )
110 delete [] _geometricTypes ;
111 if ( _count != NULL )
113 if ( _nodal != NULL )
115 if ( _descending != NULL )
117 if ( _reverseNodalConnectivity != NULL )
118 delete _reverseNodalConnectivity ;
119 if ( _reverseDescendingConnectivity != NULL )
120 delete _reverseDescendingConnectivity ;
121 if ( _constituent != NULL )
122 delete _constituent ;
126 //------------------------------------------------------------------------------------------//
127 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
128 //------------------------------------------------------------------------------------------//
130 MESSAGE("CONNECTIVITY::calculateConnectivity");
132 // a temporary limitation due to calculteDescendingConnectivity function !!!
133 if ((_entityDimension==3) & (Entity==MED_EDGE))
134 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
137 if (ConnectivityType==MED_NODAL)
138 calculateNodalConnectivity() ;
140 if (Entity==MED_CELL)
141 calculateDescendingConnectivity() ;
143 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
144 if (Entity!=_entity) {
145 calculateDescendingConnectivity() ;
146 _constituent->calculateConnectivity(ConnectivityType,Entity) ;
150 /*! Give, in full or no interlace mode (for nodal connectivity),
151 descending or nodal connectivity.
153 You must give a <medEntityMesh> (ie:MED_EDGE)
154 and a <medGeometryElement> (ie:MED_SEG3). */
156 //------------------------------------------------------------//
157 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
158 //------------------------------------------------------------//
160 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
163 int numberOfFamilies = myFamilies.size();
164 if (numberOfFamilies == 0 ) {
165 MESSAGE(LOC<<"No family") ;
168 // does we do an update ?
169 if ((_constituent != NULL)&(_descending != NULL)) {
170 MESSAGE(LOC<<"Constituent is already defined") ;
174 if ((_constituent != NULL)&(_descending == NULL)) {
175 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
176 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing") ;
180 for(int i=0; i<numberOfFamilies; i++) {
181 FAMILY * myFamily = myFamilies[i] ;
182 MESSAGE(LOC<<"updating the family (BEGIN) : " << *myFamily);
185 // well we could go !
186 CONNECTIVITY * oldConstituent = _constituent ;
187 _constituent = (CONNECTIVITY *)NULL ;
188 // for instance we must have nodal connectivity in constituent :
189 if (oldConstituent->_nodal == NULL)
191 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
192 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell")) ;
194 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf() ;
195 int * oldConstituentValue = oldConstituent->_nodal->getValue() ;
196 int * oldConstituentIndex = oldConstituent->_nodal->getIndex() ;
198 SCRUTE(oldNumberOfFace);
200 calculateDescendingConnectivity() ;
202 // if (oldConstituent->_nodal != NULL) {
203 int newNumberOfFace = _constituent->_nodal->getNumberOf() ;
204 int * newConstituentValue = _constituent->_nodal->getValue() ;
205 int * newConstituentIndex = _constituent->_nodal->getIndex() ;
207 SCRUTE(newNumberOfFace);
209 int * newReverseDescendingIndex =
210 _reverseDescendingConnectivity->getIndex();
211 int * newReverseDescendingValue =
212 _reverseDescendingConnectivity->getValue();
214 int * newDescendingIndex = _descending->getIndex();
215 int * newDescendingValue = _descending->getValue();
217 // loop on all family,
218 // for all constituent in family, we get it's old connectivity
219 // (with oldCconstituentValue and oldConstituentIndex)
220 // and search the constituent in newConstituentValue with class
223 // when a new face is found, replace old constituent
224 // number in family with new one
225 // If face is not rigth oriented, we must change _descending attribute
226 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
228 // Voila a toi de jouer Nadir :-)
230 // First we get the renumbering from the oldCconstituentValue and
231 // oldConstituentIndex in the the new one, newConstituentValue and
232 // newConstituentIndex with a negative sign if the face is not
235 int * renumberingFromOldToNew = new int [oldNumberOfFace];
239 _constituent->calculateReverseNodalConnectivity() ;
241 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
245 renumberingFromOldToNew[iOldFace] = iOldFace+1;
247 // renumberingFromOldToNew[iOldFace] = 999999;
249 int face_it_beginOld = oldConstituentIndex[iOldFace];
250 int face_it_endOld = oldConstituentIndex[iOldFace+1];
251 int face_size_itOld = face_it_endOld - face_it_beginOld;
253 int* NodesLists = oldConstituentValue + (face_it_beginOld-1) ;
256 int * reverseFaceNodal = _constituent->getReverseNodalConnectivity() ;
257 int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex() ;
259 // set an array wich contains faces numbers arround first node
260 int BeginIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]-1] ;
261 int EndIndexFaceArrayFirstNode=reverseFaceNodalIndex[NodesLists[0]] ;
262 int NumberOfFacesInList=EndIndexFaceArrayFirstNode-BeginIndexFaceArrayFirstNode;
264 int * FacesList = new int[NumberOfFacesInList] ;
265 for (int l=BeginIndexFaceArrayFirstNode; l<EndIndexFaceArrayFirstNode; l++)
266 FacesList[l-BeginIndexFaceArrayFirstNode]=reverseFaceNodal[l-1] ;
268 // foreach node in sub cell, we search elements which are in common
269 // at the end, we must have only one !
271 for (int nodeFaceOld=1; nodeFaceOld<face_size_itOld; nodeFaceOld++)
274 int NewNumberOfFacesInList = 0 ;
275 int * NewFacesList = new int[NumberOfFacesInList] ;
277 for (int l1=0; l1<NumberOfFacesInList; l1++) {
278 for (int l2=reverseFaceNodalIndex[NodesLists[nodeFaceOld]-1]; l2<reverseFaceNodalIndex[NodesLists[nodeFaceOld]]; l2++) {
279 if (FacesList[l1]<reverseFaceNodal[l2-1])
280 // increasing order : FacesList[l1] are not in elements list of node l
282 if (FacesList[l1]==reverseFaceNodal[l2-1]) {
284 NewFacesList[NewNumberOfFacesInList]=FacesList[l1];
285 NewNumberOfFacesInList++;
290 NumberOfFacesInList = NewNumberOfFacesInList;
291 delete [] FacesList ;
292 FacesList = NewFacesList;
296 if (!NumberOfFacesInList==0)
298 if (NumberOfFacesInList>1)
299 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one face found ("<<NumberOfFacesInList<<") ! " <<FacesList[0]<<" "<<FacesList[1] ));
301 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,NodesLists);
303 SCRUTE(NumberOfFacesInList);
305 SCRUTE(newConstituentIndex);
307 int face_it_beginNew = newConstituentIndex[FacesList[0]-1];
308 int face_it_endNew = newConstituentIndex[FacesList[0]];
309 face_size_itNew = face_it_endNew - face_it_beginNew;
311 int * newNodesLists = newConstituentValue+newConstituentIndex[FacesList[0]-1]-1 ;
312 MEDMODULUSARRAY modulusArrayNew(face_size_itOld,newNodesLists);
314 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
316 // Real new face found
318 if(retCompareNewOld == 1)
320 renumberingFromOldToNew[iOldFace] = FacesList[0] ;
322 MESSAGE("Renumbering index " << iOldFace << " val "<< FacesList[0]);
328 // Reverse new face found
330 if(retCompareNewOld == -1)
332 renumberingFromOldToNew[iOldFace] = FacesList[0];
334 MESSAGE("Renumbering index " << iOldFace << " val "<< FacesList[0]);
339 int face_it_begin = newReverseDescendingIndex[FacesList[0]-1];
340 int face_it_end = newReverseDescendingIndex[FacesList[0]];
341 int face_size_it = face_it_end - face_it_begin;
343 if (face_size_it == 1)
344 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
346 if (face_size_it > 2)
347 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
349 // we have always 2 neighbourings
350 int cell1 = newReverseDescendingValue[face_it_begin-1];
351 int cell2 = newReverseDescendingValue[face_it_begin];
353 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
355 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
357 if (cell2 != 0) { // we are not on border !!!!
359 newReverseDescendingValue[face_it_begin-1] = cell2;
360 // Updating _constituent->_nodal because of reversity
361 int * oldArray = oldConstituentValue+face_it_beginOld-1;
362 int * newArray = newConstituentValue+face_it_beginNew-1;
363 for(int iarray=0;iarray<face_size_itNew;iarray++)
364 newArray[iarray] = oldArray[iarray] ;
367 // Updating _reverseDescendingConnectivity
370 newReverseDescendingValue[face_it_begin] = cell1;
372 // Updating _descending for cell1 and cell2
373 for(int iface=newDescendingIndex[cell1-1];iface<newDescendingIndex[cell1];iface++)
374 if (newDescendingValue[iface-1]==FacesList[0])
375 newDescendingValue[iface-1]=-FacesList[0] ;
376 else if (newDescendingValue[iface-1]==-FacesList[0])
377 newDescendingValue[iface-1]=FacesList[0] ;
379 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
380 if (newDescendingValue[iface-1]==FacesList[0])
381 newDescendingValue[iface-1]=-FacesList[0] ;
382 else if (newDescendingValue[iface-1]==-FacesList[0])
383 newDescendingValue[iface-1]=FacesList[0] ;
384 } else {// else we are on border and we do nothing !!!!!!!!
385 INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
386 INFOS(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !") ;
387 INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
393 INFOS(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
394 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
399 MESSAGE(LOC<<"The Renumbering is finished and the status is");
403 // Updating the Family
405 for(int i=0; i<numberOfFamilies; i++) {
406 FAMILY * myFamily = myFamilies[i] ;
408 int length_skyline = myFamily->getnumber()->getLength();
409 int * value_skyline = myFamily->getnumber()->getValue();
411 for (int i=0;i<length_skyline;i++) {
412 MESSAGE("OLD " << value_skyline[i] << " NEW " << renumberingFromOldToNew[value_skyline[i]-1]);
413 value_skyline[i] = renumberingFromOldToNew[value_skyline[i]-1];
415 MESSAGE(LOC<<"updating the family (END) : " << *myFamily);
423 //------------------------------------------------------------------------------------------------------------------//
424 med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
425 //------------------------------------------------------------------------------------------------------------------//
427 const char * LOC = "CONNECTIVITY::getConnectivity" ;
430 MEDSKYLINEARRAY * Connectivity ;
431 if (Entity==_entity) {
433 if (ConnectivityType==MED_NODAL)
435 calculateNodalConnectivity();
440 calculateDescendingConnectivity() ;
441 Connectivity=_descending;
444 if (Connectivity!=NULL)
445 if (Type==MED_ALL_ELEMENTS)
446 return Connectivity->getValue() ;
448 for (med_int i=0; i<_numberOfTypes; i++)
449 if (_geometricTypes[i]==Type)
450 return Connectivity->getI(_count[i]) ;
451 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !")) ;
454 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !")) ;
456 if (_constituent != NULL)
457 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
459 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
462 /*! Give morse index array to use with
463 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
465 Each value give start index for corresponding entity in connectivity array.
467 Example : i-th element, j-th node of it :
468 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
469 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
470 //-----------------------------------------------------------------------------------------------//
471 med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
472 //------------------------------------------------------------------------------------------------//
474 const char * LOC = "CONNECTIVITY::getConnectivityIndex" ;
477 MEDSKYLINEARRAY * Connectivity ;
478 if (Entity==_entity) {
480 if (ConnectivityType==MED_NODAL)
483 Connectivity=_descending;
485 if (Connectivity!=NULL)
486 return Connectivity->getIndex() ;
488 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !")) ;
491 if (_constituent != NULL)
492 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
494 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
498 //--------------------------------------------------------------//
499 CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
500 //--------------------------------------------------------------//
502 const char * LOC = "CONNECTIVITY::getType" ;
505 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
506 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !")) ;
507 for (med_int i=0; i<_numberOfTypes; i++)
508 if (_geometricTypes[i]==Type)
510 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !")) ;
513 /*! Returns the number of elements of type <med geometrie element>.
514 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
515 //------------------------------------------------------------------------//
516 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
517 //------------------------------------------------------------------------//
519 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType" ;
522 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
523 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!")) ;
524 for (med_int i=0; i<_numberOfTypes; i++)
525 if (_geometricTypes[i]==Type)
526 return _type[i].getNumberOfNodes() ;
527 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !")) ;
530 /*! Returns the number of geometric sub cells of <med geometrie element> type.
531 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
532 //------------------------------------------------------------------------//
533 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
534 //------------------------------------------------------------------------//
536 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
537 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!") ;
538 for (med_int i=0; i<_numberOfTypes; i++)
539 if (_geometricTypes[i]==Type)
540 return _type[i].getNumberOfConstituents(1) ;
541 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !") ;
544 /*! Returns the number of elements of type <med geometrie element>.
547 - Implemented for MED_ALL_ELEMENTS
548 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
549 - Not implemented for MED_NONE */
550 //-----------------------------------------------------------------------------------//
551 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
552 //-----------------------------------------------------------------------------------//
554 const char * LOC = "CONNECTIVITY::getNumberOf" ;
557 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
559 if (Entity==_entity) {
561 return 0 ; // not defined !
562 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
563 if (Type==MED_ALL_ELEMENTS)
564 return _count[_numberOfTypes]-1;
565 for (med_int i=0; i<_numberOfTypes; i++)
566 if (_geometricTypes[i]==Type)
567 return (_count[i+1] - _count[i]);
568 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
570 if (_constituent != NULL)
571 return _constituent->getNumberOf(Entity,Type);
573 return 0 ; // valid if they are nothing else !
574 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
578 //--------------------------------------------------------------//
579 med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
580 medGeometryElement Type)
581 //--------------------------------------------------------------//
583 if (TypeConnectivity == MED_NODAL)
585 calculateNodalConnectivity();
586 if (Type==MED_ALL_ELEMENTS)
587 return _nodal->getValue();
588 for (med_int i=0; i<_numberOfTypes; i++)
589 if (_geometricTypes[i]==Type)
590 return _nodal->getI(_count[i]);
594 calculateDescendingConnectivity();
595 if (Type==MED_ALL_ELEMENTS)
596 return _descending->getValue();
597 for (med_int i=0; i<_numberOfTypes; i++)
598 if (_geometricTypes[i]==Type)
599 return _descending->getI(Type);
601 throw MEDEXCEPTION("Not found");
605 //---------------------------------------------------------------------//
606 med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
607 //---------------------------------------------------------------------//
609 if (TypeConnectivity == MED_NODAL)
611 calculateNodalConnectivity();
612 return _nodal->getIndex();
616 calculateDescendingConnectivity();
617 return _descending->getIndex();
621 /*! Not yet implemented */
622 //----------------------------------------------//
623 med_int* CONNECTIVITY:: getNeighbourhood() const
624 //----------------------------------------------//
626 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
629 /*! Returns an array which contains, for each node, all cells
631 //-------------------------------------------------//
632 med_int* CONNECTIVITY::getReverseNodalConnectivity()
633 //-------------------------------------------------//
635 calculateReverseNodalConnectivity();
636 return _reverseNodalConnectivity->getValue();
639 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
640 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
641 //-------------------------------------------------------//
642 med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
643 //-------------------------------------------------------//
645 calculateReverseNodalConnectivity();
646 return _reverseNodalConnectivity->getIndex();
649 /*! Returns an array which contains, for each face (or edge),
650 the 2 cells of each side. First is cell which face normal is outgoing.
652 //------------------------------------------------------//
653 med_int* CONNECTIVITY::getReverseDescendingConnectivity()
654 //------------------------------------------------------//
656 // it is in _constituent connectivity only if we are in MED_CELL
657 // (we could not for instance calculate face-edge connectivity !)
658 if (_entity!=MED_CELL)
659 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
661 // we need descending connectivity
662 calculateDescendingConnectivity();
663 return _reverseDescendingConnectivity->getValue();
666 /*! calculate the reverse descending Connectivity
667 and returns the index ( A DOCUMENTER MIEUX)*/
668 //-----------------------------------------------------------//
669 med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
670 //-----------------------------------------------------------//
672 // it is in _constituent connectivity only if we are in MED_CELL
673 if (_entity!=MED_CELL)
674 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
676 // we need descending connectivity
677 calculateDescendingConnectivity();
678 return _reverseDescendingConnectivity->getIndex();
681 /*! A DOCUMENTER (et a finir ???) */
682 //--------------------------------------------//
683 void CONNECTIVITY::calculateNodalConnectivity()
684 //--------------------------------------------//
688 if (_descending==NULL)
689 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
690 // calculate _nodal from _descending
694 /*! If not yet done, calculate the nodal Connectivity
695 and the reverse nodal Connectivity*/
696 //---------------------------------------------------//
697 void CONNECTIVITY::calculateReverseNodalConnectivity()
698 //---------------------------------------------------//
700 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : " ;
704 calculateNodalConnectivity() ;
706 MESSAGE(LOC<<"Number of nodes = "<<_numberOfNodes);
708 if(_reverseNodalConnectivity==NULL) {
710 med_int node_number = 0;
711 vector <vector <med_int> > reverse_connectivity;
712 reverse_connectivity.resize(_numberOfNodes+1);
714 // Treat all cells types
716 for (med_int j = 0; j < _numberOfTypes; j++)
718 // node number of the cell type
719 node_number = _type[j].getNumberOfNodes();
720 // treat all cells of a particular type
721 for (med_int k = _count[j]; k < _count[j+1]; k++)
722 // treat all nodes of the cell type
723 for (med_int local_node_number = 1 ; local_node_number < node_number+1; local_node_number++)
725 med_int global_node = _nodal->getIJ(k,local_node_number) ;
726 reverse_connectivity[global_node].push_back(k);
730 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
732 //calculate size of reverse_nodal_connectivity array
733 med_int size_reverse_nodal_connectivity = 0;
734 for (med_int i = 1; i < _numberOfNodes+1; i++)
735 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
737 MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity) ;
738 // reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
739 // reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
740 med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue() ;
741 med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex() ;
743 reverse_nodal_connectivity_index[0] = 1;
744 for (med_int i = 1; i < _numberOfNodes+1; i++)
746 med_int size = reverse_connectivity[i].size();
747 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size ;
748 for (med_int j = 0; j < size; j++)
749 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
752 _reverseNodalConnectivity = ReverseConnectivity ;
757 /*! If not yet done, calculate the Descending Connectivity */
758 //-------------------------------------------------//
759 void CONNECTIVITY::calculateDescendingConnectivity()
760 //-------------------------------------------------//
762 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
765 if (_descending==NULL) {
767 MESSAGE(LOC<<"No connectivity found !");
768 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
770 // calcul _descending from _nodal
771 // we need CONNECTIVITY for constituent
773 if (_constituent != NULL)
774 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
775 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
777 if (_entityDimension == 3)
778 _constituent = new CONNECTIVITY(MED_FACE) ;
779 else if (_entityDimension == 2)
780 _constituent = new CONNECTIVITY(MED_EDGE) ;
782 MESSAGE(LOC<<"We are in 1D");
785 _constituent->_typeConnectivity = MED_DESCENDING ;
786 _constituent->_numberOfNodes = _numberOfNodes ;
787 // foreach cells, we built array of constituent
788 int DescendingSize = 0 ;
789 for(int i=0; i<_numberOfTypes; i++)
790 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1) ;
791 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize) ;
792 int * descend_connectivity = _descending->getValue() ;
793 for (int i=0; i<DescendingSize; i++)
794 descend_connectivity[i]=0;
795 int * descend_connectivity_index = _descending->getIndex() ;
796 descend_connectivity_index[0]=1;
797 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
798 ConstituentsTypes[0]=MED_NONE ;
799 ConstituentsTypes[1]=MED_NONE ;
800 int * NumberOfConstituentsForeachType = new int[2];
801 NumberOfConstituentsForeachType[0]=0;
802 NumberOfConstituentsForeachType[1]=0;
803 for(int i=0; i<_numberOfTypes; i++) {
804 // initialize descend_connectivity_index array :
805 int NumberOfConstituents = _type[i].getNumberOfConstituents(1) ;
806 for (int j=_count[i];j<_count[i+1];j++){
807 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents ;
808 // compute number of constituent of all cell for each type
809 for(int k=1;k<NumberOfConstituents+1;k++) {
810 medGeometryElement MEDType = _type[i].getConstituentType(1,k) ;
811 if (ConstituentsTypes[0]==MED_NONE) {
812 ConstituentsTypes[0]=MEDType;
813 NumberOfConstituentsForeachType[0]++ ;
814 } else if (ConstituentsTypes[0]==MEDType)
815 NumberOfConstituentsForeachType[0]++ ;
816 else if (ConstituentsTypes[1]==MED_NONE) {
817 ConstituentsTypes[1]=MEDType;
818 NumberOfConstituentsForeachType[1]++ ;
819 } else if (ConstituentsTypes[1]==MEDType)
820 NumberOfConstituentsForeachType[1]++ ;
822 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
827 // we could built _constituent !
828 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1] ;
829 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100) ;
830 _constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes) ;
832 // we use _constituent->_nodal
833 int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
834 int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
835 ConstituentNodalConnectivityIndex[0]=1;
837 _constituent->_entityDimension=ConstituentsTypes[0]/100;
838 if (ConstituentsTypes[1]==MED_NONE)
839 _constituent->_numberOfTypes = 1 ;
841 _constituent->_numberOfTypes = 2 ;
842 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes] ;
843 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes] ;
844 _constituent->_count = new med_int[_constituent->_numberOfTypes+1] ;
845 _constituent->_count[0]=1;
846 int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1] ;
847 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
848 for (int i=0; i<_constituent->_numberOfTypes;i++) {
849 _constituent->_geometricTypes[i]=ConstituentsTypes[i] ;
850 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i] ;
851 CELLMODEL Type(ConstituentsTypes[i]);
852 _constituent->_type[i]=Type;
853 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i] ;
854 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
855 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100) ;
857 delete[] ConstituentsTypes;
858 delete[] NumberOfConstituentsForeachType ;
860 // we need reverse nodal connectivity
861 if (! _reverseNodalConnectivity)
862 calculateReverseNodalConnectivity() ;
863 int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
864 int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
866 // array to keep reverse descending connectivity
867 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
869 int TotalNumberOfSubCell = 0 ;
870 for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
871 CELLMODEL Type = _type[i] ;
872 int NumberOfNodesPerCell = Type.getNumberOfNodes() ;
873 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
874 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
875 for (int k=1 ; k<=NumberOfConstituentPerCell; k++) { // we loop on all sub cell of it
876 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
877 // it is a new sub cell !
878 // TotalNumberOfSubCell++;
880 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
881 tmp_NumberOfConstituentsForeachType[0]++;
882 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
884 tmp_NumberOfConstituentsForeachType[1]++;
885 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
887 //we have maximum two types
889 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
890 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j ;
891 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100 ;
893 int * NodesLists = new int[NumberOfNodesPerConstituent] ;
894 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
895 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
896 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
898 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
900 // all elements which contains first node of sub cell :
901 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1] ;
902 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]] ;
903 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0 ;
905 if (NumberOfCellsInList > 0) { // we could have no element !
906 int * CellsList = new int[NumberOfCellsInList] ;
907 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
908 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1] ;
910 // foreach node in sub cell, we search elements which are in common
911 // at the end, we must have only one !
913 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
914 int NewNumberOfCellsInList = 0 ;
915 int * NewCellsList = new int[NumberOfCellsInList] ;
916 for (int l1=0; l1<NumberOfCellsInList; l1++)
917 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
918 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
919 // increasing order : CellsList[l1] are not in elements list of node l
921 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
923 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
924 NewNumberOfCellsInList++;
928 NumberOfCellsInList = NewNumberOfCellsInList;
930 delete [] CellsList ;
931 CellsList = NewCellsList;
934 if (NumberOfCellsInList > 0) { // We have found some elements !
935 int CellNumber = CellsList[0] ;
936 delete [] CellsList ;
937 if (NumberOfCellsInList>1)
938 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
940 if (NumberOfCellsInList==1) {
941 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber ;
943 // we search sub cell number in this cell to not calculate it another time
946 for (int l=0; l<_numberOfTypes; l++)
947 if (CellNumber < _count[l+1]) {
951 //int sub_cell_count2 = Type2.get_entities_count(1) ;
952 //int nodes_cell_count2 = Type2.get_nodes_count() ;
954 for (int l=1; l<=Type2.getNumberOfConstituents(1) ;l++) { // on all sub cell
956 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
957 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
958 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
961 if (counter==Type.getConstituentType(1,k)%100) {
962 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
969 INFOS(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!") ; // exception ?
972 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0 ;
975 delete[] NodesLists ;
979 // we adjust _constituent
980 int NumberOfConstituent=0 ;
981 int SizeOfConstituentNodal=0 ;
982 for (int i=0;i<_constituent->_numberOfTypes; i++) {
983 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
984 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
986 // we built new _nodal attribute in _constituent
987 MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
988 int *ConstituentNodalValue = ConstituentNodal->getValue();
989 int *ConstituentNodalIndex = ConstituentNodal->getIndex();
990 ConstituentNodalIndex[0]=1;
991 // we build _reverseDescendingConnectivity
992 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent) ;
993 int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
994 int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
995 reverseDescendingConnectivityIndex[0]=1;
997 // first constituent type
998 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
999 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1000 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1001 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1003 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1004 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1005 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1008 // second type if any
1009 if (_constituent->_numberOfTypes==2) {
1010 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1 ;
1011 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1012 int offset2=offset*2 ;
1013 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1014 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
1015 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1016 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
1017 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1019 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1020 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
1021 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1024 _constituent->_count[2]=NumberOfConstituent+1 ;
1025 // we correct _descending to adjust face number
1026 for(int j=0;j<DescendingSize;j++)
1027 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
1028 descend_connectivity[j]-=offset;
1031 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1032 delete _constituent->_nodal ;
1033 _constituent->_nodal = ConstituentNodal ;
1035 delete[] ReverseDescendingConnectivityValue ;
1040 //-----------------------------------------------------------------------------------------//
1041 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
1042 //-----------------------------------------------------------------------------------------//
1044 // if (_entity==MED_CELL)
1045 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1046 // Connectivity : we could not do that with MED_CELL entity !");
1047 // if (myConnectivity->getEntity()!=MED_CELL)
1048 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
1049 // Connectivity : we need MED_CELL connectivity !");
1053 /*! Not implemented yet */
1054 //--------------------------------------------------------------------//
1055 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1056 //--------------------------------------------------------------------//
1064 Give, for one element number of a specified entity the geometric type
1067 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35) ;
1069 //--------------------------------------------------------------------//
1070 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int Number) const
1071 //--------------------------------------------------------------------//
1073 if (_entity==Entity) {
1074 for(int i=1; i<=_numberOfTypes;i++)
1075 if (Number<_count[i])
1076 return _geometricTypes[i-1] ;
1078 else if (_constituent!=NULL)
1079 return _constituent->getElementType(Entity,Number) ;
1081 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !")) ;
1082 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !")) ;