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 _coordinate = (COORDINATE *) NULL;
151 _connectivity = (CONNECTIVITY *) NULL;
153 _spaceDimension = MED_INVALID; // 0 ?!?
154 _meshDimension = MED_INVALID;
155 _numberOfNodes = MED_INVALID;
159 _arePresentOptionnalNodesNumbers = 0;
164 /*! Create an empty MESH. */
165 MESH::MESH():_coordinate(NULL),_connectivity(NULL), _isAGrid(false) {
172 _isAGrid = m._isAGrid;
174 if (m._coordinate != NULL)
175 _coordinate = new COORDINATE(* m._coordinate);
177 _coordinate = (COORDINATE *) NULL;
179 if (m._connectivity != NULL)
180 _connectivity = new CONNECTIVITY(* m._connectivity);
182 _connectivity = (CONNECTIVITY *) NULL;
184 _spaceDimension = m._spaceDimension;
185 _meshDimension = m._meshDimension;
186 _numberOfNodes = m._numberOfNodes;
188 _familyNode = m._familyNode;
189 for (int i=0; i<m._familyNode.size(); i++)
191 _familyNode[i] = new FAMILY(* m._familyNode[i]);
192 _familyNode[i]->setMesh(this);
195 _familyCell = m._familyCell;
196 for (int i=0; i<m._familyCell.size(); i++)
198 _familyCell[i] = new FAMILY(* m._familyCell[i]);
199 _familyCell[i]->setMesh(this);
202 _familyFace = m._familyFace;
203 for (int i=0; i<m._familyFace.size(); i++)
205 _familyFace[i] = new FAMILY(* m._familyFace[i]);
206 _familyFace[i]->setMesh(this);
209 _familyEdge = m._familyEdge;
210 for (int i=0; i<m._familyEdge.size(); i++)
212 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
213 _familyEdge[i]->setMesh(this);
216 _groupNode = m._groupNode;
217 for (int i=0; i<m._groupNode.size(); i++)
219 _groupNode[i] = new GROUP(* m._groupNode[i]);
220 _groupNode[i]->setMesh(this);
223 _groupCell = m._groupCell;
224 for (int i=0; i<m._groupCell.size(); i++)
226 _groupCell[i] = new GROUP(* m._groupCell[i]);
227 _groupCell[i]->setMesh(this);
230 _groupFace = m._groupFace;
231 for (int i=0; i<m._groupFace.size(); i++)
233 _groupFace[i] = new GROUP(* m._groupFace[i]);
234 _groupFace[i]->setMesh(this);
237 _groupEdge = m._groupEdge;
238 for (int i=0; i<m._groupEdge.size(); i++)
240 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
241 _groupEdge[i]->setMesh(this);
244 //_drivers = m._drivers; //Recopie des drivers?
249 MESSAGE("MESH::~MESH() : Destroying the Mesh");
250 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
251 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
253 size = _familyNode.size() ;
254 for (int i=0;i<size;i++)
255 delete _familyNode[i] ;
256 size = _familyCell.size() ;
257 for (int i=0;i<size;i++)
258 delete _familyCell[i] ;
259 size = _familyFace.size() ;
260 for (int i=0;i<size;i++)
261 delete _familyFace[i] ;
262 size = _familyEdge.size() ;
263 for (int i=0;i<size;i++)
264 delete _familyEdge[i] ;
265 size = _groupNode.size() ;
266 for (int i=0;i<size;i++)
267 delete _groupNode[i] ;
268 size = _groupCell.size() ;
269 for (int i=0;i<size;i++)
270 delete _groupCell[i] ;
271 size = _groupFace.size() ;
272 for (int i=0;i<size;i++)
273 delete _groupFace[i] ;
274 size = _groupEdge.size() ;
275 for (int i=0;i<size;i++)
276 delete _groupEdge[i] ;
278 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
280 for (unsigned int index=0; index < _drivers.size(); index++ )
282 SCRUTE(_drivers[index]);
283 if ( _drivers[index] != NULL) delete _drivers[index];
288 MESH & MESH::operator=(const MESH &m)
290 const char * LOC = "MESH & MESH::operator=(const MESH &m) : ";
293 MESSAGE(LOC <<"Not yet implemented, operating on the object " << m);
296 // ATTENTION CET OPERATEUR DE RECOPIE EST DANGEREUX POUR LES
297 // POINTEURS : ex : nodal_connectivity ???? EXPRES ????
298 // _drivers = m._drivers;
300 // space_dimension=m.space_dimension;
301 // mesh_dimension=m.mesh_dimension;
303 // nodes_count=m.nodes_count;
304 // coordinates=m.coordinates;
305 // coordinates_name=m.coordinates_name ;
306 // coordinates_unit=m.coordinates_unit ;
307 // nodes_numbers=m.nodes_numbers ;
309 // cells_types_count=m.cells_types_count;
310 // cells_type=m.cells_type;
311 // cells_count=m.cells_count;
312 // nodal_connectivity=m.nodal_connectivity;
314 // nodes_families_count=m.nodes_families_count;
315 // nodes_Families=m.nodes_Families;
317 // cells_families_count=m.cells_families_count;
318 // cells_Families=m.cells_Families;
320 // maximum_cell_number_by_node = m.maximum_cell_number_by_node;
321 // if (maximum_cell_number_by_node > 0)
323 // reverse_nodal_connectivity = m.reverse_nodal_connectivity;
324 // reverse_nodal_connectivity_index = m.reverse_nodal_connectivity_index ;
331 /*! Create a %MESH object using a %MESH driver of type %driverTypes (MED_DRIVER, ....) associated with file fileName.
332 The meshname driverName must already exists in the file.*/
333 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) throw (MEDEXCEPTION)
335 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
346 MED_MESH_RDONLY_DRIVER myDriver(fileName,this);
347 myDriver.setMeshName(driverName);
348 current = addDriver(myDriver);
353 GIBI_MESH_RDONLY_DRIVER myDriver(fileName,this);
354 current = addDriver(myDriver);
357 case PORFLOW_DRIVER :
359 PORFLOW_MESH_RDONLY_DRIVER myDriver(fileName,this);
360 current = addDriver(myDriver);
365 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
370 // current = addDriver(driverType,fileName,driverName);
371 // switch(_drivers[current]->getAccessMode() ) {
372 // case MED_WRONLY : {
373 // MESSAGE("MESH::MESH(driverTypes driverType, .....) : driverType must not be a MED_WRONLY accessMode");
374 // rmDriver(current);
379 _drivers[current]->open();
380 _drivers[current]->read();
381 _drivers[current]->close();
384 // ((GRID *) this)->fillMeshAfterRead();
390 // Node MESH::Donne_Barycentre(const Maille &m) const
392 // Node N=node[m[0]];
393 // for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
394 // N=N*(1.0/m.donne_nbr_sommets());
398 ostream & MEDMEM::operator<<(ostream &os, const MESH &myMesh)
400 int spacedimension = myMesh.getSpaceDimension();
401 int meshdimension = myMesh.getMeshDimension();
402 int numberofnodes = myMesh.getNumberOfNodes();
404 os << "Space Dimension : " << spacedimension << endl << endl;
406 os << "Mesh Dimension : " << meshdimension << endl << endl;
408 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
409 os << "SHOW NODES COORDINATES : " << endl;
411 os << "Name :" << endl;
412 const string * coordinatesnames = myMesh.getCoordinatesNames();
413 for(int i=0; i<spacedimension ; i++)
415 os << " - " << coordinatesnames[i] << endl;
417 os << "Unit :" << endl;
418 const string * coordinatesunits = myMesh.getCoordinatesUnits();
419 for(int i=0; i<spacedimension ; i++)
421 os << " - " << coordinatesunits[i] << endl;
423 for(int i=0; i<numberofnodes ; i++)
425 os << "Nodes " << i+1 << " : ";
426 for (int j=0; j<spacedimension ; j++)
427 os << coordinates[i*spacedimension+j] << " ";
431 os << endl << "SHOW CONNECTIVITY :" << endl;
432 os << *myMesh._connectivity << endl;
434 medEntityMesh entity;
435 os << endl << "SHOW FAMILIES :" << endl << endl;
436 for (int k=1; k<=4; k++)
438 if (k==1) entity = MED_NODE;
439 if (k==2) entity = MED_CELL;
440 if (k==3) entity = MED_FACE;
441 if (k==4) entity = MED_EDGE;
442 int numberoffamilies = myMesh.getNumberOfFamilies(entity);
443 using namespace MED_FR;
444 os << "NumberOfFamilies on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberoffamilies<<endl;
445 using namespace MED_EN;
446 for (int i=1; i<numberoffamilies+1;i++)
448 os << * myMesh.getFamily(entity,i) << endl;
452 os << endl << "SHOW GROUPS :" << endl << endl;
453 for (int k=1; k<=4; k++)
455 if (k==1) entity = MED_NODE;
456 if (k==2) entity = MED_CELL;
457 if (k==3) entity = MED_FACE;
458 if (k==4) entity = MED_EDGE;
459 int numberofgroups = myMesh.getNumberOfGroups(entity);
460 using namespace MED_FR;
461 os << "NumberOfGroups on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberofgroups<<endl;
462 using namespace MED_EN;
463 for (int i=1; i<numberofgroups+1;i++)
465 os << * myMesh.getGroup(entity,i) << endl;
473 Get global number of element which have same connectivity than connectivity argument.
475 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
477 Return -1 if not found.
479 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
481 const char* LOC="MESH::getElementNumber " ;
484 int numberOfValue ; // size of connectivity array
485 CELLMODEL myType(Type) ;
486 if (ConnectivityType==MED_DESCENDING)
487 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
489 numberOfValue = myType.getNumberOfNodes() ; // nodes
491 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
492 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
494 // First node or face/edge
495 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
496 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
498 // list to put cells numbers
499 list<int> cellsList ;
500 list<int>::iterator itList ;
501 for (int i=indexBegin; i<indexEnd; i++)
502 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
504 // Foreach node or face/edge in cell, we search cells which are in common.
505 // At the end, we must have only one !
507 for (int i=1; i<numberOfValue; i++) {
508 int connectivity_i = connectivity[i] ;
509 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
511 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
512 if ((*itList)==myReverseConnectivityValue[j-1]) {
518 itList=cellsList.erase(itList);
519 itList--; // well : rigth if itList = cellsList.begin() ??
524 if (cellsList.size()>1) // we have more than one cell
525 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
527 if (cellsList.size()==0)
532 return cellsList.front() ;
537 Return a support which reference all elements on the boundary of mesh.
539 For instance, we could get only face in 3D and edge in 2D.
541 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
544 const char * LOC = "MESH::getBoundaryElements : " ;
547 // actually we could only get face (in 3D) and edge (in 2D)
548 if (_spaceDimension == 3)
549 if (Entity != MED_FACE)
550 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
551 if (_spaceDimension == 2)
552 if (Entity != MED_EDGE)
553 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
556 SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
557 //mySupport.setDescription("boundary of type ...");
558 mySupport->setAll(false);
561 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
562 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
563 int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
564 list<int> myElementsList ;
566 for (int i=0 ; i<numberOf; i++)
567 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
568 myElementsList.push_back(i+1) ;
571 // Well, we must know how many geometric type we have found
572 int * myListArray = new int[size] ;
574 list<int>::iterator myElementsListIt ;
575 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
576 myListArray[id]=(*myElementsListIt) ;
580 int numberOfGeometricType ;
581 medGeometryElement* geometricType ;
582 int * numberOfGaussPoint ;
583 int * geometricTypeNumber ;
584 int * numberOfElements ;
585 //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
586 int * mySkyLineArrayIndex ;
588 int numberOfType = getNumberOfTypes(Entity) ;
589 if (numberOfType == 1) { // wonderfull : it's easy !
590 numberOfGeometricType = 1 ;
591 geometricType = new medGeometryElement[1] ;
592 const medGeometryElement * allType = getTypes(Entity);
593 geometricType[0] = allType[0] ;
594 numberOfGaussPoint = new int[1] ;
595 numberOfGaussPoint[0] = 1 ;
596 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
597 geometricTypeNumber[0] = 0 ;
598 numberOfElements = new int[1] ;
599 numberOfElements[0] = size ;
600 mySkyLineArrayIndex = new int[2] ;
601 mySkyLineArrayIndex[0]=1 ;
602 mySkyLineArrayIndex[1]=1+size ;
605 map<medGeometryElement,int> theType ;
606 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
607 medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
608 if (theType.find(myType) != theType.end() )
613 numberOfGeometricType = theType.size() ;
614 geometricType = new medGeometryElement[numberOfGeometricType] ;
615 //const medGeometryElement * allType = getTypes(Entity); !! UNUSZED VARIABLE !!
616 numberOfGaussPoint = new int[numberOfGeometricType] ;
617 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
618 numberOfElements = new int[numberOfGeometricType] ;
619 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
621 mySkyLineArrayIndex[0]=1 ;
622 map<medGeometryElement,int>::iterator theTypeIt ;
623 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
624 geometricType[index] = (*theTypeIt).first ;
625 numberOfGaussPoint[index] = 1 ;
626 geometricTypeNumber[index] = 0 ;
627 numberOfElements[index] = (*theTypeIt).second ;
628 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
632 //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
633 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
635 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
636 mySupport->setGeometricType(geometricType) ;
637 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
638 mySupport->setNumberOfElements(numberOfElements) ;
639 mySupport->setTotalNumberOfElements(size) ;
640 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
641 mySupport->setNumber(mySkyLineArray) ;
643 delete[] numberOfElements;
644 delete[] geometricTypeNumber;
645 delete[] numberOfGaussPoint;
646 delete[] geometricType;
647 delete[] mySkyLineArrayIndex;
648 delete[] myListArray;
649 // delete mySkyLineArray;
655 FIELD<double>* MESH::getVolume(const SUPPORT * Support) const throw (MEDEXCEPTION)
657 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
660 // Support must be on 3D elements
662 // Make sure that the MESH class is the same as the MESH class attribut
663 // in the class Support
664 MESH* myMesh = Support->getMesh();
666 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
668 int dim_space = getSpaceDimension();
669 medEntityMesh support_entity = Support->getEntity();
670 bool onAll = Support->isOnAllElements();
672 int nb_type, length_values;
673 const medGeometryElement* types;
675 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
676 const int* global_connectivity;
680 // nb_type = myMesh->getNumberOfTypes(support_entity);
681 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
682 // types = getTypes(support_entity);
686 nb_type = Support->getNumberOfTypes();
687 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
688 types = Support->getTypes();
692 FIELD<double>* Volume = new FIELD<double>::FIELD(Support,1);
693 // double *volume = new double [length_values];
694 Volume->setName("VOLUME");
695 Volume->setDescription("cells volume");
696 Volume->setComponentName(1,"volume");
697 Volume->setComponentDescription(1,"desc-comp");
699 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
701 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
703 Volume->setMEDComponentUnit(1,MEDComponentUnit);
705 Volume->setValueType(MED_REEL64);
707 Volume->setIterationNumber(0);
708 Volume->setOrderNumber(0);
709 Volume->setTime(0.0);
711 //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
712 MEDARRAY<double> *volume = Volume->getvalue();
715 const double * coord = getCoordinates(MED_FULL_INTERLACE);
717 for (int i=0;i<nb_type;i++)
719 medGeometryElement type = types[i] ;
724 nb_entity_type = getNumberOfElements(support_entity,type);
725 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
729 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all"));
731 nb_entity_type = Support->getNumberOfElements(type);
733 const int * supp_number = Support->getNumber(type);
734 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
735 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
736 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
738 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
739 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
740 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
743 global_connectivity = global_connectivity_tmp ;
748 case MED_TETRA4 : case MED_TETRA10 :
750 for (int tetra=0;tetra<nb_entity_type;tetra++)
752 int tetra_index = (type%100)*tetra;
754 int N1 = global_connectivity[tetra_index]-1;
755 int N2 = global_connectivity[tetra_index+1]-1;
756 int N3 = global_connectivity[tetra_index+2]-1;
757 int N4 = global_connectivity[tetra_index+3]-1;
759 double x1 = coord[dim_space*N1];
760 double x2 = coord[dim_space*N2];
761 double x3 = coord[dim_space*N3];
762 double x4 = coord[dim_space*N4];
764 double y1 = coord[dim_space*N1+1];
765 double y2 = coord[dim_space*N2+1];
766 double y3 = coord[dim_space*N3+1];
767 double y4 = coord[dim_space*N4+1];
769 double z1 = coord[dim_space*N1+2];
770 double z2 = coord[dim_space*N2+2];
771 double z3 = coord[dim_space*N3+2];
772 double z4 = coord[dim_space*N4+2];
774 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
775 (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
776 (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
778 //volume[index] = xvolume ;
779 volume->setIJ(index,1,xvolume) ;
784 case MED_PYRA5 : case MED_PYRA13 :
786 for (int pyra=0;pyra<nb_entity_type;pyra++)
788 int pyra_index = (type%100)*pyra;
790 int N1 = global_connectivity[pyra_index]-1;
791 int N2 = global_connectivity[pyra_index+1]-1;
792 int N3 = global_connectivity[pyra_index+2]-1;
793 int N4 = global_connectivity[pyra_index+3]-1;
794 int N5 = global_connectivity[pyra_index+4]-1;
796 double x1 = coord[dim_space*N1];
797 double x2 = coord[dim_space*N2];
798 double x3 = coord[dim_space*N3];
799 double x4 = coord[dim_space*N4];
800 double x5 = coord[dim_space*N5];
802 double y1 = coord[dim_space*N1+1];
803 double y2 = coord[dim_space*N2+1];
804 double y3 = coord[dim_space*N3+1];
805 double y4 = coord[dim_space*N4+1];
806 double y5 = coord[dim_space*N5+1];
808 double z1 = coord[dim_space*N1+2];
809 double z2 = coord[dim_space*N2+2];
810 double z3 = coord[dim_space*N3+2];
811 double z4 = coord[dim_space*N4+2];
812 double z5 = coord[dim_space*N5+2];
814 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
815 (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
816 (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
817 ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
818 (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
819 (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
822 //volume[index] = xvolume ;
823 volume->setIJ(index,1,xvolume) ;
828 case MED_PENTA6 : case MED_PENTA15 :
830 for (int penta=0;penta<nb_entity_type;penta++)
832 int penta_index = (type%100)*penta;
834 int N1 = global_connectivity[penta_index]-1;
835 int N2 = global_connectivity[penta_index+1]-1;
836 int N3 = global_connectivity[penta_index+2]-1;
837 int N4 = global_connectivity[penta_index+3]-1;
838 int N5 = global_connectivity[penta_index+4]-1;
839 int N6 = global_connectivity[penta_index+5]-1;
841 double x1 = coord[dim_space*N1];
842 double x2 = coord[dim_space*N2];
843 double x3 = coord[dim_space*N3];
844 double x4 = coord[dim_space*N4];
845 double x5 = coord[dim_space*N5];
846 double x6 = coord[dim_space*N6];
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];
855 double z1 = coord[dim_space*N1+2];
856 double z2 = coord[dim_space*N2+2];
857 double z3 = coord[dim_space*N3+2];
858 double z4 = coord[dim_space*N4+2];
859 double z5 = coord[dim_space*N5+2];
860 double z6 = coord[dim_space*N6+2];
862 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
863 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
864 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
865 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
866 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
867 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
868 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
870 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
872 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
874 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
875 (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
876 (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
877 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
879 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
881 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
882 (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
883 (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
884 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
886 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
888 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
889 (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
890 (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
892 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
894 //volume[index] = xvolume ;
895 volume->setIJ(index,1,xvolume) ;
900 case MED_HEXA8 : case MED_HEXA20 :
902 for (int hexa=0;hexa<nb_entity_type;hexa++)
904 int hexa_index = (type%100)*hexa;
906 int N1 = global_connectivity[hexa_index]-1;
907 int N2 = global_connectivity[hexa_index+1]-1;
908 int N3 = global_connectivity[hexa_index+2]-1;
909 int N4 = global_connectivity[hexa_index+3]-1;
910 int N5 = global_connectivity[hexa_index+4]-1;
911 int N6 = global_connectivity[hexa_index+5]-1;
912 int N7 = global_connectivity[hexa_index+6]-1;
913 int N8 = global_connectivity[hexa_index+7]-1;
915 double x1 = coord[dim_space*N1];
916 double x2 = coord[dim_space*N2];
917 double x3 = coord[dim_space*N3];
918 double x4 = coord[dim_space*N4];
919 double x5 = coord[dim_space*N5];
920 double x6 = coord[dim_space*N6];
921 double x7 = coord[dim_space*N7];
922 double x8 = coord[dim_space*N8];
924 double y1 = coord[dim_space*N1+1];
925 double y2 = coord[dim_space*N2+1];
926 double y3 = coord[dim_space*N3+1];
927 double y4 = coord[dim_space*N4+1];
928 double y5 = coord[dim_space*N5+1];
929 double y6 = coord[dim_space*N6+1];
930 double y7 = coord[dim_space*N7+1];
931 double y8 = coord[dim_space*N8+1];
933 double z1 = coord[dim_space*N1+2];
934 double z2 = coord[dim_space*N2+2];
935 double z3 = coord[dim_space*N3+2];
936 double z4 = coord[dim_space*N4+2];
937 double z5 = coord[dim_space*N5+2];
938 double z6 = coord[dim_space*N6+2];
939 double z7 = coord[dim_space*N7+2];
940 double z8 = coord[dim_space*N8+2];
942 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
943 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
944 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
945 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
946 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
947 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
948 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
949 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
950 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
951 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
952 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
953 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
955 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
957 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
959 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
960 (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
961 (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
962 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
964 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
966 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
967 (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
968 (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
969 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
970 (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
971 (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
972 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
973 (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
974 (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
975 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
976 ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
977 ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
978 ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
979 ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
980 ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
981 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
983 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
985 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
986 (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
987 (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
988 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
990 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
992 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
993 (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
994 (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
995 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
996 (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
997 (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
998 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
999 (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
1000 (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
1001 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
1002 ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
1003 ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
1004 ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
1005 ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
1006 ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
1007 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
1008 (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
1009 (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
1010 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
1011 (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
1012 (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
1013 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
1014 ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
1015 ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
1016 ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
1017 ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
1018 ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
1019 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
1020 (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
1021 (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
1022 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
1023 (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
1024 (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
1025 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
1026 ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
1027 ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
1028 ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
1029 ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
1030 ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
1031 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
1032 (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
1033 (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
1034 (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
1035 (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
1036 (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
1037 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
1038 (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
1039 (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
1040 (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
1041 (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
1042 (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
1043 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
1044 (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
1045 ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
1046 (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
1047 ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
1048 (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
1049 ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
1050 (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
1051 ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
1052 (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
1053 ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
1054 (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
1056 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
1057 4.0*(C + F + G + H + L + O + P + Q + S + T +
1058 V + W) + 2.0*(I + R + U + X + Y + Z) +
1061 //volume[index] = xvolume ;
1062 volume->setIJ(index,1,xvolume) ;
1068 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
1072 if (!onAll) delete [] global_connectivity ;
1078 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1080 const char * LOC = "MESH::getArea(SUPPORT*) : ";
1083 // Support must be on 2D elements
1085 // Make sure that the MESH class is the same as the MESH class attribut
1086 // in the class Support
1087 MESH* myMesh = Support->getMesh();
1089 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1091 int dim_space = getSpaceDimension();
1092 medEntityMesh support_entity = Support->getEntity();
1093 bool onAll = Support->isOnAllElements();
1095 int nb_type, length_values;
1096 const medGeometryElement* types;
1098 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1099 const int* global_connectivity;
1103 // nb_type = myMesh->getNumberOfTypes(support_entity);
1104 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1105 // types = getTypes(support_entity);
1109 nb_type = Support->getNumberOfTypes();
1110 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1111 types = Support->getTypes();
1115 FIELD<double>* Area;
1117 Area = new FIELD<double>::FIELD(Support,1);
1118 Area->setName("AREA");
1119 Area->setDescription("cells or faces area");
1121 Area->setComponentName(1,"area");
1122 Area->setComponentDescription(1,"desc-comp");
1124 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1126 string MEDComponentUnit="trucmuch";
1128 Area->setMEDComponentUnit(1,MEDComponentUnit);
1130 Area->setValueType(MED_REEL64);
1132 Area->setIterationNumber(0);
1133 Area->setOrderNumber(0);
1136 double *area = new double[length_values];
1137 //double *area = Area->getValue(MED_FULL_INTERLACE);
1139 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1142 for (int i=0;i<nb_type;i++)
1144 medGeometryElement type = types[i] ;
1149 nb_entity_type = getNumberOfElements(support_entity,type);
1150 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1154 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1156 nb_entity_type = Support->getNumberOfElements(type);
1158 const int * supp_number = Support->getNumber(type);
1159 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1160 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1161 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1163 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1164 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1165 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1169 global_connectivity = global_connectivity_tmp ;
1175 case MED_TRIA3 : case MED_TRIA6 :
1177 for (int tria=0;tria<nb_entity_type;tria++)
1179 int tria_index = (type%100)*tria;
1181 int N1 = global_connectivity[tria_index]-1;
1182 int N2 = global_connectivity[tria_index+1]-1;
1183 int N3 = global_connectivity[tria_index+2]-1;
1185 double x1 = coord[dim_space*N1];
1186 double x2 = coord[dim_space*N2];
1187 double x3 = coord[dim_space*N3];
1189 double y1 = coord[dim_space*N1+1];
1190 double y2 = coord[dim_space*N2+1];
1191 double y3 = coord[dim_space*N3+1];
1195 xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1199 double z1 = coord[dim_space*N1+2];
1200 double z2 = coord[dim_space*N2+2];
1201 double z3 = coord[dim_space*N3+2];
1203 xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1204 ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1205 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1206 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1207 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1208 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1211 area[index] = xarea ;
1216 case MED_QUAD4 : case MED_QUAD8 :
1218 for (int quad=0;quad<nb_entity_type;quad++)
1220 int quad_index = (type%100)*quad;
1222 int N1 = global_connectivity[quad_index]-1;
1223 int N2 = global_connectivity[quad_index+1]-1;
1224 int N3 = global_connectivity[quad_index+2]-1;
1225 int N4 = global_connectivity[quad_index+3]-1;
1227 double x1 = coord[dim_space*N1];
1228 double x2 = coord[dim_space*N2];
1229 double x3 = coord[dim_space*N3];
1230 double x4 = coord[dim_space*N4];
1232 double y1 = coord[dim_space*N1+1];
1233 double y2 = coord[dim_space*N2+1];
1234 double y3 = coord[dim_space*N3+1];
1235 double y4 = coord[dim_space*N4+1];
1239 double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1240 double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1241 double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1242 double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1244 xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1245 d1*b2 + a1*d2 - d1*a2);
1249 double z1 = coord[dim_space*N1+2];
1250 double z2 = coord[dim_space*N2+2];
1251 double z3 = coord[dim_space*N3+2];
1252 double z4 = coord[dim_space*N4+2];
1254 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1255 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1256 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1257 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1258 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1259 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1260 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1261 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1262 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1263 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1264 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1265 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1268 area[index] = xarea ;
1274 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1278 if (!onAll) delete [] global_connectivity ;
1281 Area->setValue(MED_FULL_INTERLACE,area);
1286 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1288 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1291 // Support must be on 1D elements
1293 // Make sure that the MESH class is the same as the MESH class attribut
1294 // in the class Support
1295 MESH* myMesh = Support->getMesh();
1297 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1299 int dim_space = getSpaceDimension();
1300 medEntityMesh support_entity = Support->getEntity();
1301 bool onAll = Support->isOnAllElements();
1303 int nb_type, length_values;
1304 const medGeometryElement* types;
1306 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1307 const int* global_connectivity;
1311 // nb_type = myMesh->getNumberOfTypes(support_entity);
1312 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1313 // types = getTypes(support_entity);
1317 nb_type = Support->getNumberOfTypes();
1318 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1319 types = Support->getTypes();
1323 FIELD<double>* Length;
1325 Length = new FIELD<double>::FIELD(Support,1);
1326 // double *length = new double [length_values];
1327 Length->setValueType(MED_REEL64);
1329 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1330 MEDARRAY<double> * length = Length->getvalue();
1332 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1335 for (int i=0;i<nb_type;i++)
1337 medGeometryElement type = types[i] ;
1342 nb_entity_type = getNumberOfElements(support_entity,type);
1343 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1347 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1349 nb_entity_type = Support->getNumberOfElements(type);
1351 const int * supp_number = Support->getNumber(type);
1352 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1353 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1354 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1356 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1357 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1358 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1361 global_connectivity = global_connectivity_tmp ;
1366 case MED_SEG2 : case MED_SEG3 :
1368 for (int edge=0;edge<nb_entity_type;edge++)
1370 int edge_index = (type%100)*edge;
1372 int N1 = global_connectivity[edge_index]-1;
1373 int N2 = global_connectivity[edge_index+1]-1;
1375 double x1 = coord[dim_space*N1];
1376 double x2 = coord[dim_space*N2];
1378 double y1 = coord[dim_space*N1+1];
1379 double y2 = coord[dim_space*N2+1];
1381 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1382 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1384 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1385 (z1 - z2)*(z1 - z2));
1387 length->setIJ(index,1,xlength) ;
1393 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1397 if (!onAll) delete [] global_connectivity ;
1403 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1405 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1408 // Support must be on 2D or 1D elements
1410 // Make sure that the MESH class is the same as the MESH class attribut
1411 // in the class Support
1412 MESH* myMesh = Support->getMesh();
1414 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1416 int dim_space = getSpaceDimension();
1417 medEntityMesh support_entity = Support->getEntity();
1418 bool onAll = Support->isOnAllElements();
1420 int nb_type, length_values;
1421 const medGeometryElement* types;
1423 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1424 const int* global_connectivity;
1428 // nb_type = myMesh->getNumberOfTypes(support_entity);
1429 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1430 // types = getTypes(support_entity);
1434 nb_type = Support->getNumberOfTypes();
1435 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1436 types = Support->getTypes();
1441 FIELD<double>* Normal = new FIELD<double>::FIELD(Support,dim_space);
1442 Normal->setName("NORMAL");
1443 Normal->setDescription("cells or faces normal");
1444 for (int k=1;k<=dim_space;k++) {
1445 Normal->setComponentName(k,"normal");
1446 Normal->setComponentDescription(k,"desc-comp");
1447 Normal->setMEDComponentUnit(k,"unit");
1450 Normal->setValueType(MED_REEL64);
1452 Normal->setIterationNumber(MED_NOPDT);
1453 Normal->setOrderNumber(MED_NONOR);
1454 Normal->setTime(0.0);
1456 double * normal = new double [dim_space*length_values];
1457 //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1459 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1462 for (int i=0;i<nb_type;i++)
1464 medGeometryElement type = types[i] ;
1466 // Make sure we trying to get Normals on
1467 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1468 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1470 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1471 (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1472 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1474 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1476 double xnormal1, xnormal2, xnormal3 ;
1480 nb_entity_type = getNumberOfElements(support_entity,type);
1481 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1485 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all for instance !"));
1486 nb_entity_type = Support->getNumberOfElements(type);
1488 const int * supp_number = Support->getNumber(type);
1489 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1490 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1491 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1493 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1494 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1495 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1499 global_connectivity = global_connectivity_tmp ;
1505 case MED_TRIA3 : case MED_TRIA6 :
1507 for (int tria=0;tria<nb_entity_type;tria++)
1509 int tria_index = (type%100)*tria;
1511 int N1 = global_connectivity[tria_index]-1;
1512 int N2 = global_connectivity[tria_index+1]-1;
1513 int N3 = global_connectivity[tria_index+2]-1;
1515 //double xarea; !! UNUSED VARIABLE !!
1516 double x1 = coord[dim_space*N1];
1517 double x2 = coord[dim_space*N2];
1518 double x3 = coord[dim_space*N3];
1520 double y1 = coord[dim_space*N1+1];
1521 double y2 = coord[dim_space*N2+1];
1522 double y3 = coord[dim_space*N3+1];
1524 double z1 = coord[dim_space*N1+2];
1525 double z2 = coord[dim_space*N2+2];
1526 double z3 = coord[dim_space*N3+2];
1528 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1529 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1530 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1532 normal[3*index] = xnormal1 ;
1533 normal[3*index+1] = xnormal2 ;
1534 normal[3*index+2] = xnormal3 ;
1540 case MED_QUAD4 : case MED_QUAD8 :
1542 for (int quad=0;quad<nb_entity_type;quad++)
1544 int quad_index = (type%100)*quad;
1546 int N1 = global_connectivity[quad_index]-1;
1547 int N2 = global_connectivity[quad_index+1]-1;
1548 int N3 = global_connectivity[quad_index+2]-1;
1549 int N4 = global_connectivity[quad_index+3]-1;
1552 double x1 = coord[dim_space*N1];
1553 double x2 = coord[dim_space*N2];
1554 double x3 = coord[dim_space*N3];
1555 double x4 = coord[dim_space*N4];
1557 double y1 = coord[dim_space*N1+1];
1558 double y2 = coord[dim_space*N2+1];
1559 double y3 = coord[dim_space*N3+1];
1560 double y4 = coord[dim_space*N4+1];
1562 double z1 = coord[dim_space*N1+2];
1563 double z2 = coord[dim_space*N2+2];
1564 double z3 = coord[dim_space*N3+2];
1565 double z4 = coord[dim_space*N4+2];
1567 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1568 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1569 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1570 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1573 xnormal1 = xnormal1/xarea;
1574 xnormal2 = xnormal2/xarea;
1575 xnormal3 = xnormal3/xarea;
1577 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1578 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1579 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1580 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1581 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1582 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1583 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1584 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1585 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1586 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1587 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1588 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1590 // sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1591 // ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1592 // ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1593 // ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1594 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1595 // ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1597 xnormal1 = xnormal1*xarea;
1598 xnormal2 = xnormal2*xarea;
1599 xnormal3 = xnormal3*xarea;
1601 normal[3*index] = xnormal1 ;
1602 normal[3*index+1] = xnormal2 ;
1603 normal[3*index+2] = xnormal3 ;
1609 case MED_SEG2 : case MED_SEG3 :
1611 for (int edge=0;edge<nb_entity_type;edge++)
1613 int edge_index = (type%100)*edge;
1615 int N1 = global_connectivity[edge_index]-1;
1616 int N2 = global_connectivity[edge_index+1]-1;
1618 double x1 = coord[dim_space*N1];
1619 double x2 = coord[dim_space*N2];
1621 double y1 = coord[dim_space*N1+1];
1622 double y2 = coord[dim_space*N2+1];
1624 xnormal1 = -(y2-y1);
1627 normal[2*index] = xnormal1 ;
1628 normal[2*index+1] = xnormal2 ;
1635 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1639 if (!onAll) delete [] global_connectivity ;
1642 Normal->setValue(MED_FULL_INTERLACE,normal);
1650 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1652 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1655 // Make sure that the MESH class is the same as the MESH class attribut
1656 // in the class Support
1657 MESH* myMesh = Support->getMesh();
1659 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1661 int dim_space = getSpaceDimension();
1662 medEntityMesh support_entity = Support->getEntity();
1663 bool onAll = Support->isOnAllElements();
1665 int nb_type, length_values;
1666 const medGeometryElement* types;
1668 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1669 const int* global_connectivity;
1673 // nb_type = myMesh->getNumberOfTypes(support_entity);
1674 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1675 // types = getTypes(support_entity);
1679 nb_type = Support->getNumberOfTypes();
1680 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1681 types = Support->getTypes();
1685 FIELD<double>* Barycenter;
1687 Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1688 Barycenter->setName("BARYCENTER");
1689 Barycenter->setDescription("cells or faces barycenter");
1691 for (int k=0;k<dim_space;k++) {
1693 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1694 Barycenter->setComponentDescription(kp1,"desc-comp");
1695 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1698 Barycenter->setValueType(MED_REEL64);
1700 Barycenter->setIterationNumber(0);
1701 Barycenter->setOrderNumber(0);
1702 Barycenter->setTime(0.0);
1704 double *barycenter = new double [dim_space*length_values];
1705 // double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1707 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1710 for (int i=0;i<nb_type;i++)
1712 medGeometryElement type = types[i] ;
1713 double xbarycenter1, xbarycenter2, xbarycenter3;
1717 nb_entity_type = getNumberOfElements(support_entity,type);
1718 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1722 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1723 nb_entity_type = Support->getNumberOfElements(type);
1725 const int * supp_number = Support->getNumber(type);
1726 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1727 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1728 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1730 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1731 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1732 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1735 global_connectivity = global_connectivity_tmp;
1740 case MED_TETRA4 : case MED_TETRA10 :
1742 for (int tetra =0;tetra<nb_entity_type;tetra++)
1744 int tetra_index = (type%100)*tetra;
1746 int N1 = global_connectivity[tetra_index]-1;
1747 int N2 = global_connectivity[tetra_index+1]-1;
1748 int N3 = global_connectivity[tetra_index+2]-1;
1749 int N4 = global_connectivity[tetra_index+3]-1;
1751 double x1 = coord[dim_space*N1];
1752 double x2 = coord[dim_space*N2];
1753 double x3 = coord[dim_space*N3];
1754 double x4 = coord[dim_space*N4];
1756 double y1 = coord[dim_space*N1+1];
1757 double y2 = coord[dim_space*N2+1];
1758 double y3 = coord[dim_space*N3+1];
1759 double y4 = coord[dim_space*N4+1];
1761 double z1 = coord[dim_space*N1+2];
1762 double z2 = coord[dim_space*N2+2];
1763 double z3 = coord[dim_space*N3+2];
1764 double z4 = coord[dim_space*N4+2];
1766 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1767 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1768 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1769 barycenter[3*index] = xbarycenter1 ;
1770 barycenter[3*index+1] = xbarycenter2 ;
1771 barycenter[3*index+2] = xbarycenter3 ;
1776 case MED_PYRA5 : case MED_PYRA13 :
1778 for (int pyra=0;pyra<nb_entity_type;pyra++)
1780 int pyra_index = (type%100)*pyra;
1782 int N1 = global_connectivity[pyra_index]-1;
1783 int N2 = global_connectivity[pyra_index+1]-1;
1784 int N3 = global_connectivity[pyra_index+2]-1;
1785 int N4 = global_connectivity[pyra_index+3]-1;
1786 int N5 = global_connectivity[pyra_index+4]-1;
1788 double x1 = coord[dim_space*N1];
1789 double x2 = coord[dim_space*N2];
1790 double x3 = coord[dim_space*N3];
1791 double x4 = coord[dim_space*N4];
1792 double x5 = coord[dim_space*N5];
1794 double y1 = coord[dim_space*N1+1];
1795 double y2 = coord[dim_space*N2+1];
1796 double y3 = coord[dim_space*N3+1];
1797 double y4 = coord[dim_space*N4+1];
1798 double y5 = coord[dim_space*N5+1];
1800 double z1 = coord[dim_space*N1+2];
1801 double z2 = coord[dim_space*N2+2];
1802 double z3 = coord[dim_space*N3+2];
1803 double z4 = coord[dim_space*N4+2];
1804 double z5 = coord[dim_space*N5+2];
1806 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1807 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1808 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1809 barycenter[3*index] = xbarycenter1 ;
1810 barycenter[3*index+1] = xbarycenter2 ;
1811 barycenter[3*index+2] = xbarycenter3 ;
1816 case MED_PENTA6 : case MED_PENTA15 :
1818 for (int penta=0;penta<nb_entity_type;penta++)
1820 int penta_index = (type%100)*penta;
1822 int N1 = global_connectivity[penta_index]-1;
1823 int N2 = global_connectivity[penta_index+1]-1;
1824 int N3 = global_connectivity[penta_index+2]-1;
1825 int N4 = global_connectivity[penta_index+3]-1;
1826 int N5 = global_connectivity[penta_index+4]-1;
1827 int N6 = global_connectivity[penta_index+5]-1;
1829 double x1 = coord[dim_space*N1];
1830 double x2 = coord[dim_space*N2];
1831 double x3 = coord[dim_space*N3];
1832 double x4 = coord[dim_space*N4];
1833 double x5 = coord[dim_space*N5];
1834 double x6 = coord[dim_space*N6];
1836 double y1 = coord[dim_space*N1+1];
1837 double y2 = coord[dim_space*N2+1];
1838 double y3 = coord[dim_space*N3+1];
1839 double y4 = coord[dim_space*N4+1];
1840 double y5 = coord[dim_space*N5+1];
1841 double y6 = coord[dim_space*N6+1];
1843 double z1 = coord[dim_space*N1+2];
1844 double z2 = coord[dim_space*N2+2];
1845 double z3 = coord[dim_space*N3+2];
1846 double z4 = coord[dim_space*N4+2];
1847 double z5 = coord[dim_space*N5+2];
1848 double z6 = coord[dim_space*N6+2];
1850 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1851 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1852 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1853 barycenter[3*index] = xbarycenter1 ;
1854 barycenter[3*index+1] = xbarycenter2 ;
1855 barycenter[3*index+2] = xbarycenter3 ;
1860 case MED_HEXA8 : case MED_HEXA20 :
1862 for (int hexa=0;hexa<nb_entity_type;hexa++)
1864 int hexa_index = (type%100)*hexa;
1866 int N1 = global_connectivity[hexa_index]-1;
1867 int N2 = global_connectivity[hexa_index+1]-1;
1868 int N3 = global_connectivity[hexa_index+2]-1;
1869 int N4 = global_connectivity[hexa_index+3]-1;
1870 int N5 = global_connectivity[hexa_index+4]-1;
1871 int N6 = global_connectivity[hexa_index+5]-1;
1872 int N7 = global_connectivity[hexa_index+6]-1;
1873 int N8 = global_connectivity[hexa_index+7]-1;
1875 double x1 = coord[dim_space*N1];
1876 double x2 = coord[dim_space*N2];
1877 double x3 = coord[dim_space*N3];
1878 double x4 = coord[dim_space*N4];
1879 double x5 = coord[dim_space*N5];
1880 double x6 = coord[dim_space*N6];
1881 double x7 = coord[dim_space*N7];
1882 double x8 = coord[dim_space*N8];
1884 double y1 = coord[dim_space*N1+1];
1885 double y2 = coord[dim_space*N2+1];
1886 double y3 = coord[dim_space*N3+1];
1887 double y4 = coord[dim_space*N4+1];
1888 double y5 = coord[dim_space*N5+1];
1889 double y6 = coord[dim_space*N6+1];
1890 double y7 = coord[dim_space*N7+1];
1891 double y8 = coord[dim_space*N8+1];
1893 double z1 = coord[dim_space*N1+2];
1894 double z2 = coord[dim_space*N2+2];
1895 double z3 = coord[dim_space*N3+2];
1896 double z4 = coord[dim_space*N4+2];
1897 double z5 = coord[dim_space*N5+2];
1898 double z6 = coord[dim_space*N6+2];
1899 double z7 = coord[dim_space*N7+2];
1900 double z8 = coord[dim_space*N8+2];
1902 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1903 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1904 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1906 barycenter[3*index] = xbarycenter1 ;
1907 barycenter[3*index+1] = xbarycenter2 ;
1908 barycenter[3*index+2] = xbarycenter3 ;
1914 case MED_TRIA3 : case MED_TRIA6 :
1916 for (int tria=0;tria<nb_entity_type;tria++)
1918 int tria_index = (type%100)*tria;
1920 int N1 = global_connectivity[tria_index]-1;
1921 int N2 = global_connectivity[tria_index+1]-1;
1922 int N3 = global_connectivity[tria_index+2]-1;
1924 double x1 = coord[dim_space*N1];
1925 double x2 = coord[dim_space*N2];
1926 double x3 = coord[dim_space*N3];
1928 double y1 = coord[dim_space*N1+1];
1929 double y2 = coord[dim_space*N2+1];
1930 double y3 = coord[dim_space*N3+1];
1932 xbarycenter1 = (x1 + x2 + x3)/3.0;
1933 xbarycenter2 = (y1 + y2 + y3)/3.0;
1937 barycenter[2*index] = xbarycenter1 ;
1938 barycenter[2*index+1] = xbarycenter2 ;
1943 coord[dim_space*N1+2];
1945 coord[dim_space*N2+2];
1947 coord[dim_space*N3+2];
1949 xbarycenter3 = (z1 + z2 + z3)/3.0;
1951 barycenter[3*index] = xbarycenter1 ;
1952 barycenter[3*index+1] = xbarycenter2 ;
1953 barycenter[3*index+2] = xbarycenter3 ;
1960 case MED_QUAD4 : case MED_QUAD8 :
1962 for (int quad=0;quad<nb_entity_type;quad++)
1964 int quad_index = (type%100)*quad;
1966 int N1 = global_connectivity[quad_index]-1;
1967 int N2 = global_connectivity[quad_index+1]-1;
1968 int N3 = global_connectivity[quad_index+2]-1;
1969 int N4 = global_connectivity[quad_index+3]-1;
1971 double x1 = coord[dim_space*N1];
1972 double x2 = coord[dim_space*N2];
1973 double x3 = coord[dim_space*N3];
1974 double x4 = coord[dim_space*N4];
1976 double y1 = coord[dim_space*N1+1];
1977 double y2 = coord[dim_space*N2+1];
1978 double y3 = coord[dim_space*N3+1];
1979 double y4 = coord[dim_space*N4+1];
1981 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1982 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1986 barycenter[2*index] = xbarycenter1 ;
1987 barycenter[2*index+1] = xbarycenter2 ;
1991 double z1 = coord[dim_space*N1+2];
1992 double z2 = coord[dim_space*N2+2];
1993 double z3 = coord[dim_space*N3+2];
1994 double z4 = coord[dim_space*N4+2];
1996 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1998 barycenter[3*index] = xbarycenter1 ;
1999 barycenter[3*index+1] = xbarycenter2 ;
2000 barycenter[3*index+2] = xbarycenter3 ;
2006 case MED_SEG2 : case MED_SEG3 :
2008 for (int edge=0;edge<nb_entity_type;edge++)
2010 int edge_index = (type%100)*edge;
2012 int N1 = global_connectivity[edge_index]-1;
2013 int N2 = global_connectivity[edge_index+1]-1;
2015 double x1 = coord[dim_space*N1];
2016 double x2 = coord[dim_space*N2];
2018 double y1 = coord[dim_space*N1+1];
2019 double y2 = coord[dim_space*N2+1];
2021 xbarycenter1 = (x1 + x2)/2.0;
2022 xbarycenter2 = (y1 + y2)/2.0;
2026 barycenter[2*index] = xbarycenter1 ;
2027 barycenter[2*index+1] = xbarycenter2 ;
2031 double z1 = coord[dim_space*N1+2];
2032 double z2 = coord[dim_space*N2+2];
2034 xbarycenter3 = (z1 + z2)/2.0;
2036 barycenter[3*index] = xbarycenter1 ;
2037 barycenter[3*index+1] = xbarycenter2 ;
2038 barycenter[3*index+2] = xbarycenter3 ;
2045 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
2049 if (!onAll) delete [] global_connectivity ;
2052 Barycenter->setValue(MED_FULL_INTERLACE,barycenter);
2054 delete[] barycenter ;
2061 //=======================================================================
2062 //function : checkGridFillCoords
2063 //purpose : if this->_isAGrid, assure that _coordinate is filled
2064 //=======================================================================
2066 // inline void MESH::checkGridFillCoords() const
2069 // ((GRID *) this)->fillCoordinates();
2072 //=======================================================================
2073 //function : checkGridFillConnectivity
2074 //purpose : if this->_isAGrid, assure that _connectivity is filled
2075 //=======================================================================
2077 // inline void MESH::checkGridFillConnectivity() const
2080 // ((GRID *) this)->fillConnectivity();
2083 bool MESH::isEmpty() const
2085 bool notempty = _name != "" || _coordinate != NULL || _connectivity != NULL ||
2086 _spaceDimension !=MED_INVALID || _meshDimension !=MED_INVALID ||
2087 _numberOfNodes !=MED_INVALID || _groupNode.size() != 0 ||
2088 _familyNode.size() != 0 || _groupCell.size() != 0 ||
2089 _familyCell.size() != 0 || _groupFace.size() != 0 ||
2090 _familyFace.size() != 0 || _groupEdge.size() != 0 ||
2091 _familyEdge.size() != 0 || _isAGrid != 0 ;
2095 void MESH::read(int index)
2097 const char * LOC = "MESH::read(int index=0) : ";
2100 if (_drivers[index]) {
2101 _drivers[index]->open();
2102 _drivers[index]->read();
2103 _drivers[index]->close();
2106 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2107 << "The index given is invalid, index must be between 0 and |"
2112 // ((GRID *) this)->fillMeshAfterRead();
2116 //=======================================================================
2117 //function : getSkin
2119 //=======================================================================
2121 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2123 const char * LOC = "MESH::getSkin : " ;
2126 if (this != Support3D->getMesh())
2127 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
2128 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2129 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2131 // well, all rigth !
2132 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2133 mySupport->setAll(false);
2135 list<int> myElementsList ;
2138 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2139 if (Support3D->isOnAllElements())
2141 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2142 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2143 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2145 int cellNb1 = myConnectivityValue [i];
2146 int cellNb2 = myConnectivityValue [i+1];
2147 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2148 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2150 myElementsList.push_back( j ) ;
2157 map<int,int> FaceNbEncounterNb;
2159 int * myConnectivityValue = const_cast <int *>
2160 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2161 MED_CELL, MED_ALL_ELEMENTS));
2162 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2163 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2164 int nbCells = Support3D->getnumber()->getLength();
2165 for (i=0; i<nbCells; ++i)
2167 int cellNb = myCellNbs [ i ];
2168 int faceFirst = myConnectivityIndex[ cellNb-1 ];
2169 int faceLast = myConnectivityIndex[ cellNb ];
2170 for (j = faceFirst; j < faceLast; ++j)
2172 int faceNb = abs( myConnectivityValue [ j-1 ] );
2173 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2174 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2175 FaceNbEncounterNb[ faceNb ] = 1;
2177 FaceNbEncounterNb[ faceNb ] ++;
2180 map<int,int>::iterator FaceNbEncounterNbItr;
2181 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2182 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2183 FaceNbEncounterNbItr ++)
2184 if ((*FaceNbEncounterNbItr).second == 1)
2186 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2190 // Well, we must know how many geometric type we have found
2191 int * myListArray = new int[size] ;
2193 list<int>::iterator myElementsListIt ;
2194 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2195 myListArray[id]=(*myElementsListIt) ;
2199 int numberOfGeometricType ;
2200 medGeometryElement* geometricType ;
2201 int * numberOfGaussPoint ;
2202 int * geometricTypeNumber ;
2203 int * numberOfEntities ;
2204 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2205 int * mySkyLineArrayIndex ;
2207 int numberOfType = getNumberOfTypes(MED_FACE) ;
2208 if (numberOfType == 1) { // wonderfull : it's easy !
2209 numberOfGeometricType = 1 ;
2210 geometricType = new medGeometryElement[1] ;
2211 const medGeometryElement * allType = getTypes(MED_FACE);
2212 geometricType[0] = allType[0] ;
2213 numberOfGaussPoint = new int[1] ;
2214 numberOfGaussPoint[0] = 1 ;
2215 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2216 geometricTypeNumber[0] = 0 ;
2217 numberOfEntities = new int[1] ;
2218 numberOfEntities[0] = size ;
2219 mySkyLineArrayIndex = new int[2] ;
2220 mySkyLineArrayIndex[0]=1 ;
2221 mySkyLineArrayIndex[1]=1+size ;
2224 map<medGeometryElement,int> theType ;
2225 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2226 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2227 if (theType.find(myType) != theType.end() )
2228 theType[myType]+=1 ;
2232 numberOfGeometricType = theType.size() ;
2233 geometricType = new medGeometryElement[numberOfGeometricType] ;
2234 //const medGeometryElement * allType = getTypes(MED_FACE); !! UNUSED VARIABLE !!
2235 numberOfGaussPoint = new int[numberOfGeometricType] ;
2236 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2237 numberOfEntities = new int[numberOfGeometricType] ;
2238 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2240 mySkyLineArrayIndex[0]=1 ;
2241 map<medGeometryElement,int>::iterator theTypeIt ;
2242 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2243 geometricType[index] = (*theTypeIt).first ;
2244 numberOfGaussPoint[index] = 1 ;
2245 geometricTypeNumber[index] = 0 ;
2246 numberOfEntities[index] = (*theTypeIt).second ;
2247 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2251 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2252 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2254 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2255 mySupport->setGeometricType(geometricType) ;
2256 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2257 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2258 mySupport->setNumberOfElements(numberOfEntities) ;
2259 mySupport->setTotalNumberOfElements(size) ;
2260 mySupport->setNumber(mySkyLineArray) ;
2262 delete[] numberOfEntities;
2263 delete[] geometricTypeNumber;
2264 delete[] numberOfGaussPoint;
2265 delete[] geometricType;
2266 delete[] mySkyLineArrayIndex;
2267 delete[] myListArray;
2268 // delete mySkyLineArray;
2276 return a SUPPORT pointer on the union of all SUPPORTs in Supports.
2277 You should delete this pointer after use to avoid memory leaks.
2279 SUPPORT * MESH::mergeSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2281 const char * LOC = "MESH:::mergeSupports(const vector<SUPPORT *> ) : " ;
2284 SUPPORT * returnedSupport;
2285 string returnedSupportName;
2286 string returnedSupportDescription;
2287 char * returnedSupportNameChar;
2288 char * returnedSupportDescriptionChar;
2289 int size = Supports.size();
2293 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2294 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2296 returnedSupport = new SUPPORT(*obj);
2298 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2299 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2301 returnedSupportNameChar = new char[lenName];
2302 returnedSupportDescriptionChar = new char[lenDescription];
2304 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2305 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2306 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2307 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2308 (Supports[0]->getDescription()).c_str());
2310 returnedSupportName = string(returnedSupportNameChar);
2311 returnedSupportDescription = string(returnedSupportDescriptionChar);
2313 returnedSupport->setName(returnedSupportName);
2314 returnedSupport->setDescription(returnedSupportDescription);
2316 delete [] returnedSupportNameChar;
2317 delete [] returnedSupportDescriptionChar;
2321 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2322 returnedSupport = new SUPPORT(*obj);
2324 int lenName = strlen((Supports[0]->getName()).c_str()) + 9 + 1;
2325 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 9 + 1;
2327 for(int i = 1;i<size;i++)
2329 obj = const_cast <SUPPORT *> (Supports[i]);
2330 returnedSupport->blending(obj);
2334 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2335 lenDescription = lenDescription + 5 +
2336 strlen((Supports[i]->getDescription()).c_str());
2340 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2341 lenDescription = lenDescription + 2 +
2342 strlen((Supports[i]->getDescription()).c_str());
2346 returnedSupportNameChar = new char[lenName];
2347 returnedSupportDescriptionChar = new char[lenDescription];
2349 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Merge of ");
2350 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Merge of ");
2352 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2353 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2354 (Supports[0]->getDescription()).c_str());
2356 for(int i = 1;i<size;i++)
2360 returnedSupportNameChar = strcat(returnedSupportNameChar," and ");
2361 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar," and ");
2363 returnedSupportNameChar = strcat(returnedSupportNameChar,
2364 (Supports[i]->getName()).c_str());
2365 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2366 (Supports[i]->getDescription()).c_str());
2370 returnedSupportNameChar = strcat(returnedSupportNameChar,", ");
2371 returnedSupportNameChar = strcat(returnedSupportNameChar,
2372 (Supports[i]->getName()).c_str());
2374 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,", ");
2375 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2376 (Supports[i]->getDescription()).c_str());
2380 returnedSupportName = string(returnedSupportNameChar);
2381 returnedSupport->setName(returnedSupportName);
2383 returnedSupportDescription = string(returnedSupportDescriptionChar);
2384 returnedSupport->setDescription(returnedSupportDescription);
2386 delete [] returnedSupportNameChar;
2387 delete [] returnedSupportDescriptionChar;
2391 return returnedSupport;
2395 return a SUPPORT pointer on the intersection of all SUPPORTs in Supports.
2396 The (SUPPORT *) NULL pointer is returned if the intersection is empty.
2397 You should delete this pointer after use to avois memory leaks.
2399 SUPPORT * MESH::intersectSupports(const vector<SUPPORT *> Supports) const throw (MEDEXCEPTION)
2401 const char * LOC = "MESH:::intersectSupports(const vector<SUPPORT *> ) : " ;
2404 SUPPORT * returnedSupport;
2405 string returnedSupportName;
2406 string returnedSupportDescription;
2407 char * returnedSupportNameChar;
2408 char * returnedSupportDescriptionChar;
2409 int size = Supports.size();
2413 MESSAGE(LOC <<" there is only one SUPPORT in the argument list, the method return a copy of this object !");
2414 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2416 returnedSupport = new SUPPORT(*obj);
2418 int lenName = strlen((Supports[0]->getName()).c_str()) + 8 + 1;
2419 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 8 + 1;
2421 returnedSupportNameChar = new char[lenName];
2422 returnedSupportDescriptionChar = new char[lenDescription];
2424 returnedSupportNameChar = strcpy(returnedSupportNameChar,"Copy of ");
2425 returnedSupportNameChar = strcat(returnedSupportNameChar,(Supports[0]->getName()).c_str());
2426 returnedSupportDescriptionChar = strcpy(returnedSupportDescriptionChar,"Copy of ");
2427 returnedSupportDescriptionChar = strcat(returnedSupportDescriptionChar,
2428 (Supports[0]->getDescription()).c_str());
2430 returnedSupportName = string(returnedSupportNameChar);
2431 returnedSupportDescription = string(returnedSupportDescriptionChar);
2433 returnedSupport->setName(returnedSupportName);
2434 returnedSupport->setDescription(returnedSupportDescription);
2436 delete [] returnedSupportNameChar;
2437 delete [] returnedSupportDescriptionChar;
2441 SUPPORT * obj = const_cast <SUPPORT *> (Supports[0]);
2442 returnedSupport = new SUPPORT(*obj);
2444 int lenName = strlen((Supports[0]->getName()).c_str()) + 16 + 1;
2445 int lenDescription = strlen((Supports[0]->getDescription()).c_str()) + 16 + 1;
2447 for(int i = 1;i<size;i++)
2449 obj = const_cast <SUPPORT *> (Supports[i]);
2450 returnedSupport->intersecting(obj);
2454 lenName = lenName + 5 + strlen((Supports[i]->getName()).c_str());
2455 lenDescription = lenDescription + 5 +
2456 strlen((Supports[i]->getDescription()).c_str());
2460 lenName = lenName + 2 + strlen((Supports[i]->getName()).c_str());
2461 lenDescription = lenDescription + 2 +
2462 strlen((Supports[i]->getDescription()).c_str());
2466 if(returnedSupport != (SUPPORT *) NULL)
2468 returnedSupportNameChar = new char[lenName];
2469 returnedSupportDescriptionChar = new char[lenDescription];
2471 returnedSupportNameChar = strcpy(returnedSupportNameChar,
2472 "Intersection of ");
2473 returnedSupportDescriptionChar =
2474 strcpy(returnedSupportDescriptionChar,"Intersection of ");
2476 returnedSupportNameChar = strcat(returnedSupportNameChar,
2477 (Supports[0]->getName()).c_str());
2478 returnedSupportDescriptionChar =
2479 strcat(returnedSupportDescriptionChar,
2480 (Supports[0]->getDescription()).c_str());
2482 for(int i = 1;i<size;i++)
2486 returnedSupportNameChar = strcat(returnedSupportNameChar,
2488 returnedSupportDescriptionChar =
2489 strcat(returnedSupportDescriptionChar," and ");
2491 returnedSupportNameChar =
2492 strcat(returnedSupportNameChar,
2493 (Supports[i]->getName()).c_str());
2494 returnedSupportDescriptionChar =
2495 strcat(returnedSupportDescriptionChar,
2496 (Supports[i]->getDescription()).c_str());
2500 returnedSupportNameChar = strcat(returnedSupportNameChar,
2502 returnedSupportNameChar =
2503 strcat(returnedSupportNameChar,
2504 (Supports[i]->getName()).c_str());
2506 returnedSupportDescriptionChar =
2507 strcat(returnedSupportDescriptionChar,", ");
2508 returnedSupportDescriptionChar =
2509 strcat(returnedSupportDescriptionChar,
2510 (Supports[i]->getDescription()).c_str());
2514 returnedSupportName = string(returnedSupportNameChar);
2515 returnedSupport->setName(returnedSupportName);
2517 returnedSupportDescription = string(returnedSupportDescriptionChar);
2518 returnedSupport->setDescription(returnedSupportDescription);
2520 delete [] returnedSupportNameChar;
2521 delete [] returnedSupportDescriptionChar;
2526 return returnedSupport;
2530 // internal helper type
2534 std::vector<int> groups;
2535 MED_EN::medGeometryElement geometricType;
2538 // Create families from groups
2539 void MESH::createFamilies()
2541 int nbFam=0; // count the families we create, used as identifier ???????????
2543 int idFamNode = 0; // identifier for node families
2544 int idFamElement = 0; // identifier for cell, face or edge families
2546 // Main loop on mesh's entities
2547 for (medEntityMesh entity=MED_CELL; entity!=MED_ALL_ENTITIES; ++entity)
2549 int numberofgroups = getNumberOfGroups(entity);
2551 continue; // no groups for this entity
2553 vector< vector<FAMILY*> > whichFamilyInGroup(numberofgroups); // this container is used to update groups at the end
2555 // make myFamilies points to the member corresponding to entity
2556 vector<FAMILY*>* myFamilies;
2560 myFamilies = & _familyCell;
2563 myFamilies = & _familyFace;
2566 myFamilies = & _familyEdge;
2569 myFamilies = & _familyNode;
2573 vector<GROUP*> myGroups=getGroups(entity); // get a copy of the groups ptr for the entity
2574 // get a copy of the (old) family ptrs before clearing
2575 vector<FAMILY*> myOldFamilies=getFamilies(entity);
2576 myFamilies->clear();
2579 // 1 - Create a vector containing for each cell (of the entity) an information structure
2580 // giving geometric type and the groups it belong to
2582 med_int numberOfTypes=getNumberOfTypes(entity);
2583 const int * index=getGlobalNumberingIndex(entity);
2584 const medGeometryElement* geometricTypes=_connectivity->getGeometricTypes(entity); // pb avec entity=MED_NODE???
2585 med_int numberOfCells=index[numberOfTypes]-1; // total number of cells for that entity
2586 SCRUTE(numberOfTypes);
2587 SCRUTE(numberOfCells);
2588 vector< _cell > tab_cell(numberOfCells);
2589 for(med_int t=0; t!=numberOfTypes; ++t)
2590 for(int n=index[t]-1; n!=index[t+1]-1; ++n)
2591 tab_cell[n].geometricType=geometricTypes[t];
2594 // 2 - Scan cells in groups and update in tab_cell the container of groups a cell belong to
2596 for (unsigned g=0; g!=myGroups.size(); ++g)
2598 // scan cells that belongs to the group
2599 const med_int* groupCells=myGroups[g]->getnumber()->getValue();
2600 int nbCells=myGroups[g]->getnumber()->getLength();
2601 for(int c=0; c!=nbCells; ++c)
2602 tab_cell[groupCells[c]-1].groups.push_back(g);
2606 // 3 - Scan the cells vector, genarate family name, and create a map associating the family names
2607 // whith the vector of contained cells
2609 map< string,vector<int> > tab_families;
2610 map< string,vector<int> >::iterator fam;
2611 for(int n=0; n!=numberOfCells; ++n)
2613 ostringstream key; // to generate the name of the family
2615 if(tab_cell[n].groups.empty()) // this cell don't belong to any group
2616 key << "_NONE" << entity;
2618 for(vector<int>::const_iterator it=tab_cell[n].groups.begin(); it!=tab_cell[n].groups.end(); ++it)
2620 string groupName=myGroups[*it]->getName();
2621 if(groupName.empty())
2624 key << "_" << groupName;
2627 tab_families[key.str()].push_back(n+1); // fill the vector of contained cells associated whith the family
2628 /* fam = tab_families.find(key.str());
2629 if( fam != tab_families.end())
2630 fam->second.push_back(n+1); // +1 : convention Fortran de MED
2633 vector<int> newfamily;
2634 newfamily.push_back(n+1); // +1 : convention Fortran de MED
2635 tab_families.insert(make_pair(key.str(),newfamily));
2640 // 4 - Scan the family map, create MED Families, check if it already exist.
2642 for( fam=tab_families.begin(); fam!=tab_families.end(); ++fam)
2644 vector<medGeometryElement> tab_types_geometriques;
2645 medGeometryElement geometrictype=MED_NONE;
2646 vector<int> tab_index_types_geometriques;
2647 vector<int> tab_nombres_elements;
2649 // scan family cells and fill the tab that are needed by the create a MED FAMILY
2650 for( int i=0; i!=fam->second.size(); ++i)
2652 int ncell=fam->second[i]-1;
2653 if(tab_cell[ncell].geometricType != geometrictype)
2655 // new geometric type -> we store it and complete index tabs
2656 if(!tab_index_types_geometriques.empty())
2657 tab_nombres_elements.push_back(i+1-tab_index_types_geometriques.back());
2658 tab_types_geometriques.push_back( (geometrictype=tab_cell[ncell].geometricType));
2659 tab_index_types_geometriques.push_back(i+1);
2662 // store and complete index tabs for the last geometric type
2663 tab_nombres_elements.push_back(fam->second.size()+1-tab_index_types_geometriques.back());
2664 tab_index_types_geometriques.push_back(fam->second.size()+1);
2666 // create a empty MED FAMILY and fill it with the tabs we constructed
2667 FAMILY* newFam = new FAMILY();
2668 newFam->setTotalNumberOfElements(fam->second.size());
2669 newFam->setName(fam->first);
2670 newFam->setMesh(this);
2671 newFam->setNumberOfGeometricType(tab_types_geometriques.size());
2672 newFam->setGeometricType(&tab_types_geometriques[0]); // we know the tab is not empy
2673 newFam->setNumberOfElements(&tab_nombres_elements[0]);
2674 newFam->setNumber(&tab_index_types_geometriques[0],&fam->second[0]);
2675 newFam->setEntity(entity);
2676 newFam->setAll(false);
2688 idFam = -idFamElement;
2692 idFam = -idFamElement;
2696 idFam = -idFamElement;
2700 newFam->setIdentifier(idFam);
2702 // Update links between families and groups
2704 int ncell1=fam->second[0]-1; // number of first cell in family
2705 int numberOfGroups=tab_cell[ncell1].groups.size(); // number of groups in the family
2708 newFam->setNumberOfGroups(numberOfGroups);
2709 string * groupNames=new string[numberOfGroups];
2711 // iterate on the groups the family belongs to
2712 vector<int>::const_iterator it=tab_cell[ncell1].groups.begin();
2713 for(int ng=0 ; it!=tab_cell[ncell1].groups.end(); ++it, ++ng)
2715 whichFamilyInGroup[*it].push_back(newFam);
2716 groupNames[ng]=myGroups[*it]->getName();
2718 newFam->setGroupsNames(groupNames);
2721 int sizeOfFamVect = myFamilies->size();
2723 MESSAGE(" MESH::createFamilies() entity " << entity << " size " << sizeOfFamVect);
2725 myFamilies->push_back(newFam);
2728 // delete old families
2729 for (int i=0;i<myOldFamilies.size();i++)
2730 delete myOldFamilies[i] ;
2732 // update references in groups
2733 for (int i=0;i<myGroups.size();i++)
2735 myGroups[i]->setNumberOfFamilies(whichFamilyInGroup[i].size());
2736 myGroups[i]->setFamilies(whichFamilyInGroup[i]);
2739 // re-scan the cells vector, and fill the family vector with cells.
2740 // creation of support, check if it already exist.