1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "MEDMEM_Connectivity.hxx"
24 #include "MEDMEM_Family.hxx"
25 #include "MEDMEM_Group.hxx"
26 #include "MEDMEM_CellModel.hxx"
28 #include "MEDMEM_SkyLineArray.hxx"
29 #include "MEDMEM_ModulusArray.hxx"
31 #include "MEDMEM_STRING.hxx"
32 #include "MEDMEM_Utilities.hxx"
36 using namespace MEDMEM;
37 using namespace MED_EN;
39 // Enlarge the vector if necessary to insert the element
40 static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
42 if (Indice >=(int) Vect.capacity())
43 Vect.reserve(Indice + 1000);
45 if (Indice >=(int) Vect.size())
46 Vect.resize(Indice+1);
48 Vect[Indice] = Element;
51 void mergeOrderedTabs(const int *tab1, int lgth1, const int *tab2, int lgth2, int *result, int& lgth);
53 void mergeOrderedTabs(const int *tab1, int lgth1, const int *tab2, int lgth2, int *result, int& lgth)
57 unsigned char switcher=0;
58 const int *tabS[2]={tab1,tab2};
59 while(cpt[0]<lgth1 && cpt[1]<lgth2)
61 if(tabS[1-switcher][cpt[1-switcher]]<tabS[switcher][cpt[switcher]])
63 else if(tabS[1-switcher][cpt[1-switcher]]>tabS[switcher][cpt[switcher]])
67 int tmp=tabS[switcher][cpt[switcher]];
68 cpt[switcher]++; cpt[1-switcher]++;
75 Default Constructor. /n
76 Default for Entity is MED_CELL and type of Connectivity is MED_NODAL */
77 //--------------------------------------------------------------//
78 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
79 //--------------------------------------------------------------//
81 _typeConnectivity(MED_NODAL),
83 _geometricTypes((medGeometryElement*)NULL),
84 _type((CELLMODEL*)NULL),
88 _nodal((MEDSKYLINEARRAY*)NULL),
89 _descending((MEDSKYLINEARRAY*)NULL),
90 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
91 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
92 _neighbourhood((MEDSKYLINEARRAY*)NULL),
93 _constituent((CONNECTIVITY*)NULL),
94 _isDescendingConnectivityPartial(false)
96 const char* LOC = "CONNECTIVITY(medEntityMesh Entity=MED_CELL)";
98 MESSAGE_MED("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
106 Default for Entity is MED_CELL */
107 //------------------------------------------------------------------------------//
108 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
109 //------------------------------------------------------------------------------//
111 _typeConnectivity(MED_NODAL),
112 _numberOfTypes(numberOfTypes),
115 _nodal((MEDSKYLINEARRAY*)NULL),
116 _descending((MEDSKYLINEARRAY*)NULL),
117 _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
118 _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
119 _neighbourhood((MEDSKYLINEARRAY*)NULL),
120 _constituent((CONNECTIVITY*)NULL),
121 _isDescendingConnectivityPartial(false)
123 MESSAGE_MED("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
124 _geometricTypes = new medGeometryElement[numberOfTypes];
125 _type = new CELLMODEL[numberOfTypes];
126 _count = new int[numberOfTypes+1];
129 _count[ numberOfTypes-1 ] = 0; // to know whether _count is set or not
135 //------------------------------------------------------//
136 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
137 //----------------------------------------------------//
139 _typeConnectivity (m._typeConnectivity),
140 _numberOfTypes (m._numberOfTypes),
141 _entityDimension (m._entityDimension),
142 _numberOfNodes (m._numberOfNodes),
143 _isDescendingConnectivityPartial(m._isDescendingConnectivityPartial)
145 if (m._geometricTypes != NULL)
147 _geometricTypes = new medGeometryElement[_numberOfTypes];
148 memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
151 _geometricTypes = (medGeometryElement *) NULL;
155 _type = new CELLMODEL[_numberOfTypes];
156 for (int i=0;i<_numberOfTypes;i++)
157 _type[i] = CELLMODEL(m._type[i]);
160 _type = (CELLMODEL *) NULL;
162 if (m._count != NULL)
164 _count = new int[_numberOfTypes+1];
165 memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
168 _count = (int *) NULL;
170 if (m._nodal != NULL)
171 _nodal = new MEDSKYLINEARRAY(* m._nodal);
173 _nodal = (MEDSKYLINEARRAY *) NULL;
175 if (m._descending != NULL)
176 _descending = new MEDSKYLINEARRAY(* m._descending);
178 _descending = (MEDSKYLINEARRAY *) NULL;
180 if (m._reverseNodalConnectivity != NULL)
181 _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
183 _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
185 if (m._reverseDescendingConnectivity != NULL)
186 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
188 _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
190 if (m._neighbourhood != NULL)
191 _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
193 _neighbourhood = (MEDSKYLINEARRAY *) NULL;
195 if (m._constituent != NULL)
196 _constituent = new CONNECTIVITY(* m._constituent);
198 _constituent = (CONNECTIVITY *) NULL;
203 desallocates existing pointers */
204 //----------------------------//
205 CONNECTIVITY::~CONNECTIVITY()
206 //----------------------------//
208 MESSAGE_MED("Destructeur de CONNECTIVITY()");
210 if (_geometricTypes != NULL)
211 delete [] _geometricTypes;
218 if (_descending != NULL)
220 if (_reverseNodalConnectivity != NULL)
221 delete _reverseNodalConnectivity;
222 if (_reverseDescendingConnectivity != NULL)
223 delete _reverseDescendingConnectivity;
224 if (_constituent != NULL)
229 set _constituent to Constituent
230 be aware desallocation of _constituent is done by CONNECTIVITY:~CONNECTIVITY
231 throws an exception if Constituent = MED_CELL
234 //----------------------------------------------------------//
235 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
237 //----------------------------------------------------------//
239 medEntityMesh Entity = Constituent->getEntity();
240 if (Entity == MED_CELL)
241 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
243 if ((Entity == MED_EDGE)&(_entityDimension == 3))
245 if (_constituent == NULL)
246 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
247 _constituent->setConstituent(Constituent);
252 _constituent = Constituent;
257 Return _constituent of given entity
259 //----------------------------------------------------------//
260 const CONNECTIVITY *CONNECTIVITY::getConstituent(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
261 //----------------------------------------------------------//
263 if(Entity==MED_FACE && _entityDimension==3)
265 if(Entity==MED_EDGE && _entityDimension==2)
267 if(Entity==MED_EDGE && _entityDimension==3)
270 return _constituent->getConstituent(MED_EDGE);
272 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getConstituent : FACE connectivity not defined!"));
274 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getConstituent : Invalid request !"));
277 /*! Duplicated Types array in CONNECTIVITY object. */
278 //--------------------------------------------------------------------//
279 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
280 const medEntityMesh Entity)
282 //--------------------------------------------------------------------//
284 if (Entity == _entity)
285 for (int i=0; i<_numberOfTypes; i++)
287 _geometricTypes[i] = Types[i];
288 _type[i] = CELLMODEL(_geometricTypes[i]);
289 if ( _type[i].getDimension() > _entityDimension )
290 _entityDimension = _type[i].getDimension();
294 if (_constituent == NULL)
295 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
296 _constituent->setGeometricTypes(Types,Entity);
300 /*! Set number of elements of each type of given entity */
301 //--------------------------------------------------------------------//
302 void CONNECTIVITY::setCount(const int * Count,
303 const medEntityMesh Entity)
305 //--------------------------------------------------------------------//
307 if (Entity == _entity)
309 // this exception was added at 1.18.2.2.2.4 for NPAL17043: "Correction of
310 // MEDMEM CPPUNIT tests. Integrated a work of V. Bergeaud and A. Geay."
311 // and removed for PAL19744: "regression in MEDMEM_Connectivity.cxx"
312 // Commenting this exception at least looks safe as the case
313 // _numberOfTypes==0 is previewed here
314 //if (_numberOfTypes==0)
315 // throw MEDEXCEPTION("Number of Types was not set before setting counts");
316 int * index = new int[Count[_numberOfTypes]];
319 for (int i=0; i<_numberOfTypes; i++)
321 _count[i+1] = Count[i+1];
322 int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
323 if ( _geometricTypes[i] != MED_POLYGON && _geometricTypes[i] != MED_POLYHEDRA )
324 for (int j=_count[i]; j<_count[i+1]; j++)
325 index[j] = index[j-1]+NumberOfNodesPerElement;
327 index[_count[_numberOfTypes]-1] = index[_count[_numberOfTypes-1]-1];
330 if (_nodal != NULL) delete _nodal;
331 if (_numberOfTypes != 0)
333 _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
334 _nodal->setIndex(index);
342 if (_constituent == NULL)
343 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
344 _constituent->setCount(Count,Entity);
348 //--------------------------------------------------------------------//
349 void CONNECTIVITY::setNodal(const int * Connectivity,
350 const medEntityMesh Entity,
351 const medGeometryElement Type,
352 const int * PolyConnectivityIndex)
354 //--------------------------------------------------------------------//
356 if (_entity == Entity)
358 // find geometric type
360 for (int i=0; i<_numberOfTypes; i++)
361 if (_geometricTypes[i] == Type)
363 if ( _count[i+1] == 0 )
364 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : "
365 "number of elements by type must be set before!"));
366 if ( _count[i+1] < _count[i] )
367 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : "
368 "invalid number of elements by type!"));
371 if ( Type == MED_POLYGON || Type == MED_POLYHEDRA )
373 if (! PolyConnectivityIndex )
374 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : PolyConnectivityIndex must be provided for poly- elements !"));
375 if ( _numberOfTypes == 1 )
378 _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
379 PolyConnectivityIndex[_count[_numberOfTypes]-1]-1,
380 PolyConnectivityIndex, Connectivity);
384 // append poly connectivity to the present one
385 int nbpoly = getNumberOf( Entity, Type );
386 int polyconnlen = PolyConnectivityIndex[nbpoly] - PolyConnectivityIndex[0];
387 int newconnlen = _nodal->getLength() + polyconnlen;
388 // create new MEDSKYLINEARRAY
389 int * newconn = new int[ newconnlen ];
390 int * newindex = new int[ _nodal->getNumberOf() + 1 ];
391 MEDSKYLINEARRAY* newnodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
392 newconnlen,newindex,newconn,
393 /*shallowCopy=*/true);
395 memcpy(newconn, _nodal->getValue(),sizeof(int)*(_nodal->getLength()));
396 memcpy(newindex,_nodal->getIndex(),sizeof(int)*(_nodal->getNumberOf() - nbpoly + 1));
398 memcpy(newconn+_nodal->getLength(),Connectivity, sizeof(int)*polyconnlen);
399 int* polyindex = newindex + _count[i] - 1;
400 int delta = polyindex[0] - PolyConnectivityIndex[0];
401 for( int j=0;j<nbpoly; j++)
402 polyindex[j+1] = PolyConnectivityIndex[j+1] + delta;
409 const int* index = _nodal->getIndex();
410 for( int j=_count[i];j<_count[i+1]; j++)
411 _nodal->setI(j,Connectivity+index[j]-index[_count[i]]);
415 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
419 if (_constituent == NULL)
420 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
421 _constituent->setNodal(Connectivity,Entity,Type,PolyConnectivityIndex);
426 //------------------------------------------------------------------------------------------//
427 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
428 //------------------------------------------------------------------------------------------//
430 MESSAGE_MED("CONNECTIVITY::calculateConnectivity");
432 // a temporary limitation due to calculteDescendingConnectivity function !!!
433 if ((_entityDimension==3) & (Entity==MED_EDGE))
434 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
438 if (ConnectivityType==MED_NODAL)
439 calculateNodalConnectivity();
441 if (Entity==MED_CELL)
442 calculateDescendingConnectivity();
444 throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
446 if (Entity!=_entity) {
447 calculateDescendingConnectivity();
448 if (_entityDimension == 2 || _entityDimension == 3)
449 _constituent->calculateConnectivity(ConnectivityType,Entity);
453 //------------------------------------------------------------//
454 void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
455 //------------------------------------------------------------//
457 const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
458 int numberOfFamilies = myFamilies.size();
459 if (numberOfFamilies == 0 || _constituent == NULL)
461 // does we do an update ?
462 if ((_constituent != NULL) && (_descending != NULL))
464 if (myFamilies[0]->getEntity() != _constituent->getEntity())
467 CONNECTIVITY * oldConstituent = _constituent;
468 _constituent = (CONNECTIVITY *)NULL;
469 if (oldConstituent->_nodal==NULL)
470 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
472 const int oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
473 const int * oldConstituentValue = oldConstituent->_nodal->getValue();
474 const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
476 int * renumberingFromOldToNew = new int [oldNumberOfFace+1];
478 calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
479 _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to get new face numbers from nodes numbers...
481 const int * reverseFaceNodal = _constituent->getReverseNodalConnectivity();
482 const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
484 CELLMODEL * aCELLMODEL = &oldConstituent->_type[0];
485 //int newNumberOfFaces = _constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS);
486 for ( int iOldFace = 0; iOldFace < oldNumberOfFace; iOldFace++ )
488 // Retrieving new number of iOldFace face...
490 const int *nodesOfCurrentFaceOld=oldConstituentValue+oldConstituentIndex[iOldFace]-1;
491 int nbOfNodesOfCurrentFaceOld=oldConstituentIndex[iOldFace+1]-oldConstituentIndex[iOldFace];
493 // get all faces (newFaceNb) around the first node of iOldFace
494 int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
495 int *newFaceNb=new int[sizeOfNewFaceNb];
496 memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
498 // loop on the rest nodes of iOldFace to find, among newFaceNb, a sole face including all nodes
499 for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
501 const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
502 int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
503 for(int i=0; i<sizeOfNewFaceNb; )
506 for(int j=0; j < sizeOfNewFacesNbForCurNode && !found; j++)
507 found = (newFacesNbForCurNode[j]==newFaceNb[i]);
512 newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb]; // set the last face at place of not found one
515 if(sizeOfNewFaceNb!=1)
516 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
517 renumberingFromOldToNew[iOldFace]=newFaceNb[0];
519 //end of retrieving a new number of a face...
521 // Compare nodes of a new face and those of an old one;
522 // for the second order elements, only vertex nodes are compared
523 int nbOfNodesOfCurrentFaceNew;
524 const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElement(MED_NODAL,_constituent->getEntity(), renumberingFromOldToNew[iOldFace], nbOfNodesOfCurrentFaceNew );
525 if ( aCELLMODEL && aCELLMODEL->getNumberOfNodes() != nbOfNodesOfCurrentFaceOld )
527 // type changed, find a corresponding CELLMODEL
528 int iType = 2; // 1-st type is already used at loop beginning
529 while ( iOldFace + 1 >= oldConstituent->_count[ iType ]) // check next type
531 if ( oldConstituent->_numberOfTypes >= iType &&
532 oldConstituent->_type[ iType-1 ].getNumberOfNodes() > 0 )
533 aCELLMODEL = &oldConstituent->_type[ iType-1 ];
537 int nbOfVertices = aCELLMODEL ? aCELLMODEL->getNumberOfVertexes() : nbOfNodesOfCurrentFaceOld;
538 MEDMODULUSARRAY modulusArrayOld(nbOfVertices,nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
539 int nbOfVerticesNew = std::min( nbOfVertices, nbOfNodesOfCurrentFaceNew );
540 MEDMODULUSARRAY modulusArrayNew(nbOfVerticesNew,nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
541 int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
542 if(retCompareNewOld==0)
543 throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
544 if(retCompareNewOld==-1)
545 invertConnectivityForAFace(renumberingFromOldToNew[iOldFace],nodesOfCurrentFaceOld);
548 // Updating the Family
549 for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
550 (*iter)->changeElementsNbs(_constituent->getEntity(), renumberingFromOldToNew);
553 if ( _constituent && !_constituent->_constituent ) {
554 _constituent->_constituent = oldConstituent->_constituent;
555 oldConstituent->_constituent = NULL;
559 delete oldConstituent;
560 delete [] renumberingFromOldToNew;
564 /*! Give descending or nodal connectivity.
566 You must give a %medEntityMesh (ie:MED_EDGE)
567 and a %medGeometryElement (ie:MED_SEG3). */
568 //---------------------------------------------------------------------------------//
569 const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity ConnectivityType,
570 medEntityMesh Entity,
571 medGeometryElement Type) const
572 //----------------------------------------------------------------------------------//
574 const char * LOC = "CONNECTIVITY::getConnectivity";
575 // BEGIN_OF_MED(LOC);
577 MEDSKYLINEARRAY * Connectivity;
578 if (Entity==_entity) {
580 if (ConnectivityType==MED_NODAL)
582 calculateNodalConnectivity();
587 calculateDescendingConnectivity();
588 Connectivity=_descending;
591 if (Connectivity!=NULL)
592 if (Type==MED_ALL_ELEMENTS)
593 return Connectivity->getValue();
595 for (int i=0; i<_numberOfTypes; i++)
596 if (_geometricTypes[i]==Type)
597 //return Connectivity->getI(i+1);
598 return Connectivity->getI(_count[i]);
599 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
602 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
604 if (_constituent != NULL)
605 return _constituent->getConnectivity(ConnectivityType,Entity,Type);
607 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
610 //--------------------------------------------------------------------------
611 int CONNECTIVITY::getConnectivityLength(medConnectivity ConnectivityType,
612 medEntityMesh Entity,
613 medGeometryElement Type) const
614 //--------------------------------------------------------------------------
616 const char * LOC = "CONNECTIVITY::getConnectivity";
619 MEDSKYLINEARRAY * Connectivity;
620 if (Entity == _entity) {
621 if (ConnectivityType == MED_NODAL)
623 calculateNodalConnectivity();
624 Connectivity = _nodal;
628 calculateDescendingConnectivity();
629 Connectivity = _descending;
632 if (Connectivity != NULL) {
633 if (Type == MED_ALL_ELEMENTS) {
634 return Connectivity->getLength();
637 for (int i=0; i<_numberOfTypes; i++)
638 if (_geometricTypes[i]==Type)
640 //return (_count[i+1]-_count[i])*getType(Type).getNumberOfNodes();
642 const int *ind=Connectivity->getIndex();
643 return ind[_count[i+1]-1]-ind[_count[i]-1];
645 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
649 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
652 if (_constituent != NULL)
653 return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
655 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
658 /*! Give morse index array to use with
659 getConnectivity(mode,entity,MED_ALL_ELEMENTS).
661 Each value give start index for corresponding entity in connectivity array.
663 Example : i-th element, j-th node of it :
664 - In C mode : Connectivity[ConnectivityIndex[i]-1+j-1]
665 - In fortran mode : Connectivity[ConnectivityIndex[i]+j] */
666 //-----------------------------------------------------------------------------------------------//
667 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType,
668 medEntityMesh Entity) const
669 //-----------------------------------------------------------------------------------------------//
671 const char * LOC = "CONNECTIVITY::getConnectivityIndex";
673 MEDSKYLINEARRAY * Connectivity;
674 if (Entity==_entity) {
676 if (ConnectivityType==MED_NODAL)
679 Connectivity=_descending;
681 if (Connectivity!=NULL)
682 return Connectivity->getIndex();
684 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
687 if (_constituent != NULL)
688 return _constituent->getConnectivityIndex(ConnectivityType,Entity);
690 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
695 * \brief Returns CELLMODEL of a given geometrical type
697 //--------------------------------------------------------------//
698 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
699 //--------------------------------------------------------------//
702 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
703 throw MEDEXCEPTION("CONNECTIVITY::getType : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !");
704 for (int i=0; i<_numberOfTypes; i++)
705 if (_geometricTypes[i]==Type)
707 throw MEDEXCEPTION("CONNECTIVITY::getType : medGeometryElement not found !");
710 /*! Returns the number of elements of type %medGeometryElement.
711 Note : not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
712 //------------------------------------------------------------------------//
713 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
714 //------------------------------------------------------------------------//
716 const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
719 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
720 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
721 for (int i=0; i<_numberOfTypes; i++)
722 if (_geometricTypes[i]==Type)
723 return _type[i].getNumberOfNodes();
724 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
727 /*! Returns the number of geometric sub cells of %medGeometryElement type.
728 not implemented for MED_ALL_ELEMENTS nor for MED_NONE */
729 //------------------------------------------------------------------------//
730 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
731 //------------------------------------------------------------------------//
733 if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
734 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
735 for (int i=0; i<_numberOfTypes; i++)
736 if (_geometricTypes[i]==Type)
737 return _type[i].getNumberOfConstituents(1);
738 throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
741 /*! Returns the number of elements of type %medGeometryElement.
744 - Implemented for MED_ALL_ELEMENTS
745 - Not implemented for MED_ALL_ENTITIES (A VERIFIER)
746 - Not implemented for MED_NONE */
747 //-----------------------------------------------------------------------------------//
748 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
749 //-----------------------------------------------------------------------------------//
751 if(Entity==MED_EN::MED_NODE)
752 return _numberOfNodes;
753 if (Entity==_entity) {
754 if (Type==MED_EN::MED_NONE)
755 return 0; // not defined !
756 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
757 if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity))
759 if (Type==MED_EN::MED_ALL_ELEMENTS)
760 return _count[_numberOfTypes]-1;
761 for (int i=0; i<_numberOfTypes; i++)
762 if (_geometricTypes[i]==Type)
763 return (_count[i+1] - _count[i]);
764 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
766 if (_constituent != NULL)
767 return _constituent->getNumberOf(Entity,Type);
769 return 0; // valid if they are nothing else !
770 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
774 //--------------------------------------------------------------//
775 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
776 medGeometryElement Type) const
777 //--------------------------------------------------------------//
779 if (TypeConnectivity == MED_NODAL)
781 calculateNodalConnectivity();
782 if (Type==MED_ALL_ELEMENTS)
783 return _nodal->getValue();
784 for (int i=0; i<_numberOfTypes; i++)
785 if (_geometricTypes[i]==Type)
786 return _nodal->getI(_count[i]);
790 calculateDescendingConnectivity();
791 if (Type==MED_ALL_ELEMENTS)
792 return _descending->getValue();
793 for (int i=0; i<_numberOfTypes; i++)
794 if (_geometricTypes[i]==Type)
795 return _descending->getI(_count[i]);
797 throw MEDEXCEPTION("Not found");
801 //---------------------------------------------------------------------//
802 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity) const
803 //---------------------------------------------------------------------//
805 if (TypeConnectivity == MED_NODAL)
807 calculateNodalConnectivity();
808 return _nodal->getIndex();
812 calculateDescendingConnectivity();
813 return _descending->getIndex();
817 /*! Not yet implemented */
818 //----------------------------------------------//
819 const int* CONNECTIVITY:: getNeighbourhood() const
820 //----------------------------------------------//
822 throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
825 /*! Returns an array which contains, for each node, all cells
827 //-------------------------------------------------//
828 const int* CONNECTIVITY::getReverseNodalConnectivity() const
829 //-------------------------------------------------//
831 calculateReverseNodalConnectivity();
832 return _reverseNodalConnectivity->getValue();
835 /*! Give index array to use with getReverseConnectivity(MED_NODAL).
836 It is unusefull with MED_DESCENDING mode, because we have allways two cells. */
837 //-------------------------------------------------------//
838 const int* CONNECTIVITY::getReverseNodalConnectivityIndex() const
839 //-------------------------------------------------------//
841 calculateReverseNodalConnectivity();
842 return _reverseNodalConnectivity->getIndex();
845 /*! Returns an array which contains, for each face (or edge),
846 the 2 cells of each side. First is cell which face normal is outgoing.
848 //------------------------------------------------------//
849 const int* CONNECTIVITY::getReverseDescendingConnectivity() const
850 //------------------------------------------------------//
852 // it is in _constituent connectivity only if we are in MED_CELL
853 // (we could not for instance calculate face-edge connectivity !)
854 if (_entity!=MED_CELL)
855 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
857 // we need descending connectivity
858 calculateDescendingConnectivity();
859 if (_reverseDescendingConnectivity==NULL)
860 _reverseDescendingConnectivity=_descending->makeReverseArray();
862 return _reverseDescendingConnectivity->getValue();
865 /*! calculate the reverse descending Connectivity
866 and returns the index ( A DOCUMENTER MIEUX)*/
867 //-----------------------------------------------------------//
868 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex() const
869 //-----------------------------------------------------------//
871 // it is in _constituent connectivity only if we are in MED_CELL
872 if (_entity!=MED_CELL)
873 throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
875 // we need descending connectivity
876 calculateDescendingConnectivity();
877 return _reverseDescendingConnectivity->getIndex();
880 /*! A DOCUMENTER (et a finir ???) */
881 //--------------------------------------------//
882 void CONNECTIVITY::calculateNodalConnectivity() const
883 //--------------------------------------------//
887 if (_descending==NULL)
888 throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
889 // calculate _nodal from _descending
893 /*! If not yet done, calculate the nodal Connectivity
894 and the reverse nodal Connectivity*/
895 //---------------------------------------------------//
896 void CONNECTIVITY::calculateReverseNodalConnectivity() const
897 //---------------------------------------------------//
899 const char* LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
903 SCRUTE_MED(_reverseNodalConnectivity);
906 calculateNodalConnectivity();
908 if(_reverseNodalConnectivity==NULL)
910 vector <vector <int> > reverse_connectivity;
911 reverse_connectivity.resize(_numberOfNodes+1);
913 // Treat all cells types
915 const int *index=_nodal->getIndex();
916 const int *conn =_nodal->getValue();
918 for (int j = 0; j < _numberOfTypes; j++)
920 for (int k = _count[j]; k < _count[j+1]; k++, cell_number++)
922 // node number of the cell
923 int nb_nodes = index[ cell_number ] - index[ cell_number-1 ];
924 if ( _geometricTypes[ j ] == MED_EN::MED_POLYHEDRA )
926 set<int> nodes( conn, conn + nb_nodes );
927 set<int>::iterator n = nodes.begin();
929 for ( ; n != nodes.end(); ++n )
930 reverse_connectivity[ *n ].push_back( cell_number );
935 for (int i = 0; i < nb_nodes; ++i, ++conn)
936 reverse_connectivity[ *conn ].push_back( cell_number );
941 // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
943 //calculate size of reverse_nodal_connectivity array
944 int size_reverse_nodal_connectivity = 0;
945 for (int i = 1; i < _numberOfNodes+1; i++)
946 size_reverse_nodal_connectivity += reverse_connectivity[i].size();
948 int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
949 int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
951 reverse_nodal_connectivity_index[0] = 1;
952 for (int i = 1; i < _numberOfNodes+1; i++)
954 int size = reverse_connectivity[i].size();
955 reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
956 for (int j = 0; j < size; j++)
957 reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
960 _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
961 reverse_nodal_connectivity_index,
962 reverse_nodal_connectivity,true);
967 /*! If not yet done, calculate the Descending Connectivity */
968 //-------------------------------------------------//
969 void CONNECTIVITY::calculateDescendingConnectivity() const
970 //-------------------------------------------------//
972 const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
975 if (_descending==NULL)
979 MESSAGE_MED(LOC<<"No connectivity found !");
980 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
982 // calcul _descending from _nodal
983 // we need CONNECTIVITY for constituent
985 if (_constituent != NULL && _constituent->_nodal != NULL)
987 const_cast<CONNECTIVITY*>(this)->calculatePartialDescendingConnectivity();
991 const_cast<CONNECTIVITY*>(this)->calculateFullDescendingConnectivity(_entity);
996 /*! If not yet done, calculate the full Descending Connectivity */
997 //-------------------------------------------------//
998 void CONNECTIVITY::calculateFullDescendingConnectivity(MED_EN::medEntityMesh Entity)
999 //-------------------------------------------------//
1001 const char * LOC = "CONNECTIVITY::calculateFullDescendingConnectivity() : ";
1003 if (_entity != Entity)
1005 if (_constituent == NULL)
1006 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entity not found !"));
1007 _constituent->calculateFullDescendingConnectivity(Entity);
1011 if ( _descending==NULL || _isDescendingConnectivityPartial )
1015 MESSAGE_MED(LOC<<"No connectivity found !");
1016 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
1018 delete _constituent;
1020 if (_entityDimension == 3)
1021 _constituent = new CONNECTIVITY(MED_FACE);
1022 else if (_entityDimension == 2)
1023 _constituent = new CONNECTIVITY(MED_EDGE);
1026 MESSAGE_MED(LOC<<"We are in 1D");
1030 bool hasPoly = ( _numberOfTypes && _geometricTypes[_numberOfTypes-1] == ( _entityDimension > 2 ? MED_POLYHEDRA : MED_POLYGON ));
1031 int numberOfClassicTypes = _numberOfTypes - hasPoly;
1033 _constituent->_typeConnectivity = MED_NODAL;
1034 _constituent->_numberOfNodes = _numberOfNodes;
1035 // foreach cells, we built array of constituent
1036 int DescendingSize = 0;
1037 for(int i=0; i<numberOfClassicTypes; i++)
1038 DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
1039 int * descend_connectivity = new int[DescendingSize];
1040 for (int i=0; i<DescendingSize; i++)
1041 descend_connectivity[i]=0;
1042 int * descend_connectivity_index = new int[_count[numberOfClassicTypes]];
1043 if(_count[numberOfClassicTypes]>0)
1044 descend_connectivity_index[0]=1;
1047 map<medGeometryElement,int> eltsCounter;
1048 medGeometryElement ConstituentsTypes[2] = { MED_NONE, MED_NONE };
1049 int NumberOfConstituentsForeachType [2] = { 0, 0 };
1050 map<medGeometryElement,int>::iterator status;
1051 for(int i=0; i<numberOfClassicTypes; i++)
1053 // initialize descend_connectivity_index array :
1054 int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
1055 for (int j=_count[i];j<_count[i+1];j++)
1057 descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
1058 // compute number of constituent of all cell for each type
1059 for(int k=1;k<NumberOfConstituents+1;k++)
1061 medGeometryElement MEDType = _type[i].getConstituentType(1,k);
1062 status=eltsCounter.find(MEDType);
1063 if(status!=eltsCounter.end())
1066 eltsCounter[MEDType]=1;
1070 if ( hasPoly && _entityDimension == 2 ) // there are constituent segments
1072 status=eltsCounter.insert(make_pair( MED_SEG2,0 )).first;
1073 //status->second++; // increment zero or that there was
1075 if(eltsCounter.size()>2)
1077 // free memory (issue 0020411: [CEA 342] Sigsegv on gibi writing of MESH coming from MED file without face computation)
1078 delete [] descend_connectivity;
1079 delete [] descend_connectivity_index;
1080 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" Descending connectivity does not support more than 2 types."));
1082 status=eltsCounter.begin();
1083 if(!eltsCounter.empty())
1085 ConstituentsTypes[0]=(*status).first; NumberOfConstituentsForeachType[0]=(*status).second;
1086 if(++status!=eltsCounter.end())
1087 ConstituentsTypes[1]=(*status).first, NumberOfConstituentsForeachType[1]=(*status).second;
1089 int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
1090 int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
1092 int * ConstituentNodalConnectivity = new int[TotalNumberOfNodes];
1093 int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
1094 ConstituentNodalConnectivityIndex[0]=1;
1096 _constituent->_entityDimension = _entityDimension-1;
1097 int nbClassicConstituentTypes = eltsCounter.size();
1098 int nbConstituentTypes = nbClassicConstituentTypes + hasPoly;
1099 _constituent->_numberOfTypes = nbConstituentTypes;
1100 _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
1101 _constituent->_type = new CELLMODEL[_constituent->_numberOfTypes];
1102 if(_constituent->_count) delete [] _constituent->_count;
1103 _constituent->_count = new int[_constituent->_numberOfTypes+1];
1104 _constituent->_count[0]=1;
1105 med_int* tmp_NumberOfConstituentsForeachType = new med_int[nbClassicConstituentTypes+1];
1106 tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
1107 for (int i=0; i<nbClassicConstituentTypes;i++)
1109 _constituent->_geometricTypes[i]=ConstituentsTypes[i];
1110 _constituent->_type[i]=CELLMODEL(ConstituentsTypes[i]);
1111 _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
1112 tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
1113 for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
1114 ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
1117 // we need reverse nodal connectivity
1118 if (! _reverseNodalConnectivity)
1119 calculateReverseNodalConnectivity();
1120 const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
1121 const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
1123 // array to keep reverse descending connectivity
1124 int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
1126 int TotalNumberOfSubCell = 0;
1127 for (int i=0; i<numberOfClassicTypes; i++)
1128 { // loop on all cell type
1129 CELLMODEL Type = _type[i];
1130 // int NumberOfNodesPerCell = Type.getNumberOfNodes();
1131 int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
1132 for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
1133 for (int k=1; k<=NumberOfConstituentPerCell; k++)
1134 { // we loop on all sub cell of it
1135 if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0)
1136 { // it is a new sub cell !
1137 // TotalNumberOfSubCell++;
1139 if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0])
1141 tmp_NumberOfConstituentsForeachType[0]++;
1142 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
1145 tmp_NumberOfConstituentsForeachType[1]++;
1146 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
1148 //we have maximum two types
1150 descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
1151 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
1153 int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
1154 int * NodesLists = new int[NumberOfNodesPerConstituent];
1155 for (int l=0; l<NumberOfNodesPerConstituent; l++)
1157 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
1158 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
1160 // we use reverse_nodal_connectivity to find the other element which contain this sub cell
1162 // all elements which contains first node of sub cell :
1163 int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
1164 int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
1165 //ReverseNodalConnectivityIndex[NodesLists[0]];
1166 int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
1168 if (NumberOfCellsInList > 0)
1169 { // we could have no element !
1170 int * CellsList = new int[NumberOfCellsInList];
1171 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
1172 CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
1174 // foreach node in sub cell, we search elements which are in common
1175 // at the end, we must have only one !
1177 for (int l=1; l<NumberOfNodesPerConstituent; l++)
1179 int NewNumberOfCellsInList = 0;
1180 int * NewCellsList = new int[NumberOfCellsInList];
1181 for (int l1=0; l1<NumberOfCellsInList; l1++)
1182 for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
1183 //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
1185 if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
1186 // increasing order : CellsList[l1] are not in elements list of node l
1188 if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
1189 // we have found one
1190 NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
1191 NewNumberOfCellsInList++;
1195 NumberOfCellsInList = NewNumberOfCellsInList;
1197 delete [] CellsList;
1198 CellsList = NewCellsList;
1201 if (NumberOfCellsInList > 0) { // We have found some elements !
1202 int CellNumber = CellsList[0];
1203 //delete [] CellsList;
1204 if (NumberOfCellsInList>1) {
1205 // free memory (issue 0020411: [CEA 342] Sigsegv on gibi writing of MESH coming from MED file without face computation)
1206 delete [] CellsList;
1207 delete [] NodesLists;
1208 delete [] ReverseDescendingConnectivityValue;
1209 delete [] tmp_NumberOfConstituentsForeachType;
1210 delete [] descend_connectivity;
1211 delete [] descend_connectivity_index;
1212 delete [] ConstituentNodalConnectivity;
1213 delete [] ConstituentNodalConnectivityIndex;
1214 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
1217 if (NumberOfCellsInList==1)
1219 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
1221 // we search sub cell number in this cell to not calculate it another time
1224 for (int l=0; l<numberOfClassicTypes; l++)
1225 if (CellNumber < _count[l+1])
1230 //int sub_cell_count2 = Type2.get_entities_count(1);
1231 //int nodes_cell_count2 = Type2.get_nodes_count();
1233 for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) // on all sub cell
1236 for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
1237 for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
1239 if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
1242 if (counter==Type.getConstituentType(1,k)%100)
1244 descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
1251 MESSAGE_MED(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
1254 ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
1256 delete [] CellsList;
1259 delete [] NodesLists;
1263 // we adjust _constituent
1264 int NumberOfConstituent=0;
1265 int SizeOfConstituentNodal=0;
1266 for (int i=0;i<nbClassicConstituentTypes; i++)
1268 NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
1269 SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
1271 // we built new _nodal attribute in _constituent
1272 //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
1273 //const int *ConstituentNodalValue = ConstituentNodal->getValue();
1274 //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
1275 int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
1276 int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
1277 ConstituentNodalIndex[0]=1;
1278 // we build _reverseDescendingConnectivity
1279 //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
1280 //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
1281 //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
1282 int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
1283 int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
1284 reverseDescendingConnectivityIndex[0]=1;
1286 // first constituent type
1287 for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++)
1289 ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
1290 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
1291 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
1293 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1294 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++)
1295 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
1297 // second type if any
1298 if (nbClassicConstituentTypes==2)
1300 int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
1301 int offset1=offset*_constituent->_type[0].getNumberOfNodes();
1302 int offset2=offset*2;
1303 int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
1304 for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++)
1306 ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
1307 for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
1308 ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
1310 reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
1311 for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++)
1312 reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
1314 _constituent->_count[2]=NumberOfConstituent+1;
1315 // we correct _descending to adjust face number
1316 for(int j=0;j<DescendingSize;j++)
1317 if (abs(descend_connectivity[j])>tmp_NumberOfConstituentsForeachType[0]) {
1318 if ( descend_connectivity[j] > 0 )
1319 descend_connectivity[j]-=offset;
1321 descend_connectivity[j]+=offset;
1325 delete [] ConstituentNodalConnectivityIndex;
1326 delete [] ConstituentNodalConnectivity;
1327 delete [] ReverseDescendingConnectivityValue;
1328 if (_constituent->_numberOfTypes > 0)
1329 _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
1330 delete [] tmp_NumberOfConstituentsForeachType;
1333 vector<int> Descending( descend_connectivity, descend_connectivity + DescendingSize );
1334 vector<int> DescendingIndex( descend_connectivity_index, descend_connectivity_index + _count[numberOfClassicTypes] );
1335 vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
1336 vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
1337 delete [] descend_connectivity;
1338 delete [] descend_connectivity_index;
1339 delete [] reverseDescendingConnectivityValue;
1340 delete [] reverseDescendingConnectivityIndex;
1343 // polygons (2D mesh)
1345 vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
1346 vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
1347 delete [] ConstituentNodalValue;
1348 delete [] ConstituentNodalIndex;
1349 int NumberOfNewSeg = 0;
1351 int i = _count[_numberOfTypes - ( hasPoly && _entityDimension == 2 )];
1352 int NumberOfPolygons = _count[_numberOfTypes] - i;
1353 for (; i < _count[_numberOfTypes]; i++) // for each polygon
1355 const int * vector_begin = _nodal->getValue() + _nodal->getIndex()[i-1] - 1;
1356 int vector_size = _nodal->getIndex()[i]-_nodal->getIndex()[i-1];
1357 vector<int> myPolygon(vector_begin, vector_begin+vector_size);
1358 myPolygon.push_back(myPolygon[0]);
1359 DescendingIndex.push_back( DescendingIndex.back() + vector_size);
1360 for (int j=0; j<(int)myPolygon.size()-1; j++) // for each segment of polygon
1362 MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
1363 int ret_compare = 0;
1365 // we search it in existing segments
1367 for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
1369 MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
1370 ret_compare = segment_poly.compare(segment);
1371 if (ret_compare) // segment_poly already exists
1373 Descending.push_back(ret_compare*(k+1)); // we had it to the connectivity
1374 insert_vector(Reversedescendingconnectivityvalue, 2*k+1, i); // add polygon i to reverse descending connectivity for segment_poly (in 2nd place)
1379 // segment_poly must be created
1384 Descending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
1385 Constituentnodalvalue.push_back(segment_poly[0]);
1386 Constituentnodalvalue.push_back(segment_poly[1]);
1387 Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
1388 insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewSeg-1), i ); // add polygon i to reverse descending connectivity for segment_poly (in 1st place)
1389 insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
1394 if ( NumberOfPolygons > 0 )
1396 NumberOfConstituent += NumberOfNewSeg;
1397 SizeOfConstituentNodal += 2*NumberOfNewSeg;
1400 // polyhedron (3D mesh)
1402 vector<int> Constituentpolygonsnodalvalue;
1403 vector<int> Constituentpolygonsnodalindex(1,1);
1404 int NumberOfNewFaces = 0; // by convention new faces are polygons
1405 //offset to switch between all types and classical types.
1406 //int offsetCell = getNumberOf(MED_CELL, MED_ALL_ELEMENTS);
1407 int *tabRes = new int[1000]; //temporay array for intersection calculation
1409 i = _count[_numberOfTypes - ( hasPoly && _entityDimension == 3 )];
1410 int NumberOfPolyhedron = _count[_numberOfTypes] - i;
1411 for (; i<_count[_numberOfTypes]; i++) // for each polyhedron
1413 const int * nodes = _nodal->getValue() + _nodal->getIndex()[i-1] - 1;
1414 const int * nodesEnd = _nodal->getValue() + _nodal->getIndex()[i ] - 1;
1415 int nbConstituentFaces = 0;
1416 while ( nodes < nodesEnd ) // loop on faces of polyhedron
1418 // look for end of face connectivity
1419 const int * faceNodesEnd = nodes;
1420 while ( ++faceNodesEnd < nodesEnd && *faceNodesEnd != -1 );
1421 int myFaceNumberOfNodes = faceNodesEnd - nodes;
1422 MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,nodes);
1423 // move pointer to next face while continue working with face_poly
1424 if (( nodes = faceNodesEnd ) < nodesEnd )
1425 nodes++; // skip face separator
1427 nbConstituentFaces++;
1429 int ret_compare = 0;
1431 // we search it in existing faces
1433 // we search first in TRIA3 and QUAD4
1434 if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
1436 int Begin = -1; // first TRIA3 or QUAD4
1437 int Number = 0; // number of TRIA3 or QUAD4
1438 for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
1439 if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
1446 for (int k=0; k<Number; k++)
1448 MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
1449 ret_compare = face_poly.compare(face);
1452 Descending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
1453 insert_vector(Reversedescendingconnectivityvalue, 2*(Begin+k)+1, i); // add polyhedra i to reverse descending connectivity for face_poly (in 2nd place)
1459 // we search last in POLYGONS
1463 const int *facePolyTab=face_poly.getArray(lgth);
1464 int nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[0]] -
1465 ReverseNodalConnectivityIndex[facePolyTab[0]-1];
1466 const int *candidatesCell = ReverseNodalConnectivityValue +
1467 ReverseNodalConnectivityIndex[facePolyTab[0]-1] - 1;
1468 memcpy(tabRes,candidatesCell,nbOfCandidatesCell*sizeof(int));
1469 int lgth2=nbOfCandidatesCell;
1470 for (int k=1;k<lgth && lgth2!=0;k++)
1472 nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[k]] -
1473 ReverseNodalConnectivityIndex[facePolyTab[k]-1];
1474 candidatesCell = ReverseNodalConnectivityValue +
1475 ReverseNodalConnectivityIndex[facePolyTab[k]-1] - 1;
1476 mergeOrderedTabs(tabRes,lgth2,candidatesCell,nbOfCandidatesCell,tabRes,lgth2);
1479 ret_compare=0;//here normally tabRes[0]==i
1480 else //> 2 should never happend : A face is shared by more than 2 polyhedrons...
1482 if (tabRes[0] == i) //as tabRes is ordered by construction tabRes[1] > tabRes[0] so the current
1483 // face is shared with an another cell whose id > current id. So let's create
1486 {//tabRes[0]<Constituentpolygonsnodalindex.size()-1 that is to say the current face has been built previously.
1487 int nbOfFacesConstitutingAlreadyBuiltPolyh;
1488 int* nbOfNodesPerFaces;
1489 int** nodesOfFacesOfAlreadyBuiltPolyh = getNodesPerFaceOfPolyhedron( tabRes[0], nbOfFacesConstitutingAlreadyBuiltPolyh, nbOfNodesPerFaces);
1490 for (int k1=0; k1<nbOfFacesConstitutingAlreadyBuiltPolyh && (ret_compare==0); k1++)
1492 int nbOfNodesForCurrentFace = nbOfNodesPerFaces[k1];
1493 MEDMODULUSARRAY face(nbOfNodesForCurrentFace,nodesOfFacesOfAlreadyBuiltPolyh[k1]);
1494 ret_compare = face_poly.compare(face);
1497 int curFaceId=Descending[ DescendingIndex[ tabRes[0]-1 ]-1 + k1 ];
1498 Descending.push_back(ret_compare*curFaceId); // we had it to the connectivity
1499 insert_vector(Reversedescendingconnectivityvalue, 2*(curFaceId-1)+1, i);
1502 delete [] nbOfNodesPerFaces;
1503 delete [] nodesOfFacesOfAlreadyBuiltPolyh;
1508 // if not found, face_poly must be created
1513 Descending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
1514 for (int k=0; k<myFaceNumberOfNodes; k++)
1515 Constituentpolygonsnodalvalue.push_back(face_poly[k]);
1516 Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
1517 insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewFaces-1), i); // add polyhedra i to reverse descending connectivity for face_poly (in 1st place)
1518 insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
1520 } // loop on faces of polyhedron
1522 DescendingIndex.push_back( DescendingIndex.back() + nbConstituentFaces );
1526 if (NumberOfPolyhedron > 0)
1528 NumberOfConstituent += NumberOfNewFaces;
1529 Constituentnodalvalue.insert( Constituentnodalvalue.end(),
1530 Constituentpolygonsnodalvalue.begin(),
1531 Constituentpolygonsnodalvalue.end());
1532 Constituentpolygonsnodalvalue.clear();
1533 Constituentnodalindex.reserve( Constituentnodalindex.size() + NumberOfNewFaces );
1534 int indexShift = Constituentnodalindex.back() - Constituentpolygonsnodalindex.front();
1535 for ( unsigned i = 1; i < Constituentpolygonsnodalindex.size(); ++i )
1536 Constituentnodalindex.push_back( Constituentpolygonsnodalindex[i] + indexShift );
1537 Constituentpolygonsnodalindex.clear();
1538 _constituent->_geometricTypes[ _constituent->_numberOfTypes-1 ] = MED_POLYGON;
1541 _constituent->_count[ _constituent->_numberOfTypes ] = NumberOfConstituent+1;
1542 _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
1543 Constituentnodalvalue.size(),
1544 &Constituentnodalindex[0],
1545 &Constituentnodalvalue[0]);
1549 _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
1551 &DescendingIndex[0],
1554 Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
1555 Reversedescendingconnectivityvalue.push_back(0);
1556 _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent+1,
1557 2*NumberOfConstituent,
1558 &Reversedescendingconnectivityindex[0],
1559 &Reversedescendingconnectivityvalue[0]);
1560 _isDescendingConnectivityPartial = false;
1566 void CONNECTIVITY::addToDescendingConnectivity(const set<int>& nodes,
1567 multimap<int,int>& descending,
1569 const CONNECTIVITY_HashMap & face_map ) const
1571 int dimension = getEntityDimension();
1572 vector<int> signature (dimension);
1573 set<int>::const_iterator iter=nodes.begin();
1574 for (int i=0; i< dimension;i++)
1581 CONNECTIVITY_HashMap::const_iterator itermap=face_map.find(signature);
1582 CONNECTIVITY_HashMap::const_iterator iterend=face_map.end();
1585 if (itermap!=iterend)
1586 descending.insert(make_pair(iglobal_cell,itermap->second));
1589 /*! This method calculates the descending connectivity without creating missing elements. It only maps the constituent elements that are described in the nodal representation.
1590 For instance, let us consider the following mesh with no MED_EDGE elements on the inner edges.
1592 +----1----+----2----+
1596 +---------+---------+
1600 +----6----+----5----+
1603 calculatePartialDescendingConnectivity()
1610 whereas calculateDescendingConnectivity()
1611 will create new edges, renumbering existing ones and return
1617 +----1----+----5----+
1621 +----3----+----7----+
1625 +----9----+----12---+
1628 void CONNECTIVITY::calculatePartialDescendingConnectivity() const
1630 ////////////////////////////////////////////////////////////////////////////
1631 // First stage : creating hash_map referencing faces with the triplet
1632 // of the lowest order nodes as a key and the global face number as a value
1634 CONNECTIVITY_HashMap face_map;
1637 int dimension = getEntityDimension();
1639 const int* conn =_constituent->getValue (MED_NODAL, MED_ALL_ELEMENTS);
1640 const int* index=_constituent->getValueIndex(MED_NODAL);
1641 int nbfaces=_constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS);
1642 for ( ; iglobal_face <= nbfaces; iglobal_face++ )
1644 const int* face_conn = conn + index[iglobal_face-1] - 1;
1645 const int nbnodes = index[iglobal_face] - index[iglobal_face-1];
1646 set<int> nodes( face_conn, face_conn + nbnodes );
1647 vector<int> signature (dimension);
1648 set<int>::iterator iter=nodes.begin();
1649 for (int i=0; i< dimension;i++)
1654 face_map.insert(make_pair(signature,iglobal_face));
1657 ////////////////////////////////////////////////////////////////////////////
1658 //Second stage : going through all the faces of the cells and
1659 // connecting them to the hash_map created in the first stage
1661 multimap<int,int> descending; //map storing the descending connectivity
1663 int nbcell_types = getNumberOfTypes(_entity);
1664 const medGeometryElement* cell_types = getGeometricTypes(_entity);
1665 int iglobal_cell = 1;
1666 for (int itype=0; itype<nbcell_types; itype++)
1668 medGeometryElement cell_type = cell_types[itype];
1669 int nbcells = getNumberOf(_entity,cell_type);
1671 const int* index = _nodal->getIndex();
1672 const int* conn = _nodal->getValue();
1673 switch ( cell_type )
1677 for (int icell=0;icell<nbcells;icell++)
1679 int nbnodes=index[icell+1]-index[icell];
1680 int nbfaces=nbnodes;
1681 for (int iface=0; iface<nbfaces;iface++)
1684 if (iface+1!=nbfaces)
1686 nodes.insert(conn[index[icell]-1+iface]);
1687 nodes.insert(conn[index[icell]-1+iface+1]);
1691 nodes.insert(conn[index[icell]-1+iface]);
1692 nodes.insert(conn[index[icell]-1]);
1694 addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
1703 for (int icell = 0; icell < nbcells; icell++)
1705 const int* face_conn = conn + index[icell ] - 1;
1706 const int* face_conn_end = conn + index[icell+1] - 1;
1707 while (face_conn < face_conn_end)
1709 set<int> nodes; // faces are divided by -1
1710 for ( ; face_conn < face_conn_end && *face_conn != -1; face_conn++ )
1711 nodes.insert( *face_conn );
1712 addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
1719 default: // classic cells
1721 CELLMODEL cellmodel=CELLMODEL_Map::retrieveCellModel(cell_type);
1722 const int* index=_nodal->getIndex();
1723 const int* conn=_nodal->getValue();
1724 for (int icell = 0; icell < nbcells; icell++, iglobal_cell++)
1726 int nbfaces=cellmodel.getNumberOfConstituents(1);
1727 for (int iface=0; iface<nbfaces;iface++)
1730 const int* local_index=cellmodel.getNodesConstituent(1,iface+1);
1731 medGeometryElement face_type = cellmodel.getConstituentType(1,iface+1);
1732 int nbnodes=face_type%100;
1733 for (int inode=0;inode<nbnodes;inode++)
1734 nodes.insert(conn[index[iglobal_cell-1]-1+local_index[inode]-1]);
1735 addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
1741 ////////////////////////////////////////////////////////////////////////////
1742 // Third stage : reorganizing the descending data to store it in a medskylinearray
1748 //the number of cells is given by the number
1749 //obtained by browsing all the cell types
1750 int nb_cells = iglobal_cell-1;
1752 //for (int icell = 0; icell < nb_cells; icell++)
1753 for (int icell = 1; icell <= nb_cells; icell++)
1755 multimap<int,int>::iterator beginning_of_range = descending.lower_bound(icell);
1756 multimap<int,int>::iterator end_of_range = descending.upper_bound(icell);
1758 for (multimap<int,int>::iterator iter2 = beginning_of_range; iter2 != end_of_range; iter2++)
1760 value.push_back(iter2->second);
1763 index.push_back(index.back()+nb);
1767 _descending = new MEDSKYLINEARRAY(index.size()-1, value.size(), &index[0], &value[0]);
1768 _isDescendingConnectivityPartial = true;
1772 /*! Not implemented yet */
1773 //--------------------------------------------------------------------//
1774 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
1775 //--------------------------------------------------------------------//
1777 const char* LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
1780 MESSAGE_MED(PREFIX_MED<<"method not yet implemented " << myConnectivity._entity);
1789 Returns the geometry of an element given by its entity type & its global number.
1791 Example : medGeometryElement myType = myConnectivity.getElementType(MED_CELL,35);
1793 //--------------------------------------------------------------------//
1794 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
1795 //--------------------------------------------------------------------//
1797 const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
1799 int globalNumberMin = 1;
1800 int globalNumberMax ;
1802 if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
1803 else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
1805 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1807 // The globalNumber must verify : 1 <= globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
1809 if ( (globalNumber < globalNumberMin) || (globalNumber > globalNumberMax-1 ) )
1810 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
1811 << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
1813 if (_entity==Entity) {
1814 for(int i=1; i<=_numberOfTypes;i++)
1815 if (globalNumber<_count[i])
1816 return _geometricTypes[i-1];
1818 else if (_constituent!=NULL)
1819 return _constituent->getElementType(Entity,globalNumber);
1821 throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
1822 throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
1827 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
1829 os << endl << "------------- Entity = ";
1844 case MED_ALL_ENTITIES:
1845 os << "MED_ALL_ENTITIES";
1851 os << " -------------\n\nMedConnectivity : ";
1852 switch (co._typeConnectivity)
1855 os << "MED_NODAL\n";
1857 case MED_DESCENDING:
1858 os << "MED_DESCENDING\n";
1863 os << "Entity dimension : " << co._entityDimension << endl;
1864 os << "Number of nodes : " << co._numberOfNodes << endl;
1865 os << "Number of types : " << co._numberOfTypes << endl;
1866 for (int i=0; i!=co._numberOfTypes ; ++i)
1867 os << " -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
1868 << co._count[i+1]-co._count[i] << " elements" << endl;
1870 if (co._typeConnectivity == MED_NODAL )
1872 for (int i=0; i<co._numberOfTypes; i++)
1874 os << endl << co._type[i].getName() << " : " << endl;
1875 int numberofelements = co._count[i+1]-co._count[i];
1876 const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, MED_EN::MED_ALL_ELEMENTS);
1877 const int * index = co.getConnectivityIndex(co._typeConnectivity, co._entity) + co._count[i]-1;
1878 for (int j=0;j<numberofelements;j++)
1880 const int * cellconn = connectivity+index[j]-1;
1881 const int * cellconnend = connectivity+index[j+1]-1;
1882 os << setw(6) << j+1 << " : ";
1883 while ( cellconn < cellconnend )
1884 os << *cellconn++ <<" ";
1889 else if (co._typeConnectivity == MED_DESCENDING)
1891 int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
1892 if (numberofelements > 0)
1894 const int *connectivity = co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
1895 const int *connectivity_index = co.getConnectivityIndex( co._typeConnectivity, co._entity );
1897 for ( int j=0; j!=numberofelements; j++ )
1899 os << "Element " << j+1 << " : ";
1900 for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
1901 os << connectivity[k-1] << " ";
1907 if (co._constituent)
1908 os << endl << *co._constituent << endl;
1914 method that adds to vector 'nodes' all the nodes of polyhedron with id 'polyhedronId'.
1915 WARNING the returned pointer should be deallocated. Returned nodes and polyhedronId are in form [1,...]
1917 int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
1919 int offsetWithClassicType=_count[ _numberOfTypes-1 ]-1; // hope that polyhedrons is the last type
1920 if (polyhedronId<=offsetWithClassicType || polyhedronId> getNumberOf (MED_CELL, MED_ALL_ELEMENTS))
1921 throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
1923 const int *nodes=getConnectivity(MED_NODAL, MED_CELL, MED_ALL_ELEMENTS);
1924 const int *index=getConnectivityIndex(MED_NODAL,MED_CELL);
1925 const int *cellnodes = nodes + index[polyhedronId-1] - 1;
1926 const int *cellnodesend = nodes + index[polyhedronId ] - 1;
1927 set<int> retInSet( cellnodes, cellnodesend );
1928 if ( *retInSet.begin() == -1 )
1929 retInSet.erase( retInSet.begin() );
1930 lgthOfTab=retInSet.size();
1931 int *ret=new int[lgthOfTab];
1932 set<int>::iterator iter=retInSet.begin();
1933 for(int i=0;iter!=retInSet.end();iter++)
1939 Idem as MESH::getNodesOfPolyhedron except that returned nodes are sorted by face. 'nbOfNodesPerFaces' is an array of size 'nbOfFaces'.
1940 Returned int** has a size of 'nbOfNodesPerFaces' too, and for each element j in int** the size is nbOfNodesPerFaces[j].
1941 Warning both returned 'nbOfNodesPerFaces' and returned value should be deallocated. Returned nodes and 'polyhedronId' are in form [1,...]
1943 int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
1945 int offsetWithClassicType=_count[ _numberOfTypes-1 ]-1; // hope that polyhedrons is the last type
1946 if (polyhedronId<=offsetWithClassicType || polyhedronId> getNumberOf(MED_CELL, MED_ALL_ELEMENTS))
1947 throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
1948 const int *nodes=getConnectivity(MED_EN::MED_NODAL, MED_CELL, MED_ALL_ELEMENTS);
1949 const int *index=getConnectivityIndex(MED_EN::MED_NODAL, MED_CELL);
1951 const int *cellnodes = nodes + index[polyhedronId-1] - 1;
1952 const int *cellnodesend = nodes + index[polyhedronId ] - 1;
1954 nbOfNodesPerFaces=new int[cellnodesend-cellnodes]; // more than enough
1955 int **ret=new int* [cellnodesend-cellnodes];
1956 while ( cellnodes < cellnodesend )
1958 ret[nbOfFaces] = (int*)cellnodes;
1959 while ( cellnodes < cellnodesend && *cellnodes != -1 ) ++cellnodes;
1960 nbOfNodesPerFaces[nbOfFaces] = cellnodes - ret[nbOfFaces];
1968 Method used in CalculateDescendingConnectivity. So it's typically a private method.
1969 The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
1970 are not in separate data structure, contrary to not reverse connectivities.
1972 int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk) const
1974 const int lastTypeIsPoly = (_geometricTypes[_numberOfTypes-1]==MED_EN::MED_POLYGON||
1975 _geometricTypes[_numberOfTypes-1]==MED_EN::MED_POLYHEDRA);
1976 if ( !lastTypeIsPoly )
1977 return reverseNodalIndex[rk+1];
1978 int nbOfLastElt=_count[_numberOfTypes-1]-1;
1979 int ret=reverseNodalIndex[rk];
1980 for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
1982 if(reverseNodalValue[i-1]<=nbOfLastElt)
1989 Method that inverts the face with faceId 'faceId' in the data structure.
1990 This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
1991 WARNING calculateDescendingConnectivity must have been called before.
1993 void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace)
1995 // we have always 2 neighbourings
1996 int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
1997 int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
1999 if (cell2 != 0) { // we are not on border, update compulsory. If we aren't on border no update necessary so WARNING because user specified a bad oriented face
2000 // Updating _reverseDescendingConnectivity
2001 _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
2002 _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
2003 // Updating _constituent->_nodal because of reversity
2004 const int *descendingNodalIndex=_constituent->_nodal->getIndex();
2005 MEDSKYLINEARRAY *currentNodal=_constituent->_nodal;
2006 int faceIdRelative=faceId;
2007 for(int iarray=1;iarray<=(descendingNodalIndex[faceIdRelative]-descendingNodalIndex[faceIdRelative-1]);iarray++)
2008 currentNodal->setIJ(faceIdRelative,iarray,nodalConnForFace[iarray-1]);
2010 // Updating _descending for cell1 and cell2
2011 const int NB_OF_CELLS_SHARING_A_FACE=2;
2012 int cellsToUpdate[NB_OF_CELLS_SHARING_A_FACE] = { cell1, cell2 };
2013 for(int curCell=0;curCell<NB_OF_CELLS_SHARING_A_FACE;curCell++)
2015 int cell=cellsToUpdate[curCell];
2016 const int *newDescendingIndex=_descending->getIndex();
2017 MEDSKYLINEARRAY *currentDescending=_descending;
2018 for(int iface=newDescendingIndex[cell-1];iface<newDescendingIndex[cell];iface++)
2020 int curValue=currentDescending->getIndexValue(iface);
2021 if (abs(curValue)==faceId)
2022 currentDescending->setIndexValue(iface,-curValue);
2029 Method with 2 output : the connectivity required and its length 'lgth'.
2030 This method gives the connectivity independently it is a polygons/polyhedrons or normal element.
2032 const int * CONNECTIVITY::getConnectivityOfAnElement(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth) const
2034 if(Entity==MED_EN::MED_NODE)
2035 throw MEDEXCEPTION("No connectivity attached to a node entity");
2038 if ( Number > getNumberOf( Entity, MED_ALL_ELEMENTS ))
2039 throw MEDEXCEPTION("Unknown number");
2040 const int * conn = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
2041 const int * index = getConnectivityIndex(ConnectivityType,Entity);
2042 lgth=index[Number]-index[Number-1];
2043 return conn+index[Number-1]-1;
2047 if(_constituent==NULL)
2048 calculateDescendingConnectivity();
2049 return _constituent->getConnectivityOfAnElement(ConnectivityType,Entity,Number,lgth);
2054 Perform a deep comparison of the 2 connectivities in NODAL mode and on all elements.
2056 bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
2058 CONNECTIVITY* temp=(CONNECTIVITY* )this;
2059 const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2060 int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2061 temp=(CONNECTIVITY* )(&other);
2062 const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2063 int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
2067 for(int i=0;i<size1 && ret;i++)
2069 ret=(conn1[i]==conn2[i]);