Back to index

salome-med  6.5.0
MEDMEM_Mesh.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 /*
00024  File Mesh.cxx
00025 */
00026 
00027 #include "MEDMEM_DriversDef.hxx"
00028 #include "MEDMEM_Field.hxx"
00029 #include "MEDMEM_Mesh.hxx"
00030 
00031 #include "MEDMEM_Support.hxx"
00032 #include "MEDMEM_Family.hxx"
00033 #include "MEDMEM_Group.hxx"
00034 #include "MEDMEM_Coordinate.hxx"
00035 #include "MEDMEM_Connectivity.hxx"
00036 #include "MEDMEM_CellModel.hxx"
00037 #include "VolSurfFormulae.hxx"
00038 #include "MEDMEM_DriverFactory.hxx"
00039 
00040 #include "PointLocator.hxx"
00041 
00042 #include <math.h>
00043 #include <map>
00044 #include <sstream>
00045 
00046 using namespace std;
00047 using namespace MEDMEM;
00048 using namespace MED_EN;
00049 
00050 #define MED_NOPDT -1
00051 #define MED_NONOR -1
00052 
00053 // Block defining groups for the MEDMEM_ug documentation
00074 void MESH::init()
00075 {
00076   const char* LOC = "MESH::init(): ";
00077   BEGIN_OF_MED(LOC);
00078 
00079   GMESH::init();
00080 
00081   delete _coordinate;
00082   _coordinate   = (COORDINATE   *) NULL;
00083   delete _connectivity;
00084   _connectivity = (CONNECTIVITY *) NULL;
00085 
00086   _numberOfNodes = MED_INVALID;
00087 
00088   _arePresentOptionnalNodesNumbers = 0;
00089 
00090   END_OF_MED(LOC);
00091 }
00092 
00094 MESH::MESH():GMESH(),_coordinate(NULL),_connectivity(NULL)
00095 {
00096   init();
00097 }
00098 
00107 MESH::MESH(MESH &m): GMESH(m)
00108 {
00109   if (m._coordinate != NULL)
00110     _coordinate = new COORDINATE(* m._coordinate);
00111   else
00112     _coordinate = (COORDINATE *) NULL;
00113 
00114   if (m._connectivity != NULL)
00115     _connectivity = new CONNECTIVITY(* m._connectivity);
00116   else
00117     _connectivity = (CONNECTIVITY *) NULL;
00118 
00119   _numberOfNodes = m._numberOfNodes;
00120 
00121   _arePresentOptionnalNodesNumbers = m._arePresentOptionnalNodesNumbers;
00122   _optionnalToCanonicNodesNumbers = m._optionnalToCanonicNodesNumbers;
00123 }
00124 
00131 MESH::~MESH()
00132 {
00133   MESSAGE_MED("MESH::~MESH() : Destroying the Mesh");
00134 
00135   if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
00136   if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
00137   _coordinate = 0;
00138   _connectivity = 0;
00139 }
00140 
00141 MESH & MESH::operator=(const MESH &m)
00142 {
00143   const char* LOC = "MESH & MESH::operator=(const MESH &m) : ";
00144   BEGIN_OF_MED(LOC);
00145 
00146   _name = m._name;
00147   _description = m._description;
00148 
00149   if ( _coordinate   ) delete _coordinate;
00150   _coordinate   = m._coordinate   ? new COORDINATE  ( *m._coordinate   ) : 0;
00151   if ( _connectivity ) delete _connectivity;
00152   _connectivity = m._connectivity ? new CONNECTIVITY( *m._connectivity ) : 0;
00153 
00154   _spaceDimension = m._spaceDimension;
00155   _numberOfNodes  = m._numberOfNodes;
00156 
00157   _arePresentOptionnalNodesNumbers = m._arePresentOptionnalNodesNumbers;
00158   _optionnalToCanonicNodesNumbers  = m._optionnalToCanonicNodesNumbers;
00159 
00160   vector<FAMILY*>*        fams[4] = { &_familyNode, &_familyCell, &_familyFace, &_familyEdge};
00161   const vector<FAMILY*>* mfams[4] = { &m._familyNode,&m._familyCell,&m._familyFace,&m._familyEdge };
00162   for ( int i = 0; i < 4; ++i )
00163   {
00164     for ( int f = 0; f < fams[i]->size(); ++f )
00165       fams[i]->at(f)->removeReference();
00166     fams[i]->clear();
00167     fams[i]->reserve( mfams[i]->size() );
00168     for ( int f = 0; f < mfams[i]->size(); ++f )
00169     {
00170       if ( mfams[i]->at(f) )
00171       {
00172         fams[i]->push_back( new FAMILY( *mfams[i]->at(f) ));
00173         fams[i]->back()->setMesh( this );
00174       }
00175     }
00176   }
00177   vector<GROUP*>*        groups[4] = { &_groupNode, &_groupCell, &_groupFace, &_groupEdge };
00178   const vector<GROUP*>* mgroups[4] = { &m._groupNode, &m._groupCell, &m._groupFace, &m._groupEdge };
00179   for ( int i = 0; i < 4; ++i )
00180   {
00181     for ( int g = 0; g < groups[i]->size(); ++g )
00182       groups[i]->at(g)->removeReference();
00183     groups[i]->clear();
00184     groups[i]->reserve( mgroups[i]->size() );
00185     for ( int g = 0; g < mgroups[i]->size(); ++g )
00186     {
00187       if ( mgroups[i]->at(g) )
00188       {
00189         groups[i]->push_back( new GROUP( *mgroups[i]->at(g) ));
00190         groups[i]->back()->setMesh( this );
00191       }
00192     }
00193   }
00194 
00195   for ( int drv = 0; drv < _drivers.size(); ++drv )
00196     delete _drivers[drv];
00197   _drivers.clear();
00198   _drivers.reserve( m._drivers.size());
00199   for ( int drv = 0; drv < m._drivers.size(); ++drv )
00200     if ( m._drivers[drv] )
00201       _drivers.push_back( m._drivers[drv]->copy() );
00202 
00203   return *this;
00204 }
00205 
00206 bool MESH::operator==(const MESH& other) const
00207 {
00208   const char* LOC = "MESH::operator==";
00209   BEGIN_OF_MED(LOC);
00210   return this==&other;
00211 }
00212 
00222 MESH::MESH(driverTypes driverType, const string &  fileName/*=""*/, const string &  driverName/*=""*/) throw (MEDEXCEPTION):_coordinate(0),_connectivity(0)
00223 {
00224   const char * LOC = "MESH::MESH(driverTypes driverType, const string &  fileName="", const string &  driverName/="") : ";
00225   BEGIN_OF_MED(LOC);
00226 
00227   int current;
00228   init();
00229   GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,RDONLY);
00230   current = addDriver(*myDriver);
00231   delete myDriver;
00232   try
00233   {
00234     _drivers[current]->open();
00235     _drivers[current]->read();
00236     _drivers[current]->close();
00237   }
00238   catch ( MED_EXCEPTION& e )
00239   {
00240     if ( _drivers[current] )
00241       _drivers[current]->close(), delete _drivers[current];
00242     _drivers[current] = 0;
00243     throw e;
00244   }
00245 
00246   END_OF_MED(LOC);
00247 }
00260 bool MESH::deepCompare(const GMESH& gother) const
00261 {
00262   if ( gother.getIsAGrid() != getIsAGrid())
00263     return false;
00264   const MESH& other = static_cast<const MESH&>( gother );
00265 
00266   int size1=getSpaceDimension()*getNumberOfNodes();
00267   int size2=other.getSpaceDimension()*other.getNumberOfNodes();
00268   if(size1!=size2)
00269     return false;
00270 
00271   const COORDINATE* CRD = other.getCoordinateptr();
00272   if( (!CRD && _coordinate) || (CRD && !(_coordinate)) ) {
00273     return false;
00274   }
00275 
00276   bool ret=true;
00277   if( _coordinate ) {
00278     const double* coord1=getCoordinates(MED_FULL_INTERLACE);
00279     const double* coord2=other.getCoordinates(MED_FULL_INTERLACE);
00280     for(int i=0;i<size1 && ret;i++) {
00281       ret=(fabs(coord1[i]-coord2[i])<1e-15);
00282     }
00283   }
00284   if(ret) {
00285     const CONNECTIVITY* CNV = other.getConnectivityptr();
00286     if( (!CNV && _connectivity) || (CNV && !(_connectivity)) ) {
00287       return false;
00288     }
00289     if(_connectivity) {
00290       return _connectivity->deepCompare(*other._connectivity);
00291     }
00292   }
00293   return ret;
00294 }
00295 
00299 ostream & ::MEDMEM::operator<<(ostream &os, const MESH &myMesh)
00300 {
00301   myMesh.printMySelf(os);
00302   return os;
00303 }
00304 void MESH::printMySelf(ostream &os) const
00305 {
00306   const MESH &myMesh = *this;
00307   int spacedimension = myMesh.getSpaceDimension();
00308   int meshdimension  = myMesh.getMeshDimension();
00309   int numberofnodes  = myMesh.getNumberOfNodes();
00310 
00311   if ( spacedimension == MED_INVALID ) {
00312     os << " Empty mesh ...";
00313     return;
00314   }
00315 
00316   os << "Space Dimension : " << spacedimension << endl << endl;
00317 
00318   os << "Mesh Dimension : " << meshdimension << endl << endl;
00319 
00320   if(myMesh.getCoordinateptr()) {
00321     const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
00322 
00323     os << "SHOW NODES COORDINATES : " << endl;
00324     os << "Name :" << endl;
00325     const string * coordinatesnames = myMesh.getCoordinatesNames();
00326     if ( coordinatesnames ) {
00327       for(int i=0; i<spacedimension ; i++)
00328       {
00329         os << " - " << coordinatesnames[i] << endl;
00330       }
00331     }
00332     os << "Unit :" << endl;
00333     const string * coordinatesunits = myMesh.getCoordinatesUnits();
00334     if ( coordinatesunits ) {
00335       for(int i=0; i<spacedimension ; i++)
00336       {
00337         os << " - " << coordinatesunits[i] << endl;
00338       }
00339     }
00340     for(int i=0; i<numberofnodes ; i++)
00341     {
00342       os << "Nodes " << i+1 << " : ";
00343       for (int j=0; j<spacedimension ; j++)
00344         os << coordinates[i*spacedimension+j] << " ";
00345       os << endl;
00346     }
00347   }
00348 
00349   if(myMesh.getConnectivityptr()) {
00350     os << endl << "SHOW CONNECTIVITY  :" << endl;
00351     os << *myMesh._connectivity << endl;
00352   }
00353 
00354   medEntityMesh entity;
00355   os << endl << "SHOW FAMILIES :" << endl << endl;
00356   for (int k=1; k<=4; k++)
00357   {
00358     if (k==1) entity = MED_NODE;
00359     if (k==2) entity = MED_CELL;
00360     if (k==3) entity = MED_FACE;
00361     if (k==4) entity = MED_EDGE;
00362     int numberoffamilies = myMesh.getNumberOfFamilies(entity);
00363     os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
00364     using namespace MED_EN;
00365     for (int i=1; i<numberoffamilies+1;i++)
00366     {
00367       os << * myMesh.getFamily(entity,i) << endl;
00368     }
00369   }
00370 
00371   os << endl << "SHOW GROUPS :" << endl << endl;
00372   for (int k=1; k<=4; k++)
00373   {
00374     if (k==1) entity = MED_NODE;
00375     if (k==2) entity = MED_CELL;
00376     if (k==3) entity = MED_FACE;
00377     if (k==4) entity = MED_EDGE;
00378     int numberofgroups = myMesh.getNumberOfGroups(entity);
00379     os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
00380     using namespace MED_EN;
00381     for (int i=1; i<numberofgroups+1;i++)
00382     {
00383       os << * myMesh.getGroup(entity,i) << endl;
00384     }
00385   }
00386 }
00387 
00394 int MESH::getMeshDimension() const
00395 {
00396   int dim = -1;
00397   if ( _connectivity )
00398     for ( int i = 0; i < _connectivity->getNumberOfTypes(MED_EN::MED_CELL); ++i )
00399       if ( _connectivity->getCellsTypes(MED_EN::MED_CELL)[i].getDimension() > dim )
00400         dim = _connectivity->getCellsTypes(MED_EN::MED_CELL)[i].getDimension();
00401   return dim;
00402 }
00412 int MESH::getElementNumber(MED_EN::medConnectivity ConnectivityType,
00413                            MED_EN::medEntityMesh Entity,
00414                            MED_EN::medGeometryElement Type,
00415                            int * connectivity) const
00416 {
00417   const char* LOC="MESH::getElementNumber " ;
00418   BEGIN_OF_MED(LOC) ;
00419 
00420   int numberOfValue ; // size of connectivity array
00421   CELLMODEL myType(Type) ;
00422   if (ConnectivityType==MED_DESCENDING)
00423     numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
00424   else
00425     numberOfValue = myType.getNumberOfNodes() ; // nodes
00426 
00427   const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
00428   const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
00429 
00430   // First node or face/edge
00431   int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
00432   int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
00433 
00434   // list to put cells numbers
00435   list<int> cellsList ;
00436   list<int>::iterator itList ;
00437   for (int i=indexBegin; i<indexEnd; i++)
00438     cellsList.push_back(myReverseConnectivityValue[i-1]) ;
00439 
00440   // Foreach node or face/edge in cell, we search cells which are in common.
00441   // At the end, we must have only one !
00442 
00443   for (int i=1; i<numberOfValue; i++) {
00444     int connectivity_i = connectivity[i] ;
00445     for (itList=cellsList.begin();itList!=cellsList.end();/*itList++*/) {
00446       bool find = false ;
00447       for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
00448         if ((*itList)==myReverseConnectivityValue[j-1]) {
00449           find=true ;
00450           break ;
00451         }
00452       }
00453       if (!find)
00454         itList=cellsList.erase(itList++);
00455       else
00456         itList++;
00457     }
00458   }
00459 
00460   if (cellsList.size()>1) // we have more than one cell
00461     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
00462 
00463   if (cellsList.size()==0)
00464     return -1;
00465 
00466   END_OF_MED(LOC);
00467 
00468   return cellsList.front() ;
00469 }
00470 
00493 SUPPORT * MESH::getBoundaryElements(MED_EN::medEntityMesh Entity) const throw (MEDEXCEPTION)
00494 {
00495   const char * LOC = "MESH::getBoundaryElements : " ;
00496   BEGIN_OF_MED(LOC) ;
00497   // some test :
00498   // actually we could only get face (in 3D) and edge (in 2D)
00499   medEntityMesh entityToParse=Entity;
00500   if(_spaceDimension == 3)
00501     if (Entity != MED_FACE)
00502       {
00503         if(Entity==MED_NODE)
00504           entityToParse=MED_FACE;
00505         else
00506           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
00507       }
00508   if(_spaceDimension == 2)
00509     if(Entity != MED_EDGE)
00510       {
00511         if(Entity==MED_NODE)
00512           entityToParse=MED_EDGE;
00513         else
00514           throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
00515       }
00516 
00517   // assure that descending connectivity is full
00518   if ( !_connectivity )
00519     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "no connectivity defined in MESH !"));
00520   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
00521 
00522   const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
00523   const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
00524   int numberOf = getNumberOfElements(entityToParse,MED_ALL_ELEMENTS) ;
00525   list<int> myElementsList;
00526 
00527   for (int i=0 ; i<numberOf; i++)
00528     if (myConnectivityValue[myConnectivityIndex[i]] == 0)
00529       myElementsList.push_back(i+1);
00530 
00531   if ( myElementsList.empty() && numberOf != 0 )
00532     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No boundary elements found by reverse descending connectivity for entity "<<Entity<<" !"));
00533 
00534   if(Entity==MED_NODE)
00535     return buildSupportOnNodeFromElementList(myElementsList,entityToParse);
00536   else
00537     return buildSupportOnElementsFromElementList(myElementsList,entityToParse);
00538 }
00546 void MESH::fillSupportOnNodeFromElementList(const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION)
00547 {
00548   MED_EN::medEntityMesh entity=supportToFill->getEntity();
00549   supportToFill->setAll(false);
00550   supportToFill->setMesh((MESH *)this);
00551 
00552   set<int> nodes;
00553   if ( entity == MED_NODE )
00554     {
00555       supportToFill->fillFromNodeList(listOfElt);
00556     }
00557   else
00558     {
00559       int lgth;
00560       for(list<int>::const_iterator iter=listOfElt.begin();iter!=listOfElt.end();iter++)
00561         {
00562           const int *conn=_connectivity->getConnectivityOfAnElement(MED_NODAL,entity,*iter,lgth);
00563           nodes.insert( conn, conn+lgth );
00564         }
00565       list<int> nodesList( nodes.begin(), nodes.end() );
00566       supportToFill->fillFromNodeList( nodesList );
00567     }
00568 }
00569 
00574 SUPPORT *MESH::buildSupportOnNodeFromElementList(const list<int>& listOfElt,MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
00575 {
00576   SUPPORT * mySupport = new SUPPORT;
00577   mySupport->setMesh((MESH *)this);
00578   mySupport->setName("Boundary");
00579   mySupport->setEntity( entity );
00580   fillSupportOnNodeFromElementList(listOfElt,mySupport);
00581   return mySupport;
00582 }
00583 
00591 FIELD<double, FullInterlace>* MESH::getVolume(const SUPPORT *Support, bool isAbs) const throw (MEDEXCEPTION)
00592 {
00593   const char * LOC = "MESH::getVolume(SUPPORT*) : ";
00594   BEGIN_OF_MED(LOC);
00595 
00596   // Support must be on 3D elements
00597 
00598   // Make sure that the MESH class is the same as the MESH class attribut
00599   // in the class Support
00600   const GMESH* myMesh = Support->getMesh();
00601   if (this != myMesh)
00602     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
00603 
00604   int dim_space = getSpaceDimension();
00605   medEntityMesh support_entity = Support->getEntity();
00606   bool onAll = Support->isOnAllElements();
00607 
00608   int nb_type, length_values;
00609   const medGeometryElement* types;
00610   int nb_entity_type;
00611   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
00612   const int* global_connectivity;
00613   nb_type = Support->getNumberOfTypes();
00614   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
00615   types = Support->getTypes();
00616   int index;
00617   FIELD<double, FullInterlace>* Volume =
00618     new FIELD<double, FullInterlace>(Support,1);
00619   Volume->setName("VOLUME");
00620   Volume->setDescription("cells volume");
00621   Volume->setComponentName(1,"volume");
00622   Volume->setComponentDescription(1,"desc-comp");
00623 
00624   string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
00625 
00626   Volume->setMEDComponentUnit(1,MEDComponentUnit);
00627 
00628   Volume->setIterationNumber(0);
00629   Volume->setOrderNumber(0);
00630   Volume->setTime(0.0);
00631 
00632   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
00633   ArrayNoGauss  *volume = Volume->getArrayNoGauss();
00634 
00635   index = 1;
00636   const double * coord = getCoordinates(MED_FULL_INTERLACE);
00637 
00638   for (int i=0;i<nb_type;i++)
00639   {
00640     medGeometryElement type = types[i] ;
00641     double xvolume;
00642     nb_entity_type = Support->getNumberOfElements(type);
00643     if(type != MED_EN::MED_POLYHEDRA)
00644     {
00645       if (onAll)
00646       {
00647         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
00648       }
00649       else
00650       {
00651         const int * supp_number = Support->getNumber(type);
00652         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
00653         const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
00654         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
00655 
00656         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
00657           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
00658             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
00659           }
00660         }
00661         global_connectivity = global_connectivity_tmp ;
00662       }
00663     }
00664 
00665     switch (type)
00666     {
00667     case MED_TETRA4 : case MED_TETRA10 :
00668       {
00669         for (int tetra=0;tetra<nb_entity_type;tetra++)
00670         {
00671           int tetra_index = (type%100)*tetra;
00672           int N1 = global_connectivity[tetra_index]-1;
00673           int N2 = global_connectivity[tetra_index+1]-1;
00674           int N3 = global_connectivity[tetra_index+2]-1;
00675           int N4 = global_connectivity[tetra_index+3]-1;
00676           xvolume=INTERP_KERNEL::calculateVolumeForTetra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4);
00677           if(isAbs)
00678             xvolume=fabs(xvolume);
00679           volume->setIJ(index,1,xvolume) ;
00680           index++;
00681         }
00682         break;
00683       }
00684     case MED_PYRA5 : case MED_PYRA13 :
00685       {
00686         for (int pyra=0;pyra<nb_entity_type;pyra++)
00687         {
00688           int pyra_index = (type%100)*pyra;
00689           int N1 = global_connectivity[pyra_index]-1;
00690           int N2 = global_connectivity[pyra_index+1]-1;
00691           int N3 = global_connectivity[pyra_index+2]-1;
00692           int N4 = global_connectivity[pyra_index+3]-1;
00693           int N5 = global_connectivity[pyra_index+4]-1;
00694           xvolume=INTERP_KERNEL::calculateVolumeForPyra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5);
00695           if(isAbs)
00696             xvolume=fabs(xvolume);
00697           volume->setIJ(index,1,xvolume) ;
00698           index++;
00699         }
00700         break;
00701       }
00702     case MED_PENTA6 : case MED_PENTA15 :
00703       {
00704         for (int penta=0;penta<nb_entity_type;penta++)
00705         {
00706           int penta_index = (type%100)*penta;
00707           int N1 = global_connectivity[penta_index]-1;
00708           int N2 = global_connectivity[penta_index+1]-1;
00709           int N3 = global_connectivity[penta_index+2]-1;
00710           int N4 = global_connectivity[penta_index+3]-1;
00711           int N5 = global_connectivity[penta_index+4]-1;
00712           int N6 = global_connectivity[penta_index+5]-1;
00713           xvolume=INTERP_KERNEL::calculateVolumeForPenta(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5,coord+dim_space*N6);
00714           if(isAbs)
00715             xvolume=fabs(xvolume);
00716           volume->setIJ(index,1,xvolume) ;
00717           index++;
00718         }
00719         break;
00720       }
00721     case MED_HEXA8 : case MED_HEXA20 :
00722       {
00723         for (int hexa=0;hexa<nb_entity_type;hexa++)
00724         {
00725           int hexa_index = (type%100)*hexa;
00726 
00727           int N1 = global_connectivity[hexa_index]-1;
00728           int N2 = global_connectivity[hexa_index+1]-1;
00729           int N3 = global_connectivity[hexa_index+2]-1;
00730           int N4 = global_connectivity[hexa_index+3]-1;
00731           int N5 = global_connectivity[hexa_index+4]-1;
00732           int N6 = global_connectivity[hexa_index+5]-1;
00733           int N7 = global_connectivity[hexa_index+6]-1;
00734           int N8 = global_connectivity[hexa_index+7]-1;
00735           xvolume=INTERP_KERNEL::calculateVolumeForHexa(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5,coord+dim_space*N6,coord+dim_space*N7,coord+dim_space*N8);
00736           if(isAbs)
00737             xvolume=fabs(xvolume);
00738           volume->setIJ(index,1,xvolume) ;
00739           index++;
00740         }
00741         break;
00742       }
00743     case MED_POLYHEDRA:
00744       {
00745         double bary[3];
00746         if(onAll)
00747         {
00748           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
00749           {
00750             int lgthNodes,iPts,iFaces,iPtsInFace;
00751             int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
00752             int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
00753             int nbOfFaces,*nbOfNodesPerFaces;
00754             int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(offsetWithClassicType+polyhs+1,nbOfFaces,nbOfNodesPerFaces);
00755             double **pts=new double * [lgthNodes];
00756             double ***pts1=new double ** [nbOfFaces];
00757             for(iPts=0;iPts<lgthNodes;iPts++)
00758               pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
00759             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
00760             {
00761               pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
00762               for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
00763                 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
00764             }
00765             delete [] nodes1;
00766             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
00767             delete [] nodes;
00768             delete [] pts;
00769             if(isAbs)
00770               xvolume=INTERP_KERNEL::calculateVolumeForPolyhAbs((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
00771             else
00772               xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
00773             delete [] nbOfNodesPerFaces;
00774             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
00775               delete [] pts1[iFaces];
00776             delete [] pts1;
00777             volume->setIJ(index,1,xvolume) ;
00778             index++;
00779           }
00780         }
00781         else
00782         {
00783           const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
00784           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
00785           {
00786             int lgthNodes,iPts,iFaces,iPtsInFace;
00787             int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
00788             int nbOfFaces,*nbOfNodesPerFaces;
00789             int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(supp_number[polyhs],nbOfFaces,nbOfNodesPerFaces);
00790             double **pts=new double * [lgthNodes];
00791             double ***pts1=new double ** [nbOfFaces];
00792             for(iPts=0;iPts<lgthNodes;iPts++)
00793               pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
00794             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
00795             {
00796               pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
00797               for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
00798                 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
00799             }
00800             delete [] nodes1;
00801             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
00802             delete [] nodes;
00803             delete [] pts;
00804             if(isAbs)
00805               xvolume=INTERP_KERNEL::calculateVolumeForPolyhAbs((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
00806             else
00807               xvolume=INTERP_KERNEL::calculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
00808             delete [] nbOfNodesPerFaces;
00809             for(iFaces=0;iFaces<nbOfFaces;iFaces++)
00810               delete [] pts1[iFaces];
00811             delete [] pts1;
00812             volume->setIJ(index,1,xvolume) ;
00813             index++;
00814           }
00815         }
00816         break;
00817       }
00818     default :
00819       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
00820       break;
00821     }
00822 
00823     if (!onAll && type!=MED_EN::MED_POLYHEDRA)
00824       delete [] global_connectivity ;
00825   }
00826 
00827   return Volume;
00828 }
00834 FIELD<double, FullInterlace>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
00835 {
00836   const char * LOC = "MESH::getArea(SUPPORT*) : ";
00837   BEGIN_OF_MED(LOC);
00838 
00839   // Support must be on 2D elements
00840 
00841   // Make sure that the MESH class is the same as the MESH class attribut
00842   // in the class Support
00843   const GMESH* myMesh = Support->getMesh();
00844   if (this != myMesh)
00845     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
00846 
00847   int dim_space = getSpaceDimension();
00848   medEntityMesh support_entity = Support->getEntity();
00849   bool onAll = Support->isOnAllElements();
00850 
00851   int nb_type, length_values;
00852   const medGeometryElement* types;
00853   int nb_entity_type;
00854   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
00855   const int* global_connectivity;
00856 
00857   nb_type = Support->getNumberOfTypes();
00858   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
00859   types = Support->getTypes();
00860 
00861   int index;
00862   FIELD<double, FullInterlace>* Area;
00863 
00864   Area = new FIELD<double, FullInterlace>(Support,1);
00865   Area->setName("AREA");
00866   Area->setDescription("cells or faces area");
00867 
00868   Area->setComponentName(1,"area");
00869   Area->setComponentDescription(1,"desc-comp");
00870 
00871   /*  string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
00872 
00873   string MEDComponentUnit="trucmuch";
00874 
00875   Area->setMEDComponentUnit(1,MEDComponentUnit);
00876 
00877   Area->setIterationNumber(0);
00878   Area->setOrderNumber(0);
00879   Area->setTime(0.0);
00880 
00881   double *area = (double *)Area->getValue();
00882 
00883   const double * coord = getCoordinates(MED_FULL_INTERLACE);
00884   index = 0;
00885 
00886   for (int i=0;i<nb_type;i++)
00887   {
00888     medGeometryElement type = types[i] ;
00889     nb_entity_type = Support->getNumberOfElements(type);
00890     const int *global_connectivityIndex = 0;
00891     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
00892     {
00893       global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
00894       if (onAll)
00895       {
00896         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
00897       }
00898       else
00899       {
00900         const int * supp_number = Support->getNumber(type);
00901         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
00902         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
00903 
00904         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
00905           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
00906             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
00907           }
00908         }
00909         global_connectivity = global_connectivity_tmp ;
00910       }
00911     }
00912     switch (type)
00913     {
00914     case MED_TRIA3 : case MED_TRIA6 :
00915       {
00916         for (int tria=0;tria<nb_entity_type;tria++)
00917         {
00918           int tria_index = (type%100)*tria;
00919 
00920           int N1 = global_connectivity[tria_index]-1;
00921           int N2 = global_connectivity[tria_index+1]-1;
00922           int N3 = global_connectivity[tria_index+2]-1;
00923 
00924           area[index]=INTERP_KERNEL::calculateAreaForTria(coord+(dim_space*N1),
00925                                                           coord+(dim_space*N2),
00926                                                           coord+(dim_space*N3),dim_space);
00927           index++;
00928         }
00929         break;
00930       }
00931     case MED_QUAD4 : case MED_QUAD8 :
00932       {
00933         for (int quad=0;quad<nb_entity_type;quad++)
00934         {
00935           int quad_index = (type%100)*quad;
00936 
00937           int N1 = global_connectivity[quad_index]-1;
00938           int N2 = global_connectivity[quad_index+1]-1;
00939           int N3 = global_connectivity[quad_index+2]-1;
00940           int N4 = global_connectivity[quad_index+3]-1;
00941 
00942           area[index]=INTERP_KERNEL::calculateAreaForQuad(coord+dim_space*N1,
00943                                                           coord+dim_space*N2,
00944                                                           coord+dim_space*N3,
00945                                                           coord+dim_space*N4,dim_space);
00946           index++;
00947         }
00948         break;
00949       }
00950     case MED_POLYGON :
00951       {
00952         if(onAll)
00953         {
00954           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
00955           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
00956           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
00957           for (int polygs=0;polygs<nb_entity_type;polygs++)
00958           {
00959             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
00960             double **pts=new double * [size];
00961             for(int iPts=0;iPts<size;iPts++)
00962               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1));
00963             area[index] = INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,size,dim_space);
00964             delete [] pts;
00965             index++;
00966           }
00967         }
00968         else
00969         {
00970           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
00971           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
00972           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
00973           for (int polygs=0;polygs<nb_entity_type;polygs++)
00974           {
00975             int size=connectivity_index[supp_number[polygs]]-connectivity_index[supp_number[polygs]-1];
00976             double **pts=new double * [size];
00977             for(int iPts=0;iPts<size;iPts++)
00978               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-1]+iPts-1]-1));
00979             area[index]=INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,size,dim_space);
00980             delete [] pts;
00981             index++;
00982           }
00983         }
00984         break;
00985       }
00986     default :
00987       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
00988       break;
00989     }
00990 
00991     if (!onAll)
00992       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
00993         delete [] global_connectivity ;
00994   }
00995   return Area;
00996 }
01000 FIELD<double, FullInterlace>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
01001 {
01002   const char * LOC = "MESH::getLength(SUPPORT*) : ";
01003   BEGIN_OF_MED(LOC);
01004 
01005   // Support must be on 1D elements
01006 
01007   // Make sure that the MESH class is the same as the MESH class attribut
01008   // in the class Support
01009   const GMESH* myMesh = Support->getMesh();
01010   if (this != myMesh)
01011     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
01012 
01013   int dim_space = getSpaceDimension();
01014   medEntityMesh support_entity = Support->getEntity();
01015   bool onAll = Support->isOnAllElements();
01016 
01017   int nb_type, length_values;
01018   const medGeometryElement* types;
01019   int nb_entity_type;
01020   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
01021   const int* global_connectivity;
01022 
01023   nb_type = Support->getNumberOfTypes();
01024   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
01025   types = Support->getTypes();
01026 
01027   int index;
01028   FIELD<double, FullInterlace>* Length;
01029 
01030   Length = new FIELD<double, FullInterlace>(Support,1);
01031 
01032   typedef  MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
01033   ArrayNoGauss * length = Length->getArrayNoGauss();
01034 
01035   const double * coord = getCoordinates(MED_FULL_INTERLACE);
01036   index = 1;
01037 
01038   for (int i=0;i<nb_type;i++)
01039   {
01040     medGeometryElement type = types[i] ;
01041     double xlength;
01042 
01043     if (onAll)
01044     {
01045       nb_entity_type = getNumberOfElements(support_entity,type);
01046       global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
01047     }
01048     else
01049     {
01050       nb_entity_type = Support->getNumberOfElements(type);
01051       const int * supp_number = Support->getNumber(type);
01052       const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01053       const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
01054       int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
01055 
01056       for (int k_type = 0; k_type<nb_entity_type; k_type++) {
01057         for (int j_ent = 0; j_ent<(type%100); j_ent++) {
01058           global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
01059         }
01060       }
01061       global_connectivity = global_connectivity_tmp ;
01062     }
01063 
01064     switch (type)
01065     {
01066     case MED_SEG2 : case MED_SEG3 :
01067       {
01068         for (int edge=0;edge<nb_entity_type;edge++)
01069         {
01070           int edge_index = (type%100)*edge;
01071 
01072           int N1 = global_connectivity[edge_index]-1;
01073           int N2 = global_connectivity[edge_index+1]-1;
01074 
01075           double x1 = coord[dim_space*N1];
01076           double x2 = coord[dim_space*N2];
01077 
01078           double y1 = coord[dim_space*N1+1];
01079           double y2 = coord[dim_space*N2+1];
01080 
01081           double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
01082             z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
01083 
01084           xlength =  sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
01085                           (z1 - z2)*(z1 - z2));
01086 
01087           length->setIJ(index,1,xlength) ;
01088           index++;
01089         }
01090         break;
01091       }
01092     default :
01093       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
01094       break;
01095     }
01096 
01097     if (!onAll) delete [] global_connectivity ;
01098   }
01099 
01100   return Length;
01101 }
01102 
01110 FIELD<double, FullInterlace>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
01111 {
01112   const char * LOC = "MESH::getNormal(SUPPORT*) : ";
01113   BEGIN_OF_MED(LOC);
01114 
01115   // Support must be on 2D or 1D elements
01116 
01117   // Make sure that the MESH class is the same as the MESH class attribut
01118   // in the class Support
01119   const GMESH* myMesh = Support->getMesh();
01120   if (this != myMesh)
01121     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
01122 
01123   int dim_space = getSpaceDimension();
01124   int mesh_dim=getMeshDimension();
01125   medEntityMesh support_entity = Support->getEntity();
01126   bool onAll = Support->isOnAllElements();
01127 
01128   if( support_entity!=MED_EDGE && (mesh_dim!=1 || support_entity!=MED_CELL) && ( mesh_dim!=2 || support_entity!=MED_CELL ) && ( mesh_dim!=3 || support_entity!=MED_FACE ))
01129     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"incompatible mesh dimension and entity"));
01130   int nb_type, length_values;
01131   const medGeometryElement* types;
01132   int nb_entity_type;
01133   // !!!! WARNING : use of nodal global numbering in the mesh !!!!
01134   const int* global_connectivity;
01135 
01136   nb_type = Support->getNumberOfTypes();
01137   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
01138   types = Support->getTypes();
01139 
01140   int index;
01141 
01142   FIELD<double, FullInterlace>* Normal =
01143     new FIELD<double, FullInterlace>(Support,dim_space);
01144   Normal->setName("NORMAL");
01145   Normal->setDescription("cells or faces normal");
01146   for (int k=1;k<=dim_space;k++) {
01147     Normal->setComponentName(k,"normal");
01148     Normal->setComponentDescription(k,"desc-comp");
01149     Normal->setMEDComponentUnit(k,"unit");
01150   }
01151 
01152   Normal->setIterationNumber(MED_NOPDT);
01153   Normal->setOrderNumber(MED_NONOR);
01154   Normal->setTime(0.0);
01155   double *normal = (double *)Normal->getValue();
01156 
01157   const double * coord = getCoordinates(MED_FULL_INTERLACE);
01158   index = 0;
01159 
01160   for (int i=0;i<nb_type;i++)
01161   {
01162     medGeometryElement type = types[i] ;
01163     nb_entity_type = Support->getNumberOfElements(type);
01164 
01165     // Make sure we trying to get Normals on
01166     // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
01167     // or on SEG2 or SEG3 (1D) cells on a 2D mesh
01168 
01169     if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
01170            (type==MED_QUAD4) || (type==MED_QUAD8) || (type==MED_POLYGON)) &&
01171           (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
01172                                 (dim_space != 2)) )
01173       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
01174     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
01175     {
01176       if (onAll)
01177       {
01178         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
01179       }
01180       else
01181       {
01182         const int * supp_number = Support->getNumber(type);
01183         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01184         const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
01185         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
01186 
01187         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
01188           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
01189             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
01190           }
01191         }
01192 
01193         global_connectivity = global_connectivity_tmp ;
01194       }
01195     }
01196 
01197     switch (type)
01198     {
01199     case MED_TRIA3 : case MED_TRIA6 :
01200       {
01201         for (int tria=0;tria<nb_entity_type;tria++)
01202         {
01203           int tria_index = (type%100)*tria;
01204           int N1 = global_connectivity[tria_index]-1;
01205           int N2 = global_connectivity[tria_index+1]-1;
01206           int N3 = global_connectivity[tria_index+2]-1;
01207           INTERP_KERNEL::calculateNormalForTria(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,normal+3*index);
01208           index++;
01209         }
01210         break;
01211       }
01212     case MED_QUAD4 : case MED_QUAD8 :
01213       {
01214         for (int quad=0;quad<nb_entity_type;quad++)
01215         {
01216           int quad_index = (type%100)*quad;
01217           int N1 = global_connectivity[quad_index]-1;
01218           int N2 = global_connectivity[quad_index+1]-1;
01219           int N3 = global_connectivity[quad_index+2]-1;
01220           int N4 = global_connectivity[quad_index+3]-1;
01221           INTERP_KERNEL::calculateNormalForQuad(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,normal+3*index);
01222           index++;
01223         }
01224         break;
01225       }
01226     case MED_SEG2 : case MED_SEG3 :
01227       {
01228         double xnormal1, xnormal2;
01229         for (int edge=0;edge<nb_entity_type;edge++)
01230         {
01231           int edge_index = (type%100)*edge;
01232 
01233           int N1 = global_connectivity[edge_index]-1;
01234           int N2 = global_connectivity[edge_index+1]-1;
01235 
01236           double x1 = coord[dim_space*N1];
01237           double x2 = coord[dim_space*N2];
01238 
01239           double y1 = coord[dim_space*N1+1];
01240           double y2 = coord[dim_space*N2+1];
01241 
01242           xnormal1 = -(y2-y1);
01243           xnormal2 = x2-x1;
01244 
01245           normal[2*index] = xnormal1 ;
01246           normal[2*index+1] = xnormal2 ;
01247 
01248           index++;
01249         }
01250         break;
01251       }
01252     case MED_POLYGON :
01253       {
01254         if(onAll)
01255         {
01256           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
01257           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01258           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
01259           for (int polygs=0;polygs<nb_entity_type;polygs++)
01260           {
01261             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
01262             double **pts=new double * [size];
01263             for(int iPts=0;iPts<size;iPts++)
01264               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1])-1);
01265             INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
01266             delete [] pts;
01267             index++;
01268           }
01269         }
01270         else
01271         {
01272           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
01273           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01274           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
01275           for (int polygs=0;polygs<nb_entity_type;polygs++)
01276           {
01277             int size=connectivity_index[supp_number[polygs]]-connectivity_index[supp_number[polygs]-1];
01278             double **pts=new double * [size];
01279             for(int iPts=0;iPts<size;iPts++)
01280               pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-1]+iPts-1])-1);
01281             INTERP_KERNEL::calculateNormalForPolyg((const double **)pts,size,normal+3*index);
01282             delete [] pts;
01283             index++;
01284           }
01285         }
01286         break;
01287       }
01288     default :
01289       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
01290       break;
01291     }
01292     if (!onAll && type!=MED_EN::MED_POLYGON)
01293       delete [] global_connectivity ;
01294   }
01295   END_OF_MED(LOC);
01296 
01297   return Normal;
01298 }
01303 FIELD<double, FullInterlace>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
01304 {
01305   const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
01306   const GMESH* myMesh = Support->getMesh();
01307   if (this != myMesh)
01308     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
01309 
01310   int dim_space = getSpaceDimension();
01311   medEntityMesh support_entity = Support->getEntity();
01312   bool onAll = Support->isOnAllElements();
01313 
01314   int nb_type, length_values;
01315   const medGeometryElement* types;
01316   int nb_entity_type;
01317   const int* global_connectivity;
01318 
01319   nb_type = Support->getNumberOfTypes();
01320   length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
01321   types = Support->getTypes();
01322 
01323   FIELD<double, FullInterlace>* Barycenter;
01324   Barycenter = new FIELD<double, FullInterlace>(Support,dim_space);
01325   Barycenter->setName("BARYCENTER");
01326   Barycenter->setDescription("cells or faces barycenter");
01327 
01328   for (int k=0;k<dim_space;k++) {
01329     int kp1 = k + 1;
01330     Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
01331     Barycenter->setComponentDescription(kp1,"desc-comp");
01332     Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
01333   }
01334   Barycenter->setIterationNumber(0);
01335   Barycenter->setOrderNumber(0);
01336   Barycenter->setTime(0.0);
01337   double *barycenter=(double *)Barycenter->getValue();
01338   const double * coord = getCoordinates(MED_FULL_INTERLACE);
01339   int index=0;
01340   for (int i=0;i<nb_type;i++)
01341   {
01342     medGeometryElement type = types[i] ;
01343     nb_entity_type = Support->getNumberOfElements(type);
01344     if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA )
01345     {
01346       if (onAll)
01347       {
01348         global_connectivity = getConnectivity(MED_NODAL,support_entity,type);
01349       }
01350       else
01351       {
01352         const int * supp_number = Support->getNumber(type);
01353         const int * connectivity = getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01354         int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
01355 
01356         for (int k_type = 0; k_type<nb_entity_type; k_type++) {
01357           for (int j_ent = 0; j_ent<(type%100); j_ent++) {
01358             const int *global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
01359             global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
01360           }
01361         }
01362         global_connectivity = global_connectivity_tmp;
01363       }
01364     }
01365 
01366     switch (type)
01367     {
01368     case MED_TETRA4 : case MED_TETRA10 :
01369       {
01370         for (int tetra =0;tetra<nb_entity_type;tetra++)
01371         {
01372           int tetra_index = (type%100)*tetra;
01373 
01374           int N1 = global_connectivity[tetra_index]-1;
01375           int N2 = global_connectivity[tetra_index+1]-1;
01376           int N3 = global_connectivity[tetra_index+2]-1;
01377           int N4 = global_connectivity[tetra_index+3]-1;
01378           double *pts[4];
01379           pts[0]=(double *)coord+dim_space*N1;
01380           pts[1]=(double *)coord+dim_space*N2;
01381           pts[2]=(double *)coord+dim_space*N3;
01382           pts[3]=(double *)coord+dim_space*N4;
01383           INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
01384           index++;
01385         }
01386         break;
01387       }
01388     case MED_PYRA5 : case MED_PYRA13 :
01389       {
01390         for (int pyra=0;pyra<nb_entity_type;pyra++)
01391         {
01392           int pyra_index = (type%100)*pyra;
01393 
01394           int N1 = global_connectivity[pyra_index]-1;
01395           int N2 = global_connectivity[pyra_index+1]-1;
01396           int N3 = global_connectivity[pyra_index+2]-1;
01397           int N4 = global_connectivity[pyra_index+3]-1;
01398           int N5 = global_connectivity[pyra_index+4]-1;
01399           double *pts[5];
01400           pts[0]=(double *)coord+dim_space*N1;
01401           pts[1]=(double *)coord+dim_space*N2;
01402           pts[2]=(double *)coord+dim_space*N3;
01403           pts[3]=(double *)coord+dim_space*N4;
01404           pts[4]=(double *)coord+dim_space*N5;
01405           INTERP_KERNEL::calculateBarycenter<5,3>((const double **)pts,barycenter+3*index);
01406           index++;
01407         }
01408         break;
01409       }
01410     case MED_PENTA6 : case MED_PENTA15 :
01411       {
01412         for (int penta=0;penta<nb_entity_type;penta++)
01413         {
01414           int penta_index = (type%100)*penta;
01415 
01416           int N1 = global_connectivity[penta_index]-1;
01417           int N2 = global_connectivity[penta_index+1]-1;
01418           int N3 = global_connectivity[penta_index+2]-1;
01419           int N4 = global_connectivity[penta_index+3]-1;
01420           int N5 = global_connectivity[penta_index+4]-1;
01421           int N6 = global_connectivity[penta_index+5]-1;
01422           double *pts[6];
01423           pts[0]=(double *)coord+dim_space*N1;
01424           pts[1]=(double *)coord+dim_space*N2;
01425           pts[2]=(double *)coord+dim_space*N3;
01426           pts[3]=(double *)coord+dim_space*N4;
01427           pts[4]=(double *)coord+dim_space*N5;
01428           pts[5]=(double *)coord+dim_space*N6;
01429           INTERP_KERNEL::calculateBarycenter<6,3>((const double **)pts,barycenter+3*index);
01430           index++;
01431         }
01432         break;
01433       }
01434     case MED_HEXA8 : case MED_HEXA20 :
01435       {
01436         for (int hexa=0;hexa<nb_entity_type;hexa++)
01437         {
01438           int hexa_index = (type%100)*hexa;
01439 
01440           int N1 = global_connectivity[hexa_index]-1;
01441           int N2 = global_connectivity[hexa_index+1]-1;
01442           int N3 = global_connectivity[hexa_index+2]-1;
01443           int N4 = global_connectivity[hexa_index+3]-1;
01444           int N5 = global_connectivity[hexa_index+4]-1;
01445           int N6 = global_connectivity[hexa_index+5]-1;
01446           int N7 = global_connectivity[hexa_index+6]-1;
01447           int N8 = global_connectivity[hexa_index+7]-1;
01448           double *pts[8];
01449           pts[0]=(double *)coord+dim_space*N1;
01450           pts[1]=(double *)coord+dim_space*N2;
01451           pts[2]=(double *)coord+dim_space*N3;
01452           pts[3]=(double *)coord+dim_space*N4;
01453           pts[4]=(double *)coord+dim_space*N5;
01454           pts[5]=(double *)coord+dim_space*N6;
01455           pts[6]=(double *)coord+dim_space*N7;
01456           pts[7]=(double *)coord+dim_space*N8;
01457           INTERP_KERNEL::calculateBarycenter<8,3>((const double **)pts,barycenter+3*index);
01458           index++;
01459         }
01460         break;
01461       }
01462     case MED_TRIA3 : case MED_TRIA6 :
01463       {
01464         for (int tria=0;tria<nb_entity_type;tria++)
01465         {
01466           int tria_index = (type%100)*tria;
01467           int N1 = global_connectivity[tria_index]-1;
01468           int N2 = global_connectivity[tria_index+1]-1;
01469           int N3 = global_connectivity[tria_index+2]-1;
01470           double *pts[3];
01471           pts[0]=(double *)coord+dim_space*N1;
01472           pts[1]=(double *)coord+dim_space*N2;
01473           pts[2]=(double *)coord+dim_space*N3;
01474           if (dim_space==2)
01475             INTERP_KERNEL::calculateBarycenter<3,2>((const double **)pts,barycenter+2*index);
01476           else
01477             INTERP_KERNEL::calculateBarycenter<3,3>((const double **)pts,barycenter+3*index);
01478           index++;
01479         }
01480         break;
01481       }
01482     case MED_QUAD4 : case MED_QUAD8 :
01483       {
01484         for (int quad=0;quad<nb_entity_type;quad++)
01485         {
01486           int quad_index = (type%100)*quad;
01487           int N1 = global_connectivity[quad_index]-1;
01488           int N2 = global_connectivity[quad_index+1]-1;
01489           int N3 = global_connectivity[quad_index+2]-1;
01490           int N4 = global_connectivity[quad_index+3]-1;
01491           double *pts[4];
01492           pts[0]=(double *)coord+dim_space*N1;
01493           pts[1]=(double *)coord+dim_space*N2;
01494           pts[2]=(double *)coord+dim_space*N3;
01495           pts[3]=(double *)coord+dim_space*N4;
01496           if (dim_space==2)
01497             INTERP_KERNEL::calculateBarycenter<4,2>((const double **)pts,barycenter+2*index);
01498           else
01499             INTERP_KERNEL::calculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
01500           index++;
01501         }
01502         break;
01503       }
01504     case MED_SEG2 : case MED_SEG3 :
01505       {
01506         for (int edge=0;edge<nb_entity_type;edge++)
01507         {
01508           int edge_index = (type%100)*edge;
01509           int N1 = global_connectivity[edge_index]-1;
01510           int N2 = global_connectivity[edge_index+1]-1;
01511           double *pts[2];
01512           pts[0]=(double *)coord+dim_space*N1;
01513           pts[1]=(double *)coord+dim_space*N2;
01514           if (dim_space==2)
01515             INTERP_KERNEL::calculateBarycenter<2,2>((const double **)pts,barycenter+2*index);
01516           else
01517             INTERP_KERNEL::calculateBarycenter<2,3>((const double **)pts,barycenter+3*index);
01518           index++;
01519         }
01520         break;
01521       }
01522     case MED_POLYGON :
01523       {
01524         if(onAll)
01525         {
01526           int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
01527           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01528           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity) + offsetWithClassicType;
01529           for (int polygs=0;polygs<nb_entity_type;polygs++)
01530           {
01531             int size=connectivity_index[polygs+1]-connectivity_index[polygs];
01532             double **pts=new double * [size];
01533             for(int iPts=0;iPts<size;iPts++)
01534               pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1);
01535             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
01536             delete [] pts;
01537             index++;
01538           }
01539         }
01540         else
01541         {
01542           const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
01543           const int * connectivity = _connectivity->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
01544           const int * connectivity_index = _connectivity->getConnectivityIndex(MED_EN::MED_NODAL,support_entity);
01545           for (int polygs=0;polygs<nb_entity_type;polygs++)
01546           {
01547             int localPolygsNbP1=supp_number[polygs];
01548             int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
01549             double **pts=new double * [size];
01550             for(int iPts=0;iPts<size;iPts++)
01551               pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1]-1);
01552             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
01553             delete [] pts;
01554             index++;
01555           }
01556         }
01557         break;
01558       }
01559     case MED_EN::MED_POLYHEDRA:
01560       {
01561         if(onAll)
01562         {
01563           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
01564           {
01565             int lgthNodes;
01566             int offsetWithClassicType=getGlobalNumberingIndex(support_entity)[i]-1;
01567             int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
01568             double **pts=new double * [lgthNodes];
01569             for(int iPts=0;iPts<lgthNodes;iPts++)
01570               pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
01571             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
01572             delete [] pts;
01573             delete [] nodes;
01574             index++;
01575           }
01576         }
01577         else
01578         {
01579           const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
01580           for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
01581           {
01582             int lgthNodes;
01583             int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
01584             double **pts=new double * [lgthNodes];
01585             for(int iPts=0;iPts<lgthNodes;iPts++)
01586               pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
01587             INTERP_KERNEL::calculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
01588             delete [] pts;
01589             delete [] nodes;
01590             index++;
01591           }
01592         }
01593         break;
01594       }
01595     default :
01596       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
01597       break;
01598     }
01599 
01600     if (!onAll)
01601       if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
01602         delete [] global_connectivity;
01603   }
01604   return Barycenter;
01605 }
01610 bool MESH::isEmpty() const
01611 {
01612   bool notempty = _name != "NOT DEFINED"                || _coordinate != NULL           || _connectivity != NULL ||
01613     _spaceDimension !=MED_INVALID || 
01614     _numberOfNodes !=MED_INVALID  || _groupNode.size() != 0   ||
01615     _familyNode.size() != 0       || _groupCell.size() != 0   ||
01616     _familyCell.size() != 0       || _groupFace.size() != 0   ||
01617     _familyFace.size() != 0       || _groupEdge.size() != 0   ||
01618     _familyEdge.size() != 0;
01619   return !notempty;
01620 }
01621 
01622 //================================================================================
01626 //================================================================================
01627 
01628 bool MESH::getIsAGrid() const
01629 {
01630   return false;
01631 }
01632 
01636 const MESH* MESH::convertInMESH() const
01637 {
01638   this->addReference();
01639   return this;
01640 }
01641 
01656 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
01657 {
01658   const char * LOC = "MESH::getSkin : " ;
01659   BEGIN_OF_MED(LOC) ;
01660   // some test :
01661   if (this != Support3D->getMesh())
01662     throw MEDEXCEPTION(STRING(LOC) <<  "no compatibility between *this and SUPPORT::_mesh !");
01663   if (getMeshDimension() != 3 || Support3D->getEntity() != MED_CELL)
01664     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
01665 
01666   // well, all rigth !
01667   SUPPORT * mySupport = new SUPPORT;
01668   mySupport->setMesh((MESH *)this);
01669   mySupport->setName("Skin");
01670   mySupport->setEntity( MED_FACE );
01671 
01672   list<int> myElementsList;
01673   int i,j, size = 0;
01674 
01675   // assure that descending connectivity is full
01676   if ( !_connectivity )
01677     throw MEDEXCEPTION(STRING(LOC) << "no connectivity defined in MESH !");
01678   _connectivity->calculateFullDescendingConnectivity(MED_CELL);
01679   if (Support3D->isOnAllElements())
01680   {
01681     const int* value = getReverseConnectivity(MED_DESCENDING);
01682     const int* index = getReverseConnectivityIndex(MED_DESCENDING);
01683 
01684     int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
01685     for (int i = 0; i < nbFaces; i++)
01686     {
01687       //a face is in skin if it is linked to one element or if one of the elements
01688       //it is linked to is "virtual"
01689       if ((index[i+1]-index[i]==1) || (value[index[i]-1]==0) || (value[index[i]]==0)) {
01690         myElementsList.push_back( i+1 );
01691         size++;
01692       }
01693     }
01694   }
01695   else
01696   {
01697     map<int,int> FaceNbEncounterNb;
01698 
01699     int * myConnectivityValue = const_cast <int *>
01700       (getConnectivity(MED_DESCENDING,MED_CELL, MED_ALL_ELEMENTS));
01701     int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
01702     int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
01703     int nbCells = Support3D->getnumber()->getLength();
01704     for (i=0; i<nbCells; ++i)
01705     {
01706       int cellNb = myCellNbs [ i ];
01707       int faceFirst = myConnectivityIndex[ cellNb-1 ];
01708       int faceLast  = myConnectivityIndex[ cellNb ];
01709       for (j = faceFirst; j < faceLast; ++j)
01710       {
01711         int faceNb = abs( myConnectivityValue [ j-1 ] );
01712         if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
01713           FaceNbEncounterNb[ faceNb ] = 1;
01714         else
01715           FaceNbEncounterNb[ faceNb ] ++;
01716       }
01717     }
01718     map<int,int>::iterator FaceNbEncounterNbItr;
01719     for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
01720          FaceNbEncounterNbItr != FaceNbEncounterNb.end();
01721          FaceNbEncounterNbItr ++)
01722       if ((*FaceNbEncounterNbItr).second == 1)
01723       {
01724         myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
01725         size++ ;
01726       }
01727   }
01728   // Well, we must know how many geometric type we have found
01729   int * myListArray = new int[size] ;
01730   int id = 0 ;
01731   list<int>::iterator myElementsListIt ;
01732   for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
01733     myListArray[id]=(*myElementsListIt) ;
01734     id ++ ;
01735   }
01736 
01737   int numberOfGeometricType ;
01738   medGeometryElement* geometricType ;
01739   int * geometricTypeNumber ;
01740   int * numberOfEntities ;
01741   int * mySkyLineArrayIndex ;
01742 
01743   int numberOfType = getNumberOfTypes(MED_FACE) ;
01744   if (numberOfType == 1) // wonderfull : it's easy !
01745     {
01746       numberOfGeometricType = 1 ;
01747       geometricType = new medGeometryElement[1] ;
01748       const medGeometryElement *  allType = getTypes(MED_FACE);
01749       geometricType[0] = allType[0] ;
01750       geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
01751       geometricTypeNumber[0] = 0 ;
01752       numberOfEntities = new int[1] ;
01753       numberOfEntities[0] = size ;
01754       mySkyLineArrayIndex = new int[2] ;
01755       mySkyLineArrayIndex[0]=1 ;
01756       mySkyLineArrayIndex[1]=1+size ;
01757     }
01758   else // hemmm
01759     {
01760       map<medGeometryElement,int> theType ;
01761       for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++)
01762         {
01763           medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
01764           if (theType.find(myType) != theType.end() )
01765             theType[myType]+=1 ;
01766           else
01767             theType[myType]=1 ;
01768         }
01769       numberOfGeometricType = theType.size() ;
01770       geometricType = new medGeometryElement[numberOfGeometricType] ;
01771       geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
01772       numberOfEntities = new int[numberOfGeometricType] ;
01773       mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
01774       int index = 0 ;
01775       mySkyLineArrayIndex[0]=1 ;
01776       map<medGeometryElement,int>::iterator theTypeIt ;
01777       for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++)
01778         {
01779           geometricType[index] = (*theTypeIt).first ;
01780           geometricTypeNumber[index] = 0 ;
01781           numberOfEntities[index] = (*theTypeIt).second ;
01782           mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
01783           index++ ;
01784         }
01785     }
01786   MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
01787 
01788   mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
01789   mySupport->setGeometricType(geometricType) ;
01790   mySupport->setNumberOfElements(numberOfEntities) ;
01791   mySupport->setNumber(mySkyLineArray) ;
01792 
01793   delete[] numberOfEntities;
01794   delete[] geometricTypeNumber;
01795   delete[] geometricType;
01796   delete[] mySkyLineArrayIndex;
01797   delete[] myListArray;
01798 
01799   END_OF_MED(LOC);
01800   return mySupport ;
01801 
01802 }
01803 
01804 
01805 
01806 int MESH::getElementContainingPoint(const double *coord)
01807 {
01808   if(_spaceDimension!=3 && _spaceDimension!=2)
01809     throw MEDEXCEPTION("MESH::getElementContainingPoint : invalid _spaceDimension must be equal to 2 or 3 !!!");
01810   PointLocator loc(*this);
01811   std::list<int> li=loc.locate(coord);
01812   if(li.empty())
01813     return -1;
01814   return li.front();
01815 }
01816 
01820 
01821 void MESH::convertToPoly()
01822 {
01823   if (getMeshDimension()!=3) return;
01824 
01825   CONNECTIVITY* newpolygonconnectivity = 0;
01826   CONNECTIVITY* newpolyhedraconnectivity = new CONNECTIVITY(MED_EN::MED_CELL);
01827 
01828   {
01830     // First step : Treating polygons connectivity
01832 
01833     int oldnbface = getNumberOfElements(MED_EN::MED_FACE,MED_EN::MED_ALL_ELEMENTS);
01834     int nbTypes = oldnbface > 0 ? 1 : 0;
01835     newpolygonconnectivity = new CONNECTIVITY(nbTypes, MED_EN::MED_FACE);
01836     if ( nbTypes > 0 )
01837     {
01838       medGeometryElement type = MED_POLYGON;
01839       newpolygonconnectivity->setGeometricTypes( &type, MED_FACE );
01840 
01841       const int count[] = { 1, oldnbface + 1 };
01842       newpolygonconnectivity->setCount( count, MED_FACE );
01843 
01844       const int* oldconn = getConnectivity(MED_EN::MED_NODAL, MED_EN::MED_FACE, MED_EN::MED_ALL_ELEMENTS);
01845       const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_FACE);
01846       newpolygonconnectivity->setNodal( oldconn, MED_FACE, type, oldconnindex );
01847 
01848       newpolygonconnectivity->setNumberOfNodes( getNumberOfNodes() );
01849       newpolygonconnectivity->setEntityDimension( 2 );
01850     }
01851   }
01853   // 2nd step : Treating polyhedra connectivity
01855   {
01856     vector<int> newconn;
01857     vector<int> newindex(1,1);
01858 
01859     int nboldtypes=getNumberOfTypes(MED_EN::MED_CELL);
01860     const MED_EN::medGeometryElement* oldtypes = getTypes(MED_EN::MED_CELL);
01861     const int* oldconn = getConnectivity(MED_EN::MED_NODAL, MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
01862     const int* oldconnindex= getConnectivityIndex(MED_EN::MED_NODAL,MED_EN::MED_CELL);
01863     int oldnbelem = getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
01864 
01865     for (int itype=0; itype<nboldtypes; itype++)
01866     {
01867       MED_EN::medGeometryElement type = oldtypes[itype];
01868       int nb_elems=getNumberOfElements(MED_EN::MED_CELL,type);
01869       if ( type == MED_POLYHEDRA )
01870       {
01871         const int* oldpolyindex = oldconnindex + getGlobalNumberingIndex( MED_CELL )[itype] - 1;
01872         int oldpolyconnsize = oldpolyindex[nb_elems] - oldpolyindex[0];
01873         newconn.insert( newconn.end(), oldconn, oldconn + oldpolyconnsize );
01874         int delta = newindex.back() - oldpolyindex[0];
01875         for (int ielem=0; ielem<nb_elems; ielem++)
01876           newindex.push_back( delta + oldpolyindex[ ielem+1 ]);
01877       }
01878       else
01879       {
01880         MEDMEM::CELLMODEL cellmodel(type);
01881         int nbfacespertype = cellmodel.getNumberOfConstituents(1);
01882         int nbnodespertype = cellmodel.getNumberOfNodes();
01883         for (int ielem=0; ielem<nb_elems; ielem++)
01884         {
01885           for (int iface=0; iface< nbfacespertype; iface++)
01886           {
01887             //local conn contains the local nodal connectivity for the iface-th face of type
01888             const int* local_conn = cellmodel.getNodesConstituent(1,iface+1);
01889             medGeometryElement facetype = cellmodel.getConstituentType(1,iface+1);
01890             int nbface_nodes=facetype%100;
01891             for ( int inode=0; inode<nbface_nodes;inode++)
01892               newconn.push_back(oldconn[local_conn[inode]-1]);
01893             if ( iface != nbfacespertype-1 )
01894               newconn.push_back(-1);
01895           }
01896           newindex.push_back( newconn.size() + 1 );
01897           oldconn += nbnodespertype;
01898         }
01899       }
01900     }
01901     int nbTypes = oldnbelem > 0 ? 1 : 0;
01902     if ( newpolyhedraconnectivity ) delete newpolyhedraconnectivity;
01903     newpolyhedraconnectivity = new CONNECTIVITY(nbTypes, MED_EN::MED_CELL);
01904     if ( nbTypes > 0 )
01905     {
01906       medGeometryElement type = MED_POLYHEDRA;
01907       newpolyhedraconnectivity->setGeometricTypes( &type, MED_CELL );
01908 
01909       const int count[] = { 1, oldnbelem + 1 };
01910       newpolyhedraconnectivity->setCount( count, MED_CELL );
01911 
01912       newpolyhedraconnectivity->setNodal( &newconn[0], MED_CELL, type, &newindex[0] );
01913 
01914       newpolyhedraconnectivity->setNumberOfNodes( getNumberOfNodes() );
01915       newpolyhedraconnectivity->setEntityDimension( 3 );
01916     }
01917   }
01918 
01919   delete _connectivity;
01920 
01921   _connectivity=newpolyhedraconnectivity;
01922   _connectivity->setConstituent(newpolygonconnectivity);
01923 
01924 }
01925 
01926 vector< vector<double> > MESH::getBoundingBox() const
01927 {
01928   const double *myCoords=_coordinate->getCoordinates(MED_EN::MED_FULL_INTERLACE);
01929   vector< vector<double> > ret(2);
01930   int i,j;
01931   ret[0].resize(_spaceDimension);
01932   ret[1].resize(_spaceDimension);
01933   for(i=0;i<_spaceDimension;i++)
01934   {
01935     ret[0][i]=1.e300;
01936     ret[1][i]=-1.e300;
01937   }
01938   for(i=0;i<_coordinate->getNumberOfNodes();i++)
01939   {
01940     for(j=0;j<_spaceDimension;j++)
01941     {
01942       double tmp=myCoords[i*_spaceDimension+j];
01943       if(tmp<ret[0][j])
01944         ret[0][j]=tmp;
01945       if(tmp>ret[1][j])
01946         ret[1][j]=tmp;
01947     }
01948   }
01949   return ret;
01950 }