1 #include "MEDMEM_Connectivity.hxx"
2 #include "MEDMEM_Family.hxx"
3 #include "MEDMEM_CellModel.hxx"
5 #include "MEDMEM_SkyLineArray.hxx"
6 #include "MEDMEM_ModulusArray.hxx"
8 #include "MEDMEM_STRING.hxx"
10 //------------------------------------------------------//
11 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity=MED_CELL):
12 //------------------------------------------------------//
14 _typeConnectivity(MED_NODAL),
16 _geometricTypes((medGeometryElement*)NULL),
17 _type((CELLMODEL*)NULL),
20 _nodal((MEDSKYLINEARRAY*)NULL),
21 _descending((MEDSKYLINEARRAY*)NULL),
22 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
23 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
24 _neighbourhood((MEDSKYLINEARRAY*)NULL),
25 _constituent((CONNECTIVITY*)NULL)
27 MESSAGE("CONNECTIVITY(medEntityMesh Entity=MED_CELL)") ;
30 //-------------------------------------------------------------------------//
31 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL):
32 //-------------------------------------------------------------------------//
34 _typeConnectivity(MED_NODAL),
35 _numberOfTypes(numberOfTypes),
37 _nodal((MEDSKYLINEARRAY*)NULL),
38 _descending((MEDSKYLINEARRAY*)NULL),
39 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
40 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
41 _neighbourhood((MEDSKYLINEARRAY*)NULL),
42 _constituent((CONNECTIVITY*)NULL)
44 MESSAGE("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)") ;
45 _geometricTypes = new medGeometryElement[numberOfTypes];
46 _type = new CELLMODEL[numberOfTypes];
47 _count = new int[numberOfTypes];
50 //----------------------------//
51 CONNECTIVITY::~CONNECTIVITY()
52 //----------------------------//
54 MESSAGE("Destructeur de CONNECTIVITY()") ;
55 if ( _geometricTypes != NULL )
56 delete [] _geometricTypes ;
61 if ( _descending != NULL )
63 if ( _reverseNodalConnectivity != NULL )
64 delete _reverseNodalConnectivity ;
65 if ( _reverseDescendingConnectivity != NULL )
66 delete _reverseDescendingConnectivity ;
67 if ( _constituent != NULL )
72 //------------------------------------------------------------------------------------------//
73 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
74 //------------------------------------------------------------------------------------------//
76 MESSAGE("CONNECTIVITY::calculateConnectivity");
78 // a temporary limitation due to calculteDescendingConnectivity function !!!
79 if ((_entityDimension==3) & (Entity==MED_EDGE))
80 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
83 if (ConnectivityType==MED_NODAL)
84 calculateNodalConnectivity() ;
87 calculateDescendingConnectivity() ;
89 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
90 if (Entity!=_entity) {
91 calculateDescendingConnectivity() ;
92 _constituent->calculateConnectivity(ConnectivityType,Entity) ;
96 /*! Give, in full or no interlace mode (for nodal connectivity),
97 descending or nodal connectivity.
99 You must give a <medEntityMesh> (ie:MED_EDGE)
100 and a <medGeometryElement> (ie:MED_SEG3). */
102 //------------------------------------------------------------//
103 void CONNECTIVITY::updateFamily(vector<FAMILY*> myFamilies)
104 //------------------------------------------------------------//
106 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
109 int numberOfFamilies = myFamilies.size();
110 if (numberOfFamilies == 0 ) {
111 MESSAGE(LOC<<"No family") ;
114 // does we do an update ?
115 if ((_constituent != NULL)&(_descending != NULL)) {
116 MESSAGE(LOC<<"Constituent is already defined") ;
120 if ((_constituent != NULL)&(_descending == NULL)) {
121 if (myFamilies[0]->getEntity() != _constituent->getEntity()) {
122 MESSAGE(LOC<<"Family and constituent entity are different. We do nothing") ;
126 // well we could go !
127 CONNECTIVITY * oldConstituent = _constituent ;
128 _constituent = (CONNECTIVITY *)NULL ;
129 // for instance we must have nodal connectivity in constituent :
130 if (oldConstituent->_nodal == NULL)
132 MESSAGE(LOC<<"We have no nodal connectivity of sub cell");
133 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell")) ;
135 int oldNumberOfFace = oldConstituent->_nodal->getNumberOf() ;
136 int * oldConstituentValue = oldConstituent->_nodal->getValue() ;
137 int * oldConstituentIndex = oldConstituent->_nodal->getIndex() ;
139 calculateDescendingConnectivity() ;
141 // if (oldConstituent->_nodal != NULL) {
142 int newNumberOfFace = _constituent->_nodal->getNumberOf() ;
143 int * newConstituentValue = _constituent->_nodal->getValue() ;
144 int * newConstituentIndex = _constituent->_nodal->getIndex() ;
146 int * newReverseDescendingIndex =
147 _reverseDescendingConnectivity->getIndex();
148 int * newReverseDescendingValue =
149 _reverseDescendingConnectivity->getValue();
151 int * newDescendingIndex = _descending->getIndex();
152 int * newDescendingValue = _descending->getValue();
154 // loop on all family,
155 // for all constituent in family, we get it's old connectivity
156 // (with oldCconstituentValue and oldConstituentIndex)
157 // and search the constituent in newConstituentValue with class
160 // when a new face is found, replace old constituent
161 // number in family with new one
162 // If face is not rigth oriented, we must change _descending attribute
163 // and _reverseDescendingConnectivity (see calculateDescendingConnectivity()).
165 // Voila a toi de jouer Nadir :-)
167 // First we get the renumbering from the oldCconstituentValue and
168 // oldConstituentIndex in the the new one, newConstituentValue and
169 // newConstituentIndex with a negative sign if the face is not
172 int * renumberingFromOldToNew = new int [oldNumberOfFace];
176 for (int iOldFace=0;iOldFace<oldNumberOfFace;iOldFace++)
180 // renumberingFromOldToNew[iOldFace] = 999999;
182 int face_it_beginOld = oldConstituentIndex[iOldFace];
183 int face_it_endOld = oldConstituentIndex[iOldFace+1];
184 int face_size_itOld = face_it_endOld - face_it_beginOld;
187 MEDMODULUSARRAY modulusArrayOld(face_size_itOld,oldConstituentValue+face_it_beginOld-1);
189 for (int iNewFace=0;iNewFace<newNumberOfFace && index == 0;
192 int face_it_beginNew = newConstituentIndex[iNewFace];
193 int face_it_endNew = newConstituentIndex[iNewFace+1];
194 face_size_itNew = face_it_endNew - face_it_beginNew;
196 if (face_size_itNew == face_size_itOld)
198 MEDMODULUSARRAY modulusArrayNew(face_size_itNew,newConstituentValue+face_it_beginNew-1);
200 int retCompareNewOld = modulusArrayNew.compare(modulusArrayOld);
202 // Real new face found
204 if(retCompareNewOld == 1)
206 renumberingFromOldToNew[iOldFace] = iNewFace+1;
211 // Reverse new face found
213 if(retCompareNewOld == -1)
215 renumberingFromOldToNew[iOldFace] = iNewFace+1;
219 int face_it_begin = newReverseDescendingIndex[iNewFace];
220 int face_it_end = newReverseDescendingIndex[iNewFace+1];
221 int face_size_it = face_it_end - face_it_begin;
223 if (face_size_it == 1)
224 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
226 if (face_size_it > 2)
227 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This face/edge should not be a (d-1) cell because it has "<<face_size_it<<" neighbouring"));
229 // we have always 2 neighbourings
230 int cell1 = newReverseDescendingValue[face_it_begin-1];
231 int cell2 = newReverseDescendingValue[face_it_begin];
233 // PROVISOIRE : en attendant que le SKYLINEARRAY de ReverseDescending soit correct (sans le zero)
235 // throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"This is a boundary face/edge, it should not be renumbered and it is wrongly oriented."));
237 if (cell2 != 0) { // we are not on border !!!!
239 newReverseDescendingValue[face_it_begin-1] = cell2;
240 // Updating _constituent->_nodal because of reversity
241 int * oldArray = oldConstituentValue+face_it_beginOld-1;
242 int * newArray = newConstituentValue+face_it_beginNew-1;
243 for(int iarray=0;iarray<face_size_itNew;iarray++)
244 newArray[iarray] = oldArray[iarray] ;
247 // Updating _reverseDescendingConnectivity
250 newReverseDescendingValue[face_it_begin] = cell1;
252 // Updating _descending for cell1 and cell2
253 for(int iface=newDescendingIndex[cell1-1];iface<newDescendingIndex[cell1];iface++)
254 if (newDescendingValue[iface-1]==iNewFace+1)
255 newDescendingValue[iface-1]=-iNewFace-1 ;
256 else if (newDescendingValue[iface-1]==-iNewFace-1)
257 newDescendingValue[iface-1]=iNewFace+1 ;
259 for(int iface=newDescendingIndex[cell2-1];iface<newDescendingIndex[cell2];iface++)
260 if (newDescendingValue[iface-1]==iNewFace+1)
261 newDescendingValue[iface-1]=-iNewFace-1 ;
262 else if (newDescendingValue[iface-1]==-iNewFace-1)
263 newDescendingValue[iface-1]=iNewFace+1 ;
264 } else {// else we are on border and we do nothing !!!!!!!!
265 INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
266 INFOS(LOC<<" Boudary FACE "<<iOldFace+1<<" are wrong oriented !") ;
267 INFOS("WARNING,WARNING,WARNING,WARNING,WARNING,WARNING");
275 INFOS(LOC<<"Renumbering problem with the Face connectivity given by the User and the new Connectivity computed");
276 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have a Face connectivity problem"));
280 MESSAGE(LOC<<"The Renumbering is finished and the status is");
284 // Updating the Family
286 for(int i=0; i<numberOfFamilies; i++) {
287 FAMILY * myFamily = myFamilies[i] ;
289 int length_skyline = myFamily->getnumber()->getLength();
290 int * value_skyline = myFamily->getnumber()->getValue();
292 for (int i=0;i<length_skyline;i++)
293 value_skyline[i] = renumberingFromOldToNew[value_skyline[i]-1];
301 //------------------------------------------------------------------------------------------------------------------//
302 med_int * CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type)
303 //------------------------------------------------------------------------------------------------------------------//
305 const char * LOC = "CONNECTIVITY::getConnectivity" ;
308 MEDSKYLINEARRAY * Connectivity ;
309 if (Entity==_entity) {
311 if (ConnectivityType==MED_NODAL)
314 Connectivity=_descending;
316 if (Connectivity!=NULL)
317 if (Type==MED_ALL_ELEMENTS)
318 return Connectivity->getValue() ;
320 for (med_int i=0; i<_numberOfTypes; i++)
321 if (_geometricTypes[i]==Type)
322 return Connectivity->getI(_count[i]) ;
323 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !")) ;
326 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !")) ;
328 if (_constituent != NULL)
329 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
331 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
334 /*! Give morse index array to use with
335 getConnectivity(MED_FULL_INTERLACE,mode,entity,MED_ALL_ELEMENTS).
337 Each value give start index for corresponding entity in connectivity array.
339 Example : i-th element, j-th node of it :
340 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
341 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
342 //-----------------------------------------------------------------------------------------------//
343 med_int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType, medEntityMesh Entity)
344 //------------------------------------------------------------------------------------------------//
346 const char * LOC = "CONNECTIVITY::getConnectivityIndex" ;
349 MEDSKYLINEARRAY * Connectivity ;
350 if (Entity==_entity) {
352 if (ConnectivityType==MED_NODAL)
355 Connectivity=_descending;
357 if (Connectivity!=NULL)
358 return Connectivity->getIndex() ;
360 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !")) ;
363 if (_constituent != NULL)
364 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
366 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
370 //--------------------------------------------------------------//
371 CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
372 //--------------------------------------------------------------//
374 const char * LOC = "CONNECTIVITY::getType" ;
377 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
378 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !")) ;
379 for (med_int i=0; i<_numberOfTypes; i++)
380 if (_geometricTypes[i]==Type)
382 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !")) ;
385 /*! Returns the number of elements of type <med geometrie element>.
386 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
387 //------------------------------------------------------------------------//
388 med_int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
389 //------------------------------------------------------------------------//
391 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType" ;
394 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
395 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!")) ;
396 for (med_int i=0; i<_numberOfTypes; i++)
397 if (_geometricTypes[i]==Type)
398 return _type[i].getNumberOfNodes() ;
399 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !")) ;
402 /*! Returns the number of geometric sub cells of <med geometrie element> type.
403 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
404 //------------------------------------------------------------------------//
405 med_int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
406 //------------------------------------------------------------------------//
408 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
409 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!") ;
410 for (med_int i=0; i<_numberOfTypes; i++)
411 if (_geometricTypes[i]==Type)
412 return _type[i].getNumberOfConstituents(1) ;
413 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !") ;
416 /*! Returns the number of elements of type <med geometrie element>.
419 - Implemented for MED_ALL_ELEMENTS
420 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
421 - Not implemented for MED_NONE */
422 //-----------------------------------------------------------------------------------//
423 med_int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
424 //-----------------------------------------------------------------------------------//
426 const char * LOC = "CONNECTIVITY::getNumberOf" ;
429 MESSAGE(LOC<<" Entity = "<< Entity << ", _entity = "<<_entity);
431 if (Entity==_entity) {
433 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
434 if (Type==MED_ALL_ELEMENTS)
435 return _count[_numberOfTypes]-1;
436 for (med_int i=0; i<_numberOfTypes; i++)
437 if (_geometricTypes[i]==Type)
438 return (_count[i+1] - _count[i]);
439 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
441 if (_constituent != NULL)
442 return _constituent->getNumberOf(Entity,Type);
444 return 0 ; // valid if they are nothing !
445 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !")) ;
449 //--------------------------------------------------------------//
450 med_int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
451 medGeometryElement Type)
452 //--------------------------------------------------------------//
454 if (TypeConnectivity == MED_NODAL)
456 calculateNodalConnectivity();
457 if (Type==MED_ALL_ELEMENTS)
458 return _nodal->getValue();
459 for (med_int i=0; i<_numberOfTypes; i++)
460 if (_geometricTypes[i]==Type)
461 return _nodal->getI(_count[i]);
465 calculateDescendingConnectivity();
466 if (Type==MED_ALL_ELEMENTS)
467 return _descending->getValue();
468 for (med_int i=0; i<_numberOfTypes; i++)
469 if (_geometricTypes[i]==Type)
470 return _descending->getI(Type);
472 throw MEDEXCEPTION("Not found");
476 //---------------------------------------------------------------------//
477 med_int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity)
478 //---------------------------------------------------------------------//
480 if (TypeConnectivity == MED_NODAL)
482 calculateNodalConnectivity();
483 return _nodal->getIndex();
487 calculateDescendingConnectivity();
488 return _descending->getIndex();
492 /*! Not yet implemented */
493 //----------------------------------------------//
494 med_int* CONNECTIVITY:: getNeighbourhood() const
495 //----------------------------------------------//
497 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
500 /*! Returns an array which contains, for each node, all cells
502 //-------------------------------------------------//
503 med_int* CONNECTIVITY::getReverseNodalConnectivity()
504 //-------------------------------------------------//
506 calculateReverseNodalConnectivity();
507 return _reverseNodalConnectivity->getValue();
510 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
511 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
512 //-------------------------------------------------------//
513 med_int* CONNECTIVITY::getReverseNodalConnectivityIndex()
514 //-------------------------------------------------------//
516 calculateReverseNodalConnectivity();
517 return _reverseNodalConnectivity->getIndex();
520 /*! Returns an array which contains, for each face (or edge),
521 the 2 cells of each side. First is cell which face normal is outgoing.
523 //------------------------------------------------------//
524 med_int* CONNECTIVITY::getReverseDescendingConnectivity()
525 //------------------------------------------------------//
527 // it is in _constituent connectivity only if we are in MED_CELL
528 if (_entity==MED_CELL) {
529 // we want descending connectivity
530 calculateDescendingConnectivity();
531 return _reverseDescendingConnectivity->getValue();
533 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
536 /*! calculate the reverse descending Connectivity
537 and returns the index ( A DOCUMENTER MIEUX)*/
538 //-----------------------------------------------------------//
539 med_int* CONNECTIVITY::getReverseDescendingConnectivityIndex()
540 //-----------------------------------------------------------//
542 // it is in _constituent connectivity only if we are in MED_CELL
543 if (_entity==MED_CELL) {
544 // we want descending connectivity
545 calculateDescendingConnectivity();
546 return _reverseDescendingConnectivity->getIndex();
548 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
551 /*! A DOCUMENTER (et a finir ???) */
552 //--------------------------------------------//
553 void CONNECTIVITY::calculateNodalConnectivity()
554 //--------------------------------------------//
558 if (_descending==NULL)
559 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
560 // calculate _nodal from _descending
564 /*! If not yet done, calculate the nodal Connectivity
565 and the reverse nodal Connectivity*/
566 //---------------------------------------------------//
567 void CONNECTIVITY::calculateReverseNodalConnectivity()
568 //---------------------------------------------------//
570 const char * LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : " ;
574 calculateNodalConnectivity() ;
576 MESSAGE(LOC<<"Number of nodes = "<<_numberOfNodes);
578 if(_reverseNodalConnectivity==NULL) {
580 med_int node_number = 0;
581 vector <vector <med_int> > reverse_connectivity;
582 reverse_connectivity.resize(_numberOfNodes+1);
584 // Treat all cells types
586 for (med_int j = 0; j < _numberOfTypes; j++)
588 // node number of the cell type
589 node_number = _type[j].getNumberOfNodes();
590 // treat all cells of a particular type
591 for (med_int k = _count[j]; k < _count[j+1]; k++)
592 // treat all nodes of the cell type
593 for (med_int local_node_number = 1 ; local_node_number < node_number+1; local_node_number++)
595 med_int global_node = _nodal->getIJ(k,local_node_number) ;
596 reverse_connectivity[global_node].push_back(k);
600 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
602 //calculate size of reverse_nodal_connectivity array
603 med_int size_reverse_nodal_connectivity = 0;
604 for (med_int i = 1; i < _numberOfNodes+1; i++)
605 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
607 MEDSKYLINEARRAY * ReverseConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity) ;
608 // reverse_nodal_connectivity_index = new (med_int)[_numberOfNodes+1];
609 // reverse_nodal_connectivity = new (med_int)[size_reverse_nodal_connectivity];
610 med_int * reverse_nodal_connectivity = ReverseConnectivity->getValue() ;
611 med_int * reverse_nodal_connectivity_index = ReverseConnectivity->getIndex() ;
613 reverse_nodal_connectivity_index[0] = 1;
614 for (med_int i = 1; i < _numberOfNodes+1; i++)
616 med_int size = reverse_connectivity[i].size();
617 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size ;
618 for (med_int j = 0; j < size; j++)
619 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
622 _reverseNodalConnectivity = ReverseConnectivity ;
627 /*! If not yet done, calculate the Descending Connectivity */
628 //-------------------------------------------------//
629 void CONNECTIVITY::calculateDescendingConnectivity()
630 //-------------------------------------------------//
632 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
635 if (_descending==NULL) {
637 MESSAGE(LOC<<"No connectivity found !");
638 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
640 // calcul _descending from _nodal
641 // we need CONNECTIVITY for constituent
643 if (_constituent != NULL)
644 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"ERROR : No descending connectivity defined, but constituent exist !"));
645 MESSAGE(LOC<<": No descending connectivity defined, but constituent exist !");
647 if (_entityDimension == 3)
648 _constituent = new CONNECTIVITY(MED_FACE) ;
649 else if (_entityDimension == 2)
650 _constituent = new CONNECTIVITY(MED_EDGE) ;
652 MESSAGE(LOC<<"We are in 1D");
655 _constituent->_typeConnectivity = MED_DESCENDING ;
656 _constituent->_numberOfNodes = _numberOfNodes ;
657 // foreach cells, we built array of constituent
658 int DescendingSize = 0 ;
659 for(int i=0; i<_numberOfTypes; i++)
660 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1) ;
661 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,DescendingSize) ;
662 int * descend_connectivity = _descending->getValue() ;
663 for (int i=0; i<DescendingSize; i++)
664 descend_connectivity[i]=0;
665 int * descend_connectivity_index = _descending->getIndex() ;
666 descend_connectivity_index[0]=1;
667 medGeometryElement* ConstituentsTypes = new medGeometryElement[2];
668 ConstituentsTypes[0]=MED_NONE ;
669 ConstituentsTypes[1]=MED_NONE ;
670 int * NumberOfConstituentsForeachType = new int[2];
671 NumberOfConstituentsForeachType[0]=0;
672 NumberOfConstituentsForeachType[1]=0;
673 for(int i=0; i<_numberOfTypes; i++) {
674 // initialize descend_connectivity_index array :
675 int NumberOfConstituents = _type[i].getNumberOfConstituents(1) ;
676 for (int j=_count[i];j<_count[i+1];j++){
677 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents ;
678 // compute number of constituent of all cell for each type
679 for(int k=1;k<NumberOfConstituents+1;k++) {
680 medGeometryElement MEDType = _type[i].getConstituentType(1,k) ;
681 if (ConstituentsTypes[0]==MED_NONE) {
682 ConstituentsTypes[0]=MEDType;
683 NumberOfConstituentsForeachType[0]++ ;
684 } else if (ConstituentsTypes[0]==MEDType)
685 NumberOfConstituentsForeachType[0]++ ;
686 else if (ConstituentsTypes[1]==MED_NONE) {
687 ConstituentsTypes[1]=MEDType;
688 NumberOfConstituentsForeachType[1]++ ;
689 } else if (ConstituentsTypes[1]==MEDType)
690 NumberOfConstituentsForeachType[1]++ ;
692 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<MEDType<<" is different of "<<ConstituentsTypes[0]<<" and "<<ConstituentsTypes[1]));
697 // we could built _constituent !
698 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1] ;
699 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100) ;
700 _constituent->_nodal = new MEDSKYLINEARRAY(TotalNumberOfConstituents,TotalNumberOfNodes) ;
702 // we use _constituent->_nodal
703 int * ConstituentNodalConnectivity = _constituent->_nodal->getValue();
704 int * ConstituentNodalConnectivityIndex = _constituent->_nodal->getIndex();
705 ConstituentNodalConnectivityIndex[0]=1;
707 _constituent->_entityDimension=ConstituentsTypes[0]/100;
708 if (ConstituentsTypes[1]==MED_NONE)
709 _constituent->_numberOfTypes = 1 ;
711 _constituent->_numberOfTypes = 2 ;
712 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes] ;
713 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes] ;
714 _constituent->_count = new med_int[_constituent->_numberOfTypes+1] ;
715 _constituent->_count[0]=1;
716 int* tmp_NumberOfConstituentsForeachType = new med_int[_constituent->_numberOfTypes+1] ;
717 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
718 for (int i=0; i<_constituent->_numberOfTypes;i++) {
719 _constituent->_geometricTypes[i]=ConstituentsTypes[i] ;
720 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i] ;
721 CELLMODEL Type(ConstituentsTypes[i]);
722 _constituent->_type[i]=Type;
723 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i] ;
724 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
725 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100) ;
727 delete[] ConstituentsTypes;
728 delete[] NumberOfConstituentsForeachType ;
730 // we need reverse nodal connectivity
731 if (! _reverseNodalConnectivity)
732 calculateReverseNodalConnectivity() ;
733 int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
734 int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
736 // array to keep reverse descending connectivity
737 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
739 int TotalNumberOfSubCell = 0 ;
740 for (int i=0; i<_numberOfTypes; i++) { // loop on all cell type
741 CELLMODEL Type = _type[i] ;
742 int NumberOfNodesPerCell = Type.getNumberOfNodes() ;
743 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
744 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
745 for (int k=1 ; k<=NumberOfConstituentPerCell; k++) { // we loop on all sub cell of it
746 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0) {
747 // it is a new sub cell !
748 // TotalNumberOfSubCell++;
750 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0]){
751 tmp_NumberOfConstituentsForeachType[0]++;
752 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
754 tmp_NumberOfConstituentsForeachType[1]++;
755 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
757 //we have maximum two types
759 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
760 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j ;
761 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100 ;
763 int * NodesLists = new int[NumberOfNodesPerConstituent] ;
764 for (int l=0; l<NumberOfNodesPerConstituent; l++) {
765 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
766 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
768 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
770 // all elements which contains first node of sub cell :
771 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1] ;
772 int ReverseNodalConnectivityIndex_1 = ReverseNodalConnectivityIndex[NodesLists[0]] ;
773 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0 ;
774 int * CellsList = new int[NumberOfCellsInList] ;
775 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
776 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1] ;
778 // foreach node in sub cell, we search elements which are in common
779 // at the end, we must have only one !
781 for (int l=1; l<NumberOfNodesPerConstituent; l++) {
782 int NewNumberOfCellsInList = 0 ;
783 int * NewCellsList = new int[NumberOfCellsInList] ;
784 for (int l1=0; l1<NumberOfCellsInList; l1++)
785 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++) {
786 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
787 // increasing order : CellsList[l1] are not in elements list of node l
789 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
791 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
792 NewNumberOfCellsInList++;
796 NumberOfCellsInList = NewNumberOfCellsInList;
797 delete [] CellsList ;
798 CellsList = NewCellsList;
801 int CellNumber = CellsList[0] ;
802 delete [] CellsList ;
803 if (NumberOfCellsInList>1)
804 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
807 // int cell_number_1=-1 ;
808 // int cell_number_2=-1 ;
809 // int cell_number_3=-1 ;
810 // bool find = false ;
811 // for (int l=ReverseNodalConnectivityIndex[NodesLists[0]-1]; l<ReverseNodalConnectivityIndex[NodesLists[0]]; l++) { // first node
812 // cell_number_1 = ReverseNodalConnectivityValue[l-1] ;
813 // if (cell_number_1 != j)
814 // for (int m=ReverseNodalConnectivityIndex[NodesLists[1]-1]; m<ReverseNodalConnectivityIndex[NodesLists[1]]; m++) { //second node
815 // cell_number_2 = ReverseNodalConnectivityValue[m-1] ;
816 // if ((cell_number_2 != j) && (cell_number_2 == cell_number_1))
817 // for (int n=ReverseNodalConnectivityIndex[NodesLists[2]-1]; n<ReverseNodalConnectivityIndex[NodesLists[2]]; n++) { //third node
818 // cell_number_3 = ReverseNodalConnectivityValue[n-1] ;
819 // if ((cell_number_3 != j) && (cell_number_3 == cell_number_1)) { // we found element which have three node in it
835 if (NumberOfCellsInList==1) {
836 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber ;
837 // we search sub cell number in this cell to not calculate it another time
840 for (int l=0; l<_numberOfTypes; l++)
841 if (CellNumber < _count[l+1]) {
845 //int sub_cell_count2 = Type2.get_entities_count(1) ;
846 //int nodes_cell_count2 = Type2.get_nodes_count() ;
848 for (int l=1; l<=Type2.getNumberOfConstituents(1) ;l++) { // on all sub cell
850 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
851 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++) {
852 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
855 if (counter==Type.getConstituentType(1,k)%100) {
856 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
863 INFOS(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!") ; // exception ?
865 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0 ;
867 delete[] NodesLists ;
871 // we adjust _constituent
872 int NumberOfConstituent=0 ;
873 int SizeOfConstituentNodal=0 ;
874 for (int i=0;i<_constituent->_numberOfTypes; i++) {
875 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
876 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
878 // we built new _nodal attribute in _constituent
879 MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
880 int *ConstituentNodalValue = ConstituentNodal->getValue();
881 int *ConstituentNodalIndex = ConstituentNodal->getIndex();
882 ConstituentNodalIndex[0]=1;
883 // we build _reverseDescendingConnectivity
884 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent) ;
885 int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
886 int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
887 reverseDescendingConnectivityIndex[0]=1;
889 // first constituent type
890 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++) {
891 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
892 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
893 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
895 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
896 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
897 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
900 // second type if any
901 if (_constituent->_numberOfTypes==2) {
902 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1 ;
903 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
904 int offset2=offset*2 ;
905 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
906 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++) {
907 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
908 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++){
909 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
911 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
912 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++){
913 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
916 _constituent->_count[2]=NumberOfConstituent+1 ;
917 // we correct _descending to adjust face number
918 for(int j=0;j<DescendingSize;j++)
919 if (descend_connectivity[j]>tmp_NumberOfConstituentsForeachType[0])
920 descend_connectivity[j]-=offset;
923 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
924 delete _constituent->_nodal ;
925 _constituent->_nodal = ConstituentNodal ;
927 delete[] ReverseDescendingConnectivityValue ;
932 //-----------------------------------------------------------------------------------------//
933 // void CONNECTIVITY::calculateReverseDescendingConnectivity(CONNECTIVITY *myConnectivity)
934 //-----------------------------------------------------------------------------------------//
936 // if (_entity==MED_CELL)
937 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
938 // Connectivity : we could not do that with MED_CELL entity !");
939 // if (myConnectivity->getEntity()!=MED_CELL)
940 // throw MEDEXCEPTION("CONNECTIVITY::calculateReverseDescending
941 // Connectivity : we need MED_CELL connectivity !");
945 /*! Not implemented yet */
946 //--------------------------------------------------------------------//
947 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
948 //--------------------------------------------------------------------//
956 Give, for one element number of a specified entity the geometric type
959 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35) ;
961 //--------------------------------------------------------------------//
962 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int Number) const
963 //--------------------------------------------------------------------//
965 if (_entity==Entity) {
966 for(int i=1; i<=_numberOfTypes;i++)
967 if (Number<_count[i])
968 return _geometricTypes[i-1] ;
970 else if (_constituent!=NULL)
971 return _constituent->getElementType(Entity,Number) ;
973 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !")) ;
974 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !")) ;