1 // Copyright (C) 2005 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License.
9 // This library is distributed in the hope that it will be useful
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.salome-platform.org/
29 #include "MEDMEM_DriversDef.hxx"
30 #include "MEDMEM_Field.hxx"
31 #include "MEDMEM_Mesh.hxx"
33 #include "MEDMEM_Support.hxx"
34 #include "MEDMEM_Family.hxx"
35 #include "MEDMEM_Group.hxx"
36 #include "MEDMEM_Coordinate.hxx"
37 #include "MEDMEM_Connectivity.hxx"
38 #include "MEDMEM_CellModel.hxx"
39 #include "MEDMEM_Formulae.hxx"
40 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
41 #include "MEDMEM_DriverFactory.hxx"
44 using namespace MEDMEM;
45 using namespace MED_EN;
50 // ------- Drivers Management Part
52 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
53 is driverName. addDriver returns an integer handler. */
54 int MESH::addDriver(driverTypes driverType,
55 const string & fileName/*="Default File Name.med"*/,
56 const string & driverName/*="Default Mesh Name"*/,
57 med_mode_acces access) {
59 const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\",MED_EN::med_mode_acces access) : ";
67 driver = DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,
70 _drivers.push_back(driver);
72 int current = _drivers.size()-1;
74 _drivers[current]->setMeshName(driverName);
81 /*! Add an existing MESH driver. */
82 int MESH::addDriver(GENDRIVER & driver) {
83 const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
86 // A faire : Vérifier que le driver est de type MESH.
87 GENDRIVER * newDriver = driver.copy() ;
89 _drivers.push_back(newDriver);
90 return _drivers.size()-1;
95 /*! Remove an existing MESH driver. */
96 void MESH::rmDriver (int index/*=0*/) {
97 const char * LOC = "MESH::rmDriver (int index=0): ";
100 if ( _drivers[index] ) {
101 //_drivers.erase(&_drivers[index]);
103 MESSAGE ("detruire");
106 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
107 << "The index given is invalid, index must be between 0 and |"
116 // ------ End of Drivers Management Part
121 const char * LOC = "MESH::init(): ";
125 _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
127 _coordinate = (COORDINATE *) NULL;
128 _connectivity = (CONNECTIVITY *) NULL;
130 _spaceDimension = MED_INVALID; // 0 ?!?
131 _meshDimension = MED_INVALID;
132 _numberOfNodes = MED_INVALID;
136 _arePresentOptionnalNodesNumbers = 0;
141 /*! Create an empty MESH. */
142 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
149 _isAGrid = m._isAGrid;
151 if (m._coordinate != NULL)
152 _coordinate = new COORDINATE(* m._coordinate);
154 _coordinate = (COORDINATE *) NULL;
156 if (m._connectivity != NULL)
157 _connectivity = new CONNECTIVITY(* m._connectivity);
159 _connectivity = (CONNECTIVITY *) NULL;
161 _spaceDimension = m._spaceDimension;
162 _meshDimension = m._meshDimension;
163 _numberOfNodes = m._numberOfNodes;
165 _familyNode = m._familyNode;
166 for (int i=0; i<m._familyNode.size(); i++)
168 _familyNode[i] = new FAMILY(* m._familyNode[i]);
169 _familyNode[i]->setMesh(this);
172 _familyCell = m._familyCell;
173 for (int i=0; i<m._familyCell.size(); i++)
175 _familyCell[i] = new FAMILY(* m._familyCell[i]);
176 _familyCell[i]->setMesh(this);
179 _familyFace = m._familyFace;
180 for (int i=0; i<m._familyFace.size(); i++)
182 _familyFace[i] = new FAMILY(* m._familyFace[i]);
183 _familyFace[i]->setMesh(this);
186 _familyEdge = m._familyEdge;
187 for (int i=0; i<m._familyEdge.size(); i++)
189 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
190 _familyEdge[i]->setMesh(this);
193 _groupNode = m._groupNode;
194 for (int i=0; i<m._groupNode.size(); i++)
196 _groupNode[i] = new GROUP(* m._groupNode[i]);
197 _groupNode[i]->setMesh(this);
200 _groupCell = m._groupCell;
201 for (int i=0; i<m._groupCell.size(); i++)
203 _groupCell[i] = new GROUP(* m._groupCell[i]);
204 _groupCell[i]->setMesh(this);
207 _groupFace = m._groupFace;
208 for (int i=0; i<m._groupFace.size(); i++)
210 _groupFace[i] = new GROUP(* m._groupFace[i]);
211 _groupFace[i]->setMesh(this);
214 _groupEdge = m._groupEdge;
215 for (int i=0; i<m._groupEdge.size(); i++)
217 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
218 _groupEdge[i]->setMesh(this);
221 //_drivers = m._drivers; //Recopie des drivers?
226 MESSAGE("MESH::~MESH() : Destroying the Mesh");
227 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
228 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
230 size = _familyNode.size() ;
231 for (int i=0;i<size;i++)
232 delete _familyNode[i] ;
233 size = _familyCell.size() ;
234 for (int i=0;i<size;i++)
235 delete _familyCell[i] ;
236 size = _familyFace.size() ;
237 for (int i=0;i<size;i++)
238 delete _familyFace[i] ;
239 size = _familyEdge.size() ;
240 for (int i=0;i<size;i++)
241 delete _familyEdge[i] ;
242 size = _groupNode.size() ;
243 for (int i=0;i<size;i++)
244 delete _groupNode[i] ;
245 size = _groupCell.size() ;
246 for (int i=0;i<size;i++)
247 delete _groupCell[i] ;
248 size = _groupFace.size() ;
249 for (int i=0;i<size;i++)
250 delete _groupFace[i] ;
251 size = _groupEdge.size() ;
252 for (int i=0;i<size;i++)
253 delete _groupEdge[i] ;
255 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
257 for (unsigned int index=0; index < _drivers.size(); index++ )
259 SCRUTE(_drivers[index]);
260 if ( _drivers[index] != NULL) delete _drivers[index];
266 Method equivalent to getNumberOfTypes except that it includes not only classical Types but polygons/polyhedra also.
268 int MESH::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
270 if(_connectivity!= NULL)
271 return _connectivity->getNumberOfTypesWithPoly(Entity);
272 throw MEDEXCEPTION(LOCALIZED("MESH::getNumberOfTypesWithPoly( medEntityMesh ) : Connectivity not defined !"));
276 Method equivalent to getTypesWithPoly except that it includes not only classical Types but polygons/polyhedra also.
277 WARNING the returned array MUST be deallocated.
279 MED_EN::medGeometryElement * MESH::getTypesWithPoly(MED_EN::medEntityMesh Entity) const
281 if (Entity == MED_EN::MED_NODE)
282 throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : No medGeometryElement with MED_NODE entity !"));
283 if (_connectivity != NULL)
284 return _connectivity->getGeometricTypesWithPoly(Entity);
285 throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : Connectivity not defined !"));
289 Method equivalent to getNumberOfElementsWithPoly except that it includes not only classical Types but polygons/polyhedra also.
291 int MESH::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
293 if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
295 int nbOfPolygs=_connectivity->getNumberOfElementOfPolyType(Entity);
298 else if(Type==MED_ALL_ELEMENTS)
300 int nbOfClassicalTypes=getNumberOfElements(Entity,MED_ALL_ELEMENTS);
301 int nbOfClassicalTypes2=_connectivity->getNumberOfElementOfPolyType(Entity);
302 return nbOfClassicalTypes+nbOfClassicalTypes2;
305 return getNumberOfElements(Entity,Type);
308 MESH & MESH::operator=(const MESH &m)
310 const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
313 MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
316 // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
317 // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
318 // _drivers = m._drivers;
320 // space_dimension=m.space_dimension;
321 // mesh_dimension=m.mesh_dimension;
323 // nodes_count=m.nodes_count;
324 // coordinates=m.coordinates;
325 // coordinates_name=m.coordinates_name ;
326 // coordinates_unit=m.coordinates_unit ;
327 // nodes_numbers=m.nodes_numbers ;
329 // cells_types_count=m.cells_types_count;
330 // cells_type=m.cells_type;
331 // cells_count=m.cells_count;
332 // nodal_connectivity=m.nodal_connectivity;
334 // nodes_families_count=m.nodes_families_count;
335 // nodes_Families=m.nodes_Families;
337 // cells_families_count=m.cells_families_count;
338 // cells_Families=m.cells_Families;
340 // maximum_cell_number_by_node = m.maximum_cell_number_by_node;
341 // if (maximum_cell_number_by_node > 0)
343 // reverse_nodal_connectivity = m.reverse_nodal_connectivity;
344 // reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
351 bool MESH::operator==(const MESH& other) const
353 BEGIN_OF("MESH::operator==");
357 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
358 The meshname driverName must already exists in the file.*/
359 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) throw (MEDEXCEPTION)
361 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
367 GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,MED_LECT);
368 current = addDriver(*myDriver);
370 _drivers[current]->open();
371 _drivers[current]->read();
372 _drivers[current]->close();
378 for a deep comparison of 2 meshes.
380 bool MESH::deepCompare(const MESH& other) const
382 int size1=getSpaceDimension()*getNumberOfNodes();
383 int size2=other.getSpaceDimension()*other.getNumberOfNodes();
386 const double* coord1=getCoordinates(MED_FULL_INTERLACE);
387 const double* coord2=other.getCoordinates(MED_FULL_INTERLACE);
389 for(int i=0;i<size1 && ret;i++)
391 ret=(fabs(coord1[i]-coord2[i])<1e-15);
395 return _connectivity->deepCompare(*other._connectivity);
400 ostream & ::MEDMEM::operator<<(ostream &os, const MESH &myMesh)
402 int spacedimension = myMesh.getSpaceDimension();
403 int meshdimension = myMesh.getMeshDimension();
404 int numberofnodes = myMesh.getNumberOfNodes();
406 os << "Space Dimension : " << spacedimension << endl << endl;
408 os << "Mesh Dimension : " << meshdimension << endl << endl;
410 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
411 os << "SHOW NODES COORDINATES : " << endl;
413 os << "Name :" << endl;
414 const string * coordinatesnames = myMesh.getCoordinatesNames();
415 for(int i=0; i<spacedimension ; i++)
417 os << " - " << coordinatesnames[i] << endl;
419 os << "Unit :" << endl;
420 const string * coordinatesunits = myMesh.getCoordinatesUnits();
421 for(int i=0; i<spacedimension ; i++)
423 os << " - " << coordinatesunits[i] << endl;
425 for(int i=0; i<numberofnodes ; i++)
427 os << "Nodes " << i+1 << " : ";
428 for (int j=0; j<spacedimension ; j++)
429 os << coordinates[i*spacedimension+j] << " ";
433 os << endl << "SHOW CONNECTIVITY :" << endl;
434 os << *myMesh._connectivity << endl;
436 medEntityMesh entity;
437 os << endl << "SHOW FAMILIES :" << endl << endl;
438 for (int k=1; k<=4; k++)
440 if (k==1) entity = MED_NODE;
441 if (k==2) entity = MED_CELL;
442 if (k==3) entity = MED_FACE;
443 if (k==4) entity = MED_EDGE;
444 int numberoffamilies = myMesh.getNumberOfFamilies(entity);
445 os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
446 using namespace MED_EN;
447 for (int i=1; i<numberoffamilies+1;i++)
449 os << * myMesh.getFamily(entity,i) << endl;
453 os << endl << "SHOW GROUPS :" << endl << endl;
454 for (int k=1; k<=4; k++)
456 if (k==1) entity = MED_NODE;
457 if (k==2) entity = MED_CELL;
458 if (k==3) entity = MED_FACE;
459 if (k==4) entity = MED_EDGE;
460 int numberofgroups = myMesh.getNumberOfGroups(entity);
461 os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
462 using namespace MED_EN;
463 for (int i=1; i<numberofgroups+1;i++)
465 os << * myMesh.getGroup(entity,i) << endl;
473 Get global number of element which have same connectivity than connectivity argument.
475 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
477 Return -1 if not found.
479 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
481 const char* LOC="MESH::getElementNumber " ;
484 int numberOfValue ; // size of connectivity array
485 CELLMODEL myType(Type) ;
486 if (ConnectivityType==MED_DESCENDING)
487 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
489 numberOfValue = myType.getNumberOfNodes() ; // nodes
491 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
492 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
494 // First node or face/edge
495 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
496 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
498 // list to put cells numbers
499 list<int> cellsList ;
500 list<int>::iterator itList ;
501 for (int i=indexBegin; i<indexEnd; i++)
502 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
504 // Foreach node or face/edge in cell, we search cells which are in common.
505 // At the end, we must have only one !
507 for (int i=1; i<numberOfValue; i++) {
508 int connectivity_i = connectivity[i] ;
509 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
511 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
512 if ((*itList)==myReverseConnectivityValue[j-1]) {
518 itList=cellsList.erase(itList);
519 itList--; // well : rigth if itList = cellsList.begin() ??
524 if (cellsList.size()>1) // we have more than one cell
525 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
527 if (cellsList.size()==0)
532 return cellsList.front() ;
536 Return a support which reference all elements on the boundary of mesh.
538 For instance, we could get only face in 3D and edge in 2D.
540 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
543 const char * LOC = "MESH::getBoundaryElements : " ;
546 // actually we could only get face (in 3D) and edge (in 2D)
547 medEntityMesh entityToParse=Entity;
548 if(_spaceDimension == 3)
549 if (Entity != MED_FACE)
551 entityToParse=MED_FACE;
553 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
554 if(_spaceDimension == 2)
555 if(Entity != MED_EDGE)
557 entityToParse=MED_EDGE;
559 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
561 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
562 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
563 int numberOf = getNumberOfElementsWithPoly(entityToParse,MED_ALL_ELEMENTS) ;
564 list<int> myElementsList;
566 for (int i=0 ; i<numberOf; i++)
567 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
568 myElementsList.push_back(i+1);
572 return buildSupportOnNodeFromElementList(myElementsList,entityToParse);
574 return buildSupportOnElementsFromElementList(myElementsList,entityToParse);
578 Method that do the same thing as buildSupportOnNodeFromElementList except that a SUPPORT is not created.
580 void MESH::fillSupportOnNodeFromElementList(const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION)
582 MED_EN::medEntityMesh entity=supportToFill->getEntity();
583 supportToFill->setAll(false);
584 supportToFill->setMesh((MESH *)this);
588 for(list<int>::const_iterator iter=listOfElt.begin();iter!=listOfElt.end();iter++)
591 const int *conn=_connectivity->getConnectivityOfAnElementWithPoly(MED_NODAL,entity,*iter,lgth);
593 nodes.insert(conn[i]);
596 for(set<int>::iterator iter2=nodes.begin();iter2!=nodes.end();iter2++)
597 nodesList.push_back(*iter2);
598 supportToFill->fillFromNodeList(nodesList);
602 Method created to factorize code. This method creates a new support on NODE (to deallocate) containing all the nodes id contained in elements 'listOfElt' of
605 SUPPORT *MESH::buildSupportOnNodeFromElementList(const list<int>& listOfElt,MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
607 SUPPORT * mySupport = new SUPPORT((MESH *)this,"Boundary",entity);
608 fillSupportOnNodeFromElementList(listOfElt,mySupport);
613 Method created to factorize code. This method creates a new support on entity 'entity' (to deallocate) containing all the entities contained in
614 elements 'listOfElt' of entity 'entity'.
616 SUPPORT *MESH::buildSupportOnElementsFromElementList(const list<int>& listOfElt, MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
618 const char * LOC = "MESH::buildSupportOnElementsFromElementList : " ;
620 SUPPORT *mySupport=new SUPPORT((MESH *)this,"Boundary",entity);
621 mySupport->fillFromElementList(listOfElt);
626 FIELD<double, FullInterlace>* MESH::getVolume(const SUPPORT *Support) const throw (MEDEXCEPTION)
628 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
631 // Support must be on 3D elements
633 // Make sure that the MESH class is the same as the MESH class attribut
634 // in the class Support
635 MESH* myMesh = Support->getMesh();
637 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
639 int dim_space = getSpaceDimension();
640 medEntityMesh support_entity = Support->getEntity();
641 bool onAll = Support->isOnAllElements();
643 int nb_type, length_values;
644 const medGeometryElement* types;
646 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
647 const int* global_connectivity;
648 nb_type = Support->getNumberOfTypes();
649 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
650 types = Support->getTypes();
652 FIELD<double, FullInterlace>* Volume =
653 new FIELD<double, FullInterlace>(Support,1);
654 // double *volume = new double [length_values];
655 Volume->setName("VOLUME");
656 Volume->setDescription("cells volume");
657 Volume->setComponentName(1,"volume");
658 Volume->setComponentDescription(1,"desc-comp");
660 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
662 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
664 Volume->setMEDComponentUnit(1,MEDComponentUnit);
666 Volume->setIterationNumber(0);
667 Volume->setOrderNumber(0);
668 Volume->setTime(0.0);
670 //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
671 typedef MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
672 ArrayNoGauss *volume = Volume->getArrayNoGauss();
675 const double * coord = getCoordinates(MED_FULL_INTERLACE);
677 for (int i=0;i<nb_type;i++)
679 medGeometryElement type = types[i] ;
681 nb_entity_type = Support->getNumberOfElements(type);
682 if(type != MED_EN::MED_POLYHEDRA)
686 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
690 const int * supp_number = Support->getNumber(type);
691 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
692 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
693 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
695 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
696 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
697 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
700 global_connectivity = global_connectivity_tmp ;
706 case MED_TETRA4 : case MED_TETRA10 :
708 for (int tetra=0;tetra<nb_entity_type;tetra++)
710 int tetra_index = (type%100)*tetra;
711 int N1 = global_connectivity[tetra_index]-1;
712 int N2 = global_connectivity[tetra_index+1]-1;
713 int N3 = global_connectivity[tetra_index+2]-1;
714 int N4 = global_connectivity[tetra_index+3]-1;
715 xvolume=CalculateVolumeForTetra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4);
716 volume->setIJ(index,1,xvolume) ;
721 case MED_PYRA5 : case MED_PYRA13 :
723 for (int pyra=0;pyra<nb_entity_type;pyra++)
725 int pyra_index = (type%100)*pyra;
726 int N1 = global_connectivity[pyra_index]-1;
727 int N2 = global_connectivity[pyra_index+1]-1;
728 int N3 = global_connectivity[pyra_index+2]-1;
729 int N4 = global_connectivity[pyra_index+3]-1;
730 int N5 = global_connectivity[pyra_index+4]-1;
731 xvolume=CalculateVolumeForPyra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5);
732 volume->setIJ(index,1,xvolume) ;
737 case MED_PENTA6 : case MED_PENTA15 :
739 for (int penta=0;penta<nb_entity_type;penta++)
741 int penta_index = (type%100)*penta;
742 int N1 = global_connectivity[penta_index]-1;
743 int N2 = global_connectivity[penta_index+1]-1;
744 int N3 = global_connectivity[penta_index+2]-1;
745 int N4 = global_connectivity[penta_index+3]-1;
746 int N5 = global_connectivity[penta_index+4]-1;
747 int N6 = global_connectivity[penta_index+5]-1;
748 xvolume=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);
749 volume->setIJ(index,1,xvolume) ;
754 case MED_HEXA8 : case MED_HEXA20 :
756 for (int hexa=0;hexa<nb_entity_type;hexa++)
758 int hexa_index = (type%100)*hexa;
760 int N1 = global_connectivity[hexa_index]-1;
761 int N2 = global_connectivity[hexa_index+1]-1;
762 int N3 = global_connectivity[hexa_index+2]-1;
763 int N4 = global_connectivity[hexa_index+3]-1;
764 int N5 = global_connectivity[hexa_index+4]-1;
765 int N6 = global_connectivity[hexa_index+5]-1;
766 int N7 = global_connectivity[hexa_index+6]-1;
767 int N8 = global_connectivity[hexa_index+7]-1;
768 xvolume=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);
769 volume->setIJ(index,1,xvolume) ;
779 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
781 int lgthNodes,iPts,iFaces,iPtsInFace;
782 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
783 int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
784 int nbOfFaces,*nbOfNodesPerFaces;
785 int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(offsetWithClassicType+polyhs+1,nbOfFaces,nbOfNodesPerFaces);
786 double **pts=new double * [lgthNodes];
787 double ***pts1=new double ** [nbOfFaces];
788 for(iPts=0;iPts<lgthNodes;iPts++)
789 pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
790 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
792 pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
793 for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
794 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
797 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
800 xvolume=CalculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
801 delete [] nbOfNodesPerFaces;
802 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
803 delete [] pts1[iFaces];
805 volume->setIJ(index,1,xvolume) ;
811 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
812 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
814 int lgthNodes,iPts,iFaces,iPtsInFace;
815 int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
816 int nbOfFaces,*nbOfNodesPerFaces;
817 int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(supp_number[polyhs],nbOfFaces,nbOfNodesPerFaces);
818 double **pts=new double * [lgthNodes];
819 double ***pts1=new double ** [nbOfFaces];
820 for(iPts=0;iPts<lgthNodes;iPts++)
821 pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
822 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
824 pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
825 for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
826 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
829 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
832 xvolume=CalculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
833 delete [] nbOfNodesPerFaces;
834 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
835 delete [] pts1[iFaces];
837 volume->setIJ(index,1,xvolume) ;
844 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
848 if (!onAll && type!=MED_EN::MED_POLYHEDRA)
849 delete [] global_connectivity ;
855 FIELD<double, FullInterlace>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
857 const char * LOC = "MESH::getArea(SUPPORT*) : ";
860 // Support must be on 2D elements
862 // Make sure that the MESH class is the same as the MESH class attribut
863 // in the class Support
864 MESH* myMesh = Support->getMesh();
866 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
868 int dim_space = getSpaceDimension();
869 medEntityMesh support_entity = Support->getEntity();
870 bool onAll = Support->isOnAllElements();
872 int nb_type, length_values;
873 const medGeometryElement* types;
875 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
876 const int* global_connectivity;
878 nb_type = Support->getNumberOfTypes();
879 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
880 types = Support->getTypes();
883 FIELD<double, FullInterlace>* Area;
885 Area = new FIELD<double, FullInterlace>(Support,1);
886 Area->setName("AREA");
887 Area->setDescription("cells or faces area");
889 Area->setComponentName(1,"area");
890 Area->setComponentDescription(1,"desc-comp");
892 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
894 string MEDComponentUnit="trucmuch";
896 Area->setMEDComponentUnit(1,MEDComponentUnit);
898 Area->setIterationNumber(0);
899 Area->setOrderNumber(0);
902 double *area = (double *)Area->getValue();
904 const double * coord = getCoordinates(MED_FULL_INTERLACE);
907 for (int i=0;i<nb_type;i++)
909 medGeometryElement type = types[i] ;
910 nb_entity_type = Support->getNumberOfElements(type);
911 const int *global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
912 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
916 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
920 const int * supp_number = Support->getNumber(type);
921 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
922 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
924 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
925 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
926 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
929 global_connectivity = global_connectivity_tmp ;
934 case MED_TRIA3 : case MED_TRIA6 :
936 for (int tria=0;tria<nb_entity_type;tria++)
938 int tria_index = (type%100)*tria;
940 int N1 = global_connectivity[tria_index]-1;
941 int N2 = global_connectivity[tria_index+1]-1;
942 int N3 = global_connectivity[tria_index+2]-1;
944 area[index]=CalculateAreaForTria(coord+(dim_space*N1),
945 coord+(dim_space*N2),
946 coord+(dim_space*N3),dim_space);
951 case MED_QUAD4 : case MED_QUAD8 :
953 for (int quad=0;quad<nb_entity_type;quad++)
955 int quad_index = (type%100)*quad;
957 int N1 = global_connectivity[quad_index]-1;
958 int N2 = global_connectivity[quad_index+1]-1;
959 int N3 = global_connectivity[quad_index+2]-1;
960 int N4 = global_connectivity[quad_index+3]-1;
962 area[index]=CalculateAreaForQuad(coord+dim_space*N1,
965 coord+dim_space*N4,dim_space);
974 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
975 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
976 for (int polygs=0;polygs<nb_entity_type;polygs++)
978 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
979 double **pts=new double * [size];
980 for(int iPts=0;iPts<size;iPts++)
981 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1));
982 area[index] = CalculateAreaForPolyg((const double **)pts,size,dim_space);
989 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
990 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
991 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
992 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
993 for (int polygs=0;polygs<nb_entity_type;polygs++)
995 int size=connectivity_index[supp_number[polygs]-offsetWithClassicType]-connectivity_index[supp_number[polygs]-offsetWithClassicType-1];
996 double **pts=new double * [size];
997 for(int iPts=0;iPts<size;iPts++)
998 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-offsetWithClassicType-1]+iPts-1]-1));
999 area[index]=CalculateAreaForPolyg((const double **)pts,size,dim_space);
1007 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1012 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1013 delete [] global_connectivity ;
1018 FIELD<double, FullInterlace>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1020 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1023 // Support must be on 1D elements
1025 // Make sure that the MESH class is the same as the MESH class attribut
1026 // in the class Support
1027 MESH* myMesh = Support->getMesh();
1029 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1031 int dim_space = getSpaceDimension();
1032 medEntityMesh support_entity = Support->getEntity();
1033 bool onAll = Support->isOnAllElements();
1035 int nb_type, length_values;
1036 const medGeometryElement* types;
1038 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1039 const int* global_connectivity;
1041 nb_type = Support->getNumberOfTypes();
1042 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1043 types = Support->getTypes();
1046 FIELD<double, FullInterlace>* Length;
1048 Length = new FIELD<double, FullInterlace>(Support,1);
1049 // double *length = new double [length_values];
1051 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1052 typedef MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array ArrayNoGauss;
1053 ArrayNoGauss * length = Length->getArrayNoGauss();
1055 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1058 for (int i=0;i<nb_type;i++)
1060 medGeometryElement type = types[i] ;
1065 nb_entity_type = getNumberOfElements(support_entity,type);
1066 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1070 nb_entity_type = Support->getNumberOfElements(type);
1071 const int * supp_number = Support->getNumber(type);
1072 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1073 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1074 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1076 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1077 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1078 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1081 global_connectivity = global_connectivity_tmp ;
1086 case MED_SEG2 : case MED_SEG3 :
1088 for (int edge=0;edge<nb_entity_type;edge++)
1090 int edge_index = (type%100)*edge;
1092 int N1 = global_connectivity[edge_index]-1;
1093 int N2 = global_connectivity[edge_index+1]-1;
1095 double x1 = coord[dim_space*N1];
1096 double x2 = coord[dim_space*N2];
1098 double y1 = coord[dim_space*N1+1];
1099 double y2 = coord[dim_space*N2+1];
1101 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1102 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1104 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1105 (z1 - z2)*(z1 - z2));
1107 length->setIJ(index,1,xlength) ;
1113 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1117 if (!onAll) delete [] global_connectivity ;
1123 FIELD<double, FullInterlace>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1125 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1128 // Support must be on 2D or 1D elements
1130 // Make sure that the MESH class is the same as the MESH class attribut
1131 // in the class Support
1132 MESH* myMesh = Support->getMesh();
1134 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1136 int dim_space = getSpaceDimension();
1137 int mesh_dim=getMeshDimension();
1138 medEntityMesh support_entity = Support->getEntity();
1139 bool onAll = Support->isOnAllElements();
1141 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 ))
1142 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"incompatible mesh dimension and entity"));
1143 int nb_type, length_values;
1144 const medGeometryElement* types;
1146 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1147 const int* global_connectivity;
1149 nb_type = Support->getNumberOfTypes();
1150 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1151 types = Support->getTypes();
1155 FIELD<double, FullInterlace>* Normal =
1156 new FIELD<double, FullInterlace>(Support,dim_space);
1157 Normal->setName("NORMAL");
1158 Normal->setDescription("cells or faces normal");
1159 for (int k=1;k<=dim_space;k++) {
1160 Normal->setComponentName(k,"normal");
1161 Normal->setComponentDescription(k,"desc-comp");
1162 Normal->setMEDComponentUnit(k,"unit");
1165 Normal->setIterationNumber(MED_NOPDT);
1166 Normal->setOrderNumber(MED_NONOR);
1167 Normal->setTime(0.0);
1168 double *normal = (double *)Normal->getValue();
1170 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1173 for (int i=0;i<nb_type;i++)
1175 medGeometryElement type = types[i] ;
1176 nb_entity_type = Support->getNumberOfElements(type);
1178 // Make sure we trying to get Normals on
1179 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1180 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1182 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1183 (type==MED_QUAD4) || (type==MED_QUAD8) || (type==MED_POLYGON)) &&
1184 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1186 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1187 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1191 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1195 const int * supp_number = Support->getNumber(type);
1196 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1197 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1198 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1200 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1201 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1202 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1206 global_connectivity = global_connectivity_tmp ;
1212 case MED_TRIA3 : case MED_TRIA6 :
1214 for (int tria=0;tria<nb_entity_type;tria++)
1216 int tria_index = (type%100)*tria;
1217 int N1 = global_connectivity[tria_index]-1;
1218 int N2 = global_connectivity[tria_index+1]-1;
1219 int N3 = global_connectivity[tria_index+2]-1;
1220 CalculateNormalForTria(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,normal+3*index);
1225 case MED_QUAD4 : case MED_QUAD8 :
1227 for (int quad=0;quad<nb_entity_type;quad++)
1229 int quad_index = (type%100)*quad;
1230 int N1 = global_connectivity[quad_index]-1;
1231 int N2 = global_connectivity[quad_index+1]-1;
1232 int N3 = global_connectivity[quad_index+2]-1;
1233 int N4 = global_connectivity[quad_index+3]-1;
1234 CalculateNormalForQuad(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,normal+3*index);
1239 case MED_SEG2 : case MED_SEG3 :
1241 double xnormal1, xnormal2;
1242 for (int edge=0;edge<nb_entity_type;edge++)
1244 int edge_index = (type%100)*edge;
1246 int N1 = global_connectivity[edge_index]-1;
1247 int N2 = global_connectivity[edge_index+1]-1;
1249 double x1 = coord[dim_space*N1];
1250 double x2 = coord[dim_space*N2];
1252 double y1 = coord[dim_space*N1+1];
1253 double y2 = coord[dim_space*N2+1];
1255 xnormal1 = -(y2-y1);
1258 normal[2*index] = xnormal1 ;
1259 normal[2*index+1] = xnormal2 ;
1269 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1270 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1271 for (int polygs=0;polygs<nb_entity_type;polygs++)
1273 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1274 double **pts=new double * [size];
1275 for(int iPts=0;iPts<size;iPts++)
1276 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1])-1);
1277 CalculateNormalForPolyg((const double **)pts,size,normal+3*index);
1284 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1285 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1286 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1287 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1288 for (int polygs=0;polygs<nb_entity_type;polygs++)
1290 int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1291 int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1292 double **pts=new double * [size];
1293 for(int iPts=0;iPts<size;iPts++)
1294 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1])-1);
1295 CalculateNormalForPolyg((const double **)pts,size,normal+3*index);
1303 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1306 if (!onAll && type!=MED_EN::MED_POLYGON)
1307 delete [] global_connectivity ;
1314 FIELD<double, FullInterlace>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1316 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1317 MESH* myMesh = Support->getMesh();
1319 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1321 int dim_space = getSpaceDimension();
1322 medEntityMesh support_entity = Support->getEntity();
1323 bool onAll = Support->isOnAllElements();
1325 int nb_type, length_values;
1326 const medGeometryElement* types;
1328 const int* global_connectivity;
1329 const int * global_connectivityIndex;
1331 nb_type = Support->getNumberOfTypes();
1332 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1333 types = Support->getTypes();
1335 FIELD<double, FullInterlace>* Barycenter;
1336 Barycenter = new FIELD<double, FullInterlace>(Support,dim_space);
1337 Barycenter->setName("BARYCENTER");
1338 Barycenter->setDescription("cells or faces barycenter");
1340 for (int k=0;k<dim_space;k++) {
1342 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1343 Barycenter->setComponentDescription(kp1,"desc-comp");
1344 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1346 Barycenter->setIterationNumber(0);
1347 Barycenter->setOrderNumber(0);
1348 Barycenter->setTime(0.0);
1349 double *barycenter=(double *)Barycenter->getValue();
1350 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1352 for (int i=0;i<nb_type;i++)
1354 medGeometryElement type = types[i] ;
1355 nb_entity_type = Support->getNumberOfElements(type);
1356 global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1357 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA )
1361 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1365 const int * supp_number = Support->getNumber(type);
1366 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1367 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1369 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1370 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1371 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1374 global_connectivity = global_connectivity_tmp;
1380 case MED_TETRA4 : case MED_TETRA10 :
1382 for (int tetra =0;tetra<nb_entity_type;tetra++)
1384 int tetra_index = (type%100)*tetra;
1386 int N1 = global_connectivity[tetra_index]-1;
1387 int N2 = global_connectivity[tetra_index+1]-1;
1388 int N3 = global_connectivity[tetra_index+2]-1;
1389 int N4 = global_connectivity[tetra_index+3]-1;
1391 pts[0]=(double *)coord+dim_space*N1;
1392 pts[1]=(double *)coord+dim_space*N2;
1393 pts[2]=(double *)coord+dim_space*N3;
1394 pts[3]=(double *)coord+dim_space*N4;
1395 CalculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1400 case MED_PYRA5 : case MED_PYRA13 :
1402 for (int pyra=0;pyra<nb_entity_type;pyra++)
1404 int pyra_index = (type%100)*pyra;
1406 int N1 = global_connectivity[pyra_index]-1;
1407 int N2 = global_connectivity[pyra_index+1]-1;
1408 int N3 = global_connectivity[pyra_index+2]-1;
1409 int N4 = global_connectivity[pyra_index+3]-1;
1410 int N5 = global_connectivity[pyra_index+4]-1;
1412 pts[0]=(double *)coord+dim_space*N1;
1413 pts[1]=(double *)coord+dim_space*N2;
1414 pts[2]=(double *)coord+dim_space*N3;
1415 pts[3]=(double *)coord+dim_space*N4;
1416 pts[4]=(double *)coord+dim_space*N5;
1417 CalculateBarycenter<5,3>((const double **)pts,barycenter+3*index);
1422 case MED_PENTA6 : case MED_PENTA15 :
1424 for (int penta=0;penta<nb_entity_type;penta++)
1426 int penta_index = (type%100)*penta;
1428 int N1 = global_connectivity[penta_index]-1;
1429 int N2 = global_connectivity[penta_index+1]-1;
1430 int N3 = global_connectivity[penta_index+2]-1;
1431 int N4 = global_connectivity[penta_index+3]-1;
1432 int N5 = global_connectivity[penta_index+4]-1;
1433 int N6 = global_connectivity[penta_index+5]-1;
1435 pts[0]=(double *)coord+dim_space*N1;
1436 pts[1]=(double *)coord+dim_space*N2;
1437 pts[2]=(double *)coord+dim_space*N3;
1438 pts[3]=(double *)coord+dim_space*N4;
1439 pts[4]=(double *)coord+dim_space*N5;
1440 pts[5]=(double *)coord+dim_space*N6;
1441 CalculateBarycenter<6,3>((const double **)pts,barycenter+3*index);
1446 case MED_HEXA8 : case MED_HEXA20 :
1448 for (int hexa=0;hexa<nb_entity_type;hexa++)
1450 int hexa_index = (type%100)*hexa;
1452 int N1 = global_connectivity[hexa_index]-1;
1453 int N2 = global_connectivity[hexa_index+1]-1;
1454 int N3 = global_connectivity[hexa_index+2]-1;
1455 int N4 = global_connectivity[hexa_index+3]-1;
1456 int N5 = global_connectivity[hexa_index+4]-1;
1457 int N6 = global_connectivity[hexa_index+5]-1;
1458 int N7 = global_connectivity[hexa_index+6]-1;
1459 int N8 = global_connectivity[hexa_index+7]-1;
1461 pts[0]=(double *)coord+dim_space*N1;
1462 pts[1]=(double *)coord+dim_space*N2;
1463 pts[2]=(double *)coord+dim_space*N3;
1464 pts[3]=(double *)coord+dim_space*N4;
1465 pts[4]=(double *)coord+dim_space*N5;
1466 pts[5]=(double *)coord+dim_space*N6;
1467 pts[6]=(double *)coord+dim_space*N7;
1468 pts[7]=(double *)coord+dim_space*N8;
1469 CalculateBarycenter<8,3>((const double **)pts,barycenter+3*index);
1474 case MED_TRIA3 : case MED_TRIA6 :
1476 for (int tria=0;tria<nb_entity_type;tria++)
1478 int tria_index = (type%100)*tria;
1479 int N1 = global_connectivity[tria_index]-1;
1480 int N2 = global_connectivity[tria_index+1]-1;
1481 int N3 = global_connectivity[tria_index+2]-1;
1483 pts[0]=(double *)coord+dim_space*N1;
1484 pts[1]=(double *)coord+dim_space*N2;
1485 pts[2]=(double *)coord+dim_space*N3;
1487 CalculateBarycenter<3,2>((const double **)pts,barycenter+2*index);
1489 CalculateBarycenter<3,3>((const double **)pts,barycenter+3*index);
1494 case MED_QUAD4 : case MED_QUAD8 :
1496 for (int quad=0;quad<nb_entity_type;quad++)
1498 int quad_index = (type%100)*quad;
1499 int N1 = global_connectivity[quad_index]-1;
1500 int N2 = global_connectivity[quad_index+1]-1;
1501 int N3 = global_connectivity[quad_index+2]-1;
1502 int N4 = global_connectivity[quad_index+3]-1;
1504 pts[0]=(double *)coord+dim_space*N1;
1505 pts[1]=(double *)coord+dim_space*N2;
1506 pts[2]=(double *)coord+dim_space*N3;
1507 pts[3]=(double *)coord+dim_space*N4;
1509 CalculateBarycenter<4,2>((const double **)pts,barycenter+2*index);
1511 CalculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1516 case MED_SEG2 : case MED_SEG3 :
1518 for (int edge=0;edge<nb_entity_type;edge++)
1520 int edge_index = (type%100)*edge;
1521 int N1 = global_connectivity[edge_index]-1;
1522 int N2 = global_connectivity[edge_index+1]-1;
1524 pts[0]=(double *)coord+dim_space*N1;
1525 pts[1]=(double *)coord+dim_space*N2;
1527 CalculateBarycenter<2,2>((const double **)pts,barycenter+2*index);
1529 CalculateBarycenter<2,3>((const double **)pts,barycenter+3*index);
1538 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1539 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1540 for (int polygs=0;polygs<nb_entity_type;polygs++)
1542 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1543 double **pts=new double * [size];
1544 for(int iPts=0;iPts<size;iPts++)
1545 pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1);
1546 CalculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1552 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1553 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1554 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1555 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1556 for (int polygs=0;polygs<nb_entity_type;polygs++)
1558 int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1559 int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1560 double **pts=new double * [size];
1561 for(int iPts=0;iPts<size;iPts++)
1562 pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1]-1);
1563 CalculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1570 case MED_EN::MED_POLYHEDRA:
1574 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1577 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1578 int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
1579 double **pts=new double * [lgthNodes];
1580 for(int iPts=0;iPts<lgthNodes;iPts++)
1581 pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1582 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1590 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
1591 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1594 int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
1595 double **pts=new double * [lgthNodes];
1596 for(int iPts=0;iPts<lgthNodes;iPts++)
1597 pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1598 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1607 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
1612 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1613 delete [] global_connectivity;
1619 bool MESH::isEmpty() const
1621 bool notempty = _name != "NOT DEFINED" || _coordinate != NULL || _connectivity != NULL ||
1622 _spaceDimension !=MED_INVALID || _meshDimension !=MED_INVALID ||
1623 _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
1624 _familyNode.size() != 0 || _groupCell.size() != 0 ||
1625 _familyCell.size() != 0 || _groupFace.size() != 0 ||
1626 _familyFace.size() != 0 || _groupEdge.size() != 0 ||
1627 _familyEdge.size() != 0 || _isAGrid != 0 ;
1631 void MESH::read(int index)
1633 const char * LOC = "MESH::read(int index=0) : ";
1636 if (_drivers[index]) {
1637 _drivers[index]->open();
1638 _drivers[index]->read();
1639 _drivers[index]->close();
1642 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1643 << "The index given is invalid, index must be between 0 and |"
1649 //=======================================================================
1650 //function : getSkin
1652 //=======================================================================
1654 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
1656 const char * LOC = "MESH::getSkin : " ;
1659 if (this != Support3D->getMesh())
1660 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
1661 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
1662 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
1664 // well, all rigth !
1665 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
1666 mySupport->setAll(false);
1668 list<int> myElementsList ;
1671 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
1672 if (Support3D->isOnAllElements())
1674 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
1675 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
1676 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
1678 int cellNb1 = myConnectivityValue [i];
1679 int cellNb2 = myConnectivityValue [i+1];
1680 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
1681 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
1683 myElementsList.push_back( j ) ;
1690 map<int,int> FaceNbEncounterNb;
1692 int * myConnectivityValue = const_cast <int *>
1693 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
1694 MED_CELL, MED_ALL_ELEMENTS));
1695 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
1696 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
1697 int nbCells = Support3D->getnumber()->getLength();
1698 for (i=0; i<nbCells; ++i)
1700 int cellNb = myCellNbs [ i ];
1701 int faceFirst = myConnectivityIndex[ cellNb-1 ];
1702 int faceLast = myConnectivityIndex[ cellNb ];
1703 for (j = faceFirst; j < faceLast; ++j)
1705 int faceNb = abs( myConnectivityValue [ j-1 ] );
1706 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
1707 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
1708 FaceNbEncounterNb[ faceNb ] = 1;
1710 FaceNbEncounterNb[ faceNb ] ++;
1713 map<int,int>::iterator FaceNbEncounterNbItr;
1714 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
1715 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
1716 FaceNbEncounterNbItr ++)
1717 if ((*FaceNbEncounterNbItr).second == 1)
1719 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
1723 // Well, we must know how many geometric type we have found
1724 int * myListArray = new int[size] ;
1726 list<int>::iterator myElementsListIt ;
1727 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
1728 myListArray[id]=(*myElementsListIt) ;
1732 int numberOfGeometricType ;
1733 medGeometryElement* geometricType ;
1734 int * geometricTypeNumber ;
1735 int * numberOfEntities ;
1736 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
1737 int * mySkyLineArrayIndex ;
1739 int numberOfType = getNumberOfTypes(MED_FACE) ;
1740 if (numberOfType == 1) { // wonderfull : it's easy !
1741 numberOfGeometricType = 1 ;
1742 geometricType = new medGeometryElement[1] ;
1743 const medGeometryElement * allType = getTypes(MED_FACE);
1744 geometricType[0] = allType[0] ;
1745 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
1746 geometricTypeNumber[0] = 0 ;
1747 numberOfEntities = new int[1] ;
1748 numberOfEntities[0] = size ;
1749 mySkyLineArrayIndex = new int[2] ;
1750 mySkyLineArrayIndex[0]=1 ;
1751 mySkyLineArrayIndex[1]=1+size ;
1754 map<medGeometryElement,int> theType ;
1755 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
1756 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
1757 if (theType.find(myType) != theType.end() )
1758 theType[myType]+=1 ;
1762 numberOfGeometricType = theType.size() ;
1763 geometricType = new medGeometryElement[numberOfGeometricType] ;
1764 //const medGeometryElement * allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
1765 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
1766 numberOfEntities = new int[numberOfGeometricType] ;
1767 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1769 mySkyLineArrayIndex[0]=1 ;
1770 map<medGeometryElement,int>::iterator theTypeIt ;
1771 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
1772 geometricType[index] = (*theTypeIt).first ;
1773 geometricTypeNumber[index] = 0 ;
1774 numberOfEntities[index] = (*theTypeIt).second ;
1775 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
1779 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
1780 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
1782 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
1783 mySupport->setGeometricType(geometricType) ;
1784 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
1785 mySupport->setNumberOfElements(numberOfEntities) ;
1786 mySupport->setTotalNumberOfElements(size) ;
1787 mySupport->setNumber(mySkyLineArray) ;
1789 delete[] numberOfEntities;
1790 delete[] geometricTypeNumber;
1791 delete[] geometricType;
1792 delete[] mySkyLineArrayIndex;
1793 delete[] myListArray;
1794 // delete mySkyLineArray;
1802 return a SUPPORT pointer on the union of all SUPPORTs in Supports.
1803 You should delete this pointer after use to avoid memory leaks.
1805 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
1807 const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
1810 SUPPORT * returnedSupport;
1811 string returnedSupportName;
1812 string returnedSupportDescription;
1813 char * returnedSupportNameChar;
1814 char * returnedSupportDescriptionChar;
1815 int size = Supports.size();
1819 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
1820 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1822 returnedSupport = new SUPPORT(*obj);
1824 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
1825 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
1827 returnedSupportNameChar = new char[lenName];
1828 returnedSupportDescriptionChar = new char[lenDescription];
1830 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
1831 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1832 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
1833 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1834 (Supports[0]->getDescription()).c_str());
1836 returnedSupportName = string(returnedSupportNameChar);
1837 returnedSupportDescription = string(returnedSupportDescriptionChar);
1839 returnedSupport->setName(returnedSupportName);
1840 returnedSupport->setDescription(returnedSupportDescription);
1842 delete [] returnedSupportNameChar;
1843 delete [] returnedSupportDescriptionChar;
1847 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1848 returnedSupport = new SUPPORT(*obj);
1850 int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
1851 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
1853 for(int i = 1;i<size;i++)
1855 obj = const_cast <SUPPORT *> (Supports[i]);
1856 returnedSupport->blending(obj);
1860 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
1861 lenDescription = lenDescription + 5 +
1862 strlen((Supports[i]->getDescription()).c_str());
1866 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
1867 lenDescription = lenDescription + 2 +
1868 strlen((Supports[i]->getDescription()).c_str());
1872 returnedSupportNameChar = new char[lenName];
1873 returnedSupportDescriptionChar = new char[lenDescription];
1875 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
1876 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
1878 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1879 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1880 (Supports[0]->getDescription()).c_str());
1882 for(int i = 1;i<size;i++)
1886 returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
1887 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
1889 returnedSupportNameChar = strcat(returnedSupportNameChar,
1890 (Supports[i]->getName()).c_str());
1891 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1892 (Supports[i]->getDescription()).c_str());
1896 returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
1897 returnedSupportNameChar = strcat(returnedSupportNameChar,
1898 (Supports[i]->getName()).c_str());
1900 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
1901 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1902 (Supports[i]->getDescription()).c_str());
1906 returnedSupportName = string(returnedSupportNameChar);
1907 returnedSupport->setName(returnedSupportName);
1909 returnedSupportDescription = string(returnedSupportDescriptionChar);
1910 returnedSupport->setDescription(returnedSupportDescription);
1912 delete [] returnedSupportNameChar;
1913 delete [] returnedSupportDescriptionChar;
1917 return returnedSupport;
1921 return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
1922 The (SUPPORT *) NULL pointer is returned if the intersection is empty.
1923 You should delete this pointer after use to avois memory leaks.
1925 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
1927 const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
1930 SUPPORT * returnedSupport;
1931 string returnedSupportName;
1932 string returnedSupportDescription;
1933 char * returnedSupportNameChar;
1934 char * returnedSupportDescriptionChar;
1935 int size = Supports.size();
1939 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
1940 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1942 returnedSupport = new SUPPORT(*obj);
1944 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
1945 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
1947 returnedSupportNameChar = new char[lenName];
1948 returnedSupportDescriptionChar = new char[lenDescription];
1950 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
1951 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1952 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
1953 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1954 (Supports[0]->getDescription()).c_str());
1956 returnedSupportName = string(returnedSupportNameChar);
1957 returnedSupportDescription = string(returnedSupportDescriptionChar);
1959 returnedSupport->setName(returnedSupportName);
1960 returnedSupport->setDescription(returnedSupportDescription);
1962 delete [] returnedSupportNameChar;
1963 delete [] returnedSupportDescriptionChar;
1967 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1968 returnedSupport = new SUPPORT(*obj);
1970 int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
1971 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
1973 for(int i = 1;i<size;i++)
1975 obj = const_cast <SUPPORT *> (Supports[i]);
1976 returnedSupport->intersecting(obj);
1980 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
1981 lenDescription = lenDescription + 5 +
1982 strlen((Supports[i]->getDescription()).c_str());
1986 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
1987 lenDescription = lenDescription + 2 +
1988 strlen((Supports[i]->getDescription()).c_str());
1991 if(returnedSupport != (SUPPORT *) NULL)
1993 returnedSupportNameChar = new char[lenName];
1994 returnedSupportDescriptionChar = new char[lenDescription];
1996 returnedSupportNameChar = strcpy(returnedSupportNameChar,
1997 "Intersection of ");
1998 returnedSupportDescriptionChar =
1999 strcpy(returnedSupportDescriptionChar,"Intersection of ");
2001 returnedSupportNameChar = strcat(returnedSupportNameChar,
2002 (Supports[0]->getName()).c_str());
2003 returnedSupportDescriptionChar =
2004 strcat(returnedSupportDescriptionChar,
2005 (Supports[0]->getDescription()).c_str());
2007 for(int i = 1;i<size;i++)
2011 returnedSupportNameChar = strcat(returnedSupportNameChar,
2013 returnedSupportDescriptionChar =
2014 strcat(returnedSupportDescriptionChar," and ");
2016 returnedSupportNameChar =
2017 strcat(returnedSupportNameChar,
2018 (Supports[i]->getName()).c_str());
2019 returnedSupportDescriptionChar =
2020 strcat(returnedSupportDescriptionChar,
2021 (Supports[i]->getDescription()).c_str());
2025 returnedSupportNameChar = strcat(returnedSupportNameChar,
2027 returnedSupportNameChar =
2028 strcat(returnedSupportNameChar,
2029 (Supports[i]->getName()).c_str());
2031 returnedSupportDescriptionChar =
2032 strcat(returnedSupportDescriptionChar,", ");
2033 returnedSupportDescriptionChar =
2034 strcat(returnedSupportDescriptionChar,
2035 (Supports[i]->getDescription()).c_str());
2039 returnedSupportName = string(returnedSupportNameChar);
2040 returnedSupport->setName(returnedSupportName);
2042 returnedSupportDescription = string(returnedSupportDescriptionChar);
2043 returnedSupport->setDescription(returnedSupportDescription);
2045 delete [] returnedSupportNameChar;
2046 delete [] returnedSupportDescriptionChar;
2051 return returnedSupport;
2055 // internal helper type
2058 std::vector<int> groups;
2059 MED_EN::medGeometryElement geometricType;
2062 // Create families from groups
2063 void MESH::createFamilies()
2065 int idFamNode = 0; // identifier for node families
2066 int idFamElement = 0; // identifier for cell, face or edge families
2068 // Main loop on mesh's entities
2069 for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2071 int numberofgroups = getNumberOfGroups(entity);
2073 continue; // no groups for this entity
2075 vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2077 // make myFamilies points to the member corresponding to entity
2078 vector<FAMILY*>* myFamilies;
2082 myFamilies = & _familyCell;
2085 myFamilies = & _familyFace;
2088 myFamilies = & _familyEdge;
2091 myFamilies = & _familyNode;
2095 vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2096 // get a copy of the (old) family ptrs before clearing
2097 vector<FAMILY*> myOldFamilies=getFamilies(entity);
2098 myFamilies->clear();
2101 // 1 - Create a vector containing for each cell (of the entity) an information structure
2102 // giving geometric type and the groups it belong to
2104 med_int numberOfTypes=getNumberOfTypes(entity);
2105 const int * index=getGlobalNumberingIndex(entity);
2106 const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2107 med_int numberOfCells=index[numberOfTypes]-1; // total number of cells for that entity
2108 SCRUTE(numberOfTypes);
2109 SCRUTE(numberOfCells);
2110 vector< _cell > tab_cell(numberOfCells);
2111 for(med_int t=0; t!=numberOfTypes; ++t)
2112 for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2113 tab_cell[n].geometricType=geometricTypes[t];
2116 // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2118 for (unsigned g=0; g!=myGroups.size(); ++g)
2120 // scan cells that belongs to the group
2121 const int* groupCells=myGroups[g]->getnumber()->getValue();
2122 int nbCells=myGroups[g]->getnumber()->getLength();
2123 for(int c=0; c!=nbCells; ++c)
2124 tab_cell[groupCells[c]-1].groups.push_back(g);
2128 // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2129 // whith the vector of contained cells
2131 map< string,vector<int> > tab_families;
2132 map< string,vector<int> >::iterator fam;
2133 for(int n=0; n!=numberOfCells; ++n)
2135 ostringstream key; // to generate the name of the family
2137 if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2138 key << "_NONE" << entity;
2140 for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2142 string groupName=myGroups[*it]->getName();
2143 if(groupName.empty())
2146 key << "_" << groupName;
2149 tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2153 // 4 - Scan the family map, create MED Families, check if it already exist.
2155 for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2157 vector<medGeometryElement> tab_types_geometriques;
2158 medGeometryElement geometrictype=MED_NONE;
2159 vector<int> tab_index_types_geometriques;
2160 vector<int> tab_nombres_elements;
2162 // scan family cells and fill the tab that are needed by the create a MED FAMILY
2163 for( int i=0; i!=fam->second.size(); ++i)
2165 int ncell=fam->second[i]-1;
2166 if(tab_cell[ncell].geometricType != geometrictype)
2168 // new geometric type -> we store it and complete index tabs
2169 if(!tab_index_types_geometriques.empty())
2170 tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2171 tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2172 tab_index_types_geometriques.push_back(i+1);
2175 // store and complete index tabs for the last geometric type
2176 tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2177 tab_index_types_geometriques.push_back(fam->second.size()+1);
2179 // create a empty MED FAMILY and fill it with the tabs we constructed
2180 FAMILY* newFam = new FAMILY();
2181 newFam->setTotalNumberOfElements(fam->second.size());
2182 newFam->setName(fam->first);
2183 newFam->setMesh(this);
2184 newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2185 newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2186 newFam->setNumberOfElements(&tab_nombres_elements[0]);
2187 newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2188 newFam->setEntity(entity);
2189 newFam->setAll(false);
2201 idFam = -idFamElement;
2205 idFam = -idFamElement;
2209 idFam = -idFamElement;
2213 newFam->setIdentifier(idFam);
2215 // Update links between families and groups
2217 int ncell1=fam->second[0]-1; // number of first cell in family
2218 int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2221 newFam->setNumberOfGroups(numberOfGroups);
2222 string * groupNames=new string[numberOfGroups];
2224 // iterate on the groups the family belongs to
2225 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2226 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2228 whichFamilyInGroup[*it].push_back(newFam);
2229 groupNames[ng]=myGroups[*it]->getName();
2231 newFam->setGroupsNames(groupNames);
2234 int sizeOfFamVect = myFamilies->size();
2236 MESSAGE(" MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2238 myFamilies->push_back(newFam);
2241 // delete old families
2242 for (int i=0;i<myOldFamilies.size();i++)
2243 delete myOldFamilies[i] ;
2245 // update references in groups
2246 for (int i=0;i<myGroups.size();i++)
2248 myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2249 myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2252 // re-scan the cells vector, and fill the family vector with cells.
2253 // creation of support, check if it already exist.
2257 int MESH::getElementContainingPoint(const double *coord)
2259 if(_spaceDimension==3)
2261 Meta_Wrapper<3> *fromWrapper=new Meta_Wrapper<3> (getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2262 const_cast<CONNECTIVITY *>(getConnectivityptr()));
2263 Meta_Wrapper<3> *toWrapper=new Meta_Wrapper<3> (1,const_cast<double *>(coord));
2264 Meta_Mapping<3> *mapping=new Meta_Mapping<3> (fromWrapper,toWrapper);
2265 mapping->Cree_Mapping(1);
2266 vector<int> vectormapping=mapping->Get_Mapping();
2267 return vectormapping[0]+1;
2269 else if(_spaceDimension==2)
2271 Meta_Wrapper<2> *fromWrapper=new Meta_Wrapper<2> (getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2272 const_cast<CONNECTIVITY *>(getConnectivityptr()));
2273 Meta_Wrapper<2> *toWrapper=new Meta_Wrapper<2> (1,const_cast<double *>(coord));
2274 Meta_Mapping<2> *mapping=new Meta_Mapping<2> (fromWrapper,toWrapper);
2275 mapping->Cree_Mapping(1);
2276 vector<int> vectormapping=mapping->Get_Mapping();
2277 return vectormapping[0]+1;
2280 throw MEDEXCEPTION("MESH::getElementContainingPoint : invalid _spaceDimension must be equal to 2 or 3 !!!");
2283 //Presently disconnected in C++
2284 void MESH::addReference() const
2288 //Presently disconnected in C++
2289 void MESH::removeReference() const