10 #include "MEDMEM_DriversDef.hxx"
11 #include "MEDMEM_Field.hxx"
12 #include "MEDMEM_Mesh.hxx"
14 #include "MEDMEM_Support.hxx"
15 #include "MEDMEM_Family.hxx"
16 #include "MEDMEM_Group.hxx"
17 #include "MEDMEM_Coordinate.hxx"
18 #include "MEDMEM_Connectivity.hxx"
19 #include "MEDMEM_CellModel.hxx"
20 #include "MEDMEM_Formulae.hxx"
21 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
22 #include "MEDMEM_DriverFactory.hxx"
25 using namespace MEDMEM;
26 using namespace MED_EN;
31 // ------- Drivers Management Part
33 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
34 is driverName. addDriver returns an integer handler. */
35 int MESH::addDriver(driverTypes driverType,
36 const string & fileName/*="Default File Name.med"*/,
37 const string & driverName/*="Default Mesh Name"*/,
38 med_mode_acces access) {
40 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) : ";
48 driver = DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,
51 _drivers.push_back(driver);
53 int current = _drivers.size()-1;
55 _drivers[current]->setMeshName(driverName);
62 /*! Add an existing MESH driver. */
63 int MESH::addDriver(GENDRIVER & driver) {
64 const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
67 // A faire : Vérifier que le driver est de type MESH.
68 GENDRIVER * newDriver = driver.copy() ;
70 _drivers.push_back(newDriver);
71 return _drivers.size()-1;
76 /*! Remove an existing MESH driver. */
77 void MESH::rmDriver (int index/*=0*/) {
78 const char * LOC = "MESH::rmDriver (int index=0): ";
81 if ( _drivers[index] ) {
82 //_drivers.erase(&_drivers[index]);
87 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
88 << "The index given is invalid, index must be between 0 and |"
97 // ------ End of Drivers Management Part
102 const char * LOC = "MESH::init(): ";
106 _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
108 _coordinate = (COORDINATE *) NULL;
109 _connectivity = (CONNECTIVITY *) NULL;
111 _spaceDimension = MED_INVALID; // 0 ?!?
112 _meshDimension = MED_INVALID;
113 _numberOfNodes = MED_INVALID;
117 _arePresentOptionnalNodesNumbers = 0;
122 /*! Create an empty MESH. */
123 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
130 _isAGrid = m._isAGrid;
132 if (m._coordinate != NULL)
133 _coordinate = new COORDINATE(* m._coordinate);
135 _coordinate = (COORDINATE *) NULL;
137 if (m._connectivity != NULL)
138 _connectivity = new CONNECTIVITY(* m._connectivity);
140 _connectivity = (CONNECTIVITY *) NULL;
142 _spaceDimension = m._spaceDimension;
143 _meshDimension = m._meshDimension;
144 _numberOfNodes = m._numberOfNodes;
146 _familyNode = m._familyNode;
147 for (int i=0; i<m._familyNode.size(); i++)
149 _familyNode[i] = new FAMILY(* m._familyNode[i]);
150 _familyNode[i]->setMesh(this);
153 _familyCell = m._familyCell;
154 for (int i=0; i<m._familyCell.size(); i++)
156 _familyCell[i] = new FAMILY(* m._familyCell[i]);
157 _familyCell[i]->setMesh(this);
160 _familyFace = m._familyFace;
161 for (int i=0; i<m._familyFace.size(); i++)
163 _familyFace[i] = new FAMILY(* m._familyFace[i]);
164 _familyFace[i]->setMesh(this);
167 _familyEdge = m._familyEdge;
168 for (int i=0; i<m._familyEdge.size(); i++)
170 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
171 _familyEdge[i]->setMesh(this);
174 _groupNode = m._groupNode;
175 for (int i=0; i<m._groupNode.size(); i++)
177 _groupNode[i] = new GROUP(* m._groupNode[i]);
178 _groupNode[i]->setMesh(this);
181 _groupCell = m._groupCell;
182 for (int i=0; i<m._groupCell.size(); i++)
184 _groupCell[i] = new GROUP(* m._groupCell[i]);
185 _groupCell[i]->setMesh(this);
188 _groupFace = m._groupFace;
189 for (int i=0; i<m._groupFace.size(); i++)
191 _groupFace[i] = new GROUP(* m._groupFace[i]);
192 _groupFace[i]->setMesh(this);
195 _groupEdge = m._groupEdge;
196 for (int i=0; i<m._groupEdge.size(); i++)
198 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
199 _groupEdge[i]->setMesh(this);
202 //_drivers = m._drivers; //Recopie des drivers?
207 MESSAGE("MESH::~MESH() : Destroying the Mesh");
208 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
209 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
211 size = _familyNode.size() ;
212 for (int i=0;i<size;i++)
213 delete _familyNode[i] ;
214 size = _familyCell.size() ;
215 for (int i=0;i<size;i++)
216 delete _familyCell[i] ;
217 size = _familyFace.size() ;
218 for (int i=0;i<size;i++)
219 delete _familyFace[i] ;
220 size = _familyEdge.size() ;
221 for (int i=0;i<size;i++)
222 delete _familyEdge[i] ;
223 size = _groupNode.size() ;
224 for (int i=0;i<size;i++)
225 delete _groupNode[i] ;
226 size = _groupCell.size() ;
227 for (int i=0;i<size;i++)
228 delete _groupCell[i] ;
229 size = _groupFace.size() ;
230 for (int i=0;i<size;i++)
231 delete _groupFace[i] ;
232 size = _groupEdge.size() ;
233 for (int i=0;i<size;i++)
234 delete _groupEdge[i] ;
236 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
238 for (unsigned int index=0; index < _drivers.size(); index++ )
240 SCRUTE(_drivers[index]);
241 if ( _drivers[index] != NULL) delete _drivers[index];
247 Method equivalent to getNumberOfTypes except that it includes not only classical Types but polygons/polyhedra also.
249 int MESH::getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const
251 if(_connectivity!= NULL)
252 return _connectivity->getNumberOfTypesWithPoly(Entity);
253 throw MEDEXCEPTION(LOCALIZED("MESH::getNumberOfTypesWithPoly( medEntityMesh ) : Connectivity not defined !"));
257 Method equivalent to getTypesWithPoly except that it includes not only classical Types but polygons/polyhedra also.
258 WARNING the returned array MUST be deallocated.
260 MED_EN::medGeometryElement * MESH::getTypesWithPoly(MED_EN::medEntityMesh Entity) const
262 if (Entity == MED_EN::MED_NODE)
263 throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : No medGeometryElement with MED_NODE entity !"));
264 if (_connectivity != NULL)
265 return _connectivity->getGeometricTypesWithPoly(Entity);
266 throw MEDEXCEPTION(LOCALIZED("MESH::getTypes( medEntityMesh ) : Connectivity not defined !"));
270 Method equivalent to getNumberOfElementsWithPoly except that it includes not only classical Types but polygons/polyhedra also.
272 int MESH::getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type) const
274 if(Type==MED_POLYGON || Type==MED_POLYHEDRA)
276 int nbOfPolygs=_connectivity->getNumberOfElementOfPolyType(Entity);
279 else if(Type==MED_ALL_ELEMENTS)
281 int nbOfClassicalTypes=getNumberOfElements(Entity,MED_ALL_ELEMENTS);
282 int nbOfClassicalTypes2=_connectivity->getNumberOfElementOfPolyType(Entity);
283 return nbOfClassicalTypes+nbOfClassicalTypes2;
286 return getNumberOfElements(Entity,Type);
289 MESH & MESH::operator=(const MESH &m)
291 const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
294 MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
297 // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
298 // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
299 // _drivers = m._drivers;
301 // space_dimension=m.space_dimension;
302 // mesh_dimension=m.mesh_dimension;
304 // nodes_count=m.nodes_count;
305 // coordinates=m.coordinates;
306 // coordinates_name=m.coordinates_name ;
307 // coordinates_unit=m.coordinates_unit ;
308 // nodes_numbers=m.nodes_numbers ;
310 // cells_types_count=m.cells_types_count;
311 // cells_type=m.cells_type;
312 // cells_count=m.cells_count;
313 // nodal_connectivity=m.nodal_connectivity;
315 // nodes_families_count=m.nodes_families_count;
316 // nodes_Families=m.nodes_Families;
318 // cells_families_count=m.cells_families_count;
319 // cells_Families=m.cells_Families;
321 // maximum_cell_number_by_node = m.maximum_cell_number_by_node;
322 // if (maximum_cell_number_by_node > 0)
324 // reverse_nodal_connectivity = m.reverse_nodal_connectivity;
325 // reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
332 bool MESH::operator==(const MESH& other) const
334 BEGIN_OF("MESH::operator==");
338 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
339 The meshname driverName must already exists in the file.*/
340 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) throw (MEDEXCEPTION)
342 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
348 GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,MED_LECT);
349 current = addDriver(*myDriver);
351 _drivers[current]->open();
352 _drivers[current]->read();
353 _drivers[current]->close();
359 for a deep comparison of 2 meshes.
361 bool MESH::deepCompare(const MESH& other) const
363 int size1=getSpaceDimension()*getNumberOfNodes();
364 int size2=other.getSpaceDimension()*other.getNumberOfNodes();
367 const double* coord1=getCoordinates(MED_FULL_INTERLACE);
368 const double* coord2=other.getCoordinates(MED_FULL_INTERLACE);
370 for(int i=0;i<size1 && ret;i++)
372 ret=(fabs(coord1[i]-coord2[i])<1e-15);
376 return _connectivity->deepCompare(*other._connectivity);
381 ostream & ::MEDMEM::operator<<(ostream &os, const MESH &myMesh)
383 int spacedimension = myMesh.getSpaceDimension();
384 int meshdimension = myMesh.getMeshDimension();
385 int numberofnodes = myMesh.getNumberOfNodes();
387 os << "Space Dimension : " << spacedimension << endl << endl;
389 os << "Mesh Dimension : " << meshdimension << endl << endl;
391 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
392 os << "SHOW NODES COORDINATES : " << endl;
394 os << "Name :" << endl;
395 const string * coordinatesnames = myMesh.getCoordinatesNames();
396 for(int i=0; i<spacedimension ; i++)
398 os << " - " << coordinatesnames[i] << endl;
400 os << "Unit :" << endl;
401 const string * coordinatesunits = myMesh.getCoordinatesUnits();
402 for(int i=0; i<spacedimension ; i++)
404 os << " - " << coordinatesunits[i] << endl;
406 for(int i=0; i<numberofnodes ; i++)
408 os << "Nodes " << i+1 << " : ";
409 for (int j=0; j<spacedimension ; j++)
410 os << coordinates[i*spacedimension+j] << " ";
414 os << endl << "SHOW CONNECTIVITY :" << endl;
415 os << *myMesh._connectivity << endl;
417 medEntityMesh entity;
418 os << endl << "SHOW FAMILIES :" << endl << endl;
419 for (int k=1; k<=4; k++)
421 if (k==1) entity = MED_NODE;
422 if (k==2) entity = MED_CELL;
423 if (k==3) entity = MED_FACE;
424 if (k==4) entity = MED_EDGE;
425 int numberoffamilies = myMesh.getNumberOfFamilies(entity);
426 os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
427 using namespace MED_EN;
428 for (int i=1; i<numberoffamilies+1;i++)
430 os << * myMesh.getFamily(entity,i) << endl;
434 os << endl << "SHOW GROUPS :" << endl << endl;
435 for (int k=1; k<=4; k++)
437 if (k==1) entity = MED_NODE;
438 if (k==2) entity = MED_CELL;
439 if (k==3) entity = MED_FACE;
440 if (k==4) entity = MED_EDGE;
441 int numberofgroups = myMesh.getNumberOfGroups(entity);
442 os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
443 using namespace MED_EN;
444 for (int i=1; i<numberofgroups+1;i++)
446 os << * myMesh.getGroup(entity,i) << endl;
454 Get global number of element which have same connectivity than connectivity argument.
456 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
458 Return -1 if not found.
460 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
462 const char* LOC="MESH::getElementNumber " ;
465 int numberOfValue ; // size of connectivity array
466 CELLMODEL myType(Type) ;
467 if (ConnectivityType==MED_DESCENDING)
468 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
470 numberOfValue = myType.getNumberOfNodes() ; // nodes
472 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
473 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
475 // First node or face/edge
476 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
477 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
479 // list to put cells numbers
480 list<int> cellsList ;
481 list<int>::iterator itList ;
482 for (int i=indexBegin; i<indexEnd; i++)
483 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
485 // Foreach node or face/edge in cell, we search cells which are in common.
486 // At the end, we must have only one !
488 for (int i=1; i<numberOfValue; i++) {
489 int connectivity_i = connectivity[i] ;
490 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
492 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
493 if ((*itList)==myReverseConnectivityValue[j-1]) {
499 itList=cellsList.erase(itList);
500 itList--; // well : rigth if itList = cellsList.begin() ??
505 if (cellsList.size()>1) // we have more than one cell
506 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
508 if (cellsList.size()==0)
513 return cellsList.front() ;
517 Return a support which reference all elements on the boundary of mesh.
519 For instance, we could get only face in 3D and edge in 2D.
521 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
524 const char * LOC = "MESH::getBoundaryElements : " ;
527 // actually we could only get face (in 3D) and edge (in 2D)
528 medEntityMesh entityToParse=Entity;
529 if(_spaceDimension == 3)
530 if (Entity != MED_FACE)
532 entityToParse=MED_FACE;
534 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
535 if(_spaceDimension == 2)
536 if(Entity != MED_EDGE)
538 entityToParse=MED_EDGE;
540 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
542 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
543 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
544 int numberOf = getNumberOfElementsWithPoly(entityToParse,MED_ALL_ELEMENTS) ;
545 list<int> myElementsList;
547 for (int i=0 ; i<numberOf; i++)
548 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
549 myElementsList.push_back(i+1);
553 return buildSupportOnNodeFromElementList(myElementsList,entityToParse);
555 return buildSupportOnElementsFromElementList(myElementsList,entityToParse);
559 Method that do the same thing as buildSupportOnNodeFromElementList except that a SUPPORT is not created.
561 void MESH::fillSupportOnNodeFromElementList(const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION)
563 MED_EN::medEntityMesh entity=supportToFill->getEntity();
564 supportToFill->setAll(false);
565 supportToFill->setMesh((MESH *)this);
569 for(list<int>::const_iterator iter=listOfElt.begin();iter!=listOfElt.end();iter++)
572 const int *conn=_connectivity->getConnectivityOfAnElementWithPoly(MED_NODAL,entity,*iter,lgth);
574 nodes.insert(conn[i]);
577 for(set<int>::iterator iter2=nodes.begin();iter2!=nodes.end();iter2++)
578 nodesList.push_back(*iter2);
579 supportToFill->fillFromNodeList(nodesList);
583 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
586 SUPPORT *MESH::buildSupportOnNodeFromElementList(const list<int>& listOfElt,MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
588 SUPPORT * mySupport = new SUPPORT((MESH *)this,"Boundary",entity);
589 fillSupportOnNodeFromElementList(listOfElt,mySupport);
594 Method created to factorize code. This method creates a new support on entity 'entity' (to deallocate) containing all the entities contained in
595 elements 'listOfElt' of entity 'entity'.
597 SUPPORT *MESH::buildSupportOnElementsFromElementList(const list<int>& listOfElt, MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION)
599 const char * LOC = "MESH::buildSupportOnElementsFromElementList : " ;
601 SUPPORT *mySupport=new SUPPORT((MESH *)this,"Boundary",entity);
602 mySupport->fillFromElementList(listOfElt);
607 FIELD<double>* MESH::getVolume(const SUPPORT *Support) const throw (MEDEXCEPTION)
609 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
612 // Support must be on 3D elements
614 // Make sure that the MESH class is the same as the MESH class attribut
615 // in the class Support
616 MESH* myMesh = Support->getMesh();
618 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
620 int dim_space = getSpaceDimension();
621 medEntityMesh support_entity = Support->getEntity();
622 bool onAll = Support->isOnAllElements();
624 int nb_type, length_values;
625 const medGeometryElement* types;
627 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
628 const int* global_connectivity;
629 nb_type = Support->getNumberOfTypes();
630 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
631 types = Support->getTypes();
633 FIELD<double>* Volume = new FIELD<double>(Support,1);
634 // double *volume = new double [length_values];
635 Volume->setName("VOLUME");
636 Volume->setDescription("cells volume");
637 Volume->setComponentName(1,"volume");
638 Volume->setComponentDescription(1,"desc-comp");
640 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
642 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
644 Volume->setMEDComponentUnit(1,MEDComponentUnit);
646 Volume->setValueType(MED_REEL64);
648 Volume->setIterationNumber(0);
649 Volume->setOrderNumber(0);
650 Volume->setTime(0.0);
651 MEDARRAY<double> *volume = Volume->getvalue();
653 const double * coord = getCoordinates(MED_FULL_INTERLACE);
655 for (int i=0;i<nb_type;i++)
657 medGeometryElement type = types[i] ;
659 nb_entity_type = Support->getNumberOfElements(type);
660 if(type != MED_EN::MED_POLYHEDRA)
664 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
668 const int * supp_number = Support->getNumber(type);
669 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
670 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
671 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
673 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
674 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
675 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
678 global_connectivity = global_connectivity_tmp ;
684 case MED_TETRA4 : case MED_TETRA10 :
686 for (int tetra=0;tetra<nb_entity_type;tetra++)
688 int tetra_index = (type%100)*tetra;
689 int N1 = global_connectivity[tetra_index]-1;
690 int N2 = global_connectivity[tetra_index+1]-1;
691 int N3 = global_connectivity[tetra_index+2]-1;
692 int N4 = global_connectivity[tetra_index+3]-1;
693 xvolume=CalculateVolumeForTetra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4);
694 volume->setIJ(index,1,xvolume) ;
699 case MED_PYRA5 : case MED_PYRA13 :
701 for (int pyra=0;pyra<nb_entity_type;pyra++)
703 int pyra_index = (type%100)*pyra;
704 int N1 = global_connectivity[pyra_index]-1;
705 int N2 = global_connectivity[pyra_index+1]-1;
706 int N3 = global_connectivity[pyra_index+2]-1;
707 int N4 = global_connectivity[pyra_index+3]-1;
708 int N5 = global_connectivity[pyra_index+4]-1;
709 xvolume=CalculateVolumeForPyra(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,coord+dim_space*N5);
710 volume->setIJ(index,1,xvolume) ;
715 case MED_PENTA6 : case MED_PENTA15 :
717 for (int penta=0;penta<nb_entity_type;penta++)
719 int penta_index = (type%100)*penta;
720 int N1 = global_connectivity[penta_index]-1;
721 int N2 = global_connectivity[penta_index+1]-1;
722 int N3 = global_connectivity[penta_index+2]-1;
723 int N4 = global_connectivity[penta_index+3]-1;
724 int N5 = global_connectivity[penta_index+4]-1;
725 int N6 = global_connectivity[penta_index+5]-1;
726 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);
727 volume->setIJ(index,1,xvolume) ;
732 case MED_HEXA8 : case MED_HEXA20 :
734 for (int hexa=0;hexa<nb_entity_type;hexa++)
736 int hexa_index = (type%100)*hexa;
738 int N1 = global_connectivity[hexa_index]-1;
739 int N2 = global_connectivity[hexa_index+1]-1;
740 int N3 = global_connectivity[hexa_index+2]-1;
741 int N4 = global_connectivity[hexa_index+3]-1;
742 int N5 = global_connectivity[hexa_index+4]-1;
743 int N6 = global_connectivity[hexa_index+5]-1;
744 int N7 = global_connectivity[hexa_index+6]-1;
745 int N8 = global_connectivity[hexa_index+7]-1;
746 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);
747 volume->setIJ(index,1,xvolume) ;
757 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
759 int lgthNodes,iPts,iFaces,iPtsInFace;
760 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
761 int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
762 int nbOfFaces,*nbOfNodesPerFaces;
763 int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(offsetWithClassicType+polyhs+1,nbOfFaces,nbOfNodesPerFaces);
764 double **pts=new double * [lgthNodes];
765 double ***pts1=new double ** [nbOfFaces];
766 for(iPts=0;iPts<lgthNodes;iPts++)
767 pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
768 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
770 pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
771 for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
772 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
775 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
778 xvolume=CalculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
779 delete [] nbOfNodesPerFaces;
780 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
781 delete [] pts1[iFaces];
783 volume->setIJ(index,1,xvolume) ;
789 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
790 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
792 int lgthNodes,iPts,iFaces,iPtsInFace;
793 int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
794 int nbOfFaces,*nbOfNodesPerFaces;
795 int **nodes1=_connectivity->getNodesPerFaceOfPolyhedron(supp_number[polyhs],nbOfFaces,nbOfNodesPerFaces);
796 double **pts=new double * [lgthNodes];
797 double ***pts1=new double ** [nbOfFaces];
798 for(iPts=0;iPts<lgthNodes;iPts++)
799 pts[iPts]=(double *)(coord+3*(nodes[iPts]-1));
800 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
802 pts1[iFaces]=new double* [nbOfNodesPerFaces[iFaces]];
803 for(iPtsInFace=0;iPtsInFace<nbOfNodesPerFaces[iFaces];iPtsInFace++)
804 pts1[iFaces][iPtsInFace]=(double *)(coord+3*(nodes1[iFaces][iPtsInFace]-1));
807 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,bary);
810 xvolume=CalculateVolumeForPolyh((const double ***)pts1,nbOfNodesPerFaces,nbOfFaces,bary);
811 delete [] nbOfNodesPerFaces;
812 for(iFaces=0;iFaces<nbOfFaces;iFaces++)
813 delete [] pts1[iFaces];
815 volume->setIJ(index,1,xvolume) ;
822 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
826 if (!onAll && type!=MED_EN::MED_POLYHEDRA)
827 delete [] global_connectivity ;
833 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
835 const char * LOC = "MESH::getArea(SUPPORT*) : ";
838 // Support must be on 2D elements
840 // Make sure that the MESH class is the same as the MESH class attribut
841 // in the class Support
842 MESH* myMesh = Support->getMesh();
844 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
846 int dim_space = getSpaceDimension();
847 medEntityMesh support_entity = Support->getEntity();
848 bool onAll = Support->isOnAllElements();
850 int nb_type, length_values;
851 const medGeometryElement* types;
853 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
854 const int* global_connectivity;
856 nb_type = Support->getNumberOfTypes();
857 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
858 types = Support->getTypes();
863 Area = new FIELD<double>(Support,1);
864 Area->setName("AREA");
865 Area->setDescription("cells or faces area");
867 Area->setComponentName(1,"area");
868 Area->setComponentDescription(1,"desc-comp");
870 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
872 string MEDComponentUnit="trucmuch";
874 Area->setMEDComponentUnit(1,MEDComponentUnit);
876 Area->setValueType(MED_REEL64);
878 Area->setIterationNumber(0);
879 Area->setOrderNumber(0);
882 double *area = (double *)Area->getValue(MED_FULL_INTERLACE);
884 const double * coord = getCoordinates(MED_FULL_INTERLACE);
887 for (int i=0;i<nb_type;i++)
889 medGeometryElement type = types[i] ;
890 nb_entity_type = Support->getNumberOfElements(type);
891 const int *global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
892 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
896 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
900 const int * supp_number = Support->getNumber(type);
901 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
902 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
904 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
905 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
906 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
909 global_connectivity = global_connectivity_tmp ;
914 case MED_TRIA3 : case MED_TRIA6 :
916 for (int tria=0;tria<nb_entity_type;tria++)
918 int tria_index = (type%100)*tria;
920 int N1 = global_connectivity[tria_index]-1;
921 int N2 = global_connectivity[tria_index+1]-1;
922 int N3 = global_connectivity[tria_index+2]-1;
924 area[index]=CalculateAreaForTria(coord+(dim_space*N1),
925 coord+(dim_space*N2),
926 coord+(dim_space*N3),dim_space);
931 case MED_QUAD4 : case MED_QUAD8 :
933 for (int quad=0;quad<nb_entity_type;quad++)
935 int quad_index = (type%100)*quad;
937 int N1 = global_connectivity[quad_index]-1;
938 int N2 = global_connectivity[quad_index+1]-1;
939 int N3 = global_connectivity[quad_index+2]-1;
940 int N4 = global_connectivity[quad_index+3]-1;
942 area[index]=CalculateAreaForQuad(coord+dim_space*N1,
945 coord+dim_space*N4,dim_space);
954 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
955 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
956 for (int polygs=0;polygs<nb_entity_type;polygs++)
958 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
959 double **pts=new double * [size];
960 for(int iPts=0;iPts<size;iPts++)
961 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1));
962 area[index] = CalculateAreaForPolyg((const double **)pts,size,dim_space);
969 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
970 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
971 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
972 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
973 for (int polygs=0;polygs<nb_entity_type;polygs++)
975 int size=connectivity_index[supp_number[polygs]-offsetWithClassicType]-connectivity_index[supp_number[polygs]-offsetWithClassicType-1];
976 double **pts=new double * [size];
977 for(int iPts=0;iPts<size;iPts++)
978 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[supp_number[polygs]-offsetWithClassicType-1]+iPts-1]-1));
979 area[index]=CalculateAreaForPolyg((const double **)pts,size,dim_space);
987 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
992 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
993 delete [] global_connectivity ;
998 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1000 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1003 // Support must be on 1D elements
1005 // Make sure that the MESH class is the same as the MESH class attribut
1006 // in the class Support
1007 MESH* myMesh = Support->getMesh();
1009 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1011 int dim_space = getSpaceDimension();
1012 medEntityMesh support_entity = Support->getEntity();
1013 bool onAll = Support->isOnAllElements();
1015 int nb_type, length_values;
1016 const medGeometryElement* types;
1018 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1019 const int* global_connectivity;
1021 nb_type = Support->getNumberOfTypes();
1022 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1023 types = Support->getTypes();
1026 FIELD<double>* Length;
1028 Length = new FIELD<double>(Support,1);
1029 // double *length = new double [length_values];
1030 Length->setValueType(MED_REEL64);
1032 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1033 MEDARRAY<double> * length = Length->getvalue();
1035 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1038 for (int i=0;i<nb_type;i++)
1040 medGeometryElement type = types[i] ;
1045 nb_entity_type = getNumberOfElements(support_entity,type);
1046 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1050 nb_entity_type = Support->getNumberOfElements(type);
1051 const int * supp_number = Support->getNumber(type);
1052 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1053 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1054 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1056 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1057 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1058 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1061 global_connectivity = global_connectivity_tmp ;
1066 case MED_SEG2 : case MED_SEG3 :
1068 for (int edge=0;edge<nb_entity_type;edge++)
1070 int edge_index = (type%100)*edge;
1072 int N1 = global_connectivity[edge_index]-1;
1073 int N2 = global_connectivity[edge_index+1]-1;
1075 double x1 = coord[dim_space*N1];
1076 double x2 = coord[dim_space*N2];
1078 double y1 = coord[dim_space*N1+1];
1079 double y2 = coord[dim_space*N2+1];
1081 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1082 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1084 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1085 (z1 - z2)*(z1 - z2));
1087 length->setIJ(index,1,xlength) ;
1093 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1097 if (!onAll) delete [] global_connectivity ;
1103 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1105 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1108 // Support must be on 2D or 1D elements
1110 // Make sure that the MESH class is the same as the MESH class attribut
1111 // in the class Support
1112 MESH* myMesh = Support->getMesh();
1114 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1116 int dim_space = getSpaceDimension();
1117 int mesh_dim=getMeshDimension();
1118 medEntityMesh support_entity = Support->getEntity();
1119 bool onAll = Support->isOnAllElements();
1121 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 ))
1122 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"incompatible mesh dimension and entity"));
1123 int nb_type, length_values;
1124 const medGeometryElement* types;
1126 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1127 const int* global_connectivity;
1129 nb_type = Support->getNumberOfTypes();
1130 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1131 types = Support->getTypes();
1135 FIELD<double>* Normal = new FIELD<double>(Support,dim_space);
1136 Normal->setName("NORMAL");
1137 Normal->setDescription("cells or faces normal");
1138 for (int k=1;k<=dim_space;k++) {
1139 Normal->setComponentName(k,"normal");
1140 Normal->setComponentDescription(k,"desc-comp");
1141 Normal->setMEDComponentUnit(k,"unit");
1144 Normal->setValueType(MED_REEL64);
1146 Normal->setIterationNumber(MED_NOPDT);
1147 Normal->setOrderNumber(MED_NONOR);
1148 Normal->setTime(0.0);
1149 double *normal = (double *)Normal->getValue(MED_FULL_INTERLACE);
1151 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1154 for (int i=0;i<nb_type;i++)
1156 medGeometryElement type = types[i] ;
1157 nb_entity_type = Support->getNumberOfElements(type);
1159 // Make sure we trying to get Normals on
1160 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1161 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1163 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1164 (type==MED_QUAD4) || (type==MED_QUAD8) || (type==MED_POLYGON)) &&
1165 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1167 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1168 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1172 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1176 const int * supp_number = Support->getNumber(type);
1177 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1178 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1179 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1181 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1182 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1183 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1187 global_connectivity = global_connectivity_tmp ;
1193 case MED_TRIA3 : case MED_TRIA6 :
1195 for (int tria=0;tria<nb_entity_type;tria++)
1197 int tria_index = (type%100)*tria;
1198 int N1 = global_connectivity[tria_index]-1;
1199 int N2 = global_connectivity[tria_index+1]-1;
1200 int N3 = global_connectivity[tria_index+2]-1;
1201 CalculateNormalForTria(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,normal+3*index);
1206 case MED_QUAD4 : case MED_QUAD8 :
1208 for (int quad=0;quad<nb_entity_type;quad++)
1210 int quad_index = (type%100)*quad;
1211 int N1 = global_connectivity[quad_index]-1;
1212 int N2 = global_connectivity[quad_index+1]-1;
1213 int N3 = global_connectivity[quad_index+2]-1;
1214 int N4 = global_connectivity[quad_index+3]-1;
1215 CalculateNormalForQuad(coord+dim_space*N1,coord+dim_space*N2,coord+dim_space*N3,coord+dim_space*N4,normal+3*index);
1220 case MED_SEG2 : case MED_SEG3 :
1222 double xnormal1, xnormal2;
1223 for (int edge=0;edge<nb_entity_type;edge++)
1225 int edge_index = (type%100)*edge;
1227 int N1 = global_connectivity[edge_index]-1;
1228 int N2 = global_connectivity[edge_index+1]-1;
1230 double x1 = coord[dim_space*N1];
1231 double x2 = coord[dim_space*N2];
1233 double y1 = coord[dim_space*N1+1];
1234 double y2 = coord[dim_space*N2+1];
1236 xnormal1 = -(y2-y1);
1239 normal[2*index] = xnormal1 ;
1240 normal[2*index+1] = xnormal2 ;
1250 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1251 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1252 for (int polygs=0;polygs<nb_entity_type;polygs++)
1254 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1255 double **pts=new double * [size];
1256 for(int iPts=0;iPts<size;iPts++)
1257 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1])-1);
1258 CalculateNormalForPolyg((const double **)pts,size,normal+3*index);
1265 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1266 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1267 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1268 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1269 for (int polygs=0;polygs<nb_entity_type;polygs++)
1271 int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1272 int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1273 double **pts=new double * [size];
1274 for(int iPts=0;iPts<size;iPts++)
1275 pts[iPts]=(double *)(coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1])-1);
1276 CalculateNormalForPolyg((const double **)pts,size,normal+3*index);
1284 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1287 if (!onAll && type!=MED_EN::MED_POLYGON)
1288 delete [] global_connectivity ;
1295 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1297 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1298 MESH* myMesh = Support->getMesh();
1300 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1302 int dim_space = getSpaceDimension();
1303 medEntityMesh support_entity = Support->getEntity();
1304 bool onAll = Support->isOnAllElements();
1306 int nb_type, length_values;
1307 const medGeometryElement* types;
1309 const int* global_connectivity;
1310 const int * global_connectivityIndex;
1312 nb_type = Support->getNumberOfTypes();
1313 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1314 types = Support->getTypes();
1316 FIELD<double>* Barycenter;
1317 Barycenter = new FIELD<double>(Support,dim_space);
1318 Barycenter->setName("BARYCENTER");
1319 Barycenter->setDescription("cells or faces barycenter");
1321 for (int k=0;k<dim_space;k++) {
1323 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1324 Barycenter->setComponentDescription(kp1,"desc-comp");
1325 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1327 Barycenter->setValueType(MED_REEL64);
1328 Barycenter->setIterationNumber(0);
1329 Barycenter->setOrderNumber(0);
1330 Barycenter->setTime(0.0);
1331 double *barycenter=(double *)Barycenter->getValue(MED_FULL_INTERLACE);
1332 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1334 for (int i=0;i<nb_type;i++)
1336 medGeometryElement type = types[i] ;
1337 nb_entity_type = Support->getNumberOfElements(type);
1338 global_connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1339 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA )
1343 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1347 const int * supp_number = Support->getNumber(type);
1348 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1349 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1351 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1352 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1353 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[global_connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1356 global_connectivity = global_connectivity_tmp;
1362 case MED_TETRA4 : case MED_TETRA10 :
1364 for (int tetra =0;tetra<nb_entity_type;tetra++)
1366 int tetra_index = (type%100)*tetra;
1368 int N1 = global_connectivity[tetra_index]-1;
1369 int N2 = global_connectivity[tetra_index+1]-1;
1370 int N3 = global_connectivity[tetra_index+2]-1;
1371 int N4 = global_connectivity[tetra_index+3]-1;
1373 pts[0]=(double *)coord+dim_space*N1;
1374 pts[1]=(double *)coord+dim_space*N2;
1375 pts[2]=(double *)coord+dim_space*N3;
1376 pts[3]=(double *)coord+dim_space*N4;
1377 CalculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1382 case MED_PYRA5 : case MED_PYRA13 :
1384 for (int pyra=0;pyra<nb_entity_type;pyra++)
1386 int pyra_index = (type%100)*pyra;
1388 int N1 = global_connectivity[pyra_index]-1;
1389 int N2 = global_connectivity[pyra_index+1]-1;
1390 int N3 = global_connectivity[pyra_index+2]-1;
1391 int N4 = global_connectivity[pyra_index+3]-1;
1392 int N5 = global_connectivity[pyra_index+4]-1;
1394 pts[0]=(double *)coord+dim_space*N1;
1395 pts[1]=(double *)coord+dim_space*N2;
1396 pts[2]=(double *)coord+dim_space*N3;
1397 pts[3]=(double *)coord+dim_space*N4;
1398 pts[4]=(double *)coord+dim_space*N5;
1399 CalculateBarycenter<5,3>((const double **)pts,barycenter+3*index);
1404 case MED_PENTA6 : case MED_PENTA15 :
1406 for (int penta=0;penta<nb_entity_type;penta++)
1408 int penta_index = (type%100)*penta;
1410 int N1 = global_connectivity[penta_index]-1;
1411 int N2 = global_connectivity[penta_index+1]-1;
1412 int N3 = global_connectivity[penta_index+2]-1;
1413 int N4 = global_connectivity[penta_index+3]-1;
1414 int N5 = global_connectivity[penta_index+4]-1;
1415 int N6 = global_connectivity[penta_index+5]-1;
1417 pts[0]=(double *)coord+dim_space*N1;
1418 pts[1]=(double *)coord+dim_space*N2;
1419 pts[2]=(double *)coord+dim_space*N3;
1420 pts[3]=(double *)coord+dim_space*N4;
1421 pts[4]=(double *)coord+dim_space*N5;
1422 pts[5]=(double *)coord+dim_space*N6;
1423 CalculateBarycenter<6,3>((const double **)pts,barycenter+3*index);
1428 case MED_HEXA8 : case MED_HEXA20 :
1430 for (int hexa=0;hexa<nb_entity_type;hexa++)
1432 int hexa_index = (type%100)*hexa;
1434 int N1 = global_connectivity[hexa_index]-1;
1435 int N2 = global_connectivity[hexa_index+1]-1;
1436 int N3 = global_connectivity[hexa_index+2]-1;
1437 int N4 = global_connectivity[hexa_index+3]-1;
1438 int N5 = global_connectivity[hexa_index+4]-1;
1439 int N6 = global_connectivity[hexa_index+5]-1;
1440 int N7 = global_connectivity[hexa_index+6]-1;
1441 int N8 = global_connectivity[hexa_index+7]-1;
1443 pts[0]=(double *)coord+dim_space*N1;
1444 pts[1]=(double *)coord+dim_space*N2;
1445 pts[2]=(double *)coord+dim_space*N3;
1446 pts[3]=(double *)coord+dim_space*N4;
1447 pts[4]=(double *)coord+dim_space*N5;
1448 pts[5]=(double *)coord+dim_space*N6;
1449 pts[6]=(double *)coord+dim_space*N7;
1450 pts[7]=(double *)coord+dim_space*N8;
1451 CalculateBarycenter<8,3>((const double **)pts,barycenter+3*index);
1456 case MED_TRIA3 : case MED_TRIA6 :
1458 for (int tria=0;tria<nb_entity_type;tria++)
1460 int tria_index = (type%100)*tria;
1461 int N1 = global_connectivity[tria_index]-1;
1462 int N2 = global_connectivity[tria_index+1]-1;
1463 int N3 = global_connectivity[tria_index+2]-1;
1465 pts[0]=(double *)coord+dim_space*N1;
1466 pts[1]=(double *)coord+dim_space*N2;
1467 pts[2]=(double *)coord+dim_space*N3;
1469 CalculateBarycenter<3,2>((const double **)pts,barycenter+2*index);
1471 CalculateBarycenter<3,3>((const double **)pts,barycenter+3*index);
1476 case MED_QUAD4 : case MED_QUAD8 :
1478 for (int quad=0;quad<nb_entity_type;quad++)
1480 int quad_index = (type%100)*quad;
1481 int N1 = global_connectivity[quad_index]-1;
1482 int N2 = global_connectivity[quad_index+1]-1;
1483 int N3 = global_connectivity[quad_index+2]-1;
1484 int N4 = global_connectivity[quad_index+3]-1;
1486 pts[0]=(double *)coord+dim_space*N1;
1487 pts[1]=(double *)coord+dim_space*N2;
1488 pts[2]=(double *)coord+dim_space*N3;
1489 pts[3]=(double *)coord+dim_space*N4;
1491 CalculateBarycenter<4,2>((const double **)pts,barycenter+2*index);
1493 CalculateBarycenter<4,3>((const double **)pts,barycenter+3*index);
1498 case MED_SEG2 : case MED_SEG3 :
1500 for (int edge=0;edge<nb_entity_type;edge++)
1502 int edge_index = (type%100)*edge;
1503 int N1 = global_connectivity[edge_index]-1;
1504 int N2 = global_connectivity[edge_index+1]-1;
1506 pts[0]=(double *)coord+dim_space*N1;
1507 pts[1]=(double *)coord+dim_space*N2;
1509 CalculateBarycenter<2,2>((const double **)pts,barycenter+2*index);
1511 CalculateBarycenter<2,3>((const double **)pts,barycenter+3*index);
1520 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1521 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1522 for (int polygs=0;polygs<nb_entity_type;polygs++)
1524 int size=connectivity_index[polygs+1]-connectivity_index[polygs];
1525 double **pts=new double * [size];
1526 for(int iPts=0;iPts<size;iPts++)
1527 pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[polygs]+iPts-1]-1);
1528 CalculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1534 const int * supp_number = Support->getNumber(MED_EN::MED_POLYGON);
1535 const int * connectivity = getPolygonsConnectivity(MED_EN::MED_NODAL,support_entity);
1536 const int * connectivity_index = getPolygonsConnectivityIndex(MED_EN::MED_NODAL,support_entity);
1537 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1538 for (int polygs=0;polygs<nb_entity_type;polygs++)
1540 int localPolygsNbP1=supp_number[polygs]-offsetWithClassicType;
1541 int size=connectivity_index[localPolygsNbP1]-connectivity_index[localPolygsNbP1-1];
1542 double **pts=new double * [size];
1543 for(int iPts=0;iPts<size;iPts++)
1544 pts[iPts]=(double *)coord+dim_space*(connectivity[connectivity_index[localPolygsNbP1-1]+iPts-1]-1);
1545 CalculateBarycenterDyn((const double **)pts,size,dim_space,barycenter+dim_space*index);
1552 case MED_EN::MED_POLYHEDRA:
1556 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1559 int offsetWithClassicType=getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1560 int *nodes=_connectivity->getNodesOfPolyhedron(offsetWithClassicType+polyhs+1,lgthNodes);
1561 double **pts=new double * [lgthNodes];
1562 for(int iPts=0;iPts<lgthNodes;iPts++)
1563 pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1564 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1572 const int * supp_number = Support->getNumber(MED_EN::MED_POLYHEDRA);
1573 for (int polyhs=0;polyhs<nb_entity_type;polyhs++)
1576 int *nodes=_connectivity->getNodesOfPolyhedron(supp_number[polyhs],lgthNodes);
1577 double **pts=new double * [lgthNodes];
1578 for(int iPts=0;iPts<lgthNodes;iPts++)
1579 pts[iPts]=(double *)coord+3*(nodes[iPts]-1);
1580 CalculateBarycenterDyn((const double **)pts,lgthNodes,3,barycenter+3*index);
1589 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
1594 if(type != MED_EN::MED_POLYGON && type != MED_EN::MED_POLYHEDRA)
1595 delete [] global_connectivity;
1601 bool MESH::isEmpty() const
1603 bool notempty = _name != "NOT DEFINED" || _coordinate != NULL || _connectivity != NULL ||
1604 _spaceDimension !=MED_INVALID || _meshDimension !=MED_INVALID ||
1605 _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
1606 _familyNode.size() != 0 || _groupCell.size() != 0 ||
1607 _familyCell.size() != 0 || _groupFace.size() != 0 ||
1608 _familyFace.size() != 0 || _groupEdge.size() != 0 ||
1609 _familyEdge.size() != 0 || _isAGrid != 0 ;
1613 void MESH::read(int index)
1615 const char * LOC = "MESH::read(int index=0) : ";
1618 if (_drivers[index]) {
1619 _drivers[index]->open();
1620 _drivers[index]->read();
1621 _drivers[index]->close();
1624 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1625 << "The index given is invalid, index must be between 0 and |"
1631 //=======================================================================
1632 //function : getSkin
1634 //=======================================================================
1636 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
1638 const char * LOC = "MESH::getSkin : " ;
1641 if (this != Support3D->getMesh())
1642 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
1643 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
1644 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
1646 // well, all rigth !
1647 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
1648 mySupport->setAll(false);
1650 list<int> myElementsList ;
1653 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
1654 if (Support3D->isOnAllElements())
1656 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
1657 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
1658 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
1660 int cellNb1 = myConnectivityValue [i];
1661 int cellNb2 = myConnectivityValue [i+1];
1662 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
1663 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
1665 myElementsList.push_back( j ) ;
1672 map<int,int> FaceNbEncounterNb;
1674 int * myConnectivityValue = const_cast <int *>
1675 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
1676 MED_CELL, MED_ALL_ELEMENTS));
1677 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
1678 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
1679 int nbCells = Support3D->getnumber()->getLength();
1680 for (i=0; i<nbCells; ++i)
1682 int cellNb = myCellNbs [ i ];
1683 int faceFirst = myConnectivityIndex[ cellNb-1 ];
1684 int faceLast = myConnectivityIndex[ cellNb ];
1685 for (j = faceFirst; j < faceLast; ++j)
1687 int faceNb = abs( myConnectivityValue [ j-1 ] );
1688 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
1689 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
1690 FaceNbEncounterNb[ faceNb ] = 1;
1692 FaceNbEncounterNb[ faceNb ] ++;
1695 map<int,int>::iterator FaceNbEncounterNbItr;
1696 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
1697 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
1698 FaceNbEncounterNbItr ++)
1699 if ((*FaceNbEncounterNbItr).second == 1)
1701 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
1705 // Well, we must know how many geometric type we have found
1706 int * myListArray = new int[size] ;
1708 list<int>::iterator myElementsListIt ;
1709 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
1710 myListArray[id]=(*myElementsListIt) ;
1714 int numberOfGeometricType ;
1715 medGeometryElement* geometricType ;
1716 int * numberOfGaussPoint ;
1717 int * geometricTypeNumber ;
1718 int * numberOfEntities ;
1719 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
1720 int * mySkyLineArrayIndex ;
1722 int numberOfType = getNumberOfTypes(MED_FACE) ;
1723 if (numberOfType == 1) { // wonderfull : it's easy !
1724 numberOfGeometricType = 1 ;
1725 geometricType = new medGeometryElement[1] ;
1726 const medGeometryElement * allType = getTypes(MED_FACE);
1727 geometricType[0] = allType[0] ;
1728 numberOfGaussPoint = new int[1] ;
1729 numberOfGaussPoint[0] = 1 ;
1730 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
1731 geometricTypeNumber[0] = 0 ;
1732 numberOfEntities = new int[1] ;
1733 numberOfEntities[0] = size ;
1734 mySkyLineArrayIndex = new int[2] ;
1735 mySkyLineArrayIndex[0]=1 ;
1736 mySkyLineArrayIndex[1]=1+size ;
1739 map<medGeometryElement,int> theType ;
1740 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
1741 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
1742 if (theType.find(myType) != theType.end() )
1743 theType[myType]+=1 ;
1747 numberOfGeometricType = theType.size() ;
1748 geometricType = new medGeometryElement[numberOfGeometricType] ;
1749 //const medGeometryElement * allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
1750 numberOfGaussPoint = new int[numberOfGeometricType] ;
1751 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
1752 numberOfEntities = new int[numberOfGeometricType] ;
1753 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
1755 mySkyLineArrayIndex[0]=1 ;
1756 map<medGeometryElement,int>::iterator theTypeIt ;
1757 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
1758 geometricType[index] = (*theTypeIt).first ;
1759 numberOfGaussPoint[index] = 1 ;
1760 geometricTypeNumber[index] = 0 ;
1761 numberOfEntities[index] = (*theTypeIt).second ;
1762 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
1766 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
1767 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
1769 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
1770 mySupport->setGeometricType(geometricType) ;
1771 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
1772 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
1773 mySupport->setNumberOfElements(numberOfEntities) ;
1774 mySupport->setTotalNumberOfElements(size) ;
1775 mySupport->setNumber(mySkyLineArray) ;
1777 delete[] numberOfEntities;
1778 delete[] geometricTypeNumber;
1779 delete[] numberOfGaussPoint;
1780 delete[] geometricType;
1781 delete[] mySkyLineArrayIndex;
1782 delete[] myListArray;
1783 // delete mySkyLineArray;
1791 return a SUPPORT pointer on the union of all SUPPORTs in Supports.
1792 You should delete this pointer after use to avoid memory leaks.
1794 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
1796 const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
1799 SUPPORT * returnedSupport;
1800 string returnedSupportName;
1801 string returnedSupportDescription;
1802 char * returnedSupportNameChar;
1803 char * returnedSupportDescriptionChar;
1804 int size = Supports.size();
1808 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
1809 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1811 returnedSupport = new SUPPORT(*obj);
1813 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
1814 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
1816 returnedSupportNameChar = new char[lenName];
1817 returnedSupportDescriptionChar = new char[lenDescription];
1819 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
1820 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1821 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
1822 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1823 (Supports[0]->getDescription()).c_str());
1825 returnedSupportName = string(returnedSupportNameChar);
1826 returnedSupportDescription = string(returnedSupportDescriptionChar);
1828 returnedSupport->setName(returnedSupportName);
1829 returnedSupport->setDescription(returnedSupportDescription);
1831 delete [] returnedSupportNameChar;
1832 delete [] returnedSupportDescriptionChar;
1836 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1837 returnedSupport = new SUPPORT(*obj);
1839 int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
1840 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
1842 for(int i = 1;i<size;i++)
1844 obj = const_cast <SUPPORT *> (Supports[i]);
1845 returnedSupport->blending(obj);
1849 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
1850 lenDescription = lenDescription + 5 +
1851 strlen((Supports[i]->getDescription()).c_str());
1855 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
1856 lenDescription = lenDescription + 2 +
1857 strlen((Supports[i]->getDescription()).c_str());
1861 returnedSupportNameChar = new char[lenName];
1862 returnedSupportDescriptionChar = new char[lenDescription];
1864 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
1865 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
1867 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1868 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1869 (Supports[0]->getDescription()).c_str());
1871 for(int i = 1;i<size;i++)
1875 returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
1876 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
1878 returnedSupportNameChar = strcat(returnedSupportNameChar,
1879 (Supports[i]->getName()).c_str());
1880 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1881 (Supports[i]->getDescription()).c_str());
1885 returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
1886 returnedSupportNameChar = strcat(returnedSupportNameChar,
1887 (Supports[i]->getName()).c_str());
1889 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
1890 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1891 (Supports[i]->getDescription()).c_str());
1895 returnedSupportName = string(returnedSupportNameChar);
1896 returnedSupport->setName(returnedSupportName);
1898 returnedSupportDescription = string(returnedSupportDescriptionChar);
1899 returnedSupport->setDescription(returnedSupportDescription);
1901 delete [] returnedSupportNameChar;
1902 delete [] returnedSupportDescriptionChar;
1906 return returnedSupport;
1910 return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
1911 The (SUPPORT *) NULL pointer is returned if the intersection is empty.
1912 You should delete this pointer after use to avois memory leaks.
1914 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION)
1916 const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
1919 SUPPORT * returnedSupport;
1920 string returnedSupportName;
1921 string returnedSupportDescription;
1922 char * returnedSupportNameChar;
1923 char * returnedSupportDescriptionChar;
1924 int size = Supports.size();
1928 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
1929 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1931 returnedSupport = new SUPPORT(*obj);
1933 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
1934 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
1936 returnedSupportNameChar = new char[lenName];
1937 returnedSupportDescriptionChar = new char[lenDescription];
1939 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
1940 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
1941 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
1942 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
1943 (Supports[0]->getDescription()).c_str());
1945 returnedSupportName = string(returnedSupportNameChar);
1946 returnedSupportDescription = string(returnedSupportDescriptionChar);
1948 returnedSupport->setName(returnedSupportName);
1949 returnedSupport->setDescription(returnedSupportDescription);
1951 delete [] returnedSupportNameChar;
1952 delete [] returnedSupportDescriptionChar;
1956 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
1957 returnedSupport = new SUPPORT(*obj);
1959 int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
1960 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
1962 for(int i = 1;i<size;i++)
1964 obj = const_cast <SUPPORT *> (Supports[i]);
1965 returnedSupport->intersecting(obj);
1969 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
1970 lenDescription = lenDescription + 5 +
1971 strlen((Supports[i]->getDescription()).c_str());
1975 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
1976 lenDescription = lenDescription + 2 +
1977 strlen((Supports[i]->getDescription()).c_str());
1980 if(returnedSupport != (SUPPORT *) NULL)
1982 returnedSupportNameChar = new char[lenName];
1983 returnedSupportDescriptionChar = new char[lenDescription];
1985 returnedSupportNameChar = strcpy(returnedSupportNameChar,
1986 "Intersection of ");
1987 returnedSupportDescriptionChar =
1988 strcpy(returnedSupportDescriptionChar,"Intersection of ");
1990 returnedSupportNameChar = strcat(returnedSupportNameChar,
1991 (Supports[0]->getName()).c_str());
1992 returnedSupportDescriptionChar =
1993 strcat(returnedSupportDescriptionChar,
1994 (Supports[0]->getDescription()).c_str());
1996 for(int i = 1;i<size;i++)
2000 returnedSupportNameChar = strcat(returnedSupportNameChar,
2002 returnedSupportDescriptionChar =
2003 strcat(returnedSupportDescriptionChar," and ");
2005 returnedSupportNameChar =
2006 strcat(returnedSupportNameChar,
2007 (Supports[i]->getName()).c_str());
2008 returnedSupportDescriptionChar =
2009 strcat(returnedSupportDescriptionChar,
2010 (Supports[i]->getDescription()).c_str());
2014 returnedSupportNameChar = strcat(returnedSupportNameChar,
2016 returnedSupportNameChar =
2017 strcat(returnedSupportNameChar,
2018 (Supports[i]->getName()).c_str());
2020 returnedSupportDescriptionChar =
2021 strcat(returnedSupportDescriptionChar,", ");
2022 returnedSupportDescriptionChar =
2023 strcat(returnedSupportDescriptionChar,
2024 (Supports[i]->getDescription()).c_str());
2028 returnedSupportName = string(returnedSupportNameChar);
2029 returnedSupport->setName(returnedSupportName);
2031 returnedSupportDescription = string(returnedSupportDescriptionChar);
2032 returnedSupport->setDescription(returnedSupportDescription);
2034 delete [] returnedSupportNameChar;
2035 delete [] returnedSupportDescriptionChar;
2040 return returnedSupport;
2044 // internal helper type
2047 std::vector<int> groups;
2048 MED_EN::medGeometryElement geometricType;
2051 // Create families from groups
2052 void MESH::createFamilies()
2054 int idFamNode = 0; // identifier for node families
2055 int idFamElement = 0; // identifier for cell, face or edge families
2057 // Main loop on mesh's entities
2058 for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2060 int numberofgroups = getNumberOfGroups(entity);
2062 continue; // no groups for this entity
2064 vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2066 // make myFamilies points to the member corresponding to entity
2067 vector<FAMILY*>* myFamilies;
2071 myFamilies = & _familyCell;
2074 myFamilies = & _familyFace;
2077 myFamilies = & _familyEdge;
2080 myFamilies = & _familyNode;
2084 vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2085 // get a copy of the (old) family ptrs before clearing
2086 vector<FAMILY*> myOldFamilies=getFamilies(entity);
2087 myFamilies->clear();
2090 // 1 - Create a vector containing for each cell (of the entity) an information structure
2091 // giving geometric type and the groups it belong to
2093 med_int numberOfTypes=getNumberOfTypes(entity);
2094 const int * index=getGlobalNumberingIndex(entity);
2095 const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2096 med_int numberOfCells=index[numberOfTypes]-1; // total number of cells for that entity
2097 SCRUTE(numberOfTypes);
2098 SCRUTE(numberOfCells);
2099 vector< _cell > tab_cell(numberOfCells);
2100 for(med_int t=0; t!=numberOfTypes; ++t)
2101 for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2102 tab_cell[n].geometricType=geometricTypes[t];
2105 // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2107 for (unsigned g=0; g!=myGroups.size(); ++g)
2109 // scan cells that belongs to the group
2110 const int* groupCells=myGroups[g]->getnumber()->getValue();
2111 int nbCells=myGroups[g]->getnumber()->getLength();
2112 for(int c=0; c!=nbCells; ++c)
2113 tab_cell[groupCells[c]-1].groups.push_back(g);
2117 // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2118 // whith the vector of contained cells
2120 map< string,vector<int> > tab_families;
2121 map< string,vector<int> >::iterator fam;
2122 for(int n=0; n!=numberOfCells; ++n)
2124 ostringstream key; // to generate the name of the family
2126 if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2127 key << "_NONE" << entity;
2129 for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2131 string groupName=myGroups[*it]->getName();
2132 if(groupName.empty())
2135 key << "_" << groupName;
2138 tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2142 // 4 - Scan the family map, create MED Families, check if it already exist.
2144 for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2146 vector<medGeometryElement> tab_types_geometriques;
2147 medGeometryElement geometrictype=MED_NONE;
2148 vector<int> tab_index_types_geometriques;
2149 vector<int> tab_nombres_elements;
2151 // scan family cells and fill the tab that are needed by the create a MED FAMILY
2152 for( int i=0; i!=fam->second.size(); ++i)
2154 int ncell=fam->second[i]-1;
2155 if(tab_cell[ncell].geometricType != geometrictype)
2157 // new geometric type -> we store it and complete index tabs
2158 if(!tab_index_types_geometriques.empty())
2159 tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2160 tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2161 tab_index_types_geometriques.push_back(i+1);
2164 // store and complete index tabs for the last geometric type
2165 tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2166 tab_index_types_geometriques.push_back(fam->second.size()+1);
2168 // create a empty MED FAMILY and fill it with the tabs we constructed
2169 FAMILY* newFam = new FAMILY();
2170 newFam->setTotalNumberOfElements(fam->second.size());
2171 newFam->setName(fam->first);
2172 newFam->setMesh(this);
2173 newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2174 newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2175 newFam->setNumberOfElements(&tab_nombres_elements[0]);
2176 newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2177 newFam->setEntity(entity);
2178 newFam->setAll(false);
2190 idFam = -idFamElement;
2194 idFam = -idFamElement;
2198 idFam = -idFamElement;
2202 newFam->setIdentifier(idFam);
2204 // Update links between families and groups
2206 int ncell1=fam->second[0]-1; // number of first cell in family
2207 int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2210 newFam->setNumberOfGroups(numberOfGroups);
2211 string * groupNames=new string[numberOfGroups];
2213 // iterate on the groups the family belongs to
2214 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2215 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2217 whichFamilyInGroup[*it].push_back(newFam);
2218 groupNames[ng]=myGroups[*it]->getName();
2220 newFam->setGroupsNames(groupNames);
2223 int sizeOfFamVect = myFamilies->size();
2225 MESSAGE(" MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2227 myFamilies->push_back(newFam);
2230 // delete old families
2231 for (int i=0;i<myOldFamilies.size();i++)
2232 delete myOldFamilies[i] ;
2234 // update references in groups
2235 for (int i=0;i<myGroups.size();i++)
2237 myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2238 myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2241 // re-scan the cells vector, and fill the family vector with cells.
2242 // creation of support, check if it already exist.
2246 int MESH::getElementContainingPoint(const double *coord)
2248 if(_spaceDimension==3)
2250 Meta_Wrapper<3> *fromWrapper=new Meta_Wrapper<3> (getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2251 const_cast<CONNECTIVITY *>(getConnectivityptr()));
2252 Meta_Wrapper<3> *toWrapper=new Meta_Wrapper<3> (1,const_cast<double *>(coord));
2253 Meta_Mapping<3> *mapping=new Meta_Mapping<3> (fromWrapper,toWrapper);
2254 mapping->Cree_Mapping(1);
2255 vector<int> vectormapping=mapping->Get_Mapping();
2256 return vectormapping[0]+1;
2258 else if(_spaceDimension==2)
2260 Meta_Wrapper<2> *fromWrapper=new Meta_Wrapper<2> (getNumberOfNodes(),const_cast<double *>(getCoordinates(MED_FULL_INTERLACE)),
2261 const_cast<CONNECTIVITY *>(getConnectivityptr()));
2262 Meta_Wrapper<2> *toWrapper=new Meta_Wrapper<2> (1,const_cast<double *>(coord));
2263 Meta_Mapping<2> *mapping=new Meta_Mapping<2> (fromWrapper,toWrapper);
2264 mapping->Cree_Mapping(1);
2265 vector<int> vectormapping=mapping->Get_Mapping();
2266 return vectormapping[0]+1;
2269 throw MEDEXCEPTION("MESH::getElementContainingPoint : invalid _spaceDimension must be equal to 2 or 3 !!!");
2272 //Presently disconnected in C++
2273 void MESH::addReference() const
2277 //Presently disconnected in C++
2278 void MESH::removeReference() const