Back to index

salome-med  6.5.0
MEDMEM_Connectivity.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "MEDMEM_Connectivity.hxx"
00024 #include "MEDMEM_Family.hxx"
00025 #include "MEDMEM_Group.hxx"
00026 #include "MEDMEM_CellModel.hxx"
00027 
00028 #include "MEDMEM_SkyLineArray.hxx"
00029 #include "MEDMEM_ModulusArray.hxx"
00030 
00031 #include "MEDMEM_STRING.hxx"
00032 #include "MEDMEM_Utilities.hxx"
00033 #include <iomanip>
00034 
00035 using namespace std;
00036 using namespace MEDMEM;
00037 using namespace MED_EN;
00038 
00039 // Enlarge the vector if necessary to insert the element
00040 static inline void insert_vector(vector<int> &Vect, int Indice, int Element)
00041 {
00042   if (Indice >=(int) Vect.capacity())
00043     Vect.reserve(Indice + 1000);
00044 
00045   if (Indice >=(int) Vect.size())
00046     Vect.resize(Indice+1);
00047 
00048   Vect[Indice] = Element;
00049 }
00050 
00051 void mergeOrderedTabs(const int *tab1, int lgth1, const int *tab2, int lgth2, int *result, int& lgth);
00052 
00053 void mergeOrderedTabs(const int *tab1, int lgth1, const int *tab2, int lgth2, int *result, int& lgth)
00054 {
00055   int cpt[2]={0,0};
00056   lgth=0;
00057   unsigned char switcher=0;
00058   const int *tabS[2]={tab1,tab2};
00059   while(cpt[0]<lgth1 && cpt[1]<lgth2)
00060   {
00061     if(tabS[1-switcher][cpt[1-switcher]]<tabS[switcher][cpt[switcher]])
00062       cpt[1-switcher]++;
00063     else if(tabS[1-switcher][cpt[1-switcher]]>tabS[switcher][cpt[switcher]])
00064       switcher=1-switcher;
00065     else
00066     {
00067       int tmp=tabS[switcher][cpt[switcher]];
00068       cpt[switcher]++; cpt[1-switcher]++;
00069       result[lgth++]=tmp;
00070     }
00071   }
00072 }
00073 
00077 //--------------------------------------------------------------//
00078 CONNECTIVITY::CONNECTIVITY(medEntityMesh Entity /* =MED_CELL */) :
00079   //--------------------------------------------------------------//
00080   _entity(Entity),
00081   _typeConnectivity(MED_NODAL),
00082   _numberOfTypes(0),
00083   _geometricTypes((medGeometryElement*)NULL),
00084   _type((CELLMODEL*)NULL),
00085   _entityDimension(0),
00086   _numberOfNodes(0),
00087   _count((int*)NULL),
00088   _nodal((MEDSKYLINEARRAY*)NULL),
00089   _descending((MEDSKYLINEARRAY*)NULL),
00090   _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
00091   _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
00092   _neighbourhood((MEDSKYLINEARRAY*)NULL),
00093   _constituent((CONNECTIVITY*)NULL),
00094   _isDescendingConnectivityPartial(false)
00095 {
00096   const char* LOC = "CONNECTIVITY(medEntityMesh Entity=MED_CELL)";
00097   BEGIN_OF_MED(LOC);
00098   MESSAGE_MED("CONNECTIVITY(medEntityMesh Entity=MED_CELL)");
00099   _count = new int[1];
00100   _count[0]=1;
00101   END_OF_MED(LOC);
00102 }
00103 
00107 //------------------------------------------------------------------------------//
00108 CONNECTIVITY::CONNECTIVITY(int numberOfTypes,medEntityMesh Entity /* =MED_CELL */):
00109   //------------------------------------------------------------------------------//
00110   _entity(Entity),
00111   _typeConnectivity(MED_NODAL),
00112   _numberOfTypes(numberOfTypes),
00113   _entityDimension(0),
00114   _numberOfNodes(0),
00115   _nodal((MEDSKYLINEARRAY*)NULL),
00116   _descending((MEDSKYLINEARRAY*)NULL),
00117   _reverseNodalConnectivity((MEDSKYLINEARRAY*)NULL),
00118   _reverseDescendingConnectivity((MEDSKYLINEARRAY*)NULL),
00119   _neighbourhood((MEDSKYLINEARRAY*)NULL),
00120   _constituent((CONNECTIVITY*)NULL),
00121   _isDescendingConnectivityPartial(false)
00122 {
00123   MESSAGE_MED("CONNECTIVITY(int numberOfTypes,medEntityMesh Entity=MED_CELL)");
00124   _geometricTypes = new medGeometryElement[numberOfTypes];
00125   _type = new CELLMODEL[numberOfTypes];
00126   _count = new int[numberOfTypes+1];
00127   _count[0]=1;
00128   if ( numberOfTypes )
00129     _count[ numberOfTypes-1 ] = 0; // to know whether _count is set or not
00130 }
00131 
00135 //------------------------------------------------------//
00136 CONNECTIVITY::CONNECTIVITY (const CONNECTIVITY & m):
00137   //----------------------------------------------------//
00138   _entity           (m._entity),
00139   _typeConnectivity (m._typeConnectivity),
00140   _numberOfTypes    (m._numberOfTypes),
00141   _entityDimension  (m._entityDimension),
00142   _numberOfNodes    (m._numberOfNodes),
00143   _isDescendingConnectivityPartial(m._isDescendingConnectivityPartial)
00144 {
00145   if (m._geometricTypes != NULL)
00146   {
00147     _geometricTypes = new medGeometryElement[_numberOfTypes];
00148     memcpy(_geometricTypes,m._geometricTypes,_numberOfTypes*sizeof(medGeometryElement));
00149   }
00150   else
00151     _geometricTypes = (medGeometryElement *) NULL;
00152 
00153   if (m._type != NULL)
00154   {
00155     _type = new CELLMODEL[_numberOfTypes];
00156     for (int i=0;i<_numberOfTypes;i++)
00157       _type[i] = CELLMODEL(m._type[i]);
00158   }
00159   else
00160     _type = (CELLMODEL *) NULL;
00161 
00162   if (m._count != NULL)
00163   {
00164     _count = new int[_numberOfTypes+1];
00165     memcpy(_count,m._count,(_numberOfTypes+1)*sizeof(int));
00166   }
00167   else
00168     _count = (int *) NULL;
00169 
00170   if (m._nodal != NULL)
00171     _nodal = new MEDSKYLINEARRAY(* m._nodal);
00172   else
00173     _nodal = (MEDSKYLINEARRAY *) NULL;
00174 
00175   if (m._descending != NULL)
00176     _descending = new MEDSKYLINEARRAY(* m._descending);
00177   else
00178     _descending = (MEDSKYLINEARRAY *) NULL;
00179 
00180   if (m._reverseNodalConnectivity != NULL)
00181     _reverseNodalConnectivity = new MEDSKYLINEARRAY(* m._reverseNodalConnectivity);
00182   else
00183     _reverseNodalConnectivity = (MEDSKYLINEARRAY *) NULL;
00184 
00185   if (m._reverseDescendingConnectivity != NULL)
00186     _reverseDescendingConnectivity = new MEDSKYLINEARRAY(* m._reverseDescendingConnectivity);
00187   else
00188     _reverseDescendingConnectivity = (MEDSKYLINEARRAY *) NULL;
00189 
00190   if (m._neighbourhood != NULL)
00191     _neighbourhood = new MEDSKYLINEARRAY(* m._neighbourhood);
00192   else
00193     _neighbourhood = (MEDSKYLINEARRAY *) NULL;
00194 
00195   if (m._constituent != NULL)
00196     _constituent = new CONNECTIVITY(* m._constituent);
00197   else
00198     _constituent = (CONNECTIVITY *) NULL;
00199 }
00200 
00204 //----------------------------//
00205 CONNECTIVITY::~CONNECTIVITY()
00206 //----------------------------//
00207 {
00208   MESSAGE_MED("Destructeur de CONNECTIVITY()");
00209 
00210   if (_geometricTypes != NULL)
00211     delete [] _geometricTypes;
00212   if (_type != NULL)
00213     delete [] _type;
00214   if (_count != NULL)
00215     delete [] _count;
00216   if (_nodal != NULL)
00217     delete _nodal;
00218   if (_descending != NULL)
00219     delete _descending;
00220   if (_reverseNodalConnectivity != NULL)
00221     delete _reverseNodalConnectivity;
00222   if (_reverseDescendingConnectivity != NULL)
00223     delete _reverseDescendingConnectivity;
00224   if (_constituent != NULL)
00225     delete _constituent;
00226 }
00227 
00234 //----------------------------------------------------------//
00235 void CONNECTIVITY::setConstituent(CONNECTIVITY * Constituent)
00236   throw (MEDEXCEPTION)
00237 //----------------------------------------------------------//
00238 {
00239   medEntityMesh Entity = Constituent->getEntity();
00240   if (Entity == MED_CELL)
00241     throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : could not set constituent on MED_CELL !"));
00242 
00243   if ((Entity == MED_EDGE)&(_entityDimension == 3))
00244   {
00245     if (_constituent == NULL)
00246       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setConstituent : Entity not found !"));
00247     _constituent->setConstituent(Constituent);
00248   }
00249   else
00250   {
00251     delete _constituent;
00252     _constituent = Constituent;
00253   }
00254 }
00255 
00259 //----------------------------------------------------------//
00260 const CONNECTIVITY *CONNECTIVITY::getConstituent(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
00261 //----------------------------------------------------------//
00262 {
00263   if(Entity==MED_FACE && _entityDimension==3)
00264     return _constituent;
00265   if(Entity==MED_EDGE && _entityDimension==2)
00266     return _constituent;
00267   if(Entity==MED_EDGE && _entityDimension==3)
00268     {
00269       if ( _constituent )
00270         return _constituent->getConstituent(MED_EDGE);
00271       else
00272         throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getConstituent : FACE connectivity not defined!"));
00273     }
00274   throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getConstituent : Invalid request !"));
00275 }
00276 
00278 //--------------------------------------------------------------------//
00279 void CONNECTIVITY::setGeometricTypes(const medGeometryElement * Types,
00280                                      const medEntityMesh Entity)
00281   throw (MEDEXCEPTION)
00282 //--------------------------------------------------------------------//
00283 {
00284   if (Entity == _entity)
00285     for (int i=0; i<_numberOfTypes; i++)
00286     {
00287       _geometricTypes[i] = Types[i];
00288       _type[i] = CELLMODEL(_geometricTypes[i]);
00289       if ( _type[i].getDimension() > _entityDimension )
00290         _entityDimension = _type[i].getDimension();
00291     }
00292   else
00293   {
00294     if (_constituent == NULL)
00295       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setGeometricTypes : Entity not found !"));
00296     _constituent->setGeometricTypes(Types,Entity);
00297   }
00298 }
00299 
00301 //--------------------------------------------------------------------//
00302 void CONNECTIVITY::setCount(const int *         Count,
00303                             const medEntityMesh Entity)
00304   throw (MEDEXCEPTION)
00305 //--------------------------------------------------------------------//
00306 {
00307   if (Entity == _entity)
00308   {
00309     // this exception was added at 1.18.2.2.2.4 for NPAL17043: "Correction of
00310     // MEDMEM CPPUNIT tests. Integrated a work of V. Bergeaud and A. Geay."
00311     // and removed for PAL19744: "regression in MEDMEM_Connectivity.cxx"
00312     // Commenting this exception at least looks safe as the case
00313     // _numberOfTypes==0 is previewed here
00314     //if (_numberOfTypes==0)
00315     //  throw MEDEXCEPTION("Number of Types was not set before setting counts");
00316     int * index = new int[Count[_numberOfTypes]];
00317     index[0]=1;
00318     _count[0]=1;
00319     for (int i=0; i<_numberOfTypes; i++)
00320     {
00321       _count[i+1] = Count[i+1];
00322       int NumberOfNodesPerElement = _type[i].getNumberOfNodes() ;
00323       if ( _geometricTypes[i] != MED_POLYGON && _geometricTypes[i] != MED_POLYHEDRA )
00324         for (int j=_count[i]; j<_count[i+1]; j++)
00325           index[j] = index[j-1]+NumberOfNodesPerElement;
00326       else
00327         index[_count[_numberOfTypes]-1] = index[_count[_numberOfTypes-1]-1];
00328     }
00329     // allocate _nodal
00330     if (_nodal != NULL) delete _nodal;
00331     if (_numberOfTypes != 0)
00332     {
00333       _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,index[_count[_numberOfTypes]-1]-1);
00334       _nodal->setIndex(index);
00335     }
00336     else
00337       _nodal=0;
00338     delete [] index;
00339   }
00340   else
00341   {
00342     if (_constituent == NULL)
00343       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setCount : Entity not found !"));
00344     _constituent->setCount(Count,Entity);
00345   }
00346 }
00347 
00348 //--------------------------------------------------------------------//
00349 void CONNECTIVITY::setNodal(const int *              Connectivity,
00350                             const medEntityMesh      Entity,
00351                             const medGeometryElement Type,
00352                             const int *              PolyConnectivityIndex)
00353   throw (MEDEXCEPTION)
00354 //--------------------------------------------------------------------//
00355 {
00356   if (_entity == Entity)
00357   {
00358     // find geometric type
00359     bool find = false;
00360     for (int i=0; i<_numberOfTypes; i++)
00361       if (_geometricTypes[i] == Type)
00362       {
00363         if ( _count[i+1] == 0 )
00364           throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : "
00365                                        "number of elements by type must be set before!"));
00366         if ( _count[i+1] < _count[i] )
00367           throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : "
00368                                        "invalid number of elements by type!"));
00369 
00370         find = true;
00371         if ( Type == MED_POLYGON || Type == MED_POLYHEDRA )
00372         {
00373           if (! PolyConnectivityIndex )
00374             throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : PolyConnectivityIndex must be provided for poly- elements !"));
00375           if ( _numberOfTypes == 1 )
00376           {
00377             delete _nodal;
00378             _nodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
00379                                          PolyConnectivityIndex[_count[_numberOfTypes]-1]-1,
00380                                          PolyConnectivityIndex, Connectivity);
00381           }
00382           else
00383           {
00384             // append poly connectivity to the present one
00385             int nbpoly = getNumberOf( Entity, Type );
00386             int polyconnlen = PolyConnectivityIndex[nbpoly] - PolyConnectivityIndex[0];
00387             int newconnlen = _nodal->getLength() + polyconnlen;
00388             // create new MEDSKYLINEARRAY
00389             int * newconn  = new int[ newconnlen ];
00390             int * newindex = new int[ _nodal->getNumberOf() + 1 ];
00391             MEDSKYLINEARRAY* newnodal = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
00392                                                             newconnlen,newindex,newconn,
00393                                                             /*shallowCopy=*/true);
00394             // copy old values
00395             memcpy(newconn, _nodal->getValue(),sizeof(int)*(_nodal->getLength()));
00396             memcpy(newindex,_nodal->getIndex(),sizeof(int)*(_nodal->getNumberOf() - nbpoly + 1));
00397             // append new ones
00398             memcpy(newconn+_nodal->getLength(),Connectivity, sizeof(int)*polyconnlen);
00399             int* polyindex = newindex + _count[i] - 1;
00400             int delta = polyindex[0] - PolyConnectivityIndex[0];
00401             for( int j=0;j<nbpoly; j++)
00402               polyindex[j+1] = PolyConnectivityIndex[j+1] + delta;
00403             delete _nodal;
00404             _nodal = newnodal;
00405           }
00406         }
00407         else
00408         {
00409           const int* index = _nodal->getIndex();
00410           for( int j=_count[i];j<_count[i+1]; j++)
00411             _nodal->setI(j,Connectivity+index[j]-index[_count[i]]);
00412         }
00413       }
00414     if (!find)
00415       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : geometric type not found"));
00416   }
00417   else
00418   {
00419     if (_constituent == NULL)
00420       throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::setNodal : Entity not found !"));
00421     _constituent->setNodal(Connectivity,Entity,Type,PolyConnectivityIndex);
00422   }
00423 }
00424 
00426 //------------------------------------------------------------------------------------------//
00427 void CONNECTIVITY::calculateConnectivity(medConnectivity ConnectivityType, medEntityMesh Entity)
00428 //------------------------------------------------------------------------------------------//
00429 {
00430   MESSAGE_MED("CONNECTIVITY::calculateConnectivity");
00431 
00432   // a temporary limitation due to calculteDescendingConnectivity function !!!
00433   if ((_entityDimension==3) & (Entity==MED_EDGE))
00434     throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build EDGE !");
00435 
00436   if (Entity==_entity)
00437   {
00438     if (ConnectivityType==MED_NODAL)
00439       calculateNodalConnectivity();
00440     else
00441       if (Entity==MED_CELL)
00442         calculateDescendingConnectivity();
00443       else
00444         throw MEDEXCEPTION("CONNECTIVITY::calculateConnectivity : Could not build DESCENDING !");
00445   }
00446   if (Entity!=_entity) {
00447     calculateDescendingConnectivity();
00448     if (_entityDimension == 2 || _entityDimension == 3)
00449       _constituent->calculateConnectivity(ConnectivityType,Entity);
00450   }
00451 }
00452 
00453 //------------------------------------------------------------//
00454 void CONNECTIVITY::updateFamily(const vector<FAMILY*>& myFamilies)
00455 //------------------------------------------------------------//
00456 {
00457   const char * LOC = "CONNECTIVITY::updateFamily(vector<FAMILY*>) ";
00458   int numberOfFamilies = myFamilies.size();
00459   if (numberOfFamilies == 0 || _constituent == NULL)
00460     return;
00461   // does we do an update ?
00462   if ((_constituent != NULL) && (_descending != NULL))
00463     return;
00464   if (myFamilies[0]->getEntity() != _constituent->getEntity())
00465     return;
00466 
00467   CONNECTIVITY * oldConstituent = _constituent;
00468   _constituent = (CONNECTIVITY *)NULL;
00469   if (oldConstituent->_nodal==NULL)
00470     throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"We have no nodal connectivity of sub cell"));
00471 
00472   const int       oldNumberOfFace = oldConstituent->_nodal->getNumberOf();
00473   const int * oldConstituentValue = oldConstituent->_nodal->getValue();
00474   const int * oldConstituentIndex = oldConstituent->_nodal->getIndex();
00475 
00476   int * renumberingFromOldToNew = new int [oldNumberOfFace+1];
00477 
00478   calculateDescendingConnectivity();//perform calculation of descending connectivity to update all connectivities after having taken into account of direction of faces...
00479   _constituent->calculateReverseNodalConnectivity();//getting d-1 nodal connectivity to get new face numbers from nodes numbers...
00480 
00481   const int * reverseFaceNodal      = _constituent->getReverseNodalConnectivity();
00482   const int * reverseFaceNodalIndex = _constituent->getReverseNodalConnectivityIndex();
00483 
00484   CELLMODEL * aCELLMODEL = &oldConstituent->_type[0];
00485   //int newNumberOfFaces = _constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS);
00486   for ( int iOldFace = 0; iOldFace < oldNumberOfFace; iOldFace++ )
00487   {
00488     // Retrieving new number of iOldFace face...
00489 
00490     const int *nodesOfCurrentFaceOld=oldConstituentValue+oldConstituentIndex[iOldFace]-1;
00491     int    nbOfNodesOfCurrentFaceOld=oldConstituentIndex[iOldFace+1]-oldConstituentIndex[iOldFace];
00492 
00493     // get all faces (newFaceNb) around the first node of iOldFace
00494     int sizeOfNewFaceNb=reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1];
00495     int *newFaceNb=new int[sizeOfNewFaceNb];
00496     memcpy(newFaceNb,reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[0]-1]-1,sizeOfNewFaceNb*sizeof(int));//all faces candidates are in newFaceNb !!!
00497 
00498     // loop on the rest nodes of iOldFace to find, among newFaceNb, a sole face including all nodes
00499     for(int curNode=1;curNode<nbOfNodesOfCurrentFaceOld && sizeOfNewFaceNb>1;curNode++)
00500     {
00501       const int *newFacesNbForCurNode=reverseFaceNodal+reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1]-1;
00502       int sizeOfNewFacesNbForCurNode=reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]]-reverseFaceNodalIndex[nodesOfCurrentFaceOld[curNode]-1];
00503       for(int i=0; i<sizeOfNewFaceNb; )
00504       {
00505         bool found=false;
00506         for(int j=0; j < sizeOfNewFacesNbForCurNode && !found; j++)
00507           found = (newFacesNbForCurNode[j]==newFaceNb[i]);
00508 
00509         if(found)
00510           i++;
00511         else
00512           newFaceNb[i]=newFaceNb[--sizeOfNewFaceNb]; // set the last face at place of not found one
00513       }
00514     }
00515     if(sizeOfNewFaceNb!=1)
00516       throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"unexisting face specified!!!"));
00517     renumberingFromOldToNew[iOldFace]=newFaceNb[0];
00518     delete [] newFaceNb;
00519     //end of retrieving a new number of a face...
00520 
00521     // Compare nodes of a new face and those of an old one;
00522     // for the second order elements, only vertex nodes are compared
00523     int nbOfNodesOfCurrentFaceNew;
00524     const int *nodesOfCurrentFaceNew=_constituent->getConnectivityOfAnElement(MED_NODAL,_constituent->getEntity(), renumberingFromOldToNew[iOldFace], nbOfNodesOfCurrentFaceNew );
00525     if ( aCELLMODEL && aCELLMODEL->getNumberOfNodes() != nbOfNodesOfCurrentFaceOld )
00526     {
00527       // type changed, find a corresponding CELLMODEL
00528       int iType = 2; // 1-st type is already used at loop beginning
00529       while ( iOldFace + 1 >= oldConstituent->_count[ iType ]) // check next type
00530         ++iType;
00531       if ( oldConstituent->_numberOfTypes >= iType &&
00532            oldConstituent->_type[ iType-1 ].getNumberOfNodes() > 0 )
00533         aCELLMODEL = &oldConstituent->_type[ iType-1 ];
00534       else
00535         aCELLMODEL = 0;
00536     }
00537     int nbOfVertices = aCELLMODEL ? aCELLMODEL->getNumberOfVertexes() : nbOfNodesOfCurrentFaceOld;
00538     MEDMODULUSARRAY modulusArrayOld(nbOfVertices,nbOfNodesOfCurrentFaceOld,nodesOfCurrentFaceOld);
00539     int nbOfVerticesNew = std::min( nbOfVertices, nbOfNodesOfCurrentFaceNew );
00540     MEDMODULUSARRAY modulusArrayNew(nbOfVerticesNew,nbOfNodesOfCurrentFaceNew,nodesOfCurrentFaceNew);
00541     int retCompareNewOld=modulusArrayNew.compare(modulusArrayOld);
00542     if(retCompareNewOld==0)
00543       throw MED_EXCEPTION(LOCALIZED(STRING(LOC)<<"Uncompatible given user face with calculated existing faces"));
00544     if(retCompareNewOld==-1)
00545       invertConnectivityForAFace(renumberingFromOldToNew[iOldFace],nodesOfCurrentFaceOld);
00546   }
00547 
00548   // Updating the Family
00549   for(vector<FAMILY*>::const_iterator iter=myFamilies.begin();iter!=myFamilies.end();iter++)
00550     (*iter)->changeElementsNbs(_constituent->getEntity(), renumberingFromOldToNew);
00551 
00552   // FIX PAL13414:
00553   if ( _constituent && !_constituent->_constituent ) {
00554     _constituent->_constituent = oldConstituent->_constituent;
00555     oldConstituent->_constituent = NULL;
00556   }
00557   // END FIX PAL13414:
00558 
00559   delete oldConstituent;
00560   delete [] renumberingFromOldToNew;
00561   return;
00562 }
00563 
00568 //---------------------------------------------------------------------------------//
00569 const int * MEDMEM::CONNECTIVITY::getConnectivity(medConnectivity    ConnectivityType,
00570                                                   medEntityMesh      Entity,
00571                                                   medGeometryElement Type) const
00572 //----------------------------------------------------------------------------------//
00573 {
00574   const char * LOC = "CONNECTIVITY::getConnectivity";
00575   //  BEGIN_OF_MED(LOC);
00576 
00577   MEDSKYLINEARRAY * Connectivity;
00578   if (Entity==_entity) {
00579 
00580     if (ConnectivityType==MED_NODAL)
00581     {
00582       calculateNodalConnectivity();
00583       Connectivity=_nodal;
00584     }
00585     else
00586     {
00587       calculateDescendingConnectivity();
00588       Connectivity=_descending;
00589     }
00590 
00591     if (Connectivity!=NULL)
00592       if (Type==MED_ALL_ELEMENTS)
00593         return Connectivity->getValue();
00594       else {
00595         for (int i=0; i<_numberOfTypes; i++)
00596           if (_geometricTypes[i]==Type)
00597             //return Connectivity->getI(i+1);
00598             return Connectivity->getI(_count[i]);
00599         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
00600       }
00601     else
00602       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
00603   } else
00604     if (_constituent != NULL)
00605       return _constituent->getConnectivity(ConnectivityType,Entity,Type);
00606 
00607   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
00608 }
00609 
00610 //--------------------------------------------------------------------------
00611 int CONNECTIVITY::getConnectivityLength(medConnectivity    ConnectivityType,
00612                                         medEntityMesh      Entity,
00613                                         medGeometryElement Type) const
00614 //--------------------------------------------------------------------------
00615 {
00616   const char * LOC = "CONNECTIVITY::getConnectivity";
00617   BEGIN_OF_MED(LOC);
00618 
00619   MEDSKYLINEARRAY * Connectivity;
00620   if (Entity == _entity) {
00621     if (ConnectivityType == MED_NODAL)
00622     {
00623       calculateNodalConnectivity();
00624       Connectivity = _nodal;
00625     }
00626     else
00627     {
00628       calculateDescendingConnectivity();
00629       Connectivity = _descending;
00630     }
00631 
00632     if (Connectivity != NULL) {
00633       if (Type == MED_ALL_ELEMENTS) {
00634         return Connectivity->getLength();
00635       }
00636       else {
00637         for (int i=0; i<_numberOfTypes; i++)
00638           if (_geometricTypes[i]==Type)
00639             // issue 19983
00640             //return (_count[i+1]-_count[i])*getType(Type).getNumberOfNodes();
00641           {
00642             const int *ind=Connectivity->getIndex();
00643             return ind[_count[i+1]-1]-ind[_count[i]-1];
00644           }
00645         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Type not found !"));
00646       }
00647     }
00648     else
00649       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
00650   }
00651   else
00652     if (_constituent != NULL)
00653       return _constituent->getConnectivityLength(ConnectivityType,Entity,Type);
00654 
00655   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
00656 }
00657 
00666 //-----------------------------------------------------------------------------------------------//
00667 const int * CONNECTIVITY::getConnectivityIndex(medConnectivity ConnectivityType,
00668                                                medEntityMesh Entity) const
00669 //-----------------------------------------------------------------------------------------------//
00670 {
00671   const char * LOC = "CONNECTIVITY::getConnectivityIndex";
00672 
00673   MEDSKYLINEARRAY * Connectivity;
00674   if (Entity==_entity) {
00675 
00676     if (ConnectivityType==MED_NODAL)
00677       Connectivity=_nodal;
00678     else
00679       Connectivity=_descending;
00680 
00681     if (Connectivity!=NULL)
00682       return Connectivity->getIndex();
00683     else
00684       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Connectivity not defined !"));
00685 
00686   } else
00687     if (_constituent != NULL)
00688       return _constituent->getConnectivityIndex(ConnectivityType,Entity);
00689 
00690   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
00691 }
00692 
00693 
00697 //--------------------------------------------------------------//
00698 const CELLMODEL & CONNECTIVITY::getType(medGeometryElement Type) const
00699 //--------------------------------------------------------------//
00700 {
00701 
00702   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
00703     throw MEDEXCEPTION("CONNECTIVITY::getType : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE !");
00704   for (int i=0; i<_numberOfTypes; i++)
00705     if (_geometricTypes[i]==Type)
00706       return _type[i];
00707   throw MEDEXCEPTION("CONNECTIVITY::getType :  medGeometryElement not found !");
00708 }
00709 
00712 //------------------------------------------------------------------------//
00713 int CONNECTIVITY::getNumberOfNodesInType(medGeometryElement Type) const
00714 //------------------------------------------------------------------------//
00715 {
00716   const char * LOC = "CONNECTIVITY::getNumberOfNodesInType";
00717   BEGIN_OF_MED(LOC);
00718 
00719   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
00720     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!"));
00721   for (int i=0; i<_numberOfTypes; i++)
00722     if (_geometricTypes[i]==Type)
00723       return _type[i].getNumberOfNodes();
00724   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
00725 }
00726 
00729 //------------------------------------------------------------------------//
00730 int CONNECTIVITY::getNumberOfSubCellInType(medGeometryElement Type) const
00731 //------------------------------------------------------------------------//
00732 {
00733   if ((Type==MED_ALL_ELEMENTS) || (Type==MED_NONE))
00734     throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement must be different of MED_ALL_ELEMENTS and MED_NONE!");
00735   for (int i=0; i<_numberOfTypes; i++)
00736     if (_geometricTypes[i]==Type)
00737       return _type[i].getNumberOfConstituents(1);
00738   throw MEDEXCEPTION("CONNECTIVITY::getNumberOfSubCell : medGeometryElement not found !");
00739 }
00740 
00747 //-----------------------------------------------------------------------------------//
00748 int CONNECTIVITY::getNumberOf(medEntityMesh Entity, medGeometryElement Type) const
00749 //-----------------------------------------------------------------------------------//
00750 {
00751   if(Entity==MED_EN::MED_NODE)
00752     return _numberOfNodes;
00753   if (Entity==_entity) {
00754     if (Type==MED_EN::MED_NONE)
00755       return 0; // not defined !
00756     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement must be different of MED_NONE"));
00757     if (!existConnectivity(MED_NODAL,Entity) && !existConnectivity(MED_DESCENDING,Entity))
00758       return 0;
00759     if (Type==MED_EN::MED_ALL_ELEMENTS)
00760       return _count[_numberOfTypes]-1;
00761     for (int i=0; i<_numberOfTypes; i++)
00762       if (_geometricTypes[i]==Type)
00763         return (_count[i+1] - _count[i]);
00764     //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : medGeometryElement not found !"));
00765   } else
00766     if (_constituent != NULL)
00767       return _constituent->getNumberOf(Entity,Type);
00768 
00769   return 0; // valid if they are nothing else !
00770   //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" : Entity not defined !"));
00771 }
00772 
00774 //--------------------------------------------------------------//
00775 const int* CONNECTIVITY::getValue(medConnectivity TypeConnectivity,
00776                                   medGeometryElement Type) const
00777 //--------------------------------------------------------------//
00778 {
00779   if (TypeConnectivity == MED_NODAL)
00780   {
00781     calculateNodalConnectivity();
00782     if (Type==MED_ALL_ELEMENTS)
00783       return _nodal->getValue();
00784     for (int i=0; i<_numberOfTypes; i++)
00785       if (_geometricTypes[i]==Type)
00786         return _nodal->getI(_count[i]);
00787   }
00788   else
00789   {
00790     calculateDescendingConnectivity();
00791     if (Type==MED_ALL_ELEMENTS)
00792       return _descending->getValue();
00793     for (int i=0; i<_numberOfTypes; i++)
00794       if (_geometricTypes[i]==Type)
00795         return _descending->getI(_count[i]);
00796   }
00797   throw MEDEXCEPTION("Not found");
00798 }
00799 
00801 //---------------------------------------------------------------------//
00802 const int* CONNECTIVITY:: getValueIndex(medConnectivity TypeConnectivity) const
00803 //---------------------------------------------------------------------//
00804 {
00805   if (TypeConnectivity == MED_NODAL)
00806   {
00807     calculateNodalConnectivity();
00808     return _nodal->getIndex();
00809   }
00810   else
00811   {
00812     calculateDescendingConnectivity();
00813     return _descending->getIndex();
00814   }
00815 }
00816 
00818 //----------------------------------------------//
00819 const int* CONNECTIVITY:: getNeighbourhood() const
00820 //----------------------------------------------//
00821 {
00822   throw MEDEXCEPTION("CONNECTIVITY::getNeighbourhood : Not implemented");
00823 }
00824 
00827 //-------------------------------------------------//
00828 const int* CONNECTIVITY::getReverseNodalConnectivity() const
00829 //-------------------------------------------------//
00830 {
00831   calculateReverseNodalConnectivity();
00832   return _reverseNodalConnectivity->getValue();
00833 }
00834 
00837 //-------------------------------------------------------//
00838 const int* CONNECTIVITY::getReverseNodalConnectivityIndex() const
00839 //-------------------------------------------------------//
00840 {
00841   calculateReverseNodalConnectivity();
00842   return _reverseNodalConnectivity->getIndex();
00843 }
00844 
00848 //------------------------------------------------------//
00849 const int* CONNECTIVITY::getReverseDescendingConnectivity() const
00850 //------------------------------------------------------//
00851 {
00852   // it is in _constituent connectivity only if we are in MED_CELL
00853   // (we could not for instance calculate face-edge connectivity !)
00854   if (_entity!=MED_CELL)
00855     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivity : Error Only in MED_CELL connectivity");
00856 
00857   // we need descending connectivity
00858   calculateDescendingConnectivity();
00859   if (_reverseDescendingConnectivity==NULL)
00860     _reverseDescendingConnectivity=_descending->makeReverseArray();
00861 
00862   return _reverseDescendingConnectivity->getValue();
00863 }
00864 
00867 //-----------------------------------------------------------//
00868 const int* CONNECTIVITY::getReverseDescendingConnectivityIndex() const
00869 //-----------------------------------------------------------//
00870 {
00871   // it is in _constituent connectivity only if we are in MED_CELL
00872   if (_entity!=MED_CELL)
00873     throw MEDEXCEPTION("CONNECTIVITY::getReverseDescendingConnectivityIndex : Error Only in MED_CELL connectivity");
00874 
00875   // we need descending connectivity
00876   calculateDescendingConnectivity();
00877   return _reverseDescendingConnectivity->getIndex();
00878 }
00879 
00881 //--------------------------------------------//
00882 void CONNECTIVITY::calculateNodalConnectivity() const
00883 //--------------------------------------------//
00884 {
00885   if (_nodal==NULL)
00886   {
00887     if (_descending==NULL)
00888       throw MEDEXCEPTION("CONNECTIVITY::calculateNodalConnectivity : No connectivity found !");
00889     // calculate _nodal from _descending
00890   }
00891 }
00892 
00895 //---------------------------------------------------//
00896 void CONNECTIVITY::calculateReverseNodalConnectivity() const
00897 //---------------------------------------------------//
00898 {
00899   const char* LOC = "CONNECTIVITY::calculateReverseNodalConnectivity : ";
00900   BEGIN_OF_MED(LOC);
00901 
00902   SCRUTE_MED(_nodal);
00903   SCRUTE_MED(_reverseNodalConnectivity);
00904 
00905   if (_nodal==NULL)
00906     calculateNodalConnectivity();
00907 
00908   if(_reverseNodalConnectivity==NULL)
00909   {
00910     vector <vector <int> > reverse_connectivity;
00911     reverse_connectivity.resize(_numberOfNodes+1);
00912 
00913     // Treat all cells types
00914 
00915     const int *index=_nodal->getIndex();
00916     const int *conn =_nodal->getValue();
00917     int cell_number=1;
00918     for (int j = 0; j < _numberOfTypes; j++)
00919     {
00920       for (int k = _count[j]; k < _count[j+1]; k++, cell_number++)
00921       {
00922         // node number of the cell
00923         int nb_nodes = index[ cell_number ] - index[ cell_number-1 ];
00924         if ( _geometricTypes[ j ] == MED_EN::MED_POLYHEDRA )
00925         {
00926           set<int> nodes( conn, conn + nb_nodes );
00927           set<int>::iterator n = nodes.begin();
00928           if ( *n == -1 ) ++n;
00929           for ( ; n != nodes.end(); ++n )
00930             reverse_connectivity[ *n ].push_back( cell_number );
00931           conn += nb_nodes;
00932         }
00933         else
00934         {
00935           for (int i = 0; i < nb_nodes; ++i, ++conn)
00936             reverse_connectivity[ *conn ].push_back( cell_number );
00937         }
00938       }
00939     }
00940 
00941     // Full reverse_nodal_connectivity and reverse_nodal_connectivity_index
00942 
00943     //calculate size of reverse_nodal_connectivity array
00944     int size_reverse_nodal_connectivity  = 0;
00945     for (int i = 1; i < _numberOfNodes+1; i++)
00946       size_reverse_nodal_connectivity += reverse_connectivity[i].size();
00947 
00948     int * reverse_nodal_connectivity_index = new int[_numberOfNodes+1];
00949     int * reverse_nodal_connectivity = new int[size_reverse_nodal_connectivity];
00950 
00951     reverse_nodal_connectivity_index[0] = 1;
00952     for (int i = 1; i < _numberOfNodes+1; i++)
00953     {
00954       int size = reverse_connectivity[i].size();
00955       reverse_nodal_connectivity_index[i] = reverse_nodal_connectivity_index[i-1] + size;
00956       for (int j = 0; j < size; j++)
00957         reverse_nodal_connectivity[reverse_nodal_connectivity_index[i-1]-1+j]= reverse_connectivity[i][j];
00958     }
00959 
00960     _reverseNodalConnectivity = new MEDSKYLINEARRAY(_numberOfNodes,size_reverse_nodal_connectivity,
00961                                                     reverse_nodal_connectivity_index,
00962                                                     reverse_nodal_connectivity,true);
00963   }
00964   END_OF_MED(LOC);
00965 }
00966 
00968 //-------------------------------------------------//
00969 void CONNECTIVITY::calculateDescendingConnectivity() const
00970 //-------------------------------------------------//
00971 {
00972   const char * LOC = "CONNECTIVITY::calculateDescendingConnectivity() : ";
00973   BEGIN_OF_MED(LOC);
00974 
00975   if (_descending==NULL)
00976   {
00977     if (_nodal==NULL)
00978     {
00979       MESSAGE_MED(LOC<<"No connectivity found !");
00980       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
00981     }
00982     // calcul _descending from _nodal
00983     // we need CONNECTIVITY for constituent
00984 
00985     if (_constituent != NULL && _constituent->_nodal != NULL)
00986     {
00987       const_cast<CONNECTIVITY*>(this)->calculatePartialDescendingConnectivity();
00988     }
00989     else
00990     {
00991       const_cast<CONNECTIVITY*>(this)->calculateFullDescendingConnectivity(_entity);
00992     }
00993   }
00994 }
00995 
00997 //-------------------------------------------------//
00998 void CONNECTIVITY::calculateFullDescendingConnectivity(MED_EN::medEntityMesh Entity)
00999 //-------------------------------------------------//
01000 {
01001   const char * LOC = "CONNECTIVITY::calculateFullDescendingConnectivity() : ";
01002   BEGIN_OF_MED(LOC);
01003   if (_entity != Entity)
01004   {
01005     if (_constituent == NULL)
01006       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Entity not found !"));
01007     _constituent->calculateFullDescendingConnectivity(Entity);
01008   }
01009   else
01010   {
01011     if ( _descending==NULL || _isDescendingConnectivityPartial )
01012     {
01013       if (_nodal==NULL )
01014       {
01015         MESSAGE_MED(LOC<<"No connectivity found !");
01016         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No connectivity found !"));
01017       }
01018       delete _constituent;
01019       _constituent=0;
01020       if (_entityDimension == 3)
01021         _constituent = new CONNECTIVITY(MED_FACE);
01022       else if (_entityDimension == 2)
01023         _constituent = new CONNECTIVITY(MED_EDGE);
01024       else
01025       {
01026         MESSAGE_MED(LOC<<"We are in 1D");
01027         return;
01028       }
01029 
01030       bool hasPoly = ( _numberOfTypes && _geometricTypes[_numberOfTypes-1] == ( _entityDimension > 2 ? MED_POLYHEDRA : MED_POLYGON ));
01031       int numberOfClassicTypes = _numberOfTypes - hasPoly;
01032 
01033       _constituent->_typeConnectivity = MED_NODAL;
01034       _constituent->_numberOfNodes = _numberOfNodes;
01035       // foreach cells, we built array of constituent
01036       int DescendingSize = 0;
01037       for(int i=0; i<numberOfClassicTypes; i++)
01038         DescendingSize+=(_count[i+1]-_count[i])*_type[i].getNumberOfConstituents(1);
01039       int * descend_connectivity = new int[DescendingSize];
01040       for (int i=0; i<DescendingSize; i++)
01041         descend_connectivity[i]=0;
01042       int * descend_connectivity_index = new int[_count[numberOfClassicTypes]];
01043       if(_count[numberOfClassicTypes]>0)
01044         descend_connectivity_index[0]=1;
01045 
01046 
01047       map<medGeometryElement,int> eltsCounter;
01048       medGeometryElement ConstituentsTypes[2] = { MED_NONE, MED_NONE };
01049       int NumberOfConstituentsForeachType [2] = { 0, 0 };
01050       map<medGeometryElement,int>::iterator status;
01051       for(int i=0; i<numberOfClassicTypes; i++)
01052       {
01053         // initialize descend_connectivity_index array :
01054         int NumberOfConstituents = _type[i].getNumberOfConstituents(1);
01055         for (int j=_count[i];j<_count[i+1];j++)
01056         {
01057           descend_connectivity_index[j]=descend_connectivity_index[j-1]+NumberOfConstituents;
01058           // compute number of constituent of all cell for each type
01059           for(int k=1;k<NumberOfConstituents+1;k++)
01060           {
01061             medGeometryElement MEDType = _type[i].getConstituentType(1,k);
01062             status=eltsCounter.find(MEDType);
01063             if(status!=eltsCounter.end())
01064               (*status).second++;
01065             else
01066               eltsCounter[MEDType]=1;
01067           }
01068         }
01069       }
01070       if ( hasPoly && _entityDimension == 2 ) // there are constituent segments
01071       {
01072         status=eltsCounter.insert(make_pair( MED_SEG2,0 )).first;
01073         //status->second++; // increment zero or that there was
01074       }
01075       if(eltsCounter.size()>2)
01076       {
01077         // free memory (issue 0020411: [CEA 342] Sigsegv on gibi writing of MESH coming from MED file without face computation)
01078         delete [] descend_connectivity;
01079         delete [] descend_connectivity_index;
01080         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<" Descending connectivity does not support more than 2 types."));
01081       }
01082       status=eltsCounter.begin();
01083       if(!eltsCounter.empty())
01084       {
01085         ConstituentsTypes[0]=(*status).first; NumberOfConstituentsForeachType[0]=(*status).second;
01086         if(++status!=eltsCounter.end())
01087           ConstituentsTypes[1]=(*status).first, NumberOfConstituentsForeachType[1]=(*status).second;
01088       }
01089       int TotalNumberOfConstituents = NumberOfConstituentsForeachType[0]+NumberOfConstituentsForeachType[1];
01090       int TotalNumberOfNodes = NumberOfConstituentsForeachType[0]*(ConstituentsTypes[0]%100)+NumberOfConstituentsForeachType[1]*(ConstituentsTypes[1]%100);
01091 
01092       int * ConstituentNodalConnectivity      = new int[TotalNumberOfNodes];
01093       int * ConstituentNodalConnectivityIndex = new int[TotalNumberOfConstituents+1];
01094       ConstituentNodalConnectivityIndex[0]=1;
01095 
01096       _constituent->_entityDimension = _entityDimension-1;
01097       int nbClassicConstituentTypes = eltsCounter.size();
01098       int nbConstituentTypes = nbClassicConstituentTypes + hasPoly;
01099       _constituent->_numberOfTypes  = nbConstituentTypes;
01100       _constituent->_geometricTypes = new medGeometryElement[_constituent->_numberOfTypes];
01101       _constituent->_type           = new CELLMODEL[_constituent->_numberOfTypes];
01102       if(_constituent->_count) delete [] _constituent->_count;
01103       _constituent->_count          = new int[_constituent->_numberOfTypes+1];
01104       _constituent->_count[0]=1;
01105       med_int* tmp_NumberOfConstituentsForeachType = new med_int[nbClassicConstituentTypes+1];
01106       tmp_NumberOfConstituentsForeachType[0]=0; // to count constituent of each type
01107       for (int i=0; i<nbClassicConstituentTypes;i++)
01108       {
01109         _constituent->_geometricTypes[i]=ConstituentsTypes[i];
01110         _constituent->_type[i]=CELLMODEL(ConstituentsTypes[i]);
01111         _constituent->_count[i+1]=_constituent->_count[i]+NumberOfConstituentsForeachType[i];
01112         tmp_NumberOfConstituentsForeachType[i+1]=NumberOfConstituentsForeachType[i];
01113         for (int j=tmp_NumberOfConstituentsForeachType[i]; j<tmp_NumberOfConstituentsForeachType[i+1]+tmp_NumberOfConstituentsForeachType[i]; j++)
01114           ConstituentNodalConnectivityIndex[j+1]=ConstituentNodalConnectivityIndex[j]+(ConstituentsTypes[i]%100);
01115       }
01116 
01117       // we need reverse nodal connectivity
01118       if (! _reverseNodalConnectivity)
01119         calculateReverseNodalConnectivity();
01120       const int * ReverseNodalConnectivityIndex = _reverseNodalConnectivity->getIndex();
01121       const int * ReverseNodalConnectivityValue = _reverseNodalConnectivity->getValue();
01122 
01123       // array to keep reverse descending connectivity
01124       int * ReverseDescendingConnectivityValue = new int[TotalNumberOfConstituents*2];
01125 
01126       int TotalNumberOfSubCell = 0;
01127       for (int i=0; i<numberOfClassicTypes; i++)
01128       { // loop on all cell type
01129         CELLMODEL Type = _type[i];
01130         //      int NumberOfNodesPerCell = Type.getNumberOfNodes();
01131         int NumberOfConstituentPerCell = Type.getNumberOfConstituents(1);
01132         for (int j=_count[i]; j<_count[i+1]; j++) // we loop on all cell of this type
01133           for (int k=1; k<=NumberOfConstituentPerCell; k++)
01134           { // we loop on all sub cell of it
01135             if (descend_connectivity[descend_connectivity_index[j-1]+k-2]==0)
01136             { // it is a new sub cell !
01137               //      TotalNumberOfSubCell++;
01138               // Which type ?
01139               if (Type.getConstituentType(1,k)==_constituent->_geometricTypes[0])
01140               {
01141                 tmp_NumberOfConstituentsForeachType[0]++;
01142                 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[0];
01143               } else
01144               {
01145                 tmp_NumberOfConstituentsForeachType[1]++;
01146                 TotalNumberOfSubCell=tmp_NumberOfConstituentsForeachType[1];
01147               }
01148               //we have maximum two types
01149 
01150               descend_connectivity[descend_connectivity_index[j-1]+k-2]=TotalNumberOfSubCell;
01151               ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2]=j;
01152 
01153               int NumberOfNodesPerConstituent = Type.getConstituentType(1,k)%100;
01154               int * NodesLists = new int[NumberOfNodesPerConstituent];
01155               for (int l=0; l<NumberOfNodesPerConstituent; l++)
01156               {
01157                 NodesLists[l]=_nodal->getIJ(j,Type.getNodeConstituent(1,k,l+1));
01158                 ConstituentNodalConnectivity[ConstituentNodalConnectivityIndex[TotalNumberOfSubCell-1]-1+l]=NodesLists[l];
01159               }
01160               // we use reverse_nodal_connectivity to find the other element which contain this sub cell
01161 
01162               // all elements which contains first node of sub cell :
01163               int ReverseNodalConnectivityIndex_0 = ReverseNodalConnectivityIndex[NodesLists[0]-1];
01164               int ReverseNodalConnectivityIndex_1 = getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[0]-1);
01165               //ReverseNodalConnectivityIndex[NodesLists[0]];
01166               int NumberOfCellsInList = ReverseNodalConnectivityIndex_1-ReverseNodalConnectivityIndex_0;
01167 
01168               if (NumberOfCellsInList > 0)
01169               { // we could have no element !
01170                 int * CellsList = new int[NumberOfCellsInList];
01171                 for (int l=ReverseNodalConnectivityIndex_0; l<ReverseNodalConnectivityIndex_1; l++)
01172                   CellsList[l-ReverseNodalConnectivityIndex_0]=ReverseNodalConnectivityValue[l-1];
01173 
01174                 // foreach node in sub cell, we search elements which are in common
01175                 // at the end, we must have only one !
01176 
01177                 for (int l=1; l<NumberOfNodesPerConstituent; l++)
01178                 {
01179                   int NewNumberOfCellsInList = 0;
01180                   int * NewCellsList = new int[NumberOfCellsInList];
01181                   for (int l1=0; l1<NumberOfCellsInList; l1++)
01182                     for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<getIndexOfEndClassicElementInReverseNodal(ReverseNodalConnectivityValue,ReverseNodalConnectivityIndex,NodesLists[l]-1); l2++)
01183                       //for (int l2=ReverseNodalConnectivityIndex[NodesLists[l]-1]; l2<ReverseNodalConnectivityIndex[NodesLists[l]]; l2++)
01184                     {
01185                       if (CellsList[l1]<ReverseNodalConnectivityValue[l2-1])
01186                         // increasing order : CellsList[l1] are not in elements list of node l
01187                         break;
01188                       if ((CellsList[l1]==ReverseNodalConnectivityValue[l2-1])&(CellsList[l1]!=j)) {
01189                         // we have found one
01190                         NewCellsList[NewNumberOfCellsInList]=CellsList[l1];
01191                         NewNumberOfCellsInList++;
01192                         break;
01193                       }
01194                     }
01195                   NumberOfCellsInList = NewNumberOfCellsInList;
01196 
01197                   delete [] CellsList;
01198                   CellsList = NewCellsList;
01199                 }
01200 
01201                 if (NumberOfCellsInList > 0) { // We have found some elements !
01202                   int CellNumber = CellsList[0];
01203                   //delete [] CellsList;
01204                   if (NumberOfCellsInList>1) {
01205                     // free memory (issue 0020411: [CEA 342] Sigsegv on gibi writing of MESH coming from MED file without face computation)
01206                     delete [] CellsList;
01207                     delete [] NodesLists;
01208                     delete [] ReverseDescendingConnectivityValue;
01209                     delete [] tmp_NumberOfConstituentsForeachType;
01210                     delete [] descend_connectivity;
01211                     delete [] descend_connectivity_index;
01212                     delete [] ConstituentNodalConnectivity;
01213                     delete [] ConstituentNodalConnectivityIndex;
01214                     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"More than one other Cell ("<<NumberOfCellsInList<<") !"));
01215                   }
01216 
01217                   if (NumberOfCellsInList==1)
01218                   {
01219                     ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=CellNumber;
01220 
01221                     // we search sub cell number in this cell to not calculate it another time
01222                     // which type ?
01223                     CELLMODEL Type2;
01224                     for (int l=0; l<numberOfClassicTypes; l++)
01225                       if (CellNumber < _count[l+1])
01226                       {
01227                         Type2=_type[l];
01228                         break;
01229                       }
01230                     //int sub_cell_count2 = Type2.get_entities_count(1);
01231                     //int nodes_cell_count2 = Type2.get_nodes_count();
01232                     bool find2 = false;
01233                     for (int l=1; l<=Type2.getNumberOfConstituents(1);l++) // on all sub cell
01234                     {
01235                       int counter = 0;
01236                       for (int m=1; m<=Type2.getConstituentType(1,l)%100; m++)
01237                         for (int n=1; n<=Type.getConstituentType(1,k)%100; n++)
01238                         {
01239                           if (_nodal->getIJ(CellNumber,Type2.getNodeConstituent(1,l,m)) == NodesLists[n-1])
01240                             counter++;
01241                         }
01242                       if (counter==Type.getConstituentType(1,k)%100)
01243                       {
01244                         descend_connectivity[descend_connectivity_index[CellNumber-1]+l-2]=-1*TotalNumberOfSubCell; // because, see it in other side !
01245                         find2 = true;
01246                       }
01247                       if (find2)
01248                         break;
01249                     }
01250                     if (!find2)
01251                       MESSAGE_MED(LOC<<"ERROR ERROR ERROR ERROR ERROR : we find any subcell !!!"); // exception ?
01252                   }
01253                 } else {
01254                   ReverseDescendingConnectivityValue[(TotalNumberOfSubCell-1)*2+1]=0;
01255                 }
01256                 delete [] CellsList;
01257               }
01258 
01259               delete [] NodesLists;
01260             }
01261           }
01262       }
01263       // we adjust _constituent
01264       int NumberOfConstituent=0;
01265       int SizeOfConstituentNodal=0;
01266       for (int i=0;i<nbClassicConstituentTypes; i++)
01267       {
01268         NumberOfConstituent += tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1;
01269         SizeOfConstituentNodal += (tmp_NumberOfConstituentsForeachType[i]-_constituent->_count[i]+1)*_constituent->_type[i].getNumberOfNodes();
01270       }
01271       // we built new _nodal attribute in _constituent
01272       //MEDSKYLINEARRAY *ConstituentNodal = new MEDSKYLINEARRAY(NumberOfConstituent,SizeOfConstituentNodal);
01273       //const int *ConstituentNodalValue = ConstituentNodal->getValue();
01274       //const int *ConstituentNodalIndex = ConstituentNodal->getIndex();
01275       int *ConstituentNodalValue = new int[SizeOfConstituentNodal];
01276       int *ConstituentNodalIndex = new int[NumberOfConstituent+1];
01277       ConstituentNodalIndex[0]=1;
01278       // we build _reverseDescendingConnectivity
01279       //_reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent,2*NumberOfConstituent);
01280       //const int *reverseDescendingConnectivityValue = _reverseDescendingConnectivity->getValue();
01281       //const int *reverseDescendingConnectivityIndex = _reverseDescendingConnectivity->getIndex();
01282       int *reverseDescendingConnectivityValue = new int[2*NumberOfConstituent];
01283       int *reverseDescendingConnectivityIndex = new int[NumberOfConstituent+1];
01284       reverseDescendingConnectivityIndex[0]=1;
01285 
01286       // first constituent type
01287       for(int j=0; j<tmp_NumberOfConstituentsForeachType[0]; j++)
01288       {
01289         ConstituentNodalIndex[j+1]=ConstituentNodalConnectivityIndex[j+1];
01290         for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
01291           ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[k-1];
01292 
01293         reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
01294         for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++)
01295           reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[k-1];
01296       }
01297       // second type if any
01298       if (nbClassicConstituentTypes==2)
01299       {
01300         int offset = _constituent->_count[1]-tmp_NumberOfConstituentsForeachType[0]-1;
01301         int offset1=offset*_constituent->_type[0].getNumberOfNodes();
01302         int offset2=offset*2;
01303         int NumberOfNodesPerConstituent = _constituent->_type[1].getNumberOfNodes();
01304         for(int j=tmp_NumberOfConstituentsForeachType[0]; j<NumberOfConstituent; j++)
01305         {
01306           ConstituentNodalIndex[j+1]=ConstituentNodalIndex[j]+NumberOfNodesPerConstituent;
01307           for(int k=ConstituentNodalIndex[j]; k<ConstituentNodalIndex[j+1]; k++)
01308             ConstituentNodalValue[k-1]=ConstituentNodalConnectivity[offset1+k-1];
01309 
01310           reverseDescendingConnectivityIndex[j+1]=reverseDescendingConnectivityIndex[j]+2;
01311           for(int k=reverseDescendingConnectivityIndex[j]; k<reverseDescendingConnectivityIndex[j+1]; k++)
01312             reverseDescendingConnectivityValue[k-1]=ReverseDescendingConnectivityValue[offset2+k-1];
01313         }
01314         _constituent->_count[2]=NumberOfConstituent+1;
01315         // we correct _descending to adjust face number
01316         for(int j=0;j<DescendingSize;j++)
01317           if (abs(descend_connectivity[j])>tmp_NumberOfConstituentsForeachType[0]) {
01318             if ( descend_connectivity[j] > 0 )
01319               descend_connectivity[j]-=offset;
01320             else
01321               descend_connectivity[j]+=offset;
01322           }
01323       }
01324 
01325       delete [] ConstituentNodalConnectivityIndex;
01326       delete [] ConstituentNodalConnectivity;
01327       delete [] ReverseDescendingConnectivityValue;
01328       if (_constituent->_numberOfTypes > 0)
01329         _constituent->_count[1]=tmp_NumberOfConstituentsForeachType[0]+1;
01330       delete [] tmp_NumberOfConstituentsForeachType;
01331 
01333       vector<int> Descending( descend_connectivity, descend_connectivity + DescendingSize );
01334       vector<int> DescendingIndex( descend_connectivity_index, descend_connectivity_index + _count[numberOfClassicTypes] );
01335       vector<int> Reversedescendingconnectivityvalue(reverseDescendingConnectivityValue,reverseDescendingConnectivityValue + 2*NumberOfConstituent);
01336       vector<int> Reversedescendingconnectivityindex(reverseDescendingConnectivityIndex,reverseDescendingConnectivityIndex + NumberOfConstituent);
01337       delete [] descend_connectivity;
01338       delete [] descend_connectivity_index;
01339       delete [] reverseDescendingConnectivityValue;
01340       delete [] reverseDescendingConnectivityIndex;
01341 
01342 
01343       // polygons (2D mesh)
01344 
01345       vector<int> Constituentnodalvalue(ConstituentNodalValue,ConstituentNodalValue + SizeOfConstituentNodal);
01346       vector<int> Constituentnodalindex(ConstituentNodalIndex,ConstituentNodalIndex + NumberOfConstituent+1);
01347       delete [] ConstituentNodalValue;
01348       delete [] ConstituentNodalIndex;
01349       int NumberOfNewSeg = 0;
01350 
01351       int i = _count[_numberOfTypes - ( hasPoly && _entityDimension == 2 )];
01352       int NumberOfPolygons = _count[_numberOfTypes] - i;
01353       for (; i < _count[_numberOfTypes]; i++) // for each polygon
01354       {
01355         const int * vector_begin = _nodal->getValue() + _nodal->getIndex()[i-1] - 1;
01356         int vector_size = _nodal->getIndex()[i]-_nodal->getIndex()[i-1];
01357         vector<int> myPolygon(vector_begin, vector_begin+vector_size);
01358         myPolygon.push_back(myPolygon[0]);
01359         DescendingIndex.push_back( DescendingIndex.back() + vector_size);
01360         for (int j=0; j<(int)myPolygon.size()-1; j++) // for each segment of polygon
01361         {
01362           MEDMODULUSARRAY segment_poly(2,&myPolygon[j]);
01363           int ret_compare = 0;
01364 
01365           // we search it in existing segments
01366 
01367           for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
01368           {
01369             MEDMODULUSARRAY segment(2,&Constituentnodalvalue[0] + Constituentnodalindex[k]-1);
01370             ret_compare = segment_poly.compare(segment);
01371             if (ret_compare) // segment_poly already exists
01372             {
01373               Descending.push_back(ret_compare*(k+1)); // we had it to the connectivity
01374               insert_vector(Reversedescendingconnectivityvalue, 2*k+1, i); // add polygon i to reverse descending connectivity for segment_poly (in 2nd place)
01375               break;
01376             }
01377           }
01378 
01379           // segment_poly must be created
01380 
01381           if (!ret_compare)
01382           {
01383             NumberOfNewSeg++;
01384             Descending.push_back(NumberOfConstituent+NumberOfNewSeg); // we had it to the connectivity
01385             Constituentnodalvalue.push_back(segment_poly[0]);
01386             Constituentnodalvalue.push_back(segment_poly[1]);
01387             Constituentnodalindex.push_back(Constituentnodalindex[Constituentnodalindex.size()-1] + 2); // we have only segments
01388             insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewSeg-1), i ); // add polygon i to reverse descending connectivity for segment_poly (in 1st place)
01389             insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewSeg-1, 2*(NumberOfConstituent+NumberOfNewSeg-1)+1); // idem with index
01390           }
01391         }
01392       }
01393 
01394       if ( NumberOfPolygons > 0 )
01395       {
01396         NumberOfConstituent += NumberOfNewSeg;
01397         SizeOfConstituentNodal += 2*NumberOfNewSeg;
01398       }
01399 
01400       // polyhedron (3D mesh)
01401 
01402       vector<int> Constituentpolygonsnodalvalue;
01403       vector<int> Constituentpolygonsnodalindex(1,1);
01404       int NumberOfNewFaces = 0; // by convention new faces are polygons
01405       //offset to switch between all types and classical types.
01406       //int offsetCell = getNumberOf(MED_CELL, MED_ALL_ELEMENTS);
01407       int *tabRes = new int[1000]; //temporay array for intersection calculation
01408 
01409       i = _count[_numberOfTypes - ( hasPoly && _entityDimension == 3 )];
01410       int NumberOfPolyhedron = _count[_numberOfTypes] - i;
01411       for (; i<_count[_numberOfTypes]; i++) // for each polyhedron
01412       {
01413         const int * nodes    = _nodal->getValue() + _nodal->getIndex()[i-1] - 1;
01414         const int * nodesEnd = _nodal->getValue() + _nodal->getIndex()[i  ] - 1;
01415         int nbConstituentFaces = 0;
01416         while ( nodes < nodesEnd ) // loop on faces of polyhedron 
01417         {
01418           // look for end of face connectivity
01419           const int * faceNodesEnd = nodes;
01420           while ( ++faceNodesEnd < nodesEnd && *faceNodesEnd != -1 );
01421           int myFaceNumberOfNodes = faceNodesEnd - nodes;
01422           MEDMODULUSARRAY face_poly(myFaceNumberOfNodes,nodes);
01423           // move pointer to next face while continue working with face_poly
01424           if (( nodes = faceNodesEnd ) < nodesEnd )
01425             nodes++; // skip face separator
01426 
01427           nbConstituentFaces++;
01428 
01429           int ret_compare = 0;
01430 
01431           // we search it in existing faces
01432 
01433           // we search first in TRIA3 and QUAD4
01434           if (myFaceNumberOfNodes == 3 || myFaceNumberOfNodes == 4)
01435           {
01436             int Begin = -1; // first TRIA3 or QUAD4
01437             int Number = 0; // number of TRIA3 or QUAD4
01438             for (int k=0; k<(int)Constituentnodalindex.size()-1; k++)
01439               if (Constituentnodalindex[k+1]-Constituentnodalindex[k] == myFaceNumberOfNodes)
01440               {
01441                 if (Begin == -1)
01442                   Begin = k;
01443                 Number++;
01444               }
01445 
01446             for (int k=0; k<Number; k++)
01447             {
01448               MEDMODULUSARRAY face(myFaceNumberOfNodes,&Constituentnodalvalue[0] + Constituentnodalindex[Begin+k]-1);
01449               ret_compare = face_poly.compare(face);
01450               if (ret_compare)
01451               {
01452                 Descending.push_back(ret_compare*(Begin+k+1)); // we had it to the connectivity
01453                 insert_vector(Reversedescendingconnectivityvalue, 2*(Begin+k)+1, i); // add polyhedra i to reverse descending connectivity for face_poly (in 2nd place)
01454                 break;
01455               }
01456             }
01457           }
01458 
01459           // we search last in POLYGONS
01460           if (!ret_compare)
01461           {
01462             int lgth;
01463             const int *facePolyTab=face_poly.getArray(lgth);
01464             int nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[0]] -
01465               ReverseNodalConnectivityIndex[facePolyTab[0]-1];
01466             const int *candidatesCell = ReverseNodalConnectivityValue +
01467               ReverseNodalConnectivityIndex[facePolyTab[0]-1] - 1;
01468             memcpy(tabRes,candidatesCell,nbOfCandidatesCell*sizeof(int));
01469             int lgth2=nbOfCandidatesCell;
01470             for (int k=1;k<lgth && lgth2!=0;k++)
01471             {
01472               nbOfCandidatesCell = ReverseNodalConnectivityIndex[facePolyTab[k]] -
01473                 ReverseNodalConnectivityIndex[facePolyTab[k]-1];
01474               candidatesCell = ReverseNodalConnectivityValue +
01475                 ReverseNodalConnectivityIndex[facePolyTab[k]-1] - 1;
01476               mergeOrderedTabs(tabRes,lgth2,candidatesCell,nbOfCandidatesCell,tabRes,lgth2);
01477             }
01478             if (lgth2<=1)
01479               ret_compare=0;//here normally tabRes[0]==i
01480             else //> 2 should never happend : A face is shared by more than 2 polyhedrons...
01481             {
01482               if (tabRes[0] == i) //as tabRes is ordered by construction tabRes[1] > tabRes[0] so the current
01483                 // face is shared with an another cell whose id > current id. So let's create
01484                 ret_compare=0;
01485               else
01486               {//tabRes[0]<Constituentpolygonsnodalindex.size()-1 that is to say the current face has been built previously.
01487                 int nbOfFacesConstitutingAlreadyBuiltPolyh;
01488                 int* nbOfNodesPerFaces;
01489                 int** nodesOfFacesOfAlreadyBuiltPolyh = getNodesPerFaceOfPolyhedron( tabRes[0], nbOfFacesConstitutingAlreadyBuiltPolyh, nbOfNodesPerFaces);
01490                 for (int k1=0; k1<nbOfFacesConstitutingAlreadyBuiltPolyh && (ret_compare==0); k1++)
01491                 {
01492                   int nbOfNodesForCurrentFace = nbOfNodesPerFaces[k1];
01493                   MEDMODULUSARRAY face(nbOfNodesForCurrentFace,nodesOfFacesOfAlreadyBuiltPolyh[k1]);
01494                   ret_compare = face_poly.compare(face);
01495                   if (ret_compare)
01496                   {
01497                     int curFaceId=Descending[ DescendingIndex[ tabRes[0]-1 ]-1 + k1 ];
01498                     Descending.push_back(ret_compare*curFaceId); // we had it to the connectivity
01499                     insert_vector(Reversedescendingconnectivityvalue, 2*(curFaceId-1)+1, i);
01500                   }
01501                 }
01502                 delete [] nbOfNodesPerFaces;
01503                 delete [] nodesOfFacesOfAlreadyBuiltPolyh;
01504               }
01505             }
01506           }
01507 
01508           // if not found, face_poly must be created
01509 
01510           if (!ret_compare)
01511           {
01512             NumberOfNewFaces++;
01513             Descending.push_back(NumberOfConstituent+NumberOfNewFaces); // we had it to the connectivity
01514             for (int k=0; k<myFaceNumberOfNodes; k++)
01515               Constituentpolygonsnodalvalue.push_back(face_poly[k]);
01516             Constituentpolygonsnodalindex.push_back(Constituentpolygonsnodalindex[Constituentpolygonsnodalindex.size()-1]+myFaceNumberOfNodes);
01517             insert_vector(Reversedescendingconnectivityvalue, 2*(NumberOfConstituent+NumberOfNewFaces-1), i); // add polyhedra i to reverse descending connectivity for face_poly (in 1st place)
01518             insert_vector(Reversedescendingconnectivityindex, NumberOfConstituent+NumberOfNewFaces-1, 2*(NumberOfConstituent+NumberOfNewFaces-1)+1); // idem with index
01519           }
01520         }  // loop on faces of polyhedron
01521 
01522         DescendingIndex.push_back( DescendingIndex.back() + nbConstituentFaces );
01523       }
01524       delete [] tabRes;
01525 
01526       if (NumberOfPolyhedron > 0)
01527       {
01528         NumberOfConstituent += NumberOfNewFaces;
01529         Constituentnodalvalue.insert( Constituentnodalvalue.end(),
01530                                       Constituentpolygonsnodalvalue.begin(),
01531                                       Constituentpolygonsnodalvalue.end());
01532         Constituentpolygonsnodalvalue.clear();
01533         Constituentnodalindex.reserve( Constituentnodalindex.size() + NumberOfNewFaces );
01534         int indexShift = Constituentnodalindex.back() - Constituentpolygonsnodalindex.front();
01535         for ( unsigned i = 1; i < Constituentpolygonsnodalindex.size(); ++i )
01536           Constituentnodalindex.push_back( Constituentpolygonsnodalindex[i] + indexShift );
01537         Constituentpolygonsnodalindex.clear();
01538         _constituent->_geometricTypes[ _constituent->_numberOfTypes-1 ] = MED_POLYGON;
01539       }
01540 
01541       _constituent->_count[ _constituent->_numberOfTypes ] = NumberOfConstituent+1;
01542       _constituent->_nodal = new MEDSKYLINEARRAY(NumberOfConstituent,
01543                                                  Constituentnodalvalue.size(),
01544                                                  &Constituentnodalindex[0],
01545                                                  &Constituentnodalvalue[0]);
01546 
01547       if(_descending)
01548         delete _descending;
01549       _descending = new MEDSKYLINEARRAY(_count[_numberOfTypes]-1,
01550                                         Descending.size(),
01551                                         &DescendingIndex[0],
01552                                         &Descending[0]);
01553 
01554       Reversedescendingconnectivityindex.push_back(Reversedescendingconnectivityindex[Reversedescendingconnectivityindex.size()-1]+2); // we complete the index
01555       Reversedescendingconnectivityvalue.push_back(0);
01556       _reverseDescendingConnectivity = new MEDSKYLINEARRAY(NumberOfConstituent+1,
01557                                                            2*NumberOfConstituent,
01558                                                            &Reversedescendingconnectivityindex[0],
01559                                                            &Reversedescendingconnectivityvalue[0]);
01560       _isDescendingConnectivityPartial = false;
01561     }
01562     END_OF_MED(LOC);
01563   }
01564 }
01565 
01566 void CONNECTIVITY::addToDescendingConnectivity(const set<int>& nodes,
01567                                                multimap<int,int>& descending,
01568                                                int iglobal_cell,
01569                                                const CONNECTIVITY_HashMap & face_map ) const
01570 {
01571   int dimension = getEntityDimension();
01572   vector<int> signature (dimension);
01573   set<int>::const_iterator iter=nodes.begin();
01574   for (int i=0; i< dimension;i++)
01575   {
01576     signature[i]=*iter;
01577     iter++;
01578   }
01579 
01580 
01581   CONNECTIVITY_HashMap::const_iterator itermap=face_map.find(signature);
01582   CONNECTIVITY_HashMap::const_iterator iterend=face_map.end();
01583 
01584 
01585   if (itermap!=iterend)
01586     descending.insert(make_pair(iglobal_cell,itermap->second));
01587 }
01588 
01628 void CONNECTIVITY::calculatePartialDescendingConnectivity() const
01629 {
01631   // First stage : creating hash_map referencing faces with the triplet
01632   // of the lowest order nodes as a key and the global face number as a value
01633 
01634   CONNECTIVITY_HashMap face_map;
01635   int iglobal_face=1;
01636 
01637   int dimension = getEntityDimension();
01638   {
01639     const int* conn =_constituent->getValue     (MED_NODAL, MED_ALL_ELEMENTS);
01640     const int* index=_constituent->getValueIndex(MED_NODAL);
01641     int nbfaces=_constituent->getNumberOf(_constituent->getEntity(),MED_ALL_ELEMENTS);
01642     for ( ; iglobal_face <= nbfaces; iglobal_face++ )
01643     {
01644       const int* face_conn = conn + index[iglobal_face-1] - 1;
01645       const int nbnodes = index[iglobal_face] - index[iglobal_face-1];
01646       set<int> nodes( face_conn, face_conn + nbnodes );
01647       vector<int> signature (dimension);
01648       set<int>::iterator iter=nodes.begin();
01649       for (int i=0; i< dimension;i++)
01650       {
01651         signature[i]=*iter;
01652         iter++;
01653       }
01654       face_map.insert(make_pair(signature,iglobal_face));
01655     }
01656   }
01658   //Second stage : going through all the faces of the cells and
01659   // connecting them to the hash_map created in the first stage
01660 
01661   multimap<int,int> descending; //map storing the descending connectivity
01662 
01663   int nbcell_types = getNumberOfTypes(_entity);
01664   const medGeometryElement* cell_types = getGeometricTypes(_entity);
01665   int iglobal_cell = 1;
01666   for (int itype=0; itype<nbcell_types; itype++)
01667   {
01668     medGeometryElement cell_type = cell_types[itype];
01669     int nbcells = getNumberOf(_entity,cell_type);
01670 
01671     const int* index = _nodal->getIndex();
01672     const int* conn  = _nodal->getValue();
01673     switch ( cell_type )
01674     {
01675     case MED_POLYGON:
01676       {
01677         for (int icell=0;icell<nbcells;icell++)
01678         {
01679           int nbnodes=index[icell+1]-index[icell];
01680           int nbfaces=nbnodes;
01681           for (int iface=0; iface<nbfaces;iface++)
01682           {
01683             set<int> nodes;
01684             if (iface+1!=nbfaces)
01685             {
01686               nodes.insert(conn[index[icell]-1+iface]);
01687               nodes.insert(conn[index[icell]-1+iface+1]);
01688             }
01689             else
01690             {
01691               nodes.insert(conn[index[icell]-1+iface]);
01692               nodes.insert(conn[index[icell]-1]);
01693             }
01694             addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
01695           }
01696           iglobal_cell++;
01697         }
01698       }
01699       break;
01700 
01701     case MED_POLYHEDRA:
01702       {
01703         for (int icell = 0; icell < nbcells; icell++)
01704         {
01705           const int* face_conn     = conn + index[icell  ] - 1;
01706           const int* face_conn_end = conn + index[icell+1] - 1;
01707           while (face_conn < face_conn_end)
01708           {
01709             set<int> nodes; // faces are divided by -1
01710             for ( ; face_conn < face_conn_end && *face_conn != -1; face_conn++ )
01711               nodes.insert( *face_conn );
01712             addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
01713           }
01714           iglobal_cell++;
01715         }
01716       }
01717       break;
01718 
01719     default: // classic cells
01720       {
01721         CELLMODEL cellmodel=CELLMODEL_Map::retrieveCellModel(cell_type);
01722         const int* index=_nodal->getIndex();
01723         const int* conn=_nodal->getValue();
01724         for (int icell = 0; icell < nbcells; icell++, iglobal_cell++)
01725         {
01726           int nbfaces=cellmodel.getNumberOfConstituents(1);
01727           for (int iface=0; iface<nbfaces;iface++)
01728           {
01729             set<int> nodes;
01730             const int* local_index=cellmodel.getNodesConstituent(1,iface+1);
01731             medGeometryElement face_type = cellmodel.getConstituentType(1,iface+1);
01732             int nbnodes=face_type%100;
01733             for (int inode=0;inode<nbnodes;inode++)
01734               nodes.insert(conn[index[iglobal_cell-1]-1+local_index[inode]-1]);
01735             addToDescendingConnectivity(nodes,descending,iglobal_cell,face_map);
01736           }
01737         }
01738       }
01739     }
01740   }
01742   // Third stage : reorganizing the descending data to store it in a medskylinearray
01743 
01744   vector<int> index;
01745   vector<int> value;
01746   index.push_back(1);
01747 
01748   //the number of cells is given by the number
01749   //obtained by browsing all the cell types
01750   int nb_cells = iglobal_cell-1;
01751 
01752   //for (int icell = 0; icell < nb_cells; icell++)
01753   for (int icell = 1; icell <= nb_cells; icell++)
01754   {
01755     multimap<int,int>::iterator beginning_of_range = descending.lower_bound(icell);
01756     multimap<int,int>::iterator end_of_range = descending.upper_bound(icell);
01757     int nb=0;
01758     for (multimap<int,int>::iterator iter2 = beginning_of_range; iter2 != end_of_range; iter2++)
01759     {
01760       value.push_back(iter2->second);
01761       nb++;
01762     }
01763     index.push_back(index.back()+nb);
01764   }
01765   if(_descending)
01766     delete _descending;
01767   _descending = new MEDSKYLINEARRAY(index.size()-1, value.size(), &index[0], &value[0]);
01768   _isDescendingConnectivityPartial = true;
01769 }
01770 
01771 
01773 //--------------------------------------------------------------------//
01774 void CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity)
01775 //--------------------------------------------------------------------//
01776 {
01777   const char* LOC = "CONNECTIVITY::calculateNeighbourhood(CONNECTIVITY &myConnectivity) : ";
01778   BEGIN_OF_MED(LOC);
01779 
01780   MESSAGE_MED(PREFIX_MED<<"method not yet implemented " << myConnectivity._entity);
01781   // Mesh dimension !
01782 
01783   END_OF_MED(LOC);
01784   return;
01785 }
01786 
01787 
01793 //--------------------------------------------------------------------//
01794 medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity,int globalNumber) const
01795 //--------------------------------------------------------------------//
01796 {
01797   const char * LOC = "medGeometryElement CONNECTIVITY::getElementType(medEntityMesh Entity, int globalNumber) const : ";
01798   BEGIN_OF_MED(LOC);
01799   int globalNumberMin = 1;
01800   int globalNumberMax ;
01801 
01802   if (_entity==Entity) globalNumberMax = _count[_numberOfTypes];
01803   else if (_constituent!=NULL) globalNumberMax = _constituent->_count[_constituent->_numberOfTypes];
01804   else
01805     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
01806 
01807   // The globalNumber must verify  : 1 <=  globalNumber < _count[_numberOfTypes] (== totalNumberOfElement+1)
01808 
01809   if ( (globalNumber < globalNumberMin) || (globalNumber >  globalNumberMax-1 )  )
01810     throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "globalNumber must be between >= |"
01811                                  << globalNumberMin <<"| and <= |" << globalNumberMax-1 ));
01812 
01813   if (_entity==Entity) {
01814     for(int i=1; i<=_numberOfTypes;i++)
01815       if (globalNumber<_count[i])
01816         return _geometricTypes[i-1];
01817   }
01818   else if (_constituent!=NULL)
01819     return _constituent->getElementType(Entity,globalNumber);
01820   else
01821     throw MEDEXCEPTION(LOCALIZED("getElementType : Entity not defined !"));
01822   throw MEDEXCEPTION(LOCALIZED("getElementType : Wrong Number !"));
01823 
01824   END_OF_MED(LOC);
01825 }
01826 
01827 ostream & MEDMEM::operator<<(ostream &os, CONNECTIVITY &co)
01828 {
01829   os << endl << "------------- Entity = ";
01830   switch (co._entity)
01831   {
01832   case MED_CELL:
01833     os << "MED_CELL";
01834     break;
01835   case MED_FACE:
01836     os << "MED_FACE";
01837     break;
01838   case MED_EDGE:
01839     os << "MED_EDGE";
01840     break;
01841   case MED_NODE:
01842     os << "MED_NODE";
01843     break;
01844   case MED_ALL_ENTITIES:
01845     os << "MED_ALL_ENTITIES";
01846     break;
01847   default:
01848     os << "Unknown";
01849     break;
01850   }
01851   os  << " -------------\n\nMedConnectivity : ";
01852   switch (co._typeConnectivity)
01853   {
01854   case MED_NODAL:
01855     os << "MED_NODAL\n";
01856     break;
01857   case MED_DESCENDING:
01858     os << "MED_DESCENDING\n";
01859     break;
01860   default:
01861     break;
01862   }
01863   os << "Entity dimension : " << co._entityDimension << endl;
01864   os << "Number of nodes : " << co._numberOfNodes << endl;
01865   os << "Number of types : " << co._numberOfTypes << endl;
01866   for (int i=0; i!=co._numberOfTypes ; ++i)
01867     os << "  -> Type " << co._geometricTypes[i] << " (" << co._type[i].getName() << ") : "
01868        << co._count[i+1]-co._count[i] << " elements" << endl;
01869 
01870   if (co._typeConnectivity == MED_NODAL )
01871   {
01872     for (int i=0; i<co._numberOfTypes; i++)
01873     {
01874       os << endl << co._type[i].getName() << " : " << endl;
01875       int numberofelements = co._count[i+1]-co._count[i];
01876       const int * connectivity = co.getConnectivity(co._typeConnectivity, co._entity, MED_EN::MED_ALL_ELEMENTS);
01877       const int * index = co.getConnectivityIndex(co._typeConnectivity, co._entity) + co._count[i]-1;
01878       for (int j=0;j<numberofelements;j++)
01879       {
01880         const int * cellconn = connectivity+index[j]-1;
01881         const int * cellconnend = connectivity+index[j+1]-1;
01882         os << setw(6) << j+1 << " : ";
01883         while ( cellconn < cellconnend )
01884           os << *cellconn++ <<" ";
01885         os << endl;
01886       }
01887     }
01888   }
01889   else if (co._typeConnectivity == MED_DESCENDING)
01890   {
01891     int numberofelements = co.getNumberOf( co._entity , MED_ALL_ELEMENTS);
01892     if (numberofelements > 0)
01893     {
01894       const int *connectivity =  co.getConnectivity( co._typeConnectivity, co._entity, MED_ALL_ELEMENTS);
01895       const int *connectivity_index =  co.getConnectivityIndex( co._typeConnectivity, co._entity );
01896 
01897       for ( int j=0; j!=numberofelements; j++ )
01898       {
01899         os << "Element " << j+1 << " : ";
01900         for ( int k=connectivity_index[j]; k!=connectivity_index[j+1]; k++ )
01901           os << connectivity[k-1] << " ";
01902         os << endl;
01903       }
01904     }
01905   }
01906 
01907   if (co._constituent)
01908     os << endl << *co._constituent << endl;
01909 
01910   return os;
01911 }
01912 
01917 int *CONNECTIVITY::getNodesOfPolyhedron(int polyhedronId, int& lgthOfTab) const
01918 {
01919   int offsetWithClassicType=_count[ _numberOfTypes-1 ]-1; // hope that polyhedrons is the last type
01920   if (polyhedronId<=offsetWithClassicType || polyhedronId> getNumberOf (MED_CELL, MED_ALL_ELEMENTS))
01921     throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
01922 
01923   const int *nodes=getConnectivity(MED_NODAL, MED_CELL, MED_ALL_ELEMENTS);
01924   const int *index=getConnectivityIndex(MED_NODAL,MED_CELL);
01925   const int *cellnodes    = nodes + index[polyhedronId-1] - 1;
01926   const int *cellnodesend = nodes + index[polyhedronId  ] - 1;
01927   set<int> retInSet( cellnodes, cellnodesend );
01928   if ( *retInSet.begin() == -1 )
01929     retInSet.erase( retInSet.begin() );
01930   lgthOfTab=retInSet.size();
01931   int *ret=new int[lgthOfTab];
01932   set<int>::iterator iter=retInSet.begin();
01933   for(int i=0;iter!=retInSet.end();iter++)
01934     ret[i++]=*iter;
01935   return ret;
01936 }
01937 
01943 int **CONNECTIVITY::getNodesPerFaceOfPolyhedron(int polyhedronId, int& nbOfFaces, int* & nbOfNodesPerFaces) const
01944 {
01945   int offsetWithClassicType=_count[ _numberOfTypes-1 ]-1; // hope that polyhedrons is the last type
01946   if (polyhedronId<=offsetWithClassicType || polyhedronId> getNumberOf(MED_CELL, MED_ALL_ELEMENTS))
01947     throw MEDEXCEPTION("Polyhedron ID does not match a polyhedron in the element range");
01948   const int *nodes=getConnectivity(MED_EN::MED_NODAL, MED_CELL, MED_ALL_ELEMENTS);
01949   const int *index=getConnectivityIndex(MED_EN::MED_NODAL, MED_CELL);
01950 
01951   const int *cellnodes    = nodes + index[polyhedronId-1] - 1;
01952   const int *cellnodesend = nodes + index[polyhedronId  ] - 1;
01953   nbOfFaces=0;
01954   nbOfNodesPerFaces=new int[cellnodesend-cellnodes]; // more than enough
01955   int **ret=new int* [cellnodesend-cellnodes];
01956   while ( cellnodes < cellnodesend )
01957   {
01958     ret[nbOfFaces] = (int*)cellnodes;
01959     while ( cellnodes < cellnodesend && *cellnodes != -1 ) ++cellnodes;
01960     nbOfNodesPerFaces[nbOfFaces] = cellnodes - ret[nbOfFaces];
01961     cellnodes++;
01962     nbOfFaces++;
01963   }
01964   return ret;
01965 }
01966 
01967 /*
01968   Method used in CalculateDescendingConnectivity. So it's typically a private method.
01969   The aim of this method is to hide to CalculateDescendingConnectivity algorithm the fact that in reverse connectivity polygons and polyhedrons
01970   are not in separate data structure, contrary to not reverse connectivities.
01971 */
01972 int CONNECTIVITY::getIndexOfEndClassicElementInReverseNodal(const int *reverseNodalValue, const int *reverseNodalIndex, int rk)  const
01973 {
01974   const int lastTypeIsPoly = (_geometricTypes[_numberOfTypes-1]==MED_EN::MED_POLYGON||
01975                               _geometricTypes[_numberOfTypes-1]==MED_EN::MED_POLYHEDRA);
01976   if ( !lastTypeIsPoly )
01977     return reverseNodalIndex[rk+1];
01978   int nbOfLastElt=_count[_numberOfTypes-1]-1;
01979   int ret=reverseNodalIndex[rk];
01980   for(int i=reverseNodalIndex[rk];i<reverseNodalIndex[rk+1];i++)
01981   {
01982     if(reverseNodalValue[i-1]<=nbOfLastElt)
01983       ret++;
01984   }
01985   return ret;
01986 }
01987 
01988 /*
01989   Method that inverts the face with faceId 'faceId' in the data structure.
01990   This method has to be applied on a instance of MEDMEM::CONNECTIVITY with _entityDimension==3.
01991   WARNING calculateDescendingConnectivity must have been called before.
01992 */
01993 void CONNECTIVITY::invertConnectivityForAFace(int faceId, const int *nodalConnForFace)
01994 {
01995   // we have always 2 neighbourings
01996   int cell1 = _reverseDescendingConnectivity->getIJ(faceId,1);
01997   int cell2 = _reverseDescendingConnectivity->getIJ(faceId,2);
01998 
01999   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
02000     // Updating _reverseDescendingConnectivity
02001     _reverseDescendingConnectivity->setIJ(faceId,1,cell2);
02002     _reverseDescendingConnectivity->setIJ(faceId,2,cell1);
02003     // Updating _constituent->_nodal because of reversity
02004     const int *descendingNodalIndex=_constituent->_nodal->getIndex();
02005     MEDSKYLINEARRAY *currentNodal=_constituent->_nodal;
02006     int faceIdRelative=faceId;
02007     for(int iarray=1;iarray<=(descendingNodalIndex[faceIdRelative]-descendingNodalIndex[faceIdRelative-1]);iarray++)
02008       currentNodal->setIJ(faceIdRelative,iarray,nodalConnForFace[iarray-1]);
02009 
02010     // Updating _descending for cell1 and cell2
02011     const int NB_OF_CELLS_SHARING_A_FACE=2;
02012     int cellsToUpdate[NB_OF_CELLS_SHARING_A_FACE] = { cell1, cell2 };
02013     for(int curCell=0;curCell<NB_OF_CELLS_SHARING_A_FACE;curCell++)
02014     {
02015       int cell=cellsToUpdate[curCell];
02016       const int *newDescendingIndex=_descending->getIndex();
02017       MEDSKYLINEARRAY *currentDescending=_descending;
02018       for(int iface=newDescendingIndex[cell-1];iface<newDescendingIndex[cell];iface++)
02019       {
02020         int curValue=currentDescending->getIndexValue(iface);
02021         if (abs(curValue)==faceId)
02022           currentDescending->setIndexValue(iface,-curValue);
02023       }
02024     }
02025   }
02026 }
02027 
02028 /*
02029   Method with 2 output : the connectivity required and its length 'lgth'.
02030   This method gives the connectivity independently it is a polygons/polyhedrons or normal element.
02031 */
02032 const int * CONNECTIVITY::getConnectivityOfAnElement(MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity, int Number, int &lgth) const
02033 {
02034   if(Entity==MED_EN::MED_NODE)
02035     throw  MEDEXCEPTION("No connectivity attached to a node entity");
02036   if(Entity==_entity)
02037   {
02038     if ( Number > getNumberOf( Entity, MED_ALL_ELEMENTS ))
02039       throw  MEDEXCEPTION("Unknown number");
02040     const int * conn  = getConnectivity(ConnectivityType,Entity,MED_ALL_ELEMENTS);
02041     const int * index = getConnectivityIndex(ConnectivityType,Entity);
02042     lgth=index[Number]-index[Number-1];
02043     return conn+index[Number-1]-1;
02044   }
02045   else
02046   {
02047     if(_constituent==NULL)
02048       calculateDescendingConnectivity();
02049     return _constituent->getConnectivityOfAnElement(ConnectivityType,Entity,Number,lgth);
02050   }
02051 }
02052 
02056 bool CONNECTIVITY::deepCompare(const CONNECTIVITY& other) const
02057 {
02058   CONNECTIVITY* temp=(CONNECTIVITY* )this;
02059   const int *conn1=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
02060   int size1=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
02061   temp=(CONNECTIVITY* )(&other);
02062   const int *conn2=temp->getConnectivity(MED_NODAL,_entity,MED_ALL_ELEMENTS);
02063   int size2=temp->getConnectivityLength(MED_NODAL,_entity,MED_ALL_ELEMENTS);
02064   if(size1!=size2)
02065     return false;
02066   bool ret=true;
02067   for(int i=0;i<size1 && ret;i++)
02068   {
02069     ret=(conn1[i]==conn2[i]);
02070   }
02071   return ret;
02072 }