12 #include "MEDMEM_DriversDef.hxx"
13 #include "MEDMEM_Field.hxx"
14 #include "MEDMEM_Mesh.hxx"
16 #include "MEDMEM_Support.hxx"
17 #include "MEDMEM_Family.hxx"
18 #include "MEDMEM_Group.hxx"
19 #include "MEDMEM_Coordinate.hxx"
20 #include "MEDMEM_Connectivity.hxx"
21 #include "MEDMEM_CellModel.hxx"
23 #include "MEDMEM_DriverFactory.hxx"
26 using namespace MEDMEM;
27 using namespace MED_EN;
29 //#include "MEDMEM_Grid.hxx" this inclision should have never be here !!!
31 //update Families with content list
32 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
34 // ------- Drivers Management Part
36 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
37 is driverName. addDriver returns an integer handler. */
38 int MESH::addDriver(driverTypes driverType,
39 const string & fileName/*="Default File Name.med"*/,
40 const string & driverName/*="Default Mesh Name"*/,
41 med_mode_acces access) {
43 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) : ";
51 driver = DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,
54 _drivers.push_back(driver);
56 int current = _drivers.size()-1;
58 _drivers[current]->setMeshName(driverName);
65 /*! Add an existing MESH driver. */
66 int MESH::addDriver(GENDRIVER & driver) {
67 const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
70 // A faire : Vérifier que le driver est de type MESH.
71 GENDRIVER * newDriver = driver.copy() ;
73 _drivers.push_back(newDriver);
74 return _drivers.size()-1;
79 /*! Remove an existing MESH driver. */
80 void MESH::rmDriver (int index/*=0*/) {
81 const char * LOC = "MESH::rmDriver (int index=0): ";
84 if ( _drivers[index] ) {
85 //_drivers.erase(&_drivers[index]);
90 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
91 << "The index given is invalid, index must be between 0 and |"
100 // ------ End of Drivers Management Part
105 const char * LOC = "MESH::init(): ";
109 _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
111 _coordinate = (COORDINATE *) NULL;
112 _connectivity = (CONNECTIVITY *) NULL;
114 _spaceDimension = MED_INVALID; // 0 ?!?
115 _meshDimension = MED_INVALID;
116 _numberOfNodes = MED_INVALID;
120 _arePresentOptionnalNodesNumbers = 0;
125 /*! Create an empty MESH. */
126 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
133 _isAGrid = m._isAGrid;
135 if (m._coordinate != NULL)
136 _coordinate = new COORDINATE(* m._coordinate);
138 _coordinate = (COORDINATE *) NULL;
140 if (m._connectivity != NULL)
141 _connectivity = new CONNECTIVITY(* m._connectivity);
143 _connectivity = (CONNECTIVITY *) NULL;
145 _spaceDimension = m._spaceDimension;
146 _meshDimension = m._meshDimension;
147 _numberOfNodes = m._numberOfNodes;
149 _familyNode = m._familyNode;
150 for (int i=0; i<m._familyNode.size(); i++)
152 _familyNode[i] = new FAMILY(* m._familyNode[i]);
153 _familyNode[i]->setMesh(this);
156 _familyCell = m._familyCell;
157 for (int i=0; i<m._familyCell.size(); i++)
159 _familyCell[i] = new FAMILY(* m._familyCell[i]);
160 _familyCell[i]->setMesh(this);
163 _familyFace = m._familyFace;
164 for (int i=0; i<m._familyFace.size(); i++)
166 _familyFace[i] = new FAMILY(* m._familyFace[i]);
167 _familyFace[i]->setMesh(this);
170 _familyEdge = m._familyEdge;
171 for (int i=0; i<m._familyEdge.size(); i++)
173 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
174 _familyEdge[i]->setMesh(this);
177 _groupNode = m._groupNode;
178 for (int i=0; i<m._groupNode.size(); i++)
180 _groupNode[i] = new GROUP(* m._groupNode[i]);
181 _groupNode[i]->setMesh(this);
184 _groupCell = m._groupCell;
185 for (int i=0; i<m._groupCell.size(); i++)
187 _groupCell[i] = new GROUP(* m._groupCell[i]);
188 _groupCell[i]->setMesh(this);
191 _groupFace = m._groupFace;
192 for (int i=0; i<m._groupFace.size(); i++)
194 _groupFace[i] = new GROUP(* m._groupFace[i]);
195 _groupFace[i]->setMesh(this);
198 _groupEdge = m._groupEdge;
199 for (int i=0; i<m._groupEdge.size(); i++)
201 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
202 _groupEdge[i]->setMesh(this);
205 //_drivers = m._drivers; //Recopie des drivers?
210 MESSAGE("MESH::~MESH() : Destroying the Mesh");
211 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
212 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
214 size = _familyNode.size() ;
215 for (int i=0;i<size;i++)
216 delete _familyNode[i] ;
217 size = _familyCell.size() ;
218 for (int i=0;i<size;i++)
219 delete _familyCell[i] ;
220 size = _familyFace.size() ;
221 for (int i=0;i<size;i++)
222 delete _familyFace[i] ;
223 size = _familyEdge.size() ;
224 for (int i=0;i<size;i++)
225 delete _familyEdge[i] ;
226 size = _groupNode.size() ;
227 for (int i=0;i<size;i++)
228 delete _groupNode[i] ;
229 size = _groupCell.size() ;
230 for (int i=0;i<size;i++)
231 delete _groupCell[i] ;
232 size = _groupFace.size() ;
233 for (int i=0;i<size;i++)
234 delete _groupFace[i] ;
235 size = _groupEdge.size() ;
236 for (int i=0;i<size;i++)
237 delete _groupEdge[i] ;
239 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
241 for (unsigned int index=0; index < _drivers.size(); index++ )
243 SCRUTE(_drivers[index]);
244 if ( _drivers[index] != NULL) delete _drivers[index];
249 MESH & MESH::operator=(const MESH &m)
251 const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
254 MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
257 // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
258 // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
259 // _drivers = m._drivers;
261 // space_dimension=m.space_dimension;
262 // mesh_dimension=m.mesh_dimension;
264 // nodes_count=m.nodes_count;
265 // coordinates=m.coordinates;
266 // coordinates_name=m.coordinates_name ;
267 // coordinates_unit=m.coordinates_unit ;
268 // nodes_numbers=m.nodes_numbers ;
270 // cells_types_count=m.cells_types_count;
271 // cells_type=m.cells_type;
272 // cells_count=m.cells_count;
273 // nodal_connectivity=m.nodal_connectivity;
275 // nodes_families_count=m.nodes_families_count;
276 // nodes_Families=m.nodes_Families;
278 // cells_families_count=m.cells_families_count;
279 // cells_Families=m.cells_Families;
281 // maximum_cell_number_by_node = m.maximum_cell_number_by_node;
282 // if (maximum_cell_number_by_node > 0)
284 // reverse_nodal_connectivity = m.reverse_nodal_connectivity;
285 // reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
292 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
293 The meshname driverName must already exists in the file.*/
294 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) throw (MEDEXCEPTION)
296 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
302 GENDRIVER *myDriver=DRIVERFACTORY::buildDriverForMesh(driverType,fileName,this,driverName,MED_LECT);
303 current = addDriver(*myDriver);
305 _drivers[current]->open();
306 _drivers[current]->read();
307 _drivers[current]->close();
310 // ((GRID *) this)->fillMeshAfterRead();
316 // Node MESH::Donne_Barycentre(const Maille &m) const
318 // Node N=node[m[0]];
319 // for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
320 // N=N*(1.0/m.donne_nbr_sommets());
324 ostream & MEDMEM::operator<<(ostream &os, const MESH &myMesh)
326 int spacedimension = myMesh.getSpaceDimension();
327 int meshdimension = myMesh.getMeshDimension();
328 int numberofnodes = myMesh.getNumberOfNodes();
330 os << "Space Dimension : " << spacedimension << endl << endl;
332 os << "Mesh Dimension : " << meshdimension << endl << endl;
334 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
335 os << "SHOW NODES COORDINATES : " << endl;
337 os << "Name :" << endl;
338 const string * coordinatesnames = myMesh.getCoordinatesNames();
339 for(int i=0; i<spacedimension ; i++)
341 os << " - " << coordinatesnames[i] << endl;
343 os << "Unit :" << endl;
344 const string * coordinatesunits = myMesh.getCoordinatesUnits();
345 for(int i=0; i<spacedimension ; i++)
347 os << " - " << coordinatesunits[i] << endl;
349 for(int i=0; i<numberofnodes ; i++)
351 os << "Nodes " << i+1 << " : ";
352 for (int j=0; j<spacedimension ; j++)
353 os << coordinates[i*spacedimension+j] << " ";
357 os << endl << "SHOW CONNECTIVITY :" << endl;
358 os << *myMesh._connectivity << endl;
360 medEntityMesh entity;
361 os << endl << "SHOW FAMILIES :" << endl << endl;
362 for (int k=1; k<=4; k++)
364 if (k==1) entity = MED_NODE;
365 if (k==2) entity = MED_CELL;
366 if (k==3) entity = MED_FACE;
367 if (k==4) entity = MED_EDGE;
368 int numberoffamilies = myMesh.getNumberOfFamilies(entity);
369 os << "NumberOfFamilies on "<< entNames[entity] <<" : "<<numberoffamilies<<endl;
370 using namespace MED_EN;
371 for (int i=1; i<numberoffamilies+1;i++)
373 os << * myMesh.getFamily(entity,i) << endl;
377 os << endl << "SHOW GROUPS :" << endl << endl;
378 for (int k=1; k<=4; k++)
380 if (k==1) entity = MED_NODE;
381 if (k==2) entity = MED_CELL;
382 if (k==3) entity = MED_FACE;
383 if (k==4) entity = MED_EDGE;
384 int numberofgroups = myMesh.getNumberOfGroups(entity);
385 os << "NumberOfGroups on "<< entNames[entity] <<" : "<<numberofgroups<<endl;
386 using namespace MED_EN;
387 for (int i=1; i<numberofgroups+1;i++)
389 os << * myMesh.getGroup(entity,i) << endl;
397 Get global number of element which have same connectivity than connectivity argument.
399 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
401 Return -1 if not found.
403 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
405 const char* LOC="MESH::getElementNumber " ;
408 int numberOfValue ; // size of connectivity array
409 CELLMODEL myType(Type) ;
410 if (ConnectivityType==MED_DESCENDING)
411 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
413 numberOfValue = myType.getNumberOfNodes() ; // nodes
415 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
416 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
418 // First node or face/edge
419 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
420 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
422 // list to put cells numbers
423 list<int> cellsList ;
424 list<int>::iterator itList ;
425 for (int i=indexBegin; i<indexEnd; i++)
426 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
428 // Foreach node or face/edge in cell, we search cells which are in common.
429 // At the end, we must have only one !
431 for (int i=1; i<numberOfValue; i++) {
432 int connectivity_i = connectivity[i] ;
433 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
435 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
436 if ((*itList)==myReverseConnectivityValue[j-1]) {
442 itList=cellsList.erase(itList);
443 itList--; // well : rigth if itList = cellsList.begin() ??
448 if (cellsList.size()>1) // we have more than one cell
449 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
451 if (cellsList.size()==0)
456 return cellsList.front() ;
461 Return a support which reference all elements on the boundary of mesh.
463 For instance, we could get only face in 3D and edge in 2D.
465 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
468 const char * LOC = "MESH::getBoundaryElements : " ;
471 // actually we could only get face (in 3D) and edge (in 2D)
472 if (_spaceDimension == 3)
473 if (Entity != MED_FACE)
474 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
475 if (_spaceDimension == 2)
476 if (Entity != MED_EDGE)
477 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
480 SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
481 //mySupport.setDescription("boundary of type ...");
482 mySupport->setAll(false);
485 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
486 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
487 int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
488 list<int> myElementsList ;
490 for (int i=0 ; i<numberOf; i++)
491 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
492 myElementsList.push_back(i+1) ;
495 // Well, we must know how many geometric type we have found
496 int * myListArray = new int[size] ;
498 list<int>::iterator myElementsListIt ;
499 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
500 myListArray[id]=(*myElementsListIt) ;
504 int numberOfGeometricType ;
505 medGeometryElement* geometricType ;
506 int * numberOfGaussPoint ;
507 int * geometricTypeNumber ;
508 int * numberOfElements ;
509 //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
510 int * mySkyLineArrayIndex ;
512 int numberOfType = getNumberOfTypes(Entity) ;
513 if (numberOfType == 1) { // wonderfull : it's easy !
514 numberOfGeometricType = 1 ;
515 geometricType = new medGeometryElement[1] ;
516 const medGeometryElement * allType = getTypes(Entity);
517 geometricType[0] = allType[0] ;
518 numberOfGaussPoint = new int[1] ;
519 numberOfGaussPoint[0] = 1 ;
520 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
521 geometricTypeNumber[0] = 0 ;
522 numberOfElements = new int[1] ;
523 numberOfElements[0] = size ;
524 mySkyLineArrayIndex = new int[2] ;
525 mySkyLineArrayIndex[0]=1 ;
526 mySkyLineArrayIndex[1]=1+size ;
529 map<medGeometryElement,int> theType ;
530 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
531 medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
532 if (theType.find(myType) != theType.end() )
537 numberOfGeometricType = theType.size() ;
538 geometricType = new medGeometryElement[numberOfGeometricType] ;
539 //const medGeometryElement * allType = getTypes(Entity); !! UNUSZED VARIABLE !!
540 numberOfGaussPoint = new int[numberOfGeometricType] ;
541 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
542 numberOfElements = new int[numberOfGeometricType] ;
543 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
545 mySkyLineArrayIndex[0]=1 ;
546 map<medGeometryElement,int>::iterator theTypeIt ;
547 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
548 geometricType[index] = (*theTypeIt).first ;
549 numberOfGaussPoint[index] = 1 ;
550 geometricTypeNumber[index] = 0 ;
551 numberOfElements[index] = (*theTypeIt).second ;
552 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
556 //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
557 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
559 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
560 mySupport->setGeometricType(geometricType) ;
561 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
562 mySupport->setNumberOfElements(numberOfElements) ;
563 mySupport->setTotalNumberOfElements(size) ;
564 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
565 mySupport->setNumber(mySkyLineArray) ;
567 delete[] numberOfElements;
568 delete[] geometricTypeNumber;
569 delete[] numberOfGaussPoint;
570 delete[] geometricType;
571 delete[] mySkyLineArrayIndex;
572 delete[] myListArray;
573 // delete mySkyLineArray;
579 FIELD<double>* MESH::getVolume(const SUPPORT * Support) const throw (MEDEXCEPTION)
581 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
584 // Support must be on 3D elements
586 // Make sure that the MESH class is the same as the MESH class attribut
587 // in the class Support
588 MESH* myMesh = Support->getMesh();
590 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
592 int dim_space = getSpaceDimension();
593 medEntityMesh support_entity = Support->getEntity();
594 bool onAll = Support->isOnAllElements();
596 int nb_type, length_values;
597 const medGeometryElement* types;
599 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
600 const int* global_connectivity;
604 // nb_type = myMesh->getNumberOfTypes(support_entity);
605 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
606 // types = getTypes(support_entity);
610 nb_type = Support->getNumberOfTypes();
611 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
612 types = Support->getTypes();
616 FIELD<double>* Volume = new FIELD<double>(Support,1);
617 // double *volume = new double [length_values];
618 Volume->setName("VOLUME");
619 Volume->setDescription("cells volume");
620 Volume->setComponentName(1,"volume");
621 Volume->setComponentDescription(1,"desc-comp");
623 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
625 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
627 Volume->setMEDComponentUnit(1,MEDComponentUnit);
629 Volume->setValueType(MED_REEL64);
631 Volume->setIterationNumber(0);
632 Volume->setOrderNumber(0);
633 Volume->setTime(0.0);
635 //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
636 MEDARRAY<double> *volume = Volume->getvalue();
639 const double * coord = getCoordinates(MED_FULL_INTERLACE);
641 for (int i=0;i<nb_type;i++)
643 medGeometryElement type = types[i] ;
648 nb_entity_type = getNumberOfElements(support_entity,type);
649 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
653 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all"));
655 nb_entity_type = Support->getNumberOfElements(type);
657 const int * supp_number = Support->getNumber(type);
658 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
659 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
660 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
662 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
663 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
664 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
667 global_connectivity = global_connectivity_tmp ;
672 case MED_TETRA4 : case MED_TETRA10 :
674 for (int tetra=0;tetra<nb_entity_type;tetra++)
676 int tetra_index = (type%100)*tetra;
678 int N1 = global_connectivity[tetra_index]-1;
679 int N2 = global_connectivity[tetra_index+1]-1;
680 int N3 = global_connectivity[tetra_index+2]-1;
681 int N4 = global_connectivity[tetra_index+3]-1;
683 double x1 = coord[dim_space*N1];
684 double x2 = coord[dim_space*N2];
685 double x3 = coord[dim_space*N3];
686 double x4 = coord[dim_space*N4];
688 double y1 = coord[dim_space*N1+1];
689 double y2 = coord[dim_space*N2+1];
690 double y3 = coord[dim_space*N3+1];
691 double y4 = coord[dim_space*N4+1];
693 double z1 = coord[dim_space*N1+2];
694 double z2 = coord[dim_space*N2+2];
695 double z3 = coord[dim_space*N3+2];
696 double z4 = coord[dim_space*N4+2];
698 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
699 (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
700 (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
702 //volume[index] = xvolume ;
703 volume->setIJ(index,1,xvolume) ;
708 case MED_PYRA5 : case MED_PYRA13 :
710 for (int pyra=0;pyra<nb_entity_type;pyra++)
712 int pyra_index = (type%100)*pyra;
714 int N1 = global_connectivity[pyra_index]-1;
715 int N2 = global_connectivity[pyra_index+1]-1;
716 int N3 = global_connectivity[pyra_index+2]-1;
717 int N4 = global_connectivity[pyra_index+3]-1;
718 int N5 = global_connectivity[pyra_index+4]-1;
720 double x1 = coord[dim_space*N1];
721 double x2 = coord[dim_space*N2];
722 double x3 = coord[dim_space*N3];
723 double x4 = coord[dim_space*N4];
724 double x5 = coord[dim_space*N5];
726 double y1 = coord[dim_space*N1+1];
727 double y2 = coord[dim_space*N2+1];
728 double y3 = coord[dim_space*N3+1];
729 double y4 = coord[dim_space*N4+1];
730 double y5 = coord[dim_space*N5+1];
732 double z1 = coord[dim_space*N1+2];
733 double z2 = coord[dim_space*N2+2];
734 double z3 = coord[dim_space*N3+2];
735 double z4 = coord[dim_space*N4+2];
736 double z5 = coord[dim_space*N5+2];
738 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
739 (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
740 (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
741 ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
742 (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
743 (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
746 //volume[index] = xvolume ;
747 volume->setIJ(index,1,xvolume) ;
752 case MED_PENTA6 : case MED_PENTA15 :
754 for (int penta=0;penta<nb_entity_type;penta++)
756 int penta_index = (type%100)*penta;
758 int N1 = global_connectivity[penta_index]-1;
759 int N2 = global_connectivity[penta_index+1]-1;
760 int N3 = global_connectivity[penta_index+2]-1;
761 int N4 = global_connectivity[penta_index+3]-1;
762 int N5 = global_connectivity[penta_index+4]-1;
763 int N6 = global_connectivity[penta_index+5]-1;
765 double x1 = coord[dim_space*N1];
766 double x2 = coord[dim_space*N2];
767 double x3 = coord[dim_space*N3];
768 double x4 = coord[dim_space*N4];
769 double x5 = coord[dim_space*N5];
770 double x6 = coord[dim_space*N6];
772 double y1 = coord[dim_space*N1+1];
773 double y2 = coord[dim_space*N2+1];
774 double y3 = coord[dim_space*N3+1];
775 double y4 = coord[dim_space*N4+1];
776 double y5 = coord[dim_space*N5+1];
777 double y6 = coord[dim_space*N6+1];
779 double z1 = coord[dim_space*N1+2];
780 double z2 = coord[dim_space*N2+2];
781 double z3 = coord[dim_space*N3+2];
782 double z4 = coord[dim_space*N4+2];
783 double z5 = coord[dim_space*N5+2];
784 double z6 = coord[dim_space*N6+2];
786 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
787 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
788 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
789 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
790 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
791 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
792 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
794 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
796 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
798 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
799 (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
800 (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
801 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
803 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
805 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
806 (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
807 (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
808 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
810 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
812 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
813 (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
814 (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
816 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
818 //volume[index] = xvolume ;
819 volume->setIJ(index,1,xvolume) ;
824 case MED_HEXA8 : case MED_HEXA20 :
826 for (int hexa=0;hexa<nb_entity_type;hexa++)
828 int hexa_index = (type%100)*hexa;
830 int N1 = global_connectivity[hexa_index]-1;
831 int N2 = global_connectivity[hexa_index+1]-1;
832 int N3 = global_connectivity[hexa_index+2]-1;
833 int N4 = global_connectivity[hexa_index+3]-1;
834 int N5 = global_connectivity[hexa_index+4]-1;
835 int N6 = global_connectivity[hexa_index+5]-1;
836 int N7 = global_connectivity[hexa_index+6]-1;
837 int N8 = global_connectivity[hexa_index+7]-1;
839 double x1 = coord[dim_space*N1];
840 double x2 = coord[dim_space*N2];
841 double x3 = coord[dim_space*N3];
842 double x4 = coord[dim_space*N4];
843 double x5 = coord[dim_space*N5];
844 double x6 = coord[dim_space*N6];
845 double x7 = coord[dim_space*N7];
846 double x8 = coord[dim_space*N8];
848 double y1 = coord[dim_space*N1+1];
849 double y2 = coord[dim_space*N2+1];
850 double y3 = coord[dim_space*N3+1];
851 double y4 = coord[dim_space*N4+1];
852 double y5 = coord[dim_space*N5+1];
853 double y6 = coord[dim_space*N6+1];
854 double y7 = coord[dim_space*N7+1];
855 double y8 = coord[dim_space*N8+1];
857 double z1 = coord[dim_space*N1+2];
858 double z2 = coord[dim_space*N2+2];
859 double z3 = coord[dim_space*N3+2];
860 double z4 = coord[dim_space*N4+2];
861 double z5 = coord[dim_space*N5+2];
862 double z6 = coord[dim_space*N6+2];
863 double z7 = coord[dim_space*N7+2];
864 double z8 = coord[dim_space*N8+2];
866 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
867 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
868 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
869 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
870 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
871 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
872 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
873 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
874 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
875 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
876 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
877 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
879 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
881 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
883 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
884 (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
885 (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
886 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
888 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
890 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
891 (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
892 (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
893 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
894 (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
895 (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
896 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
897 (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
898 (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
899 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
900 ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
901 ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
902 ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
903 ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
904 ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
905 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
907 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
909 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
910 (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
911 (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
912 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
914 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
916 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
917 (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
918 (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
919 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
920 (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
921 (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
922 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
923 (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
924 (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
925 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
926 ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
927 ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
928 ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
929 ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
930 ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
931 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
932 (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
933 (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
934 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
935 (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
936 (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
937 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
938 ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
939 ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
940 ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
941 ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
942 ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
943 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
944 (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
945 (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
946 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
947 (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
948 (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
949 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
950 ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
951 ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
952 ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
953 ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
954 ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
955 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
956 (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
957 (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
958 (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
959 (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
960 (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
961 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
962 (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
963 (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
964 (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
965 (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
966 (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
967 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
968 (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
969 ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
970 (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
971 ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
972 (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
973 ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
974 (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
975 ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
976 (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
977 ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
978 (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
980 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
981 4.0*(C + F + G + H + L + O + P + Q + S + T +
982 V + W) + 2.0*(I + R + U + X + Y + Z) +
985 //volume[index] = xvolume ;
986 volume->setIJ(index,1,xvolume) ;
992 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
996 if (!onAll) delete [] global_connectivity ;
1002 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1004 const char * LOC = "MESH::getArea(SUPPORT*) : ";
1007 // Support must be on 2D elements
1009 // Make sure that the MESH class is the same as the MESH class attribut
1010 // in the class Support
1011 MESH* myMesh = Support->getMesh();
1013 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1015 int dim_space = getSpaceDimension();
1016 medEntityMesh support_entity = Support->getEntity();
1017 bool onAll = Support->isOnAllElements();
1019 int nb_type, length_values;
1020 const medGeometryElement* types;
1022 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1023 const int* global_connectivity;
1027 // nb_type = myMesh->getNumberOfTypes(support_entity);
1028 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1029 // types = getTypes(support_entity);
1033 nb_type = Support->getNumberOfTypes();
1034 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1035 types = Support->getTypes();
1039 FIELD<double>* Area;
1041 Area = new FIELD<double>(Support,1);
1042 Area->setName("AREA");
1043 Area->setDescription("cells or faces area");
1045 Area->setComponentName(1,"area");
1046 Area->setComponentDescription(1,"desc-comp");
1048 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1050 string MEDComponentUnit="trucmuch";
1052 Area->setMEDComponentUnit(1,MEDComponentUnit);
1054 Area->setValueType(MED_REEL64);
1056 Area->setIterationNumber(0);
1057 Area->setOrderNumber(0);
1060 double *area = new double[length_values];
1061 //double *area = Area->getValue(MED_FULL_INTERLACE);
1063 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1066 for (int i=0;i<nb_type;i++)
1068 medGeometryElement type = types[i] ;
1073 nb_entity_type = getNumberOfElements(support_entity,type);
1074 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1078 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1080 nb_entity_type = Support->getNumberOfElements(type);
1082 const int * supp_number = Support->getNumber(type);
1083 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1084 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1085 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1087 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1088 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1089 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1093 global_connectivity = global_connectivity_tmp ;
1099 case MED_TRIA3 : case MED_TRIA6 :
1101 for (int tria=0;tria<nb_entity_type;tria++)
1103 int tria_index = (type%100)*tria;
1105 int N1 = global_connectivity[tria_index]-1;
1106 int N2 = global_connectivity[tria_index+1]-1;
1107 int N3 = global_connectivity[tria_index+2]-1;
1109 double x1 = coord[dim_space*N1];
1110 double x2 = coord[dim_space*N2];
1111 double x3 = coord[dim_space*N3];
1113 double y1 = coord[dim_space*N1+1];
1114 double y2 = coord[dim_space*N2+1];
1115 double y3 = coord[dim_space*N3+1];
1119 xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1123 double z1 = coord[dim_space*N1+2];
1124 double z2 = coord[dim_space*N2+2];
1125 double z3 = coord[dim_space*N3+2];
1127 xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1128 ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1129 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1130 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1131 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1132 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1135 area[index] = xarea ;
1140 case MED_QUAD4 : case MED_QUAD8 :
1142 for (int quad=0;quad<nb_entity_type;quad++)
1144 int quad_index = (type%100)*quad;
1146 int N1 = global_connectivity[quad_index]-1;
1147 int N2 = global_connectivity[quad_index+1]-1;
1148 int N3 = global_connectivity[quad_index+2]-1;
1149 int N4 = global_connectivity[quad_index+3]-1;
1151 double x1 = coord[dim_space*N1];
1152 double x2 = coord[dim_space*N2];
1153 double x3 = coord[dim_space*N3];
1154 double x4 = coord[dim_space*N4];
1156 double y1 = coord[dim_space*N1+1];
1157 double y2 = coord[dim_space*N2+1];
1158 double y3 = coord[dim_space*N3+1];
1159 double y4 = coord[dim_space*N4+1];
1163 double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1164 double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1165 double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1166 double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1168 xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1169 d1*b2 + a1*d2 - d1*a2);
1173 double z1 = coord[dim_space*N1+2];
1174 double z2 = coord[dim_space*N2+2];
1175 double z3 = coord[dim_space*N3+2];
1176 double z4 = coord[dim_space*N4+2];
1178 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1179 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1180 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1181 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1182 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1183 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1184 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1185 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1186 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1187 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1188 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1189 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1192 area[index] = xarea ;
1198 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1202 if (!onAll) delete [] global_connectivity ;
1205 Area->setValue(MED_FULL_INTERLACE,area);
1210 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1212 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1215 // Support must be on 1D elements
1217 // Make sure that the MESH class is the same as the MESH class attribut
1218 // in the class Support
1219 MESH* myMesh = Support->getMesh();
1221 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1223 int dim_space = getSpaceDimension();
1224 medEntityMesh support_entity = Support->getEntity();
1225 bool onAll = Support->isOnAllElements();
1227 int nb_type, length_values;
1228 const medGeometryElement* types;
1230 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1231 const int* global_connectivity;
1235 // nb_type = myMesh->getNumberOfTypes(support_entity);
1236 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1237 // types = getTypes(support_entity);
1241 nb_type = Support->getNumberOfTypes();
1242 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1243 types = Support->getTypes();
1247 FIELD<double>* Length;
1249 Length = new FIELD<double>(Support,1);
1250 // double *length = new double [length_values];
1251 Length->setValueType(MED_REEL64);
1253 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1254 MEDARRAY<double> * length = Length->getvalue();
1256 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1259 for (int i=0;i<nb_type;i++)
1261 medGeometryElement type = types[i] ;
1266 nb_entity_type = getNumberOfElements(support_entity,type);
1267 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1271 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1273 nb_entity_type = Support->getNumberOfElements(type);
1275 const int * supp_number = Support->getNumber(type);
1276 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1277 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1278 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1280 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1281 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1282 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1285 global_connectivity = global_connectivity_tmp ;
1290 case MED_SEG2 : case MED_SEG3 :
1292 for (int edge=0;edge<nb_entity_type;edge++)
1294 int edge_index = (type%100)*edge;
1296 int N1 = global_connectivity[edge_index]-1;
1297 int N2 = global_connectivity[edge_index+1]-1;
1299 double x1 = coord[dim_space*N1];
1300 double x2 = coord[dim_space*N2];
1302 double y1 = coord[dim_space*N1+1];
1303 double y2 = coord[dim_space*N2+1];
1305 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1306 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1308 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1309 (z1 - z2)*(z1 - z2));
1311 length->setIJ(index,1,xlength) ;
1317 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1321 if (!onAll) delete [] global_connectivity ;
1327 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1329 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1332 // Support must be on 2D or 1D elements
1334 // Make sure that the MESH class is the same as the MESH class attribut
1335 // in the class Support
1336 MESH* myMesh = Support->getMesh();
1338 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1340 int dim_space = getSpaceDimension();
1341 medEntityMesh support_entity = Support->getEntity();
1342 bool onAll = Support->isOnAllElements();
1344 int nb_type, length_values;
1345 const medGeometryElement* types;
1347 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1348 const int* global_connectivity;
1352 // nb_type = myMesh->getNumberOfTypes(support_entity);
1353 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1354 // types = getTypes(support_entity);
1358 nb_type = Support->getNumberOfTypes();
1359 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1360 types = Support->getTypes();
1365 FIELD<double>* Normal = new FIELD<double>(Support,dim_space);
1366 Normal->setName("NORMAL");
1367 Normal->setDescription("cells or faces normal");
1368 for (int k=1;k<=dim_space;k++) {
1369 Normal->setComponentName(k,"normal");
1370 Normal->setComponentDescription(k,"desc-comp");
1371 Normal->setMEDComponentUnit(k,"unit");
1374 Normal->setValueType(MED_REEL64);
1376 Normal->setIterationNumber(MED_NOPDT);
1377 Normal->setOrderNumber(MED_NONOR);
1378 Normal->setTime(0.0);
1380 double * normal = new double [dim_space*length_values];
1381 //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1383 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1386 for (int i=0;i<nb_type;i++)
1388 medGeometryElement type = types[i] ;
1390 // Make sure we trying to get Normals on
1391 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1392 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1394 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1395 (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1396 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1398 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1400 double xnormal1, xnormal2, xnormal3 ;
1404 nb_entity_type = getNumberOfElements(support_entity,type);
1405 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1409 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all for instance !"));
1410 nb_entity_type = Support->getNumberOfElements(type);
1412 const int * supp_number = Support->getNumber(type);
1413 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1414 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1415 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1417 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1418 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1419 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1423 global_connectivity = global_connectivity_tmp ;
1429 case MED_TRIA3 : case MED_TRIA6 :
1431 for (int tria=0;tria<nb_entity_type;tria++)
1433 int tria_index = (type%100)*tria;
1435 int N1 = global_connectivity[tria_index]-1;
1436 int N2 = global_connectivity[tria_index+1]-1;
1437 int N3 = global_connectivity[tria_index+2]-1;
1439 //double xarea; !! UNUSED VARIABLE !!
1440 double x1 = coord[dim_space*N1];
1441 double x2 = coord[dim_space*N2];
1442 double x3 = coord[dim_space*N3];
1444 double y1 = coord[dim_space*N1+1];
1445 double y2 = coord[dim_space*N2+1];
1446 double y3 = coord[dim_space*N3+1];
1448 double z1 = coord[dim_space*N1+2];
1449 double z2 = coord[dim_space*N2+2];
1450 double z3 = coord[dim_space*N3+2];
1452 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1453 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1454 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1456 normal[3*index] = xnormal1 ;
1457 normal[3*index+1] = xnormal2 ;
1458 normal[3*index+2] = xnormal3 ;
1464 case MED_QUAD4 : case MED_QUAD8 :
1466 for (int quad=0;quad<nb_entity_type;quad++)
1468 int quad_index = (type%100)*quad;
1470 int N1 = global_connectivity[quad_index]-1;
1471 int N2 = global_connectivity[quad_index+1]-1;
1472 int N3 = global_connectivity[quad_index+2]-1;
1473 int N4 = global_connectivity[quad_index+3]-1;
1476 double x1 = coord[dim_space*N1];
1477 double x2 = coord[dim_space*N2];
1478 double x3 = coord[dim_space*N3];
1479 double x4 = coord[dim_space*N4];
1481 double y1 = coord[dim_space*N1+1];
1482 double y2 = coord[dim_space*N2+1];
1483 double y3 = coord[dim_space*N3+1];
1484 double y4 = coord[dim_space*N4+1];
1486 double z1 = coord[dim_space*N1+2];
1487 double z2 = coord[dim_space*N2+2];
1488 double z3 = coord[dim_space*N3+2];
1489 double z4 = coord[dim_space*N4+2];
1491 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1492 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1493 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1494 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1497 xnormal1 = xnormal1/xarea;
1498 xnormal2 = xnormal2/xarea;
1499 xnormal3 = xnormal3/xarea;
1501 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1502 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1503 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1504 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1505 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1506 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1507 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1508 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1509 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1510 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1511 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1512 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1514 // sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1515 // ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1516 // ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1517 // ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1518 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1519 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1521 xnormal1 = xnormal1*xarea;
1522 xnormal2 = xnormal2*xarea;
1523 xnormal3 = xnormal3*xarea;
1525 normal[3*index] = xnormal1 ;
1526 normal[3*index+1] = xnormal2 ;
1527 normal[3*index+2] = xnormal3 ;
1533 case MED_SEG2 : case MED_SEG3 :
1535 for (int edge=0;edge<nb_entity_type;edge++)
1537 int edge_index = (type%100)*edge;
1539 int N1 = global_connectivity[edge_index]-1;
1540 int N2 = global_connectivity[edge_index+1]-1;
1542 double x1 = coord[dim_space*N1];
1543 double x2 = coord[dim_space*N2];
1545 double y1 = coord[dim_space*N1+1];
1546 double y2 = coord[dim_space*N2+1];
1548 xnormal1 = -(y2-y1);
1551 normal[2*index] = xnormal1 ;
1552 normal[2*index+1] = xnormal2 ;
1559 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1563 if (!onAll) delete [] global_connectivity ;
1566 Normal->setValue(MED_FULL_INTERLACE,normal);
1574 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1576 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1579 // Make sure that the MESH class is the same as the MESH class attribut
1580 // in the class Support
1581 MESH* myMesh = Support->getMesh();
1583 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1585 int dim_space = getSpaceDimension();
1586 medEntityMesh support_entity = Support->getEntity();
1587 bool onAll = Support->isOnAllElements();
1589 int nb_type, length_values;
1590 const medGeometryElement* types;
1592 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1593 const int* global_connectivity;
1597 // nb_type = myMesh->getNumberOfTypes(support_entity);
1598 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1599 // types = getTypes(support_entity);
1603 nb_type = Support->getNumberOfTypes();
1604 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1605 types = Support->getTypes();
1609 FIELD<double>* Barycenter;
1611 Barycenter = new FIELD<double>(Support,dim_space);
1612 Barycenter->setName("BARYCENTER");
1613 Barycenter->setDescription("cells or faces barycenter");
1615 for (int k=0;k<dim_space;k++) {
1617 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1618 Barycenter->setComponentDescription(kp1,"desc-comp");
1619 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1622 Barycenter->setValueType(MED_REEL64);
1624 Barycenter->setIterationNumber(0);
1625 Barycenter->setOrderNumber(0);
1626 Barycenter->setTime(0.0);
1628 double *barycenter = new double [dim_space*length_values];
1629 // double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1631 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1634 for (int i=0;i<nb_type;i++)
1636 medGeometryElement type = types[i] ;
1637 double xbarycenter1, xbarycenter2, xbarycenter3;
1641 nb_entity_type = getNumberOfElements(support_entity,type);
1642 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1646 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1647 nb_entity_type = Support->getNumberOfElements(type);
1649 const int * supp_number = Support->getNumber(type);
1650 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1651 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1652 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1654 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1655 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1656 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1659 global_connectivity = global_connectivity_tmp;
1664 case MED_TETRA4 : case MED_TETRA10 :
1666 for (int tetra =0;tetra<nb_entity_type;tetra++)
1668 int tetra_index = (type%100)*tetra;
1670 int N1 = global_connectivity[tetra_index]-1;
1671 int N2 = global_connectivity[tetra_index+1]-1;
1672 int N3 = global_connectivity[tetra_index+2]-1;
1673 int N4 = global_connectivity[tetra_index+3]-1;
1675 double x1 = coord[dim_space*N1];
1676 double x2 = coord[dim_space*N2];
1677 double x3 = coord[dim_space*N3];
1678 double x4 = coord[dim_space*N4];
1680 double y1 = coord[dim_space*N1+1];
1681 double y2 = coord[dim_space*N2+1];
1682 double y3 = coord[dim_space*N3+1];
1683 double y4 = coord[dim_space*N4+1];
1685 double z1 = coord[dim_space*N1+2];
1686 double z2 = coord[dim_space*N2+2];
1687 double z3 = coord[dim_space*N3+2];
1688 double z4 = coord[dim_space*N4+2];
1690 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1691 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1692 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1693 barycenter[3*index] = xbarycenter1 ;
1694 barycenter[3*index+1] = xbarycenter2 ;
1695 barycenter[3*index+2] = xbarycenter3 ;
1700 case MED_PYRA5 : case MED_PYRA13 :
1702 for (int pyra=0;pyra<nb_entity_type;pyra++)
1704 int pyra_index = (type%100)*pyra;
1706 int N1 = global_connectivity[pyra_index]-1;
1707 int N2 = global_connectivity[pyra_index+1]-1;
1708 int N3 = global_connectivity[pyra_index+2]-1;
1709 int N4 = global_connectivity[pyra_index+3]-1;
1710 int N5 = global_connectivity[pyra_index+4]-1;
1712 double x1 = coord[dim_space*N1];
1713 double x2 = coord[dim_space*N2];
1714 double x3 = coord[dim_space*N3];
1715 double x4 = coord[dim_space*N4];
1716 double x5 = coord[dim_space*N5];
1718 double y1 = coord[dim_space*N1+1];
1719 double y2 = coord[dim_space*N2+1];
1720 double y3 = coord[dim_space*N3+1];
1721 double y4 = coord[dim_space*N4+1];
1722 double y5 = coord[dim_space*N5+1];
1724 double z1 = coord[dim_space*N1+2];
1725 double z2 = coord[dim_space*N2+2];
1726 double z3 = coord[dim_space*N3+2];
1727 double z4 = coord[dim_space*N4+2];
1728 double z5 = coord[dim_space*N5+2];
1730 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1731 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1732 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1733 barycenter[3*index] = xbarycenter1 ;
1734 barycenter[3*index+1] = xbarycenter2 ;
1735 barycenter[3*index+2] = xbarycenter3 ;
1740 case MED_PENTA6 : case MED_PENTA15 :
1742 for (int penta=0;penta<nb_entity_type;penta++)
1744 int penta_index = (type%100)*penta;
1746 int N1 = global_connectivity[penta_index]-1;
1747 int N2 = global_connectivity[penta_index+1]-1;
1748 int N3 = global_connectivity[penta_index+2]-1;
1749 int N4 = global_connectivity[penta_index+3]-1;
1750 int N5 = global_connectivity[penta_index+4]-1;
1751 int N6 = global_connectivity[penta_index+5]-1;
1753 double x1 = coord[dim_space*N1];
1754 double x2 = coord[dim_space*N2];
1755 double x3 = coord[dim_space*N3];
1756 double x4 = coord[dim_space*N4];
1757 double x5 = coord[dim_space*N5];
1758 double x6 = coord[dim_space*N6];
1760 double y1 = coord[dim_space*N1+1];
1761 double y2 = coord[dim_space*N2+1];
1762 double y3 = coord[dim_space*N3+1];
1763 double y4 = coord[dim_space*N4+1];
1764 double y5 = coord[dim_space*N5+1];
1765 double y6 = coord[dim_space*N6+1];
1767 double z1 = coord[dim_space*N1+2];
1768 double z2 = coord[dim_space*N2+2];
1769 double z3 = coord[dim_space*N3+2];
1770 double z4 = coord[dim_space*N4+2];
1771 double z5 = coord[dim_space*N5+2];
1772 double z6 = coord[dim_space*N6+2];
1774 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1775 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1776 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1777 barycenter[3*index] = xbarycenter1 ;
1778 barycenter[3*index+1] = xbarycenter2 ;
1779 barycenter[3*index+2] = xbarycenter3 ;
1784 case MED_HEXA8 : case MED_HEXA20 :
1786 for (int hexa=0;hexa<nb_entity_type;hexa++)
1788 int hexa_index = (type%100)*hexa;
1790 int N1 = global_connectivity[hexa_index]-1;
1791 int N2 = global_connectivity[hexa_index+1]-1;
1792 int N3 = global_connectivity[hexa_index+2]-1;
1793 int N4 = global_connectivity[hexa_index+3]-1;
1794 int N5 = global_connectivity[hexa_index+4]-1;
1795 int N6 = global_connectivity[hexa_index+5]-1;
1796 int N7 = global_connectivity[hexa_index+6]-1;
1797 int N8 = global_connectivity[hexa_index+7]-1;
1799 double x1 = coord[dim_space*N1];
1800 double x2 = coord[dim_space*N2];
1801 double x3 = coord[dim_space*N3];
1802 double x4 = coord[dim_space*N4];
1803 double x5 = coord[dim_space*N5];
1804 double x6 = coord[dim_space*N6];
1805 double x7 = coord[dim_space*N7];
1806 double x8 = coord[dim_space*N8];
1808 double y1 = coord[dim_space*N1+1];
1809 double y2 = coord[dim_space*N2+1];
1810 double y3 = coord[dim_space*N3+1];
1811 double y4 = coord[dim_space*N4+1];
1812 double y5 = coord[dim_space*N5+1];
1813 double y6 = coord[dim_space*N6+1];
1814 double y7 = coord[dim_space*N7+1];
1815 double y8 = coord[dim_space*N8+1];
1817 double z1 = coord[dim_space*N1+2];
1818 double z2 = coord[dim_space*N2+2];
1819 double z3 = coord[dim_space*N3+2];
1820 double z4 = coord[dim_space*N4+2];
1821 double z5 = coord[dim_space*N5+2];
1822 double z6 = coord[dim_space*N6+2];
1823 double z7 = coord[dim_space*N7+2];
1824 double z8 = coord[dim_space*N8+2];
1826 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1827 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1828 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1830 barycenter[3*index] = xbarycenter1 ;
1831 barycenter[3*index+1] = xbarycenter2 ;
1832 barycenter[3*index+2] = xbarycenter3 ;
1838 case MED_TRIA3 : case MED_TRIA6 :
1840 for (int tria=0;tria<nb_entity_type;tria++)
1842 int tria_index = (type%100)*tria;
1844 int N1 = global_connectivity[tria_index]-1;
1845 int N2 = global_connectivity[tria_index+1]-1;
1846 int N3 = global_connectivity[tria_index+2]-1;
1848 double x1 = coord[dim_space*N1];
1849 double x2 = coord[dim_space*N2];
1850 double x3 = coord[dim_space*N3];
1852 double y1 = coord[dim_space*N1+1];
1853 double y2 = coord[dim_space*N2+1];
1854 double y3 = coord[dim_space*N3+1];
1856 xbarycenter1 = (x1 + x2 + x3)/3.0;
1857 xbarycenter2 = (y1 + y2 + y3)/3.0;
1861 barycenter[2*index] = xbarycenter1 ;
1862 barycenter[2*index+1] = xbarycenter2 ;
1867 coord[dim_space*N1+2];
1869 coord[dim_space*N2+2];
1871 coord[dim_space*N3+2];
1873 xbarycenter3 = (z1 + z2 + z3)/3.0;
1875 barycenter[3*index] = xbarycenter1 ;
1876 barycenter[3*index+1] = xbarycenter2 ;
1877 barycenter[3*index+2] = xbarycenter3 ;
1884 case MED_QUAD4 : case MED_QUAD8 :
1886 for (int quad=0;quad<nb_entity_type;quad++)
1888 int quad_index = (type%100)*quad;
1890 int N1 = global_connectivity[quad_index]-1;
1891 int N2 = global_connectivity[quad_index+1]-1;
1892 int N3 = global_connectivity[quad_index+2]-1;
1893 int N4 = global_connectivity[quad_index+3]-1;
1895 double x1 = coord[dim_space*N1];
1896 double x2 = coord[dim_space*N2];
1897 double x3 = coord[dim_space*N3];
1898 double x4 = coord[dim_space*N4];
1900 double y1 = coord[dim_space*N1+1];
1901 double y2 = coord[dim_space*N2+1];
1902 double y3 = coord[dim_space*N3+1];
1903 double y4 = coord[dim_space*N4+1];
1905 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1906 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1910 barycenter[2*index] = xbarycenter1 ;
1911 barycenter[2*index+1] = xbarycenter2 ;
1915 double z1 = coord[dim_space*N1+2];
1916 double z2 = coord[dim_space*N2+2];
1917 double z3 = coord[dim_space*N3+2];
1918 double z4 = coord[dim_space*N4+2];
1920 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1922 barycenter[3*index] = xbarycenter1 ;
1923 barycenter[3*index+1] = xbarycenter2 ;
1924 barycenter[3*index+2] = xbarycenter3 ;
1930 case MED_SEG2 : case MED_SEG3 :
1932 for (int edge=0;edge<nb_entity_type;edge++)
1934 int edge_index = (type%100)*edge;
1936 int N1 = global_connectivity[edge_index]-1;
1937 int N2 = global_connectivity[edge_index+1]-1;
1939 double x1 = coord[dim_space*N1];
1940 double x2 = coord[dim_space*N2];
1942 double y1 = coord[dim_space*N1+1];
1943 double y2 = coord[dim_space*N2+1];
1945 xbarycenter1 = (x1 + x2)/2.0;
1946 xbarycenter2 = (y1 + y2)/2.0;
1950 barycenter[2*index] = xbarycenter1 ;
1951 barycenter[2*index+1] = xbarycenter2 ;
1955 double z1 = coord[dim_space*N1+2];
1956 double z2 = coord[dim_space*N2+2];
1958 xbarycenter3 = (z1 + z2)/2.0;
1960 barycenter[3*index] = xbarycenter1 ;
1961 barycenter[3*index+1] = xbarycenter2 ;
1962 barycenter[3*index+2] = xbarycenter3 ;
1969 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
1973 if (!onAll) delete [] global_connectivity ;
1976 Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
1978 delete[] barycenter ;
1985 //=======================================================================
1986 //function : checkGridFillCoords
1987 //purpose : if this->_isAGrid, assure that _coordinate is filled
1988 //=======================================================================
1990 // inline void MESH::checkGridFillCoords() const
1993 // ((GRID *) this)->fillCoordinates();
1996 //=======================================================================
1997 //function : checkGridFillConnectivity
1998 //purpose : if this->_isAGrid, assure that _connectivity is filled
1999 //=======================================================================
2001 // inline void MESH::checkGridFillConnectivity() const
2004 // ((GRID *) this)->fillConnectivity();
2007 bool MESH::isEmpty() const
2009 bool notempty = _name != "NOT DEFINED" || _coordinate != NULL || _connectivity != NULL ||
2010 _spaceDimension !=MED_INVALID || _meshDimension !=MED_INVALID ||
2011 _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
2012 _familyNode.size() != 0 || _groupCell.size() != 0 ||
2013 _familyCell.size() != 0 || _groupFace.size() != 0 ||
2014 _familyFace.size() != 0 || _groupEdge.size() != 0 ||
2015 _familyEdge.size() != 0 || _isAGrid != 0 ;
2019 void MESH::read(int index)
2021 const char * LOC = "MESH::read(int index=0) : ";
2024 if (_drivers[index]) {
2025 _drivers[index]->open();
2026 _drivers[index]->read();
2027 _drivers[index]->close();
2030 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2031 << "The index given is invalid, index must be between 0 and |"
2036 // ((GRID *) this)->fillMeshAfterRead();
2040 //=======================================================================
2041 //function : getSkin
2043 //=======================================================================
2045 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2047 const char * LOC = "MESH::getSkin : " ;
2050 if (this != Support3D->getMesh())
2051 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
2052 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2053 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2055 // well, all rigth !
2056 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2057 mySupport->setAll(false);
2059 list<int> myElementsList ;
2062 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2063 if (Support3D->isOnAllElements())
2065 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2066 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2067 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2069 int cellNb1 = myConnectivityValue [i];
2070 int cellNb2 = myConnectivityValue [i+1];
2071 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2072 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2074 myElementsList.push_back( j ) ;
2081 map<int,int> FaceNbEncounterNb;
2083 int * myConnectivityValue = const_cast <int *>
2084 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2085 MED_CELL, MED_ALL_ELEMENTS));
2086 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2087 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2088 int nbCells = Support3D->getnumber()->getLength();
2089 for (i=0; i<nbCells; ++i)
2091 int cellNb = myCellNbs [ i ];
2092 int faceFirst = myConnectivityIndex[ cellNb-1 ];
2093 int faceLast = myConnectivityIndex[ cellNb ];
2094 for (j = faceFirst; j < faceLast; ++j)
2096 int faceNb = abs( myConnectivityValue [ j-1 ] );
2097 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2098 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2099 FaceNbEncounterNb[ faceNb ] = 1;
2101 FaceNbEncounterNb[ faceNb ] ++;
2104 map<int,int>::iterator FaceNbEncounterNbItr;
2105 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2106 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2107 FaceNbEncounterNbItr ++)
2108 if ((*FaceNbEncounterNbItr).second == 1)
2110 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2114 // Well, we must know how many geometric type we have found
2115 int * myListArray = new int[size] ;
2117 list<int>::iterator myElementsListIt ;
2118 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2119 myListArray[id]=(*myElementsListIt) ;
2123 int numberOfGeometricType ;
2124 medGeometryElement* geometricType ;
2125 int * numberOfGaussPoint ;
2126 int * geometricTypeNumber ;
2127 int * numberOfEntities ;
2128 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2129 int * mySkyLineArrayIndex ;
2131 int numberOfType = getNumberOfTypes(MED_FACE) ;
2132 if (numberOfType == 1) { // wonderfull : it's easy !
2133 numberOfGeometricType = 1 ;
2134 geometricType = new medGeometryElement[1] ;
2135 const medGeometryElement * allType = getTypes(MED_FACE);
2136 geometricType[0] = allType[0] ;
2137 numberOfGaussPoint = new int[1] ;
2138 numberOfGaussPoint[0] = 1 ;
2139 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2140 geometricTypeNumber[0] = 0 ;
2141 numberOfEntities = new int[1] ;
2142 numberOfEntities[0] = size ;
2143 mySkyLineArrayIndex = new int[2] ;
2144 mySkyLineArrayIndex[0]=1 ;
2145 mySkyLineArrayIndex[1]=1+size ;
2148 map<medGeometryElement,int> theType ;
2149 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2150 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2151 if (theType.find(myType) != theType.end() )
2152 theType[myType]+=1 ;
2156 numberOfGeometricType = theType.size() ;
2157 geometricType = new medGeometryElement[numberOfGeometricType] ;
2158 //const medGeometryElement * allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
2159 numberOfGaussPoint = new int[numberOfGeometricType] ;
2160 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2161 numberOfEntities = new int[numberOfGeometricType] ;
2162 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2164 mySkyLineArrayIndex[0]=1 ;
2165 map<medGeometryElement,int>::iterator theTypeIt ;
2166 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2167 geometricType[index] = (*theTypeIt).first ;
2168 numberOfGaussPoint[index] = 1 ;
2169 geometricTypeNumber[index] = 0 ;
2170 numberOfEntities[index] = (*theTypeIt).second ;
2171 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2175 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2176 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2178 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2179 mySupport->setGeometricType(geometricType) ;
2180 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2181 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2182 mySupport->setNumberOfElements(numberOfEntities) ;
2183 mySupport->setTotalNumberOfElements(size) ;
2184 mySupport->setNumber(mySkyLineArray) ;
2186 delete[] numberOfEntities;
2187 delete[] geometricTypeNumber;
2188 delete[] numberOfGaussPoint;
2189 delete[] geometricType;
2190 delete[] mySkyLineArrayIndex;
2191 delete[] myListArray;
2192 // delete mySkyLineArray;
2200 return a SUPPORT pointer on the union of all SUPPORTs in Supports.
2201 You should delete this pointer after use to avoid memory leaks.
2203 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2205 const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
2208 SUPPORT * returnedSupport;
2209 string returnedSupportName;
2210 string returnedSupportDescription;
2211 char * returnedSupportNameChar;
2212 char * returnedSupportDescriptionChar;
2213 int size = Supports.size();
2217 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2218 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2220 returnedSupport = new SUPPORT(*obj);
2222 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2223 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2225 returnedSupportNameChar = new char[lenName];
2226 returnedSupportDescriptionChar = new char[lenDescription];
2228 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2229 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2230 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2231 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2232 (Supports[0]->getDescription()).c_str());
2234 returnedSupportName = string(returnedSupportNameChar);
2235 returnedSupportDescription = string(returnedSupportDescriptionChar);
2237 returnedSupport->setName(returnedSupportName);
2238 returnedSupport->setDescription(returnedSupportDescription);
2240 delete [] returnedSupportNameChar;
2241 delete [] returnedSupportDescriptionChar;
2245 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2246 returnedSupport = new SUPPORT(*obj);
2248 int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
2249 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
2251 for(int i = 1;i<size;i++)
2253 obj = const_cast <SUPPORT *> (Supports[i]);
2254 returnedSupport->blending(obj);
2258 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2259 lenDescription = lenDescription + 5 +
2260 strlen((Supports[i]->getDescription()).c_str());
2264 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2265 lenDescription = lenDescription + 2 +
2266 strlen((Supports[i]->getDescription()).c_str());
2270 returnedSupportNameChar = new char[lenName];
2271 returnedSupportDescriptionChar = new char[lenDescription];
2273 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
2274 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
2276 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2277 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2278 (Supports[0]->getDescription()).c_str());
2280 for(int i = 1;i<size;i++)
2284 returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
2285 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
2287 returnedSupportNameChar = strcat(returnedSupportNameChar,
2288 (Supports[i]->getName()).c_str());
2289 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2290 (Supports[i]->getDescription()).c_str());
2294 returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
2295 returnedSupportNameChar = strcat(returnedSupportNameChar,
2296 (Supports[i]->getName()).c_str());
2298 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
2299 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2300 (Supports[i]->getDescription()).c_str());
2304 returnedSupportName = string(returnedSupportNameChar);
2305 returnedSupport->setName(returnedSupportName);
2307 returnedSupportDescription = string(returnedSupportDescriptionChar);
2308 returnedSupport->setDescription(returnedSupportDescription);
2310 delete [] returnedSupportNameChar;
2311 delete [] returnedSupportDescriptionChar;
2315 return returnedSupport;
2319 return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
2320 The (SUPPORT *) NULL pointer is returned if the intersection is empty.
2321 You should delete this pointer after use to avois memory leaks.
2323 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2325 const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
2328 SUPPORT * returnedSupport;
2329 string returnedSupportName;
2330 string returnedSupportDescription;
2331 char * returnedSupportNameChar;
2332 char * returnedSupportDescriptionChar;
2333 int size = Supports.size();
2337 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2338 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2340 returnedSupport = new SUPPORT(*obj);
2342 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2343 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2345 returnedSupportNameChar = new char[lenName];
2346 returnedSupportDescriptionChar = new char[lenDescription];
2348 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2349 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2350 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2351 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2352 (Supports[0]->getDescription()).c_str());
2354 returnedSupportName = string(returnedSupportNameChar);
2355 returnedSupportDescription = string(returnedSupportDescriptionChar);
2357 returnedSupport->setName(returnedSupportName);
2358 returnedSupport->setDescription(returnedSupportDescription);
2360 delete [] returnedSupportNameChar;
2361 delete [] returnedSupportDescriptionChar;
2365 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2366 returnedSupport = new SUPPORT(*obj);
2368 int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
2369 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
2371 for(int i = 1;i<size;i++)
2373 obj = const_cast <SUPPORT *> (Supports[i]);
2374 returnedSupport->intersecting(obj);
2378 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2379 lenDescription = lenDescription + 5 +
2380 strlen((Supports[i]->getDescription()).c_str());
2384 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2385 lenDescription = lenDescription + 2 +
2386 strlen((Supports[i]->getDescription()).c_str());
2390 if(returnedSupport != (SUPPORT *) NULL)
2392 returnedSupportNameChar = new char[lenName];
2393 returnedSupportDescriptionChar = new char[lenDescription];
2395 returnedSupportNameChar = strcpy(returnedSupportNameChar,
2396 "Intersection of ");
2397 returnedSupportDescriptionChar =
2398 strcpy(returnedSupportDescriptionChar,"Intersection of ");
2400 returnedSupportNameChar = strcat(returnedSupportNameChar,
2401 (Supports[0]->getName()).c_str());
2402 returnedSupportDescriptionChar =
2403 strcat(returnedSupportDescriptionChar,
2404 (Supports[0]->getDescription()).c_str());
2406 for(int i = 1;i<size;i++)
2410 returnedSupportNameChar = strcat(returnedSupportNameChar,
2412 returnedSupportDescriptionChar =
2413 strcat(returnedSupportDescriptionChar," and ");
2415 returnedSupportNameChar =
2416 strcat(returnedSupportNameChar,
2417 (Supports[i]->getName()).c_str());
2418 returnedSupportDescriptionChar =
2419 strcat(returnedSupportDescriptionChar,
2420 (Supports[i]->getDescription()).c_str());
2424 returnedSupportNameChar = strcat(returnedSupportNameChar,
2426 returnedSupportNameChar =
2427 strcat(returnedSupportNameChar,
2428 (Supports[i]->getName()).c_str());
2430 returnedSupportDescriptionChar =
2431 strcat(returnedSupportDescriptionChar,", ");
2432 returnedSupportDescriptionChar =
2433 strcat(returnedSupportDescriptionChar,
2434 (Supports[i]->getDescription()).c_str());
2438 returnedSupportName = string(returnedSupportNameChar);
2439 returnedSupport->setName(returnedSupportName);
2441 returnedSupportDescription = string(returnedSupportDescriptionChar);
2442 returnedSupport->setDescription(returnedSupportDescription);
2444 delete [] returnedSupportNameChar;
2445 delete [] returnedSupportDescriptionChar;
2450 return returnedSupport;
2454 // internal helper type
2458 std::vector<int> groups;
2459 MED_EN::medGeometryElement geometricType;
2462 // Create families from groups
2463 void MESH::createFamilies()
2465 int nbFam=0; // count the families we create, used as identifier ???????????
2467 int idFamNode = 0; // identifier for node families
2468 int idFamElement = 0; // identifier for cell, face or edge families
2470 // Main loop on mesh's entities
2471 for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2473 int numberofgroups = getNumberOfGroups(entity);
2475 continue; // no groups for this entity
2477 vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2479 // make myFamilies points to the member corresponding to entity
2480 vector<FAMILY*>* myFamilies;
2484 myFamilies = & _familyCell;
2487 myFamilies = & _familyFace;
2490 myFamilies = & _familyEdge;
2493 myFamilies = & _familyNode;
2497 vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2498 // get a copy of the (old) family ptrs before clearing
2499 vector<FAMILY*> myOldFamilies=getFamilies(entity);
2500 myFamilies->clear();
2503 // 1 - Create a vector containing for each cell (of the entity) an information structure
2504 // giving geometric type and the groups it belong to
2506 med_int numberOfTypes=getNumberOfTypes(entity);
2507 const int * index=getGlobalNumberingIndex(entity);
2508 const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2509 med_int numberOfCells=index[numberOfTypes]-1; // total number of cells for that entity
2510 SCRUTE(numberOfTypes);
2511 SCRUTE(numberOfCells);
2512 vector< _cell > tab_cell(numberOfCells);
2513 for(med_int t=0; t!=numberOfTypes; ++t)
2514 for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2515 tab_cell[n].geometricType=geometricTypes[t];
2518 // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2520 for (unsigned g=0; g!=myGroups.size(); ++g)
2522 // scan cells that belongs to the group
2523 const int* groupCells=myGroups[g]->getnumber()->getValue();
2524 int nbCells=myGroups[g]->getnumber()->getLength();
2525 for(int c=0; c!=nbCells; ++c)
2526 tab_cell[groupCells[c]-1].groups.push_back(g);
2530 // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2531 // whith the vector of contained cells
2533 map< string,vector<int> > tab_families;
2534 map< string,vector<int> >::iterator fam;
2535 for(int n=0; n!=numberOfCells; ++n)
2537 ostringstream key; // to generate the name of the family
2539 if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2540 key << "_NONE" << entity;
2542 for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2544 string groupName=myGroups[*it]->getName();
2545 if(groupName.empty())
2548 key << "_" << groupName;
2551 tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2552 /* fam = tab_families.find(key.str());
2553 if( fam != tab_families.end())
2554 fam->second.push_back(n+1); // +1 : convention Fortran de MED
2557 vector<int> newfamily;
2558 newfamily.push_back(n+1); // +1 : convention Fortran de MED
2559 tab_families.insert(make_pair(key.str(),newfamily));
2564 // 4 - Scan the family map, create MED Families, check if it already exist.
2566 for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2568 vector<medGeometryElement> tab_types_geometriques;
2569 medGeometryElement geometrictype=MED_NONE;
2570 vector<int> tab_index_types_geometriques;
2571 vector<int> tab_nombres_elements;
2573 // scan family cells and fill the tab that are needed by the create a MED FAMILY
2574 for( int i=0; i!=fam->second.size(); ++i)
2576 int ncell=fam->second[i]-1;
2577 if(tab_cell[ncell].geometricType != geometrictype)
2579 // new geometric type -> we store it and complete index tabs
2580 if(!tab_index_types_geometriques.empty())
2581 tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2582 tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2583 tab_index_types_geometriques.push_back(i+1);
2586 // store and complete index tabs for the last geometric type
2587 tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2588 tab_index_types_geometriques.push_back(fam->second.size()+1);
2590 // create a empty MED FAMILY and fill it with the tabs we constructed
2591 FAMILY* newFam = new FAMILY();
2592 newFam->setTotalNumberOfElements(fam->second.size());
2593 newFam->setName(fam->first);
2594 newFam->setMesh(this);
2595 newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2596 newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2597 newFam->setNumberOfElements(&tab_nombres_elements[0]);
2598 newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2599 newFam->setEntity(entity);
2600 newFam->setAll(false);
2612 idFam = -idFamElement;
2616 idFam = -idFamElement;
2620 idFam = -idFamElement;
2624 newFam->setIdentifier(idFam);
2626 // Update links between families and groups
2628 int ncell1=fam->second[0]-1; // number of first cell in family
2629 int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2632 newFam->setNumberOfGroups(numberOfGroups);
2633 string * groupNames=new string[numberOfGroups];
2635 // iterate on the groups the family belongs to
2636 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2637 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2639 whichFamilyInGroup[*it].push_back(newFam);
2640 groupNames[ng]=myGroups[*it]->getName();
2642 newFam->setGroupsNames(groupNames);
2645 int sizeOfFamVect = myFamilies->size();
2647 MESSAGE(" MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2649 myFamilies->push_back(newFam);
2652 // delete old families
2653 for (int i=0;i<myOldFamilies.size();i++)
2654 delete myOldFamilies[i] ;
2656 // update references in groups
2657 for (int i=0;i<myGroups.size();i++)
2659 myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2660 myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2663 // re-scan the cells vector, and fill the family vector with cells.
2664 // creation of support, check if it already exist.