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"
22 using namespace MEDMEM;
23 //#include "MEDMEM_Grid.hxx" this inclision should have never be here !!!
25 //update Families with content list
26 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
28 // ------- Drivers Management Part
30 // MESH::INSTANCE_DE<MED_MESH_RDONLY_DRIVER> MESH::inst_med_rdonly ;
31 //const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med_rdonly , &MESH::inst_med_rdwr } ;
33 // Add a similar line for your personnal driver (step 3)
35 MESH::INSTANCE_DE<MED_MESH_RDWR_DRIVER> MESH::inst_med ;
36 MESH::INSTANCE_DE<GIBI_MESH_RDWR_DRIVER> MESH::inst_gibi ;
37 MESH::INSTANCE_DE<PORFLOW_MESH_RDWR_DRIVER> MESH::inst_porflow ;
38 MESH::INSTANCE_DE<VTK_MESH_DRIVER> MESH::inst_vtk;
40 // Add your own driver in the driver list (step 4)
41 // Note the list must be coherent with the driver type list defined in MEDMEM_DRIVER.hxx.
42 const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med, &MESH::inst_gibi, &MESH::inst_porflow, &MESH::inst_vtk } ;
44 /*! Add a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName. The meshname used in the file
45 is driverName. addDriver returns an integer handler. */
46 int MESH::addDriver(driverTypes driverType,
47 const string & fileName/*="Default File Name.med"*/,const string & driverName/*="Default Mesh Name"*/) {
49 const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\") : ";
53 int itDriver = (int) NO_DRIVER;
59 SCRUTE(instances[driverType]);
64 itDriver = (int) driverType ;
69 itDriver = (int) driverType ;
73 case PORFLOW_DRIVER : {
74 itDriver = (int) driverType ;
84 throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "NO_DRIVER has been specified to the method which is not allowed"));
88 if (itDriver == ((int) NO_DRIVER))
89 throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "othe driver than MED_DRIVER GIBI_DRIVER PORFLOW_DRIVER and VT_DRIVER has been specified to the method which is not allowed"));
91 driver = instances[itDriver]->run(fileName, this) ;
93 _drivers.push_back(driver);
95 current = _drivers.size()-1;
97 _drivers[current]->setMeshName(driverName);
104 /*! Add an existing MESH driver. */
105 int MESH::addDriver(GENDRIVER & driver) {
106 const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
109 // A faire : Vérifier que le driver est de type MESH.
110 GENDRIVER * newDriver = driver.copy() ;
112 _drivers.push_back(newDriver);
113 return _drivers.size()-1;
118 /*! Remove an existing MESH driver. */
119 void MESH::rmDriver (int index/*=0*/) {
120 const char * LOC = "MESH::rmDriver (int index=0): ";
123 if ( _drivers[index] ) {
124 //_drivers.erase(&_drivers[index]);
126 MESSAGE ("detruire");
129 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
130 << "The index given is invalid, index must be between 0 and |"
139 // ------ End of Drivers Management Part
144 const char * LOC = "MESH::init(): ";
148 string _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
150 string _decription = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
152 _coordinate = (COORDINATE *) NULL;
153 _connectivity = (CONNECTIVITY *) NULL;
155 _spaceDimension = MED_INVALID; // 0 ?!?
156 _meshDimension = MED_INVALID;
157 _numberOfNodes = MED_INVALID;
161 _arePresentOptionnalNodesNumbers = 0;
166 /*! Create an empty MESH. */
167 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
174 _description = m._description;
175 _isAGrid = m._isAGrid;
177 if (m._coordinate != NULL)
178 _coordinate = new COORDINATE(* m._coordinate);
180 _coordinate = (COORDINATE *) NULL;
182 if (m._connectivity != NULL)
183 _connectivity = new CONNECTIVITY(* m._connectivity);
185 _connectivity = (CONNECTIVITY *) NULL;
187 _spaceDimension = m._spaceDimension;
188 _meshDimension = m._meshDimension;
189 _numberOfNodes = m._numberOfNodes;
191 _familyNode = m._familyNode;
192 for (int i=0; i<m._familyNode.size(); i++)
194 _familyNode[i] = new FAMILY(* m._familyNode[i]);
195 _familyNode[i]->setMesh(this);
198 _familyCell = m._familyCell;
199 for (int i=0; i<m._familyCell.size(); i++)
201 _familyCell[i] = new FAMILY(* m._familyCell[i]);
202 _familyCell[i]->setMesh(this);
205 _familyFace = m._familyFace;
206 for (int i=0; i<m._familyFace.size(); i++)
208 _familyFace[i] = new FAMILY(* m._familyFace[i]);
209 _familyFace[i]->setMesh(this);
212 _familyEdge = m._familyEdge;
213 for (int i=0; i<m._familyEdge.size(); i++)
215 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
216 _familyEdge[i]->setMesh(this);
219 _groupNode = m._groupNode;
220 for (int i=0; i<m._groupNode.size(); i++)
222 _groupNode[i] = new GROUP(* m._groupNode[i]);
223 _groupNode[i]->setMesh(this);
226 _groupCell = m._groupCell;
227 for (int i=0; i<m._groupCell.size(); i++)
229 _groupCell[i] = new GROUP(* m._groupCell[i]);
230 _groupCell[i]->setMesh(this);
233 _groupFace = m._groupFace;
234 for (int i=0; i<m._groupFace.size(); i++)
236 _groupFace[i] = new GROUP(* m._groupFace[i]);
237 _groupFace[i]->setMesh(this);
240 _groupEdge = m._groupEdge;
241 for (int i=0; i<m._groupEdge.size(); i++)
243 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
244 _groupEdge[i]->setMesh(this);
247 //_drivers = m._drivers; //Recopie des drivers?
252 MESSAGE("MESH::~MESH() : Destroying the Mesh");
253 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
254 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
256 size = _familyNode.size() ;
257 for (int i=0;i<size;i++)
258 delete _familyNode[i] ;
259 size = _familyCell.size() ;
260 for (int i=0;i<size;i++)
261 delete _familyCell[i] ;
262 size = _familyFace.size() ;
263 for (int i=0;i<size;i++)
264 delete _familyFace[i] ;
265 size = _familyEdge.size() ;
266 for (int i=0;i<size;i++)
267 delete _familyEdge[i] ;
268 size = _groupNode.size() ;
269 for (int i=0;i<size;i++)
270 delete _groupNode[i] ;
271 size = _groupCell.size() ;
272 for (int i=0;i<size;i++)
273 delete _groupCell[i] ;
274 size = _groupFace.size() ;
275 for (int i=0;i<size;i++)
276 delete _groupFace[i] ;
277 size = _groupEdge.size() ;
278 for (int i=0;i<size;i++)
279 delete _groupEdge[i] ;
281 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
283 for (unsigned int index=0; index < _drivers.size(); index++ )
285 SCRUTE(_drivers[index]);
286 if ( _drivers[index] != NULL) delete _drivers[index];
291 MESH & MESH::operator=(const MESH &m)
293 const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
296 MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
299 // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
300 // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
301 // _drivers = m._drivers;
303 // space_dimension=m.space_dimension;
304 // mesh_dimension=m.mesh_dimension;
306 // nodes_count=m.nodes_count;
307 // coordinates=m.coordinates;
308 // coordinates_name=m.coordinates_name ;
309 // coordinates_unit=m.coordinates_unit ;
310 // nodes_numbers=m.nodes_numbers ;
312 // cells_types_count=m.cells_types_count;
313 // cells_type=m.cells_type;
314 // cells_count=m.cells_count;
315 // nodal_connectivity=m.nodal_connectivity;
317 // nodes_families_count=m.nodes_families_count;
318 // nodes_Families=m.nodes_Families;
320 // cells_families_count=m.cells_families_count;
321 // cells_Families=m.cells_Families;
323 // maximum_cell_number_by_node = m.maximum_cell_number_by_node;
324 // if (maximum_cell_number_by_node > 0)
326 // reverse_nodal_connectivity = m.reverse_nodal_connectivity;
327 // reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
334 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
335 The meshname driverName must already exists in the file.*/
336 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) throw (MEDEXCEPTION)
338 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
349 MED_MESH_RDONLY_DRIVER myDriver(fileName,this);
350 myDriver.setMeshName(driverName);
351 current = addDriver(myDriver);
356 GIBI_MESH_RDONLY_DRIVER myDriver(fileName,this);
357 current = addDriver(myDriver);
360 case PORFLOW_DRIVER :
362 PORFLOW_MESH_RDONLY_DRIVER myDriver(fileName,this);
363 current = addDriver(myDriver);
368 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
373 // current = addDriver(driverType,fileName,driverName);
374 // switch(_drivers[current]->getAccessMode() ) {
375 // case MED_WRONLY : {
376 // MESSAGE("MESH::MESH(driverTypes driverType, .....) : driverType must not be a MED_WRONLY accessMode");
377 // rmDriver(current);
382 _drivers[current]->open();
383 _drivers[current]->read();
384 _drivers[current]->close();
387 // ((GRID *) this)->fillMeshAfterRead();
393 // Node MESH::Donne_Barycentre(const Maille &m) const
395 // Node N=node[m[0]];
396 // for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
397 // N=N*(1.0/m.donne_nbr_sommets());
401 ostream & MEDMEM::operator<<(ostream &os, const MESH &myMesh)
403 int spacedimension = myMesh.getSpaceDimension();
404 int meshdimension = myMesh.getMeshDimension();
405 int numberofnodes = myMesh.getNumberOfNodes();
407 os << "Space Dimension : " << spacedimension << endl << endl;
409 os << "Mesh Dimension : " << meshdimension << endl << endl;
411 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
412 os << "SHOW NODES COORDINATES : " << endl;
414 os << "Name :" << endl;
415 const string * coordinatesnames = myMesh.getCoordinatesNames();
416 for(int i=0; i<spacedimension ; i++)
418 os << " - " << coordinatesnames[i] << endl;
420 os << "Unit :" << endl;
421 const string * coordinatesunits = myMesh.getCoordinatesUnits();
422 for(int i=0; i<spacedimension ; i++)
424 os << " - " << coordinatesunits[i] << endl;
426 for(int i=0; i<numberofnodes ; i++)
428 os << "Nodes " << i+1 << " : ";
429 for (int j=0; j<spacedimension ; j++)
430 os << coordinates[i*spacedimension+j] << " ";
434 os << endl << "SHOW CONNECTIVITY :" << endl;
435 os << *myMesh._connectivity << endl;
437 medEntityMesh entity;
438 os << endl << "SHOW FAMILIES :" << endl << endl;
439 for (int k=1; k<=4; k++)
441 if (k==1) entity = MED_NODE;
442 if (k==2) entity = MED_CELL;
443 if (k==3) entity = MED_FACE;
444 if (k==4) entity = MED_EDGE;
445 int numberoffamilies = myMesh.getNumberOfFamilies(entity);
446 using namespace MED_FR;
447 os << "NumberOfFamilies on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberoffamilies<<endl;
448 using namespace MED_EN;
449 for (int i=1; i<numberoffamilies+1;i++)
451 os << * myMesh.getFamily(entity,i) << endl;
455 os << endl << "SHOW GROUPS :" << endl << endl;
456 for (int k=1; k<=4; k++)
458 if (k==1) entity = MED_NODE;
459 if (k==2) entity = MED_CELL;
460 if (k==3) entity = MED_FACE;
461 if (k==4) entity = MED_EDGE;
462 int numberofgroups = myMesh.getNumberOfGroups(entity);
463 using namespace MED_FR;
464 os << "NumberOfGroups on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberofgroups<<endl;
465 using namespace MED_EN;
466 for (int i=1; i<numberofgroups+1;i++)
468 os << * myMesh.getGroup(entity,i) << endl;
476 Get global number of element which have same connectivity than connectivity argument.
478 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
480 Return -1 if not found.
482 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
484 const char* LOC="MESH::getElementNumber " ;
487 int numberOfValue ; // size of connectivity array
488 CELLMODEL myType(Type) ;
489 if (ConnectivityType==MED_DESCENDING)
490 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
492 numberOfValue = myType.getNumberOfNodes() ; // nodes
494 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
495 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
497 // First node or face/edge
498 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
499 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
501 // list to put cells numbers
502 list<int> cellsList ;
503 list<int>::iterator itList ;
504 for (int i=indexBegin; i<indexEnd; i++)
505 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
507 // Foreach node or face/edge in cell, we search cells which are in common.
508 // At the end, we must have only one !
510 for (int i=1; i<numberOfValue; i++) {
511 int connectivity_i = connectivity[i] ;
512 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
514 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
515 if ((*itList)==myReverseConnectivityValue[j-1]) {
521 itList=cellsList.erase(itList);
522 itList--; // well : rigth if itList = cellsList.begin() ??
527 if (cellsList.size()>1) // we have more than one cell
528 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
530 if (cellsList.size()==0)
535 return cellsList.front() ;
540 Return a support which reference all elements on the boundary of mesh.
542 For instance, we could get only face in 3D and edge in 2D.
544 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
547 const char * LOC = "MESH::getBoundaryElements : " ;
550 // actually we could only get face (in 3D) and edge (in 2D)
551 if (_spaceDimension == 3)
552 if (Entity != MED_FACE)
553 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
554 if (_spaceDimension == 2)
555 if (Entity != MED_EDGE)
556 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
559 SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
560 //mySupport.setDescription("boundary of type ...");
561 mySupport->setAll(false);
564 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
565 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
566 int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
567 list<int> myElementsList ;
569 for (int i=0 ; i<numberOf; i++)
570 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
571 myElementsList.push_back(i+1) ;
574 // Well, we must know how many geometric type we have found
575 int * myListArray = new int[size] ;
577 list<int>::iterator myElementsListIt ;
578 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
579 myListArray[id]=(*myElementsListIt) ;
583 int numberOfGeometricType ;
584 medGeometryElement* geometricType ;
585 int * numberOfGaussPoint ;
586 int * geometricTypeNumber ;
587 int * numberOfElements ;
588 //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
589 int * mySkyLineArrayIndex ;
591 int numberOfType = getNumberOfTypes(Entity) ;
592 if (numberOfType == 1) { // wonderfull : it's easy !
593 numberOfGeometricType = 1 ;
594 geometricType = new medGeometryElement[1] ;
595 const medGeometryElement * allType = getTypes(Entity);
596 geometricType[0] = allType[0] ;
597 numberOfGaussPoint = new int[1] ;
598 numberOfGaussPoint[0] = 1 ;
599 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
600 geometricTypeNumber[0] = 0 ;
601 numberOfElements = new int[1] ;
602 numberOfElements[0] = size ;
603 mySkyLineArrayIndex = new int[2] ;
604 mySkyLineArrayIndex[0]=1 ;
605 mySkyLineArrayIndex[1]=1+size ;
608 map<medGeometryElement,int> theType ;
609 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
610 medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
611 if (theType.find(myType) != theType.end() )
616 numberOfGeometricType = theType.size() ;
617 geometricType = new medGeometryElement[numberOfGeometricType] ;
618 //const medGeometryElement * allType = getTypes(Entity); !! UNUSZED VARIABLE !!
619 numberOfGaussPoint = new int[numberOfGeometricType] ;
620 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
621 numberOfElements = new int[numberOfGeometricType] ;
622 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
624 mySkyLineArrayIndex[0]=1 ;
625 map<medGeometryElement,int>::iterator theTypeIt ;
626 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
627 geometricType[index] = (*theTypeIt).first ;
628 numberOfGaussPoint[index] = 1 ;
629 geometricTypeNumber[index] = 0 ;
630 numberOfElements[index] = (*theTypeIt).second ;
631 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
635 //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
636 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
638 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
639 mySupport->setGeometricType(geometricType) ;
640 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
641 mySupport->setNumberOfElements(numberOfElements) ;
642 mySupport->setTotalNumberOfElements(size) ;
643 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
644 mySupport->setNumber(mySkyLineArray) ;
646 delete[] numberOfElements;
647 delete[] geometricTypeNumber;
648 delete[] numberOfGaussPoint;
649 delete[] geometricType;
650 delete[] mySkyLineArrayIndex;
651 delete[] myListArray;
652 // delete mySkyLineArray;
658 FIELD<double>* MESH::getVolume(const SUPPORT * Support) const throw (MEDEXCEPTION)
660 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
663 // Support must be on 3D elements
665 // Make sure that the MESH class is the same as the MESH class attribut
666 // in the class Support
667 MESH* myMesh = Support->getMesh();
669 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
671 int dim_space = getSpaceDimension();
672 medEntityMesh support_entity = Support->getEntity();
673 bool onAll = Support->isOnAllElements();
675 int nb_type, length_values;
676 const medGeometryElement* types;
678 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
679 const int* global_connectivity;
683 // nb_type = myMesh->getNumberOfTypes(support_entity);
684 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
685 // types = getTypes(support_entity);
689 nb_type = Support->getNumberOfTypes();
690 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
691 types = Support->getTypes();
695 FIELD<double>* Volume = new FIELD<double>::FIELD(Support,1);
696 // double *volume = new double [length_values];
697 Volume->setName("VOLUME");
698 Volume->setDescription("cells volume");
699 Volume->setComponentName(1,"volume");
700 Volume->setComponentDescription(1,"desc-comp");
702 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
704 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
706 Volume->setMEDComponentUnit(1,MEDComponentUnit);
708 Volume->setValueType(MED_REEL64);
710 Volume->setIterationNumber(0);
711 Volume->setOrderNumber(0);
712 Volume->setTime(0.0);
714 //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
715 MEDARRAY<double> *volume = Volume->getvalue();
718 const double * coord = getCoordinates(MED_FULL_INTERLACE);
720 for (int i=0;i<nb_type;i++)
722 medGeometryElement type = types[i] ;
727 nb_entity_type = getNumberOfElements(support_entity,type);
728 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
732 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all"));
734 nb_entity_type = Support->getNumberOfElements(type);
736 const int * supp_number = Support->getNumber(type);
737 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
738 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
739 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
741 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
742 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
743 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
746 global_connectivity = global_connectivity_tmp ;
751 case MED_TETRA4 : case MED_TETRA10 :
753 for (int tetra=0;tetra<nb_entity_type;tetra++)
755 int tetra_index = (type%100)*tetra;
757 int N1 = global_connectivity[tetra_index]-1;
758 int N2 = global_connectivity[tetra_index+1]-1;
759 int N3 = global_connectivity[tetra_index+2]-1;
760 int N4 = global_connectivity[tetra_index+3]-1;
762 double x1 = coord[dim_space*N1];
763 double x2 = coord[dim_space*N2];
764 double x3 = coord[dim_space*N3];
765 double x4 = coord[dim_space*N4];
767 double y1 = coord[dim_space*N1+1];
768 double y2 = coord[dim_space*N2+1];
769 double y3 = coord[dim_space*N3+1];
770 double y4 = coord[dim_space*N4+1];
772 double z1 = coord[dim_space*N1+2];
773 double z2 = coord[dim_space*N2+2];
774 double z3 = coord[dim_space*N3+2];
775 double z4 = coord[dim_space*N4+2];
777 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
778 (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
779 (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
781 //volume[index] = xvolume ;
782 volume->setIJ(index,1,xvolume) ;
787 case MED_PYRA5 : case MED_PYRA13 :
789 for (int pyra=0;pyra<nb_entity_type;pyra++)
791 int pyra_index = (type%100)*pyra;
793 int N1 = global_connectivity[pyra_index]-1;
794 int N2 = global_connectivity[pyra_index+1]-1;
795 int N3 = global_connectivity[pyra_index+2]-1;
796 int N4 = global_connectivity[pyra_index+3]-1;
797 int N5 = global_connectivity[pyra_index+4]-1;
799 double x1 = coord[dim_space*N1];
800 double x2 = coord[dim_space*N2];
801 double x3 = coord[dim_space*N3];
802 double x4 = coord[dim_space*N4];
803 double x5 = coord[dim_space*N5];
805 double y1 = coord[dim_space*N1+1];
806 double y2 = coord[dim_space*N2+1];
807 double y3 = coord[dim_space*N3+1];
808 double y4 = coord[dim_space*N4+1];
809 double y5 = coord[dim_space*N5+1];
811 double z1 = coord[dim_space*N1+2];
812 double z2 = coord[dim_space*N2+2];
813 double z3 = coord[dim_space*N3+2];
814 double z4 = coord[dim_space*N4+2];
815 double z5 = coord[dim_space*N5+2];
817 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
818 (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
819 (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
820 ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
821 (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
822 (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
825 //volume[index] = xvolume ;
826 volume->setIJ(index,1,xvolume) ;
831 case MED_PENTA6 : case MED_PENTA15 :
833 for (int penta=0;penta<nb_entity_type;penta++)
835 int penta_index = (type%100)*penta;
837 int N1 = global_connectivity[penta_index]-1;
838 int N2 = global_connectivity[penta_index+1]-1;
839 int N3 = global_connectivity[penta_index+2]-1;
840 int N4 = global_connectivity[penta_index+3]-1;
841 int N5 = global_connectivity[penta_index+4]-1;
842 int N6 = global_connectivity[penta_index+5]-1;
844 double x1 = coord[dim_space*N1];
845 double x2 = coord[dim_space*N2];
846 double x3 = coord[dim_space*N3];
847 double x4 = coord[dim_space*N4];
848 double x5 = coord[dim_space*N5];
849 double x6 = coord[dim_space*N6];
851 double y1 = coord[dim_space*N1+1];
852 double y2 = coord[dim_space*N2+1];
853 double y3 = coord[dim_space*N3+1];
854 double y4 = coord[dim_space*N4+1];
855 double y5 = coord[dim_space*N5+1];
856 double y6 = coord[dim_space*N6+1];
858 double z1 = coord[dim_space*N1+2];
859 double z2 = coord[dim_space*N2+2];
860 double z3 = coord[dim_space*N3+2];
861 double z4 = coord[dim_space*N4+2];
862 double z5 = coord[dim_space*N5+2];
863 double z6 = coord[dim_space*N6+2];
865 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
866 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
867 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
868 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
869 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
870 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
871 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
873 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
875 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
877 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
878 (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
879 (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
880 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
882 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
884 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
885 (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
886 (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
887 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
889 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
891 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
892 (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
893 (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
895 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
897 //volume[index] = xvolume ;
898 volume->setIJ(index,1,xvolume) ;
903 case MED_HEXA8 : case MED_HEXA20 :
905 for (int hexa=0;hexa<nb_entity_type;hexa++)
907 int hexa_index = (type%100)*hexa;
909 int N1 = global_connectivity[hexa_index]-1;
910 int N2 = global_connectivity[hexa_index+1]-1;
911 int N3 = global_connectivity[hexa_index+2]-1;
912 int N4 = global_connectivity[hexa_index+3]-1;
913 int N5 = global_connectivity[hexa_index+4]-1;
914 int N6 = global_connectivity[hexa_index+5]-1;
915 int N7 = global_connectivity[hexa_index+6]-1;
916 int N8 = global_connectivity[hexa_index+7]-1;
918 double x1 = coord[dim_space*N1];
919 double x2 = coord[dim_space*N2];
920 double x3 = coord[dim_space*N3];
921 double x4 = coord[dim_space*N4];
922 double x5 = coord[dim_space*N5];
923 double x6 = coord[dim_space*N6];
924 double x7 = coord[dim_space*N7];
925 double x8 = coord[dim_space*N8];
927 double y1 = coord[dim_space*N1+1];
928 double y2 = coord[dim_space*N2+1];
929 double y3 = coord[dim_space*N3+1];
930 double y4 = coord[dim_space*N4+1];
931 double y5 = coord[dim_space*N5+1];
932 double y6 = coord[dim_space*N6+1];
933 double y7 = coord[dim_space*N7+1];
934 double y8 = coord[dim_space*N8+1];
936 double z1 = coord[dim_space*N1+2];
937 double z2 = coord[dim_space*N2+2];
938 double z3 = coord[dim_space*N3+2];
939 double z4 = coord[dim_space*N4+2];
940 double z5 = coord[dim_space*N5+2];
941 double z6 = coord[dim_space*N6+2];
942 double z7 = coord[dim_space*N7+2];
943 double z8 = coord[dim_space*N8+2];
945 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
946 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
947 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
948 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
949 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
950 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
951 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
952 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
953 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
954 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
955 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
956 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
958 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
960 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
962 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
963 (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
964 (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
965 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
967 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
969 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
970 (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
971 (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
972 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
973 (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
974 (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
975 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
976 (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
977 (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
978 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
979 ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
980 ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
981 ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
982 ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
983 ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
984 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
986 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
988 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
989 (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
990 (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
991 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
993 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
995 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
996 (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
997 (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
998 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
999 (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
1000 (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
1001 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
1002 (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
1003 (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
1004 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
1005 ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
1006 ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
1007 ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
1008 ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
1009 ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
1010 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
1011 (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
1012 (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
1013 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
1014 (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
1015 (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
1016 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
1017 ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
1018 ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
1019 ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
1020 ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
1021 ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
1022 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
1023 (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
1024 (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
1025 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
1026 (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
1027 (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
1028 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
1029 ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
1030 ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
1031 ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
1032 ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
1033 ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
1034 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
1035 (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
1036 (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
1037 (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
1038 (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
1039 (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
1040 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
1041 (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
1042 (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
1043 (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
1044 (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
1045 (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
1046 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
1047 (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
1048 ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
1049 (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
1050 ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
1051 (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
1052 ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
1053 (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
1054 ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
1055 (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
1056 ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
1057 (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
1059 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
1060 4.0*(C + F + G + H + L + O + P + Q + S + T +
1061 V + W) + 2.0*(I + R + U + X + Y + Z) +
1064 //volume[index] = xvolume ;
1065 volume->setIJ(index,1,xvolume) ;
1071 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
1075 if (!onAll) delete [] global_connectivity ;
1081 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1083 const char * LOC = "MESH::getArea(SUPPORT*) : ";
1086 // Support must be on 2D elements
1088 // Make sure that the MESH class is the same as the MESH class attribut
1089 // in the class Support
1090 MESH* myMesh = Support->getMesh();
1092 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1094 int dim_space = getSpaceDimension();
1095 medEntityMesh support_entity = Support->getEntity();
1096 bool onAll = Support->isOnAllElements();
1098 int nb_type, length_values;
1099 const medGeometryElement* types;
1101 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1102 const int* global_connectivity;
1106 // nb_type = myMesh->getNumberOfTypes(support_entity);
1107 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1108 // types = getTypes(support_entity);
1112 nb_type = Support->getNumberOfTypes();
1113 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1114 types = Support->getTypes();
1118 FIELD<double>* Area;
1120 Area = new FIELD<double>::FIELD(Support,1);
1121 Area->setName("AREA");
1122 Area->setDescription("cells or faces area");
1124 Area->setComponentName(1,"area");
1125 Area->setComponentDescription(1,"desc-comp");
1127 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1129 string MEDComponentUnit="trucmuch";
1131 Area->setMEDComponentUnit(1,MEDComponentUnit);
1133 Area->setValueType(MED_REEL64);
1135 Area->setIterationNumber(0);
1136 Area->setOrderNumber(0);
1139 double *area = new double[length_values];
1140 //double *area = Area->getValue(MED_FULL_INTERLACE);
1142 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1145 for (int i=0;i<nb_type;i++)
1147 medGeometryElement type = types[i] ;
1152 nb_entity_type = getNumberOfElements(support_entity,type);
1153 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1157 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1159 nb_entity_type = Support->getNumberOfElements(type);
1161 const int * supp_number = Support->getNumber(type);
1162 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1163 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1164 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1166 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1167 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1168 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1172 global_connectivity = global_connectivity_tmp ;
1178 case MED_TRIA3 : case MED_TRIA6 :
1180 for (int tria=0;tria<nb_entity_type;tria++)
1182 int tria_index = (type%100)*tria;
1184 int N1 = global_connectivity[tria_index]-1;
1185 int N2 = global_connectivity[tria_index+1]-1;
1186 int N3 = global_connectivity[tria_index+2]-1;
1188 double x1 = coord[dim_space*N1];
1189 double x2 = coord[dim_space*N2];
1190 double x3 = coord[dim_space*N3];
1192 double y1 = coord[dim_space*N1+1];
1193 double y2 = coord[dim_space*N2+1];
1194 double y3 = coord[dim_space*N3+1];
1198 xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1202 double z1 = coord[dim_space*N1+2];
1203 double z2 = coord[dim_space*N2+2];
1204 double z3 = coord[dim_space*N3+2];
1206 xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1207 ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1208 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1209 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1210 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1211 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1214 area[index] = xarea ;
1219 case MED_QUAD4 : case MED_QUAD8 :
1221 for (int quad=0;quad<nb_entity_type;quad++)
1223 int quad_index = (type%100)*quad;
1225 int N1 = global_connectivity[quad_index]-1;
1226 int N2 = global_connectivity[quad_index+1]-1;
1227 int N3 = global_connectivity[quad_index+2]-1;
1228 int N4 = global_connectivity[quad_index+3]-1;
1230 double x1 = coord[dim_space*N1];
1231 double x2 = coord[dim_space*N2];
1232 double x3 = coord[dim_space*N3];
1233 double x4 = coord[dim_space*N4];
1235 double y1 = coord[dim_space*N1+1];
1236 double y2 = coord[dim_space*N2+1];
1237 double y3 = coord[dim_space*N3+1];
1238 double y4 = coord[dim_space*N4+1];
1242 double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1243 double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1244 double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1245 double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1247 xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1248 d1*b2 + a1*d2 - d1*a2);
1252 double z1 = coord[dim_space*N1+2];
1253 double z2 = coord[dim_space*N2+2];
1254 double z3 = coord[dim_space*N3+2];
1255 double z4 = coord[dim_space*N4+2];
1257 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1258 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1259 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1260 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1261 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1262 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1263 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1264 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1265 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1266 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1267 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1268 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1271 area[index] = xarea ;
1277 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1281 if (!onAll) delete [] global_connectivity ;
1284 Area->setValue(MED_FULL_INTERLACE,area);
1289 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1291 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1294 // Support must be on 1D elements
1296 // Make sure that the MESH class is the same as the MESH class attribut
1297 // in the class 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 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1310 const int* global_connectivity;
1314 // nb_type = myMesh->getNumberOfTypes(support_entity);
1315 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1316 // types = getTypes(support_entity);
1320 nb_type = Support->getNumberOfTypes();
1321 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1322 types = Support->getTypes();
1326 FIELD<double>* Length;
1328 Length = new FIELD<double>::FIELD(Support,1);
1329 // double *length = new double [length_values];
1330 Length->setValueType(MED_REEL64);
1332 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1333 MEDARRAY<double> * length = Length->getvalue();
1335 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1338 for (int i=0;i<nb_type;i++)
1340 medGeometryElement type = types[i] ;
1345 nb_entity_type = getNumberOfElements(support_entity,type);
1346 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1350 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1352 nb_entity_type = Support->getNumberOfElements(type);
1354 const int * supp_number = Support->getNumber(type);
1355 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1356 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1357 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1359 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1360 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1361 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1364 global_connectivity = global_connectivity_tmp ;
1369 case MED_SEG2 : case MED_SEG3 :
1371 for (int edge=0;edge<nb_entity_type;edge++)
1373 int edge_index = (type%100)*edge;
1375 int N1 = global_connectivity[edge_index]-1;
1376 int N2 = global_connectivity[edge_index+1]-1;
1378 double x1 = coord[dim_space*N1];
1379 double x2 = coord[dim_space*N2];
1381 double y1 = coord[dim_space*N1+1];
1382 double y2 = coord[dim_space*N2+1];
1384 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1385 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1387 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1388 (z1 - z2)*(z1 - z2));
1390 length->setIJ(index,1,xlength) ;
1396 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1400 if (!onAll) delete [] global_connectivity ;
1406 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1408 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1411 // Support must be on 2D or 1D elements
1413 // Make sure that the MESH class is the same as the MESH class attribut
1414 // in the class Support
1415 MESH* myMesh = Support->getMesh();
1417 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1419 int dim_space = getSpaceDimension();
1420 medEntityMesh support_entity = Support->getEntity();
1421 bool onAll = Support->isOnAllElements();
1423 int nb_type, length_values;
1424 const medGeometryElement* types;
1426 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1427 const int* global_connectivity;
1431 // nb_type = myMesh->getNumberOfTypes(support_entity);
1432 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1433 // types = getTypes(support_entity);
1437 nb_type = Support->getNumberOfTypes();
1438 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1439 types = Support->getTypes();
1444 FIELD<double>* Normal = new FIELD<double>::FIELD(Support,dim_space);
1445 Normal->setName("NORMAL");
1446 Normal->setDescription("cells or faces normal");
1447 for (int k=1;k<=dim_space;k++) {
1448 Normal->setComponentName(k,"normal");
1449 Normal->setComponentDescription(k,"desc-comp");
1450 Normal->setMEDComponentUnit(k,"unit");
1453 Normal->setValueType(MED_REEL64);
1455 Normal->setIterationNumber(MED_NOPDT);
1456 Normal->setOrderNumber(MED_NONOR);
1457 Normal->setTime(0.0);
1459 double * normal = new double [dim_space*length_values];
1460 //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1462 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1465 for (int i=0;i<nb_type;i++)
1467 medGeometryElement type = types[i] ;
1469 // Make sure we trying to get Normals on
1470 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1471 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1473 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1474 (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1475 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1477 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1479 double xnormal1, xnormal2, xnormal3 ;
1483 nb_entity_type = getNumberOfElements(support_entity,type);
1484 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1488 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all for instance !"));
1489 nb_entity_type = Support->getNumberOfElements(type);
1491 const int * supp_number = Support->getNumber(type);
1492 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1493 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1494 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1496 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1497 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1498 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1502 global_connectivity = global_connectivity_tmp ;
1508 case MED_TRIA3 : case MED_TRIA6 :
1510 for (int tria=0;tria<nb_entity_type;tria++)
1512 int tria_index = (type%100)*tria;
1514 int N1 = global_connectivity[tria_index]-1;
1515 int N2 = global_connectivity[tria_index+1]-1;
1516 int N3 = global_connectivity[tria_index+2]-1;
1518 //double xarea; !! UNUSED VARIABLE !!
1519 double x1 = coord[dim_space*N1];
1520 double x2 = coord[dim_space*N2];
1521 double x3 = coord[dim_space*N3];
1523 double y1 = coord[dim_space*N1+1];
1524 double y2 = coord[dim_space*N2+1];
1525 double y3 = coord[dim_space*N3+1];
1527 double z1 = coord[dim_space*N1+2];
1528 double z2 = coord[dim_space*N2+2];
1529 double z3 = coord[dim_space*N3+2];
1531 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1532 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1533 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1535 normal[3*index] = xnormal1 ;
1536 normal[3*index+1] = xnormal2 ;
1537 normal[3*index+2] = xnormal3 ;
1543 case MED_QUAD4 : case MED_QUAD8 :
1545 for (int quad=0;quad<nb_entity_type;quad++)
1547 int quad_index = (type%100)*quad;
1549 int N1 = global_connectivity[quad_index]-1;
1550 int N2 = global_connectivity[quad_index+1]-1;
1551 int N3 = global_connectivity[quad_index+2]-1;
1552 int N4 = global_connectivity[quad_index+3]-1;
1555 double x1 = coord[dim_space*N1];
1556 double x2 = coord[dim_space*N2];
1557 double x3 = coord[dim_space*N3];
1558 double x4 = coord[dim_space*N4];
1560 double y1 = coord[dim_space*N1+1];
1561 double y2 = coord[dim_space*N2+1];
1562 double y3 = coord[dim_space*N3+1];
1563 double y4 = coord[dim_space*N4+1];
1565 double z1 = coord[dim_space*N1+2];
1566 double z2 = coord[dim_space*N2+2];
1567 double z3 = coord[dim_space*N3+2];
1568 double z4 = coord[dim_space*N4+2];
1570 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1571 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1572 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1573 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1576 xnormal1 = xnormal1/xarea;
1577 xnormal2 = xnormal2/xarea;
1578 xnormal3 = xnormal3/xarea;
1580 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1581 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1582 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1583 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1584 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1585 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1586 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1587 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1588 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1589 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1590 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1591 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1593 // sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1594 // ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1595 // ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1596 // ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1597 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1598 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1600 xnormal1 = xnormal1*xarea;
1601 xnormal2 = xnormal2*xarea;
1602 xnormal3 = xnormal3*xarea;
1604 normal[3*index] = xnormal1 ;
1605 normal[3*index+1] = xnormal2 ;
1606 normal[3*index+2] = xnormal3 ;
1612 case MED_SEG2 : case MED_SEG3 :
1614 for (int edge=0;edge<nb_entity_type;edge++)
1616 int edge_index = (type%100)*edge;
1618 int N1 = global_connectivity[edge_index]-1;
1619 int N2 = global_connectivity[edge_index+1]-1;
1621 double x1 = coord[dim_space*N1];
1622 double x2 = coord[dim_space*N2];
1624 double y1 = coord[dim_space*N1+1];
1625 double y2 = coord[dim_space*N2+1];
1627 xnormal1 = -(y2-y1);
1630 normal[2*index] = xnormal1 ;
1631 normal[2*index+1] = xnormal2 ;
1638 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1642 if (!onAll) delete [] global_connectivity ;
1645 Normal->setValue(MED_FULL_INTERLACE,normal);
1653 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1655 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1658 // Make sure that the MESH class is the same as the MESH class attribut
1659 // in the class Support
1660 MESH* myMesh = Support->getMesh();
1662 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1664 int dim_space = getSpaceDimension();
1665 medEntityMesh support_entity = Support->getEntity();
1666 bool onAll = Support->isOnAllElements();
1668 int nb_type, length_values;
1669 const medGeometryElement* types;
1671 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1672 const int* global_connectivity;
1676 // nb_type = myMesh->getNumberOfTypes(support_entity);
1677 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1678 // types = getTypes(support_entity);
1682 nb_type = Support->getNumberOfTypes();
1683 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1684 types = Support->getTypes();
1688 FIELD<double>* Barycenter;
1690 Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1691 Barycenter->setName("BARYCENTER");
1692 Barycenter->setDescription("cells or faces barycenter");
1694 for (int k=0;k<dim_space;k++) {
1696 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1697 Barycenter->setComponentDescription(kp1,"desc-comp");
1698 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1701 Barycenter->setValueType(MED_REEL64);
1703 Barycenter->setIterationNumber(0);
1704 Barycenter->setOrderNumber(0);
1705 Barycenter->setTime(0.0);
1707 double *barycenter = new double [dim_space*length_values];
1708 // double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1710 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1713 for (int i=0;i<nb_type;i++)
1715 medGeometryElement type = types[i] ;
1716 double xbarycenter1, xbarycenter2, xbarycenter3;
1720 nb_entity_type = getNumberOfElements(support_entity,type);
1721 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1725 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1726 nb_entity_type = Support->getNumberOfElements(type);
1728 const int * supp_number = Support->getNumber(type);
1729 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1730 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1731 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1733 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1734 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1735 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1738 global_connectivity = global_connectivity_tmp;
1743 case MED_TETRA4 : case MED_TETRA10 :
1745 for (int tetra =0;tetra<nb_entity_type;tetra++)
1747 int tetra_index = (type%100)*tetra;
1749 int N1 = global_connectivity[tetra_index]-1;
1750 int N2 = global_connectivity[tetra_index+1]-1;
1751 int N3 = global_connectivity[tetra_index+2]-1;
1752 int N4 = global_connectivity[tetra_index+3]-1;
1754 double x1 = coord[dim_space*N1];
1755 double x2 = coord[dim_space*N2];
1756 double x3 = coord[dim_space*N3];
1757 double x4 = coord[dim_space*N4];
1759 double y1 = coord[dim_space*N1+1];
1760 double y2 = coord[dim_space*N2+1];
1761 double y3 = coord[dim_space*N3+1];
1762 double y4 = coord[dim_space*N4+1];
1764 double z1 = coord[dim_space*N1+2];
1765 double z2 = coord[dim_space*N2+2];
1766 double z3 = coord[dim_space*N3+2];
1767 double z4 = coord[dim_space*N4+2];
1769 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1770 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1771 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1772 barycenter[3*index] = xbarycenter1 ;
1773 barycenter[3*index+1] = xbarycenter2 ;
1774 barycenter[3*index+2] = xbarycenter3 ;
1779 case MED_PYRA5 : case MED_PYRA13 :
1781 for (int pyra=0;pyra<nb_entity_type;pyra++)
1783 int pyra_index = (type%100)*pyra;
1785 int N1 = global_connectivity[pyra_index]-1;
1786 int N2 = global_connectivity[pyra_index+1]-1;
1787 int N3 = global_connectivity[pyra_index+2]-1;
1788 int N4 = global_connectivity[pyra_index+3]-1;
1789 int N5 = global_connectivity[pyra_index+4]-1;
1791 double x1 = coord[dim_space*N1];
1792 double x2 = coord[dim_space*N2];
1793 double x3 = coord[dim_space*N3];
1794 double x4 = coord[dim_space*N4];
1795 double x5 = coord[dim_space*N5];
1797 double y1 = coord[dim_space*N1+1];
1798 double y2 = coord[dim_space*N2+1];
1799 double y3 = coord[dim_space*N3+1];
1800 double y4 = coord[dim_space*N4+1];
1801 double y5 = coord[dim_space*N5+1];
1803 double z1 = coord[dim_space*N1+2];
1804 double z2 = coord[dim_space*N2+2];
1805 double z3 = coord[dim_space*N3+2];
1806 double z4 = coord[dim_space*N4+2];
1807 double z5 = coord[dim_space*N5+2];
1809 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1810 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1811 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1812 barycenter[3*index] = xbarycenter1 ;
1813 barycenter[3*index+1] = xbarycenter2 ;
1814 barycenter[3*index+2] = xbarycenter3 ;
1819 case MED_PENTA6 : case MED_PENTA15 :
1821 for (int penta=0;penta<nb_entity_type;penta++)
1823 int penta_index = (type%100)*penta;
1825 int N1 = global_connectivity[penta_index]-1;
1826 int N2 = global_connectivity[penta_index+1]-1;
1827 int N3 = global_connectivity[penta_index+2]-1;
1828 int N4 = global_connectivity[penta_index+3]-1;
1829 int N5 = global_connectivity[penta_index+4]-1;
1830 int N6 = global_connectivity[penta_index+5]-1;
1832 double x1 = coord[dim_space*N1];
1833 double x2 = coord[dim_space*N2];
1834 double x3 = coord[dim_space*N3];
1835 double x4 = coord[dim_space*N4];
1836 double x5 = coord[dim_space*N5];
1837 double x6 = coord[dim_space*N6];
1839 double y1 = coord[dim_space*N1+1];
1840 double y2 = coord[dim_space*N2+1];
1841 double y3 = coord[dim_space*N3+1];
1842 double y4 = coord[dim_space*N4+1];
1843 double y5 = coord[dim_space*N5+1];
1844 double y6 = coord[dim_space*N6+1];
1846 double z1 = coord[dim_space*N1+2];
1847 double z2 = coord[dim_space*N2+2];
1848 double z3 = coord[dim_space*N3+2];
1849 double z4 = coord[dim_space*N4+2];
1850 double z5 = coord[dim_space*N5+2];
1851 double z6 = coord[dim_space*N6+2];
1853 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1854 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1855 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1856 barycenter[3*index] = xbarycenter1 ;
1857 barycenter[3*index+1] = xbarycenter2 ;
1858 barycenter[3*index+2] = xbarycenter3 ;
1863 case MED_HEXA8 : case MED_HEXA20 :
1865 for (int hexa=0;hexa<nb_entity_type;hexa++)
1867 int hexa_index = (type%100)*hexa;
1869 int N1 = global_connectivity[hexa_index]-1;
1870 int N2 = global_connectivity[hexa_index+1]-1;
1871 int N3 = global_connectivity[hexa_index+2]-1;
1872 int N4 = global_connectivity[hexa_index+3]-1;
1873 int N5 = global_connectivity[hexa_index+4]-1;
1874 int N6 = global_connectivity[hexa_index+5]-1;
1875 int N7 = global_connectivity[hexa_index+6]-1;
1876 int N8 = global_connectivity[hexa_index+7]-1;
1878 double x1 = coord[dim_space*N1];
1879 double x2 = coord[dim_space*N2];
1880 double x3 = coord[dim_space*N3];
1881 double x4 = coord[dim_space*N4];
1882 double x5 = coord[dim_space*N5];
1883 double x6 = coord[dim_space*N6];
1884 double x7 = coord[dim_space*N7];
1885 double x8 = coord[dim_space*N8];
1887 double y1 = coord[dim_space*N1+1];
1888 double y2 = coord[dim_space*N2+1];
1889 double y3 = coord[dim_space*N3+1];
1890 double y4 = coord[dim_space*N4+1];
1891 double y5 = coord[dim_space*N5+1];
1892 double y6 = coord[dim_space*N6+1];
1893 double y7 = coord[dim_space*N7+1];
1894 double y8 = coord[dim_space*N8+1];
1896 double z1 = coord[dim_space*N1+2];
1897 double z2 = coord[dim_space*N2+2];
1898 double z3 = coord[dim_space*N3+2];
1899 double z4 = coord[dim_space*N4+2];
1900 double z5 = coord[dim_space*N5+2];
1901 double z6 = coord[dim_space*N6+2];
1902 double z7 = coord[dim_space*N7+2];
1903 double z8 = coord[dim_space*N8+2];
1905 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1906 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1907 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1909 barycenter[3*index] = xbarycenter1 ;
1910 barycenter[3*index+1] = xbarycenter2 ;
1911 barycenter[3*index+2] = xbarycenter3 ;
1917 case MED_TRIA3 : case MED_TRIA6 :
1919 for (int tria=0;tria<nb_entity_type;tria++)
1921 int tria_index = (type%100)*tria;
1923 int N1 = global_connectivity[tria_index]-1;
1924 int N2 = global_connectivity[tria_index+1]-1;
1925 int N3 = global_connectivity[tria_index+2]-1;
1927 double x1 = coord[dim_space*N1];
1928 double x2 = coord[dim_space*N2];
1929 double x3 = coord[dim_space*N3];
1931 double y1 = coord[dim_space*N1+1];
1932 double y2 = coord[dim_space*N2+1];
1933 double y3 = coord[dim_space*N3+1];
1935 xbarycenter1 = (x1 + x2 + x3)/3.0;
1936 xbarycenter2 = (y1 + y2 + y3)/3.0;
1940 barycenter[2*index] = xbarycenter1 ;
1941 barycenter[2*index+1] = xbarycenter2 ;
1946 coord[dim_space*N1+2];
1948 coord[dim_space*N2+2];
1950 coord[dim_space*N3+2];
1952 xbarycenter3 = (z1 + z2 + z3)/3.0;
1954 barycenter[3*index] = xbarycenter1 ;
1955 barycenter[3*index+1] = xbarycenter2 ;
1956 barycenter[3*index+2] = xbarycenter3 ;
1963 case MED_QUAD4 : case MED_QUAD8 :
1965 for (int quad=0;quad<nb_entity_type;quad++)
1967 int quad_index = (type%100)*quad;
1969 int N1 = global_connectivity[quad_index]-1;
1970 int N2 = global_connectivity[quad_index+1]-1;
1971 int N3 = global_connectivity[quad_index+2]-1;
1972 int N4 = global_connectivity[quad_index+3]-1;
1974 double x1 = coord[dim_space*N1];
1975 double x2 = coord[dim_space*N2];
1976 double x3 = coord[dim_space*N3];
1977 double x4 = coord[dim_space*N4];
1979 double y1 = coord[dim_space*N1+1];
1980 double y2 = coord[dim_space*N2+1];
1981 double y3 = coord[dim_space*N3+1];
1982 double y4 = coord[dim_space*N4+1];
1984 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1985 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1989 barycenter[2*index] = xbarycenter1 ;
1990 barycenter[2*index+1] = xbarycenter2 ;
1994 double z1 = coord[dim_space*N1+2];
1995 double z2 = coord[dim_space*N2+2];
1996 double z3 = coord[dim_space*N3+2];
1997 double z4 = coord[dim_space*N4+2];
1999 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
2001 barycenter[3*index] = xbarycenter1 ;
2002 barycenter[3*index+1] = xbarycenter2 ;
2003 barycenter[3*index+2] = xbarycenter3 ;
2009 case MED_SEG2 : case MED_SEG3 :
2011 for (int edge=0;edge<nb_entity_type;edge++)
2013 int edge_index = (type%100)*edge;
2015 int N1 = global_connectivity[edge_index]-1;
2016 int N2 = global_connectivity[edge_index+1]-1;
2018 double x1 = coord[dim_space*N1];
2019 double x2 = coord[dim_space*N2];
2021 double y1 = coord[dim_space*N1+1];
2022 double y2 = coord[dim_space*N2+1];
2024 xbarycenter1 = (x1 + x2)/2.0;
2025 xbarycenter2 = (y1 + y2)/2.0;
2029 barycenter[2*index] = xbarycenter1 ;
2030 barycenter[2*index+1] = xbarycenter2 ;
2034 double z1 = coord[dim_space*N1+2];
2035 double z2 = coord[dim_space*N2+2];
2037 xbarycenter3 = (z1 + z2)/2.0;
2039 barycenter[3*index] = xbarycenter1 ;
2040 barycenter[3*index+1] = xbarycenter2 ;
2041 barycenter[3*index+2] = xbarycenter3 ;
2048 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
2052 if (!onAll) delete [] global_connectivity ;
2055 Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
2057 delete[] barycenter ;
2064 //=======================================================================
2065 //function : checkGridFillCoords
2066 //purpose : if this->_isAGrid, assure that _coordinate is filled
2067 //=======================================================================
2069 // inline void MESH::checkGridFillCoords() const
2072 // ((GRID *) this)->fillCoordinates();
2075 //=======================================================================
2076 //function : checkGridFillConnectivity
2077 //purpose : if this->_isAGrid, assure that _connectivity is filled
2078 //=======================================================================
2080 // inline void MESH::checkGridFillConnectivity() const
2083 // ((GRID *) this)->fillConnectivity();
2086 bool MESH::isEmpty() const
2088 bool notempty = _name != "" || _description != "" ||
2089 _coordinate != NULL || _connectivity != NULL ||
2090 _spaceDimension !=MED_INVALID ||
2091 _meshDimension !=MED_INVALID ||
2092 _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
2093 _familyNode.size() != 0 || _groupCell.size() != 0 ||
2094 _familyCell.size() != 0 || _groupFace.size() != 0 ||
2095 _familyFace.size() != 0 || _groupEdge.size() != 0 ||
2096 _familyEdge.size() != 0 || _isAGrid != 0 ;
2100 void MESH::read(int index)
2102 const char * LOC = "MESH::read(int index=0) : ";
2105 if (_drivers[index]) {
2106 _drivers[index]->open();
2107 _drivers[index]->read();
2108 _drivers[index]->close();
2111 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2112 << "The index given is invalid, index must be between 0 and |"
2117 // ((GRID *) this)->fillMeshAfterRead();
2121 //=======================================================================
2122 //function : getSkin
2124 //=======================================================================
2126 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2128 const char * LOC = "MESH::getSkin : " ;
2131 if (this != Support3D->getMesh())
2132 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
2133 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2134 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2136 // well, all rigth !
2137 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2138 mySupport->setAll(false);
2140 list<int> myElementsList ;
2143 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2144 if (Support3D->isOnAllElements())
2146 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2147 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2148 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2150 int cellNb1 = myConnectivityValue [i];
2151 int cellNb2 = myConnectivityValue [i+1];
2152 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2153 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2155 myElementsList.push_back( j ) ;
2162 map<int,int> FaceNbEncounterNb;
2164 int * myConnectivityValue = const_cast <int *>
2165 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2166 MED_CELL, MED_ALL_ELEMENTS));
2167 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2168 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2169 int nbCells = Support3D->getnumber()->getLength();
2170 for (i=0; i<nbCells; ++i)
2172 int cellNb = myCellNbs [ i ];
2173 int faceFirst = myConnectivityIndex[ cellNb-1 ];
2174 int faceLast = myConnectivityIndex[ cellNb ];
2175 for (j = faceFirst; j < faceLast; ++j)
2177 int faceNb = abs( myConnectivityValue [ j-1 ] );
2178 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2179 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2180 FaceNbEncounterNb[ faceNb ] = 1;
2182 FaceNbEncounterNb[ faceNb ] ++;
2185 map<int,int>::iterator FaceNbEncounterNbItr;
2186 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2187 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2188 FaceNbEncounterNbItr ++)
2189 if ((*FaceNbEncounterNbItr).second == 1)
2191 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2195 // Well, we must know how many geometric type we have found
2196 int * myListArray = new int[size] ;
2198 list<int>::iterator myElementsListIt ;
2199 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2200 myListArray[id]=(*myElementsListIt) ;
2204 int numberOfGeometricType ;
2205 medGeometryElement* geometricType ;
2206 int * numberOfGaussPoint ;
2207 int * geometricTypeNumber ;
2208 int * numberOfEntities ;
2209 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2210 int * mySkyLineArrayIndex ;
2212 int numberOfType = getNumberOfTypes(MED_FACE) ;
2213 if (numberOfType == 1) { // wonderfull : it's easy !
2214 numberOfGeometricType = 1 ;
2215 geometricType = new medGeometryElement[1] ;
2216 const medGeometryElement * allType = getTypes(MED_FACE);
2217 geometricType[0] = allType[0] ;
2218 numberOfGaussPoint = new int[1] ;
2219 numberOfGaussPoint[0] = 1 ;
2220 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2221 geometricTypeNumber[0] = 0 ;
2222 numberOfEntities = new int[1] ;
2223 numberOfEntities[0] = size ;
2224 mySkyLineArrayIndex = new int[2] ;
2225 mySkyLineArrayIndex[0]=1 ;
2226 mySkyLineArrayIndex[1]=1+size ;
2229 map<medGeometryElement,int> theType ;
2230 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2231 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2232 if (theType.find(myType) != theType.end() )
2233 theType[myType]+=1 ;
2237 numberOfGeometricType = theType.size() ;
2238 geometricType = new medGeometryElement[numberOfGeometricType] ;
2239 //const medGeometryElement * allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
2240 numberOfGaussPoint = new int[numberOfGeometricType] ;
2241 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2242 numberOfEntities = new int[numberOfGeometricType] ;
2243 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2245 mySkyLineArrayIndex[0]=1 ;
2246 map<medGeometryElement,int>::iterator theTypeIt ;
2247 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2248 geometricType[index] = (*theTypeIt).first ;
2249 numberOfGaussPoint[index] = 1 ;
2250 geometricTypeNumber[index] = 0 ;
2251 numberOfEntities[index] = (*theTypeIt).second ;
2252 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2256 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2257 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2259 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2260 mySupport->setGeometricType(geometricType) ;
2261 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2262 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2263 mySupport->setNumberOfElements(numberOfEntities) ;
2264 mySupport->setTotalNumberOfElements(size) ;
2265 mySupport->setNumber(mySkyLineArray) ;
2267 delete[] numberOfEntities;
2268 delete[] geometricTypeNumber;
2269 delete[] numberOfGaussPoint;
2270 delete[] geometricType;
2271 delete[] mySkyLineArrayIndex;
2272 delete[] myListArray;
2273 // delete mySkyLineArray;
2281 return a SUPPORT pointer on the union of all SUPPORTs in Supports.
2282 You should delete this pointer after use to avoid memory leaks.
2284 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2286 const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
2289 SUPPORT * returnedSupport;
2290 string returnedSupportName;
2291 string returnedSupportDescription;
2292 char * returnedSupportNameChar;
2293 char * returnedSupportDescriptionChar;
2294 int size = Supports.size();
2298 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2299 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2301 returnedSupport = new SUPPORT(*obj);
2303 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2304 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2306 returnedSupportNameChar = new char[lenName];
2307 returnedSupportDescriptionChar = new char[lenDescription];
2309 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2310 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2311 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2312 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2313 (Supports[0]->getDescription()).c_str());
2315 returnedSupportName = string(returnedSupportNameChar);
2316 returnedSupportDescription = string(returnedSupportDescriptionChar);
2318 returnedSupport->setName(returnedSupportName);
2319 returnedSupport->setDescription(returnedSupportDescription);
2321 delete [] returnedSupportNameChar;
2322 delete [] returnedSupportDescriptionChar;
2326 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2327 returnedSupport = new SUPPORT(*obj);
2329 int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
2330 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
2332 for(int i = 1;i<size;i++)
2334 obj = const_cast <SUPPORT *> (Supports[i]);
2335 returnedSupport->blending(obj);
2339 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2340 lenDescription = lenDescription + 5 +
2341 strlen((Supports[i]->getDescription()).c_str());
2345 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2346 lenDescription = lenDescription + 2 +
2347 strlen((Supports[i]->getDescription()).c_str());
2351 returnedSupportNameChar = new char[lenName];
2352 returnedSupportDescriptionChar = new char[lenDescription];
2354 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
2355 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
2357 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2358 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2359 (Supports[0]->getDescription()).c_str());
2361 for(int i = 1;i<size;i++)
2365 returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
2366 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
2368 returnedSupportNameChar = strcat(returnedSupportNameChar,
2369 (Supports[i]->getName()).c_str());
2370 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2371 (Supports[i]->getDescription()).c_str());
2375 returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
2376 returnedSupportNameChar = strcat(returnedSupportNameChar,
2377 (Supports[i]->getName()).c_str());
2379 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
2380 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2381 (Supports[i]->getDescription()).c_str());
2385 returnedSupportName = string(returnedSupportNameChar);
2386 returnedSupport->setName(returnedSupportName);
2388 returnedSupportDescription = string(returnedSupportDescriptionChar);
2389 returnedSupport->setDescription(returnedSupportDescription);
2391 delete [] returnedSupportNameChar;
2392 delete [] returnedSupportDescriptionChar;
2396 return returnedSupport;
2400 return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
2401 The (SUPPORT *) NULL pointer is returned if the intersection is empty.
2402 You should delete this pointer after use to avois memory leaks.
2404 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2406 const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
2409 SUPPORT * returnedSupport;
2410 string returnedSupportName;
2411 string returnedSupportDescription;
2412 char * returnedSupportNameChar;
2413 char * returnedSupportDescriptionChar;
2414 int size = Supports.size();
2418 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2419 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2421 returnedSupport = new SUPPORT(*obj);
2423 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2424 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2426 returnedSupportNameChar = new char[lenName];
2427 returnedSupportDescriptionChar = new char[lenDescription];
2429 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2430 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2431 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2432 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2433 (Supports[0]->getDescription()).c_str());
2435 returnedSupportName = string(returnedSupportNameChar);
2436 returnedSupportDescription = string(returnedSupportDescriptionChar);
2438 returnedSupport->setName(returnedSupportName);
2439 returnedSupport->setDescription(returnedSupportDescription);
2441 delete [] returnedSupportNameChar;
2442 delete [] returnedSupportDescriptionChar;
2446 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2447 returnedSupport = new SUPPORT(*obj);
2449 int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
2450 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
2452 for(int i = 1;i<size;i++)
2454 obj = const_cast <SUPPORT *> (Supports[i]);
2455 returnedSupport->intersecting(obj);
2459 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2460 lenDescription = lenDescription + 5 +
2461 strlen((Supports[i]->getDescription()).c_str());
2465 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2466 lenDescription = lenDescription + 2 +
2467 strlen((Supports[i]->getDescription()).c_str());
2471 if(returnedSupport != (SUPPORT *) NULL)
2473 returnedSupportNameChar = new char[lenName];
2474 returnedSupportDescriptionChar = new char[lenDescription];
2476 returnedSupportNameChar = strcpy(returnedSupportNameChar,
2477 "Intersection of ");
2478 returnedSupportDescriptionChar =
2479 strcpy(returnedSupportDescriptionChar,"Intersection of ");
2481 returnedSupportNameChar = strcat(returnedSupportNameChar,
2482 (Supports[0]->getName()).c_str());
2483 returnedSupportDescriptionChar =
2484 strcat(returnedSupportDescriptionChar,
2485 (Supports[0]->getDescription()).c_str());
2487 for(int i = 1;i<size;i++)
2491 returnedSupportNameChar = strcat(returnedSupportNameChar,
2493 returnedSupportDescriptionChar =
2494 strcat(returnedSupportDescriptionChar," and ");
2496 returnedSupportNameChar =
2497 strcat(returnedSupportNameChar,
2498 (Supports[i]->getName()).c_str());
2499 returnedSupportDescriptionChar =
2500 strcat(returnedSupportDescriptionChar,
2501 (Supports[i]->getDescription()).c_str());
2505 returnedSupportNameChar = strcat(returnedSupportNameChar,
2507 returnedSupportNameChar =
2508 strcat(returnedSupportNameChar,
2509 (Supports[i]->getName()).c_str());
2511 returnedSupportDescriptionChar =
2512 strcat(returnedSupportDescriptionChar,", ");
2513 returnedSupportDescriptionChar =
2514 strcat(returnedSupportDescriptionChar,
2515 (Supports[i]->getDescription()).c_str());
2519 returnedSupportName = string(returnedSupportNameChar);
2520 returnedSupport->setName(returnedSupportName);
2522 returnedSupportDescription = string(returnedSupportDescriptionChar);
2523 returnedSupport->setDescription(returnedSupportDescription);
2525 delete [] returnedSupportNameChar;
2526 delete [] returnedSupportDescriptionChar;
2531 return returnedSupport;
2535 // internal helper type
2539 std::vector<int> groups;
2540 MED_EN::medGeometryElement geometricType;
2543 // Create families from groups
2544 void MESH::createFamilies()
2546 int nbFam=0; // count the families we create, used as identifier ???????????
2548 int idFamNode = 0; // identifier for node families
2549 int idFamElement = 0; // identifier for cell, face or edge families
2551 // Main loop on mesh's entities
2552 for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2554 int numberofgroups = getNumberOfGroups(entity);
2556 continue; // no groups for this entity
2558 vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2560 // make myFamilies points to the member corresponding to entity
2561 vector<FAMILY*>* myFamilies;
2565 myFamilies = & _familyCell;
2568 myFamilies = & _familyFace;
2571 myFamilies = & _familyEdge;
2574 myFamilies = & _familyNode;
2578 vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2579 // get a copy of the (old) family ptrs before clearing
2580 vector<FAMILY*> myOldFamilies=getFamilies(entity);
2581 myFamilies->clear();
2584 // 1 - Create a vector containing for each cell (of the entity) an information structure
2585 // giving geometric type and the groups it belong to
2587 med_int numberOfTypes=getNumberOfTypes(entity);
2588 const int * index=getGlobalNumberingIndex(entity);
2589 const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2590 med_int numberOfCells=index[numberOfTypes]-1; // total number of cells for that entity
2591 SCRUTE(numberOfTypes);
2592 SCRUTE(numberOfCells);
2593 vector< _cell > tab_cell(numberOfCells);
2594 for(med_int t=0; t!=numberOfTypes; ++t)
2595 for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2596 tab_cell[n].geometricType=geometricTypes[t];
2599 // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2601 for (unsigned g=0; g!=myGroups.size(); ++g)
2603 // scan cells that belongs to the group
2604 const med_int* groupCells=myGroups[g]->getnumber()->getValue();
2605 int nbCells=myGroups[g]->getnumber()->getLength();
2606 for(int c=0; c!=nbCells; ++c)
2607 tab_cell[groupCells[c]-1].groups.push_back(g);
2611 // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2612 // whith the vector of contained cells
2614 map< string,vector<int> > tab_families;
2615 map< string,vector<int> >::iterator fam;
2616 for(int n=0; n!=numberOfCells; ++n)
2618 ostringstream key; // to generate the name of the family
2620 if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2621 key << "_NONE" << entity;
2623 for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2625 string groupName=myGroups[*it]->getName();
2626 if(groupName.empty())
2629 key << "_" << groupName;
2632 tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2633 /* fam = tab_families.find(key.str());
2634 if( fam != tab_families.end())
2635 fam->second.push_back(n+1); // +1 : convention Fortran de MED
2638 vector<int> newfamily;
2639 newfamily.push_back(n+1); // +1 : convention Fortran de MED
2640 tab_families.insert(make_pair(key.str(),newfamily));
2645 // 4 - Scan the family map, create MED Families, check if it already exist.
2647 for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2649 vector<medGeometryElement> tab_types_geometriques;
2650 medGeometryElement geometrictype=MED_NONE;
2651 vector<int> tab_index_types_geometriques;
2652 vector<int> tab_nombres_elements;
2654 // scan family cells and fill the tab that are needed by the create a MED FAMILY
2655 for( int i=0; i!=fam->second.size(); ++i)
2657 int ncell=fam->second[i]-1;
2658 if(tab_cell[ncell].geometricType != geometrictype)
2660 // new geometric type -> we store it and complete index tabs
2661 if(!tab_index_types_geometriques.empty())
2662 tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2663 tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2664 tab_index_types_geometriques.push_back(i+1);
2667 // store and complete index tabs for the last geometric type
2668 tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2669 tab_index_types_geometriques.push_back(fam->second.size()+1);
2671 // create a empty MED FAMILY and fill it with the tabs we constructed
2672 FAMILY* newFam = new FAMILY();
2673 newFam->setTotalNumberOfElements(fam->second.size());
2674 newFam->setName(fam->first);
2675 newFam->setMesh(this);
2676 newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2677 newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2678 newFam->setNumberOfElements(&tab_nombres_elements[0]);
2679 newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2680 newFam->setEntity(entity);
2681 newFam->setAll(false);
2693 idFam = -idFamElement;
2697 idFam = -idFamElement;
2701 idFam = -idFamElement;
2705 newFam->setIdentifier(idFam);
2707 // Update links between families and groups
2709 int ncell1=fam->second[0]-1; // number of first cell in family
2710 int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2713 newFam->setNumberOfGroups(numberOfGroups);
2714 string * groupNames=new string[numberOfGroups];
2716 // iterate on the groups the family belongs to
2717 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2718 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2720 whichFamilyInGroup[*it].push_back(newFam);
2721 groupNames[ng]=myGroups[*it]->getName();
2723 newFam->setGroupsNames(groupNames);
2726 int sizeOfFamVect = myFamilies->size();
2728 MESSAGE(" MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2730 myFamilies->push_back(newFam);
2733 // delete old families
2734 for (int i=0;i<myOldFamilies.size();i++)
2735 delete myOldFamilies[i] ;
2737 // update references in groups
2738 for (int i=0;i<myGroups.size();i++)
2740 myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2741 myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2744 // re-scan the cells vector, and fill the family vector with cells.
2745 // creation of support, check if it already exist.