1 // MED MEDMEM : MED files in memory
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : MEDMEM_Mesh.cxx
38 #include "MEDMEM_DriversDef.hxx"
39 #include "MEDMEM_Field.hxx"
40 #include "MEDMEM_Mesh.hxx"
42 #include "MEDMEM_Support.hxx"
43 #include "MEDMEM_Family.hxx"
44 #include "MEDMEM_Group.hxx"
45 #include "MEDMEM_Coordinate.hxx"
46 #include "MEDMEM_Connectivity.hxx"
47 #include "MEDMEM_CellModel.hxx"
48 #include "MEDMEM_Grid.hxx"
50 //update Families with content list
51 //int family_count(int family_number, int count, int * entities_number, int * entities_list) ;
53 // ------- Drivers Management Part
55 // MESH::INSTANCE_DE<MED_MESH_RDONLY_DRIVER> MESH::inst_med_rdonly ;
56 //const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med_rdonly , &MESH::inst_med_rdwr } ;
58 // Add a similar line for your personnal driver (step 3)
59 MESH::INSTANCE_DE<MED_MESH_RDWR_DRIVER> MESH::inst_med ;
60 MESH::INSTANCE_DE<GIBI_MESH_RDWR_DRIVER> MESH::inst_gibi ;
61 // Add your own driver in the driver list (step 4)
62 // Note the list must be coherent with the driver type list defined in MEDMEM_DRIVER.hxx.
63 const MESH::INSTANCE * const MESH::instances[] = { &MESH::inst_med, &MESH::inst_gibi } ;
65 /*! Add a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>. The meshname used in the file
66 is <driverName>. addDriver returns an int handler. */
67 int MESH::addDriver(driverTypes driverType,
68 const string & fileName/*="Default File Name.med"*/,const string & driverName/*="Default Mesh Name"*/) {
70 const char * LOC = "MESH::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Mesh Name\") : ";
77 driver = instances[driverType]->run(fileName, this) ;
78 _drivers.push_back(driver);
79 current = _drivers.size()-1;
81 _drivers[current]->setMeshName(driverName);
88 /*! Add an existing MESH driver. */
89 int MESH::addDriver(GENDRIVER & driver) {
90 const char * LOC = "MESH::addDriver(GENDRIVER &) : ";
93 // A faire : Vérifier que le driver est de type MESH.
94 GENDRIVER * newDriver = driver.copy() ;
96 _drivers.push_back(newDriver);
97 return _drivers.size()-1;
102 /*! Remove an existing MESH driver. */
103 void MESH::rmDriver (int index/*=0*/) {
104 const char * LOC = "MESH::rmDriver (int index=0): ";
107 if ( _drivers[index] ) {
108 //_drivers.erase(&_drivers[index]);
110 MESSAGE ("detruire");
113 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
114 << "The index given is invalid, index must be between 0 and |"
123 // ------ End of Drivers Management Part
128 const char * LOC = "MESH::init(): ";
132 string _name = "NOT DEFINED"; // A POSITIONNER EN FCT DES IOS ?
134 COORDINATE * _coordinate = (COORDINATE *) NULL;
135 CONNECTIVITY * _connectivity = (CONNECTIVITY *) NULL;
137 _spaceDimension = MED_INVALID; // 0 ?!?
138 _meshDimension = MED_INVALID;
139 _numberOfNodes = MED_INVALID;
141 _numberOfNodesFamilies = 0; // MED_INVALID ?!?
142 _numberOfCellsFamilies = 0;
143 _numberOfFacesFamilies = 0;
144 _numberOfEdgesFamilies = 0;
146 _numberOfNodesGroups = 0; // MED_INVALID ?!?
147 _numberOfCellsGroups = 0;
148 _numberOfFacesGroups = 0;
149 _numberOfEdgesGroups = 0;
154 /*! Create an empty MESH. */
155 MESH::MESH():_coordinate(NULL),_connectivity(NULL) {
162 _isAGrid = m._isAGrid;
164 if (m._coordinate != NULL)
165 _coordinate = new COORDINATE(* m._coordinate);
167 _coordinate = (COORDINATE *) NULL;
169 if (m._connectivity != NULL)
170 _connectivity = new CONNECTIVITY(* m._connectivity);
172 _connectivity = (CONNECTIVITY *) NULL;
174 _spaceDimension = m._spaceDimension;
175 _meshDimension = m._meshDimension;
176 _numberOfNodes = m._numberOfNodes;
178 _numberOfNodesFamilies = m._numberOfNodesFamilies;
179 _familyNode = m._familyNode;
180 for (int i=0; i<m._numberOfNodesFamilies; i++)
182 _familyNode[i] = new FAMILY(* m._familyNode[i]);
183 _familyNode[i]->setMesh(this);
186 _numberOfCellsFamilies = m._numberOfCellsFamilies;
187 _familyCell = m._familyCell;
188 for (int i=0; i<m._numberOfCellsFamilies; i++)
190 _familyCell[i] = new FAMILY(* m._familyCell[i]);
191 _familyCell[i]->setMesh(this);
194 _numberOfFacesFamilies = m._numberOfFacesFamilies;
195 _familyFace = m._familyFace;
196 for (int i=0; i<m._numberOfFacesFamilies; i++)
198 _familyFace[i] = new FAMILY(* m._familyFace[i]);
199 _familyFace[i]->setMesh(this);
202 _numberOfEdgesFamilies = m._numberOfEdgesFamilies;
203 _familyEdge = m._familyEdge;
204 for (int i=0; i<m._numberOfEdgesFamilies; i++)
206 _familyEdge[i] = new FAMILY(* m._familyEdge[i]);
207 _familyEdge[i]->setMesh(this);
210 _numberOfNodesGroups = m._numberOfNodesGroups;
211 _groupNode = m._groupNode;
212 for (int i=0; i<m._numberOfNodesGroups; i++)
214 _groupNode[i] = new GROUP(* m._groupNode[i]);
215 _groupNode[i]->setMesh(this);
218 _numberOfCellsGroups = m._numberOfCellsGroups;
219 _groupCell = m._groupCell;
220 for (int i=0; i<m._numberOfCellsGroups; i++)
222 _groupCell[i] = new GROUP(* m._groupCell[i]);
223 _groupCell[i]->setMesh(this);
226 _numberOfFacesGroups = m._numberOfFacesGroups;
227 _groupFace = m._groupFace;
228 for (int i=0; i<m._numberOfFacesGroups; i++)
230 _groupFace[i] = new GROUP(* m._groupFace[i]);
231 _groupFace[i]->setMesh(this);
234 _numberOfEdgesGroups = m._numberOfEdgesGroups;
235 _groupEdge = m._groupEdge;
236 for (int i=0; i<m._numberOfEdgesGroups; i++)
238 _groupEdge[i] = new GROUP(* m._groupEdge[i]);
239 _groupEdge[i]->setMesh(this);
242 //_drivers = m._drivers; //Recopie des drivers?
247 MESSAGE("MESH::~MESH() : Destroying the Mesh");
248 if (_coordinate != ((COORDINATE *) NULL)) delete _coordinate ;
249 if (_connectivity != ((CONNECTIVITY *) NULL)) delete _connectivity ;
251 size = _familyNode.size() ;
252 for (int i=0;i<size;i++)
253 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
254 //delete _familyNode[i] ;
255 size = _familyCell.size() ;
256 for (int i=0;i<size;i++)
257 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
258 //delete _familyCell[i] ;
259 size = _familyFace.size() ;
260 for (int i=0;i<size;i++)
261 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
262 //delete _familyFace[i] ;
263 size = _familyEdge.size() ;
264 for (int i=0;i<size;i++)
265 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
266 //delete _familyEdge[i] ;
267 size = _groupNode.size() ;
268 for (int i=0;i<size;i++)
269 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
270 //delete _groupNode[i] ;
271 size = _groupCell.size() ;
272 for (int i=0;i<size;i++)
273 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
274 //delete _groupCell[i] ;
275 size = _groupFace.size() ;
276 for (int i=0;i<size;i++)
277 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
278 //delete _groupFace[i] ;
279 size = _groupEdge.size() ;
280 for (int i=0;i<size;i++)
281 ;// mpv such destructors call is problemmatic for ALLIANCES algorithms
282 //delete _groupEdge[i] ;
284 MESSAGE("In this object MESH there is(are) " << _drivers.size() << " driver(s)");
286 for (int index=0; index < _drivers.size(); index++ )
288 SCRUTE(_drivers[index]);
289 if ( _drivers[index] != NULL) delete _drivers[index];
294 MESH & MESH::operator=(const MESH &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 ;
332 /*! Create a MESH object using a MESH driver of type (MED_DRIVER, ....) associated with file <fileName>.
333 The meshname <driverName> must already exists in the file.*/
334 MESH::MESH(driverTypes driverType, const string & fileName/*=""*/, const string & driverName/*=""*/) {
335 const char * LOC ="MESH::MESH(driverTypes driverType, const string & fileName="", const string & driverName="") : ";
343 MED_MESH_RDONLY_DRIVER myDriver(fileName,this) ;
344 myDriver.setMeshName(driverName);
345 current = addDriver(myDriver);
346 // current = addDriver(driverType,fileName,driverName);
347 // switch(_drivers[current]->getAccessMode() ) {
348 // case MED_WRONLY : {
349 // MESSAGE("MESH::MESH(driverTypes driverType, .....) : driverType must not be a MED_WRONLY accessMode");
350 // rmDriver(current);
355 _drivers[current]->open();
356 _drivers[current]->read();
357 _drivers[current]->close();
360 ((GRID *) this)->fillMeshAfterRead();
366 // Node MESH::Donne_Barycentre(const Maille &m) const
368 // Node N=node[m[0]];
369 // for (int i=1;i<m.donne_nbr_sommets();i++) N=N+node[m[i]];
370 // N=N*(1.0/m.donne_nbr_sommets());
374 ostream & operator<<(ostream &os, MESH &myMesh)
376 int spacedimension = myMesh.getSpaceDimension();
377 int meshdimension = myMesh.getMeshDimension();
378 int numberofnodes = myMesh.getNumberOfNodes();
380 os << "Space Dimension : " << spacedimension << endl << endl;
382 os << "Mesh Dimension : " << meshdimension << endl << endl;
384 const double * coordinates = myMesh.getCoordinates(MED_FULL_INTERLACE);
385 os << "SHOW NODES COORDINATES : " << endl;
387 os << "Name :" << endl;
388 const string * coordinatesnames = myMesh.getCoordinatesNames();
389 for(int i=0; i<spacedimension ; i++)
391 os << " - " << coordinatesnames[i] << endl;
393 os << "Unit :" << endl;
394 const string * coordinatesunits = myMesh.getCoordinatesUnits();
395 for(int i=0; i<spacedimension ; i++)
397 os << " - " << coordinatesunits[i] << endl;
399 for(int i=0; i<numberofnodes ; i++)
401 os << "Nodes " << i+1 << " : ";
402 for (int j=0; j<spacedimension ; j++)
403 os << coordinates[i*spacedimension+j] << " ";
407 const CONNECTIVITY * myConnectivity = myMesh.getConnectivityptr();
408 if (!myConnectivity->existConnectivity(MED_NODAL,MED_CELL) && myConnectivity->existConnectivity(MED_DESCENDING,MED_CELL))
410 os << endl << "SHOW CONNECTIVITY (DESCENDING) :" << endl;
411 int numberofelements;
412 const int * connectivity;
413 const int * connectivity_index;
414 myMesh.calculateConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL);
416 numberofelements = myMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
417 connectivity = myMesh.getConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
418 connectivity_index = myMesh.getConnectivityIndex(MED_DESCENDING,MED_CELL);
420 catch (MEDEXCEPTION m) {
421 os << m.what() << endl;
424 for (int j=0;j<numberofelements;j++) {
425 os << "Element "<<j+1<<" : ";
426 for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
427 os << connectivity[k-1]<<" ";
433 int numberoftypes = myMesh.getNumberOfTypes(MED_CELL);
434 const medGeometryElement * types = myMesh.getTypes(MED_CELL);
435 os << endl << "SHOW CONNECTIVITY (NODAL) :" << endl;
436 for (int i=0; i<numberoftypes; i++) {
437 os << "For type " << types[i] << " : " << endl;
438 int numberofelements = myMesh.getNumberOfElements(MED_CELL,types[i]);
439 const int * connectivity = myMesh.getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,types[i]);
440 int numberofnodespercell = types[i]%100;
441 for (int j=0;j<numberofelements;j++){
442 os << "Element "<< j+1 <<" : ";
443 for (int k=0;k<numberofnodespercell;k++)
444 os << connectivity[j*numberofnodespercell+k]<<" ";
451 medEntityMesh entity;
452 os << endl << "SHOW FAMILIES :" << 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 numberoffamilies = myMesh.getNumberOfFamilies(entity);
460 using namespace MED_FR;
461 os << "NumberOfFamilies on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberoffamilies<<endl;
462 using namespace MED_EN;
463 for (int i=1; i<numberoffamilies+1;i++)
465 os << * myMesh.getFamily(entity,i) << endl;
469 os << endl << "SHOW GROUPS :" << endl << endl;
470 for (int k=1; k<=4; k++)
472 if (k==1) entity = MED_NODE;
473 if (k==2) entity = MED_CELL;
474 if (k==3) entity = MED_FACE;
475 if (k==4) entity = MED_EDGE;
476 int numberofgroups = myMesh.getNumberOfGroups(entity);
477 using namespace MED_FR;
478 os << "NumberOfGroups on "<< entNames[(MED_FR::med_entite_maillage) entity] <<" : "<<numberofgroups<<endl;
479 using namespace MED_EN;
480 for (int i=1; i<numberofgroups+1;i++)
482 os << * myMesh.getGroup(entity,i) << endl;
490 Get global number of element which have same connectivity than connectivity argument.
492 It do not take care of connectivity order (3,4,7,10 is same as 7,3,10,4).
494 Return -1 if not found.
496 int MESH::getElementNumber(medConnectivity ConnectivityType, medEntityMesh Entity, medGeometryElement Type, int * connectivity) const
498 const char* LOC="MESH::getElementNumber " ;
501 int numberOfValue ; // size of connectivity array
502 CELLMODEL myType(Type) ;
503 if (ConnectivityType==MED_DESCENDING)
504 numberOfValue = myType.getNumberOfConstituents(1) ; // faces (3D) or edges (2D)
506 numberOfValue = myType.getNumberOfNodes() ; // nodes
508 const int * myReverseConnectivityValue = getReverseConnectivity(ConnectivityType,Entity) ;
509 const int * myReverseConnectivityIndex = getReverseConnectivityIndex(ConnectivityType,Entity) ;
511 // First node or face/edge
512 int indexBegin = myReverseConnectivityIndex[connectivity[0]-1] ;
513 int indexEnd = myReverseConnectivityIndex[connectivity[0]] ;
515 // list to put cells numbers
516 list<int> cellsList ;
517 list<int>::iterator itList ;
518 for (int i=indexBegin; i<indexEnd; i++)
519 cellsList.push_back(myReverseConnectivityValue[i-1]) ;
521 // Foreach node or face/edge in cell, we search cells which are in common.
522 // At the end, we must have only one !
524 for (int i=1; i<numberOfValue; i++) {
525 int connectivity_i = connectivity[i] ;
526 for (itList=cellsList.begin();itList!=cellsList.end();itList++) {
528 for (int j=myReverseConnectivityIndex[connectivity_i-1]; j<myReverseConnectivityIndex[connectivity_i]; j++) {
529 if ((*itList)==myReverseConnectivityValue[j-1]) {
535 itList=cellsList.erase(itList);
536 itList--; // well : rigth if itList = cellsList.begin() ??
541 if (cellsList.size()>1) // we have more than one cell
542 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Found more than one element !")) ;
544 if (cellsList.size()==0)
549 return cellsList.front() ;
554 Return a support which reference all elements on the boundary of mesh.
556 For instance, we could get only face in 3D and edge in 2D.
558 SUPPORT * MESH::getBoundaryElements(medEntityMesh Entity)
561 const char * LOC = "MESH::getBoundaryElements : " ;
564 // actually we could only get face (in 3D) and edge (in 2D)
565 if (_spaceDimension == 3)
566 if (Entity != MED_FACE)
567 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 3D mesh for entity "<<Entity<<" !"));
568 if (_spaceDimension == 2)
569 if (Entity != MED_EDGE)
570 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not defined in 2D mesh for entity "<<Entity<<" !"));
573 SUPPORT * mySupport = new SUPPORT(this,"Boundary",Entity);
574 //mySupport.setDescription("boundary of type ...");
575 mySupport->setAll(false);
578 const int * myConnectivityValue = getReverseConnectivity(MED_DESCENDING) ;
579 const int * myConnectivityIndex = getReverseConnectivityIndex(MED_DESCENDING) ;
580 int numberOf = getNumberOfElements(Entity,MED_ALL_ELEMENTS) ;
581 list<int> myElementsList ;
583 for (int i=0 ; i<numberOf; i++)
584 if (myConnectivityValue[myConnectivityIndex[i]] == 0) {
585 myElementsList.push_back(i+1) ;
588 // Well, we must know how many geometric type we have found
589 int * myListArray = new int[size] ;
591 list<int>::iterator myElementsListIt ;
592 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
593 myListArray[id]=(*myElementsListIt) ;
597 int numberOfGeometricType ;
598 medGeometryElement* geometricType ;
599 int * numberOfGaussPoint ;
600 int * geometricTypeNumber ;
601 int * numberOfElements ;
602 //MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
603 int * mySkyLineArrayIndex ;
605 int numberOfType = getNumberOfTypes(Entity) ;
606 if (numberOfType == 1) { // wonderfull : it's easy !
607 numberOfGeometricType = 1 ;
608 geometricType = new medGeometryElement[1] ;
609 const medGeometryElement * allType = getTypes(Entity);
610 geometricType[0] = allType[0] ;
611 numberOfGaussPoint = new int[1] ;
612 numberOfGaussPoint[0] = 1 ;
613 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
614 geometricTypeNumber[0] = 0 ;
615 numberOfElements = new int[1] ;
616 numberOfElements[0] = size ;
617 mySkyLineArrayIndex = new int[2] ;
618 mySkyLineArrayIndex[0]=1 ;
619 mySkyLineArrayIndex[1]=1+size ;
622 map<medGeometryElement,int> theType ;
623 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
624 medGeometryElement myType = getElementType(Entity,*myElementsListIt) ;
625 if (theType.find(myType) != theType.end() )
630 numberOfGeometricType = theType.size() ;
631 geometricType = new medGeometryElement[numberOfGeometricType] ;
632 const medGeometryElement * allType = getTypes(Entity);
633 numberOfGaussPoint = new int[numberOfGeometricType] ;
634 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
635 numberOfElements = new int[numberOfGeometricType] ;
636 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
638 mySkyLineArrayIndex[0]=1 ;
639 map<medGeometryElement,int>::iterator theTypeIt ;
640 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
641 geometricType[index] = (*theTypeIt).first ;
642 numberOfGaussPoint[index] = 1 ;
643 geometricTypeNumber[index] = 0 ;
644 numberOfElements[index] = (*theTypeIt).second ;
645 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfElements[index] ;
649 //mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
650 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
652 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
653 mySupport->setGeometricType(geometricType) ;
654 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
655 mySupport->setNumberOfElements(numberOfElements) ;
656 mySupport->setTotalNumberOfElements(size) ;
657 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
658 mySupport->setNumber(mySkyLineArray) ;
660 delete[] numberOfElements;
661 delete[] geometricTypeNumber;
662 delete[] numberOfGaussPoint;
663 delete[] geometricType;
664 delete[] mySkyLineArrayIndex;
665 delete[] myListArray;
666 // delete mySkyLineArray;
672 FIELD<double>* MESH::getVolume(const SUPPORT * Support) const throw (MEDEXCEPTION)
674 const char * LOC = "MESH::getVolume(SUPPORT*) : ";
677 // Support must be on 3D elements
679 // Make sure that the MESH class is the same as the MESH class attribut
680 // in the class Support
681 MESH* myMesh = Support->getMesh();
683 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
685 int dim_space = getSpaceDimension();
686 medEntityMesh support_entity = Support->getEntity();
687 bool onAll = Support->isOnAllElements();
689 int nb_type, length_values;
690 const medGeometryElement* types;
692 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
693 const int* global_connectivity;
697 // nb_type = myMesh->getNumberOfTypes(support_entity);
698 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
699 // types = getTypes(support_entity);
703 nb_type = Support->getNumberOfTypes();
704 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
705 types = Support->getTypes();
709 FIELD<double>* Volume = new FIELD<double>::FIELD(Support,1);
710 // double *volume = new double [length_values];
711 Volume->setName("VOLUME");
712 Volume->setDescription("cells volume");
713 Volume->setComponentName(1,"volume");
714 Volume->setComponentDescription(1,"desc-comp");
716 /* string MEDComponentUnit(MED_TAILLE_PNOM,' ');*/
718 string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1]+"*"+myMesh->getCoordinatesUnits()[2];
720 Volume->setMEDComponentUnit(1,MEDComponentUnit);
722 Volume->setValueType(MED_REEL64);
724 Volume->setIterationNumber(0);
725 Volume->setOrderNumber(0);
726 Volume->setTime(0.0);
728 //const double *volume = Volume->getValue(MED_FULL_INTERLACE);
729 MEDARRAY<double> *volume = Volume->getvalue();
732 const double * coord = getCoordinates(MED_FULL_INTERLACE);
734 for (int i=0;i<nb_type;i++)
736 medGeometryElement type = types[i] ;
741 nb_entity_type = getNumberOfElements(support_entity,type);
742 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
746 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all"));
748 nb_entity_type = Support->getNumberOfElements(type);
750 const int * supp_number = Support->getNumber(type);
751 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
752 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
753 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
755 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
756 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
757 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
760 global_connectivity = global_connectivity_tmp ;
765 case MED_TETRA4 : case MED_TETRA10 :
767 for (int tetra=0;tetra<nb_entity_type;tetra++)
769 int tetra_index = (type%100)*tetra;
771 int N1 = global_connectivity[tetra_index]-1;
772 int N2 = global_connectivity[tetra_index+1]-1;
773 int N3 = global_connectivity[tetra_index+2]-1;
774 int N4 = global_connectivity[tetra_index+3]-1;
776 double x1 = coord[dim_space*N1];
777 double x2 = coord[dim_space*N2];
778 double x3 = coord[dim_space*N3];
779 double x4 = coord[dim_space*N4];
781 double y1 = coord[dim_space*N1+1];
782 double y2 = coord[dim_space*N2+1];
783 double y3 = coord[dim_space*N3+1];
784 double y4 = coord[dim_space*N4+1];
786 double z1 = coord[dim_space*N1+2];
787 double z2 = coord[dim_space*N2+2];
788 double z3 = coord[dim_space*N3+2];
789 double z4 = coord[dim_space*N4+2];
791 xvolume = ((x3-x1)*((y2-y1)*(z4-z1) - (z2-z1)*(y4-y1)) -
792 (x2-x1)*((y3-y1)*(z4-z1) - (z3-z1)*(y4-y1)) +
793 (x4-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1)))/6.0;
795 //volume[index] = xvolume ;
796 volume->setIJ(index,1,xvolume) ;
801 case MED_PYRA5 : case MED_PYRA13 :
803 for (int pyra=0;pyra<nb_entity_type;pyra++)
805 int pyra_index = (type%100)*pyra;
807 int N1 = global_connectivity[pyra_index]-1;
808 int N2 = global_connectivity[pyra_index+1]-1;
809 int N3 = global_connectivity[pyra_index+2]-1;
810 int N4 = global_connectivity[pyra_index+3]-1;
811 int N5 = global_connectivity[pyra_index+4]-1;
813 double x1 = coord[dim_space*N1];
814 double x2 = coord[dim_space*N2];
815 double x3 = coord[dim_space*N3];
816 double x4 = coord[dim_space*N4];
817 double x5 = coord[dim_space*N5];
819 double y1 = coord[dim_space*N1+1];
820 double y2 = coord[dim_space*N2+1];
821 double y3 = coord[dim_space*N3+1];
822 double y4 = coord[dim_space*N4+1];
823 double y5 = coord[dim_space*N5+1];
825 double z1 = coord[dim_space*N1+2];
826 double z2 = coord[dim_space*N2+2];
827 double z3 = coord[dim_space*N3+2];
828 double z4 = coord[dim_space*N4+2];
829 double z5 = coord[dim_space*N5+2];
831 xvolume = (((x3-x1)*((y2-y1)*(z5-z1) - (z2-z1)*(y5-y1)) -
832 (x2-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) +
833 (x5-x1)*((y3-y1)*(z2-z1) - (z3-z1)*(y2-y1))) +
834 ((x4-x1)*((y3-y1)*(z5-z1) - (z3-z1)*(y5-y1)) -
835 (x3-x1)*((y4-y1)*(z5-z1) - (z4-z1)*(y5-y1)) +
836 (x5-x1)*((y4-y1)*(z3-z1) - (z4-z1)*(y3-y1)))
839 //volume[index] = xvolume ;
840 volume->setIJ(index,1,xvolume) ;
845 case MED_PENTA6 : case MED_PENTA15 :
847 for (int penta=0;penta<nb_entity_type;penta++)
849 int penta_index = (type%100)*penta;
851 int N1 = global_connectivity[penta_index]-1;
852 int N2 = global_connectivity[penta_index+1]-1;
853 int N3 = global_connectivity[penta_index+2]-1;
854 int N4 = global_connectivity[penta_index+3]-1;
855 int N5 = global_connectivity[penta_index+4]-1;
856 int N6 = global_connectivity[penta_index+5]-1;
858 double x1 = coord[dim_space*N1];
859 double x2 = coord[dim_space*N2];
860 double x3 = coord[dim_space*N3];
861 double x4 = coord[dim_space*N4];
862 double x5 = coord[dim_space*N5];
863 double x6 = coord[dim_space*N6];
865 double y1 = coord[dim_space*N1+1];
866 double y2 = coord[dim_space*N2+1];
867 double y3 = coord[dim_space*N3+1];
868 double y4 = coord[dim_space*N4+1];
869 double y5 = coord[dim_space*N5+1];
870 double y6 = coord[dim_space*N6+1];
872 double z1 = coord[dim_space*N1+2];
873 double z2 = coord[dim_space*N2+2];
874 double z3 = coord[dim_space*N3+2];
875 double z4 = coord[dim_space*N4+2];
876 double z5 = coord[dim_space*N5+2];
877 double z6 = coord[dim_space*N6+2];
879 double a1 = (x2-x3)/2.0, a2 = (y2-y3)/2.0, a3 = (z2-z3)/2.0;
880 double b1 = (x5-x6)/2.0, b2 = (y5-y6)/2.0, b3 = (z5-z6)/2.0;
881 double c1 = (x4-x1)/2.0, c2 = (y4-y1)/2.0, c3 = (z4-z1)/2.0;
882 double d1 = (x5-x2)/2.0, d2 = (y5-y2)/2.0, d3 = (z5-z2)/2.0;
883 double e1 = (x6-x3)/2.0, e2 = (y6-y3)/2.0, e3 = (z6-z3)/2.0;
884 double f1 = (x1-x3)/2.0, f2 = (y1-y3)/2.0, f3 = (z1-z3)/2.0;
885 double h1 = (x4-x6)/2.0, h2 = (y4-y6)/2.0, h3 = (z4-z6)/2.0;
887 double A = a1*c2*f3 - a1*c3*f2 - a2*c1*f3 + a2*c3*f1 +
889 double B = b1*c2*h3 - b1*c3*h2 - b2*c1*h3 + b2*c3*h1 +
891 double C = (a1*c2*h3 + b1*c2*f3) - (a1*c3*h2 + b1*c3*f2) -
892 (a2*c1*h3 + b2*c1*f3) + (a2*c3*h1 + b2*c3*f1) +
893 (a3*c1*h2 + b3*c1*f2) - (a3*c2*h1 + b3*c2*f1);
894 double D = a1*d2*f3 - a1*d3*f2 - a2*d1*f3 + a2*d3*f1 +
896 double E = b1*d2*h3 - b1*d3*h2 - b2*d1*h3 + b2*d3*h1 +
898 double F = (a1*d2*h3 + b1*d2*f3) - (a1*d3*h2 + b1*d3*f2) -
899 (a2*d1*h3 + b2*d1*f3) + (a2*d3*h1 + b2*d3*f1) +
900 (a3*d1*h2 + b3*d1*f2) - (a3*d2*h1 + b3*d2*f1);
901 double G = a1*e2*f3 - a1*e3*f2 - a2*e1*f3 + a2*e3*f1 +
903 double H = b1*e2*h3 - b1*e3*h2 - b2*e1*h3 + b2*e3*h1 +
905 double P = (a1*e2*h3 + b1*e2*f3) - (a1*e3*h2 + b1*e3*f2) -
906 (a2*e1*h3 + b2*e1*f3) + (a2*e3*h1 + b2*e3*f1) +
907 (a3*e1*h2 + b3*e1*f2) - (a3*e2*h1 + b3*e2*f1);
909 xvolume = -2.0*(2.0*(A + B + D + E + G + H) + C + F + P)/9.0;
911 //volume[index] = xvolume ;
912 volume->setIJ(index,1,xvolume) ;
917 case MED_HEXA8 : case MED_HEXA20 :
919 for (int hexa=0;hexa<nb_entity_type;hexa++)
921 int hexa_index = (type%100)*hexa;
923 int N1 = global_connectivity[hexa_index]-1;
924 int N2 = global_connectivity[hexa_index+1]-1;
925 int N3 = global_connectivity[hexa_index+2]-1;
926 int N4 = global_connectivity[hexa_index+3]-1;
927 int N5 = global_connectivity[hexa_index+4]-1;
928 int N6 = global_connectivity[hexa_index+5]-1;
929 int N7 = global_connectivity[hexa_index+6]-1;
930 int N8 = global_connectivity[hexa_index+7]-1;
932 double x1 = coord[dim_space*N1];
933 double x2 = coord[dim_space*N2];
934 double x3 = coord[dim_space*N3];
935 double x4 = coord[dim_space*N4];
936 double x5 = coord[dim_space*N5];
937 double x6 = coord[dim_space*N6];
938 double x7 = coord[dim_space*N7];
939 double x8 = coord[dim_space*N8];
941 double y1 = coord[dim_space*N1+1];
942 double y2 = coord[dim_space*N2+1];
943 double y3 = coord[dim_space*N3+1];
944 double y4 = coord[dim_space*N4+1];
945 double y5 = coord[dim_space*N5+1];
946 double y6 = coord[dim_space*N6+1];
947 double y7 = coord[dim_space*N7+1];
948 double y8 = coord[dim_space*N8+1];
950 double z1 = coord[dim_space*N1+2];
951 double z2 = coord[dim_space*N2+2];
952 double z3 = coord[dim_space*N3+2];
953 double z4 = coord[dim_space*N4+2];
954 double z5 = coord[dim_space*N5+2];
955 double z6 = coord[dim_space*N6+2];
956 double z7 = coord[dim_space*N7+2];
957 double z8 = coord[dim_space*N8+2];
959 double a1 = (x3-x4)/8.0, a2 = (y3-y4)/8.0, a3 = (z3-z4)/8.0;
960 double b1 = (x2-x1)/8.0, b2 = (y2-y1)/8.0, b3 = (z2-z1)/8.0;
961 double c1 = (x7-x8)/8.0, c2 = (y7-y8)/8.0, c3 = (z7-z8)/8.0;
962 double d1 = (x6-x5)/8.0, d2 = (y6-y5)/8.0, d3 = (z6-z5)/8.0;
963 double e1 = (x3-x2)/8.0, e2 = (y3-y2)/8.0, e3 = (z3-z2)/8.0;
964 double f1 = (x4-x1)/8.0, f2 = (y4-y1)/8.0, f3 = (z4-z1)/8.0;
965 double h1 = (x7-x6)/8.0, h2 = (y7-y6)/8.0, h3 = (z7-z6)/8.0;
966 double p1 = (x8-x5)/8.0, p2 = (y8-y5)/8.0, p3 = (z8-z5)/8.0;
967 double q1 = (x3-x7)/8.0, q2 = (y3-y7)/8.0, q3 = (z3-z7)/8.0;
968 double r1 = (x4-x8)/8.0, r2 = (y4-y8)/8.0, r3 = (z4-z8)/8.0;
969 double s1 = (x2-x6)/8.0, s2 = (y2-y6)/8.0, s3 = (z2-z6)/8.0;
970 double t1 = (x1-x5)/8.0, t2 = (y1-y5)/8.0, t3 = (z1-z5)/8.0;
972 double A = a1*e2*q3 - a1*e3*q2 - a2*e1*q3 + a2*e3*q1 +
974 double B = c1*h2*q3 - c1*h3*q2 - c2*h1*q3 + c2*h3*q1 +
976 double C = (a1*h2 + c1*e2)*q3 - (a1*h3 + c1*e3)*q2 -
977 (a2*h1 + c2*e1)*q3 + (a2*h3 + c2*e3)*q1 +
978 (a3*h1 + c3*e1)*q2 - (a3*h2 + c3*e2)*q1;
979 double D = b1*e2*s3 - b1*e3*s2 - b2*e1*s3 + b2*e3*s1 +
981 double E = d1*h2*s3 - d1*h3*s2 - d2*h1*s3 + d2*h3*s1 +
983 double F = (b1*h2 + d1*e2)*s3 - (b1*h3 + d1*e3)*s2 -
984 (b2*h1 + d2*e1)*s3 + (b2*h3 + d2*e3)*s1 +
985 (b3*h1 + d3*e1)*s2 - (b3*h2 + d3*e2)*s1;
986 double G = (a1*e2*s3 + b1*e2*q3) - (a1*e3*s2 + b1*e3*q2) -
987 (a2*e1*s3 + b2*e1*q3) + (a2*e3*s1 + b2*e3*q1) +
988 (a3*e1*s2 + b3*e1*q2) - (a3*e2*s1 + b3*e2*q1);
989 double H = (c1*h2*s3 + d1*h2*q3) - (c1*h3*s2 + d1*h3*q2) -
990 (c2*h1*s3 + d2*h1*q3) + (c2*h3*s1 + d2*h3*q1) +
991 (c3*h1*s2 + d3*h1*q2) - (c3*h2*s1 + d3*h2*q1);
992 double I = ((a1*h2 + c1*e2)*s3 + (b1*h2 + d1*e2)*q3) -
993 ((a1*h3 + c1*e3)*s2 + (b1*h3 + d1*e3)*q2) -
994 ((a2*h1 + c2*e1)*s3 + (b2*h1 + d2*e1)*q3) +
995 ((a2*h3 + c2*e3)*s1 + (b2*h3 + d2*e3)*q1) +
996 ((a3*h1 + c3*e1)*s2 + (b3*h1 + d3*e1)*q2) -
997 ((a3*h2 + c3*e2)*s1 + (b3*h2 + d3*e2)*q1);
998 double J = a1*f2*r3 - a1*f3*r2 - a2*f1*r3 + a2*f3*r1 +
1000 double K = c1*p2*r3 - c1*p3*r2 - c2*p1*r3 + c2*p3*r1 +
1001 c3*p1*r2 - c3*p2*r1;
1002 double L = (a1*p2 + c1*f2)*r3 - (a1*p3 + c1*f3)*r2 -
1003 (a2*p1 + c2*f1)*r3 + (a2*p3 + c2*f3)*r1 +
1004 (a3*p1 + c3*f1)*r2 - (a3*p2 + c3*f2)*r1;
1005 double M = b1*f2*t3 - b1*f3*t2 - b2*f1*t3 + b2*f3*t1 +
1006 b3*f1*t2 - b3*f2*t1;
1007 double N = d1*p2*t3 - d1*p3*t2 - d2*p1*t3 + d2*p3*t1 +
1008 d3*p1*t2 - d3*p2*t1;
1009 double O = (b1*p2 + d1*f2)*t3 - (b1*p3 + d1*f3)*t2 -
1010 (b2*p1 + d2*f1)*t3 + (b2*p3 + d2*f3)*t1 +
1011 (b3*p1 + d3*f1)*t2 - (b3*p2 + d3*f2)*t1;
1012 double P = (a1*f2*t3 + b1*f2*r3) - (a1*f3*t2 + b1*f3*r2) -
1013 (a2*f1*t3 + b2*f1*r3) + (a2*f3*t1 + b2*f3*r1) +
1014 (a3*f1*t2 + b3*f1*r2) - (a3*f2*t1 + b3*f2*r1);
1015 double Q = (c1*p2*t3 + d1*p2*r3) - (c1*p3*t2 + d1*p3*r2) -
1016 (c2*p1*t3 + d2*p1*r3) + (c2*p3*t1 + d2*p3*r1) +
1017 (c3*p1*t2 + d3*p1*r2) - (c3*p2*t1 + d3*p2*r1);
1018 double R = ((a1*p2 + c1*f2)*t3 + (b1*p2 + d1*f2)*r3) -
1019 ((a1*p3 + c1*f3)*t2 + (b1*p3 + d1*f3)*r2) -
1020 ((a2*p1 + c2*f1)*t3 + (b2*p1 + d2*f1)*r3) +
1021 ((a2*p3 + c2*f3)*t1 + (b2*p3 + d2*f3)*r1) +
1022 ((a3*p1 + c3*f1)*t2 + (b3*p1 + d3*f1)*r2) -
1023 ((a3*p2 + c3*f2)*t1 + (b3*p2 + d3*f2)*r1);
1024 double S = (a1*e2*r3 + a1*f2*q3) - (a1*e3*r2 + a1*f3*q2) -
1025 (a2*e1*r3 + a2*f1*q3) + (a2*e3*r1 + a2*f3*q1) +
1026 (a3*e1*r2 + a3*f1*q2) - (a3*e2*r1 + a3*f2*q1);
1027 double T = (c1*h2*r3 + c1*p2*q3) - (c1*h3*r2 + c1*p3*q2) -
1028 (c2*h1*r3 + c2*p1*q3) + (c2*h3*r1 + c2*p3*q1) +
1029 (c3*h1*r2 + c3*p1*q2) - (c3*h2*r1 + c3*p2*q1);
1030 double U = ((a1*h2 + c1*e2)*r3 + (a1*p2 + c1*f2)*q3) -
1031 ((a1*h3 + c1*e3)*r2 + (a1*p3 + c1*f3)*q2) -
1032 ((a2*h1 + c2*e1)*r3 + (a2*p1 + c2*f1)*q3) +
1033 ((a2*h3 + c2*e3)*r1 + (a2*p3 + c2*f3)*q1) +
1034 ((a3*h1 + c3*e1)*r2 + (a3*p1 + c3*f1)*q2) -
1035 ((a3*h2 + c3*e2)*r1 + (a3*p2 + c3*f2)*q1);
1036 double V = (b1*e2*t3 + b1*f2*s3) - (b1*e3*t2 + b1*f3*s2) -
1037 (b2*e1*t3 + b2*f1*s3) + (b2*e3*t1 + b2*f3*s1) +
1038 (b3*e1*t2 + b3*f1*s2) - (b3*e2*t1 + b3*f2*s1);
1039 double W = (d1*h2*t3 + d1*p2*s3) - (d1*h3*t2 + d1*p3*s2) -
1040 (d2*h1*t3 + d2*p1*s3) + (d2*h3*t1 + d2*p3*s1) +
1041 (d3*h1*t2 + d3*p1*s2) - (d3*h2*t1 + d3*p2*s1);
1042 double X = ((b1*h2 + d1*e2)*t3 + (b1*p2 + d1*f2)*s3) -
1043 ((b1*h3 + d1*e3)*t2 + (b1*p3 + d1*f3)*s2) -
1044 ((b2*h1 + d2*e1)*t3 + (b2*p1 + d2*f1)*s3) +
1045 ((b2*h3 + d2*e3)*t1 + (b2*p3 + d2*f3)*s1) +
1046 ((b3*h1 + d3*e1)*t2 + (b3*p1 + d3*f1)*s2) -
1047 ((b3*h2 + d3*e2)*t1 + (b3*p2 + d3*f2)*s1);
1048 double Y = (a1*e2*t3 + a1*f2*s3 + b1*e2*r3 + b1*f2*q3) -
1049 (a1*e3*t2 + a1*f3*s2 + b1*e3*r2 + b1*f3*q2) -
1050 (a2*e1*t3 + a2*f1*s3 + b2*e1*r3 + b2*f1*q3) +
1051 (a2*e3*t1 + a2*f3*s1 + b2*e3*r1 + b2*f3*q1) +
1052 (a3*e1*t2 + a3*f1*s2 + b3*e1*r2 + b3*f1*q2) -
1053 (a3*e2*t1 + a3*f2*s1 + b3*e2*r1 + b3*f2*q1);
1054 double Z = (c1*h2*t3 + c1*p2*s3 + d1*h2*r3 + d1*p2*q3) -
1055 (c1*h3*t2 + c1*p3*s2 + d1*h3*r2 + d1*p3*q2) -
1056 (c2*h1*t3 + c2*p1*s3 + d2*h1*r3 + d2*p1*q3) +
1057 (c2*h3*t1 + c2*p3*s1 + d2*h3*r1 + d2*p3*q1) +
1058 (c3*h1*t2 + c3*p1*s2 + d3*h1*r2 + d3*p1*q2) -
1059 (c3*h2*t1 + c3*p2*s1 + d3*h2*r1 + d3*p2*q1);
1060 double AA = ((a1*h2 + c1*e2)*t3 + (a1*p2 + c1*f2)*s3 +
1061 (b1*h2 + d1*e2)*r3 + (b1*p2 + d1*f2)*q3) -
1062 ((a1*h3 + c1*e3)*t2 + (a1*p3 + c1*f3)*s2 +
1063 (b1*h3 + d1*e3)*r2 + (b1*p3 + d1*f3)*q2) -
1064 ((a2*h1 + c2*e1)*t3 + (a2*p1 + c2*f1)*s3 +
1065 (b2*h1 + d2*e1)*r3 + (b2*p1 + d2*f1)*q3) +
1066 ((a2*h3 + c2*e3)*t1 + (a2*p3 + c2*f3)*s1 +
1067 (b2*h3 + d2*e3)*r1 + (b2*p3 + d2*f3)*q1) +
1068 ((a3*h1 + c3*e1)*t2 + (a3*p1 + c3*f1)*s2 +
1069 (b3*h1 + d3*e1)*r2 + (b3*p1 + d3*f1)*q2) -
1070 ((a3*h2 + c3*e2)*t1 + (a3*p2 + c3*f2)*s1 +
1071 (b3*h2 + d3*e2)*r1 + (b3*p2 + d3*f2)*q1);
1073 xvolume = 64.0*(8.0*(A + B + D + E + J + K + M + N) +
1074 4.0*(C + F + G + H + L + O + P + Q + S + T +
1075 V + W) + 2.0*(I + R + U + X + Y + Z) +
1078 //volume[index] = xvolume ;
1079 volume->setIJ(index,1,xvolume) ;
1085 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Volumes on it !"));
1093 FIELD<double>* MESH::getArea(const SUPPORT * Support) const throw (MEDEXCEPTION)
1095 const char * LOC = "MESH::getArea(SUPPORT*) : ";
1098 // Support must be on 2D elements
1100 // Make sure that the MESH class is the same as the MESH class attribut
1101 // in the class Support
1102 MESH* myMesh = Support->getMesh();
1104 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1106 int dim_space = getSpaceDimension();
1107 medEntityMesh support_entity = Support->getEntity();
1108 bool onAll = Support->isOnAllElements();
1110 int nb_type, length_values;
1111 const medGeometryElement* types;
1113 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1114 const int* global_connectivity;
1118 // nb_type = myMesh->getNumberOfTypes(support_entity);
1119 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1120 // types = getTypes(support_entity);
1124 nb_type = Support->getNumberOfTypes();
1125 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1126 types = Support->getTypes();
1130 FIELD<double>* Area;
1132 Area = new FIELD<double>::FIELD(Support,1);
1133 Area->setName("AREA");
1134 Area->setDescription("cells or faces area");
1136 Area->setComponentName(1,"area");
1137 Area->setComponentDescription(1,"desc-comp");
1139 /* string MEDComponentUnit = myMesh->getCoordinatesUnits()[0]+"*"+myMesh->getCoordinatesUnits()[1];*/
1141 string MEDComponentUnit="trucmuch";
1143 Area->setMEDComponentUnit(1,MEDComponentUnit);
1145 Area->setValueType(MED_REEL64);
1147 Area->setIterationNumber(0);
1148 Area->setOrderNumber(0);
1151 double *area = new double[length_values];
1152 //double *area = Area->getValue(MED_FULL_INTERLACE);
1154 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1157 for (int i=0;i<nb_type;i++)
1159 medGeometryElement type = types[i] ;
1164 nb_entity_type = getNumberOfElements(support_entity,type);
1165 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1169 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1171 nb_entity_type = Support->getNumberOfElements(type);
1173 const int * supp_number = Support->getNumber(type);
1174 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1175 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1176 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1178 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1179 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1180 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1184 global_connectivity = global_connectivity_tmp ;
1190 case MED_TRIA3 : case MED_TRIA6 :
1192 for (int tria=0;tria<nb_entity_type;tria++)
1194 int tria_index = (type%100)*tria;
1196 int N1 = global_connectivity[tria_index]-1;
1197 int N2 = global_connectivity[tria_index+1]-1;
1198 int N3 = global_connectivity[tria_index+2]-1;
1200 double x1 = coord[dim_space*N1];
1201 double x2 = coord[dim_space*N2];
1202 double x3 = coord[dim_space*N3];
1204 double y1 = coord[dim_space*N1+1];
1205 double y2 = coord[dim_space*N2+1];
1206 double y3 = coord[dim_space*N3+1];
1210 xarea = - ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1214 double z1 = coord[dim_space*N1+2];
1215 double z2 = coord[dim_space*N2+2];
1216 double z3 = coord[dim_space*N3+2];
1218 xarea = sqrt(((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))*
1219 ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)) +
1220 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))*
1221 ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1)) +
1222 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))*
1223 ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)))/2.0;
1226 area[index] = xarea ;
1231 case MED_QUAD4 : case MED_QUAD8 :
1233 for (int quad=0;quad<nb_entity_type;quad++)
1235 int quad_index = (type%100)*quad;
1237 int N1 = global_connectivity[quad_index]-1;
1238 int N2 = global_connectivity[quad_index+1]-1;
1239 int N3 = global_connectivity[quad_index+2]-1;
1240 int N4 = global_connectivity[quad_index+3]-1;
1242 double x1 = coord[dim_space*N1];
1243 double x2 = coord[dim_space*N2];
1244 double x3 = coord[dim_space*N3];
1245 double x4 = coord[dim_space*N4];
1247 double y1 = coord[dim_space*N1+1];
1248 double y2 = coord[dim_space*N2+1];
1249 double y3 = coord[dim_space*N3+1];
1250 double y4 = coord[dim_space*N4+1];
1254 double a1 = (x2-x1)/4.0, a2 = (y2-y1)/4.0;
1255 double b1 = (x3-x4)/4.0, b2 = (y3-y4)/4.0;
1256 double c1 = (x3-x2)/4.0, c2 = (y3-y2)/4.0;
1257 double d1 = (x4-x1)/4.0, d2 = (y4-y1)/4.0;
1259 xarea = - 4.0*(b1*c2 - c1*b2 + a1*c2 - c1*a2 + b1*d2 -
1260 d1*b2 + a1*d2 - d1*a2);
1264 double z1 = coord[dim_space*N1+2];
1265 double z2 = coord[dim_space*N2+2];
1266 double z3 = coord[dim_space*N3+2];
1267 double z4 = coord[dim_space*N4+2];
1269 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1270 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1271 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1272 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1273 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1274 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1275 sqrt(((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3))*
1276 ((y4-y3)*(z2-z3) - (y2-y3)*(z4-z3)) +
1277 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3))*
1278 ((x2-x3)*(z4-z3) - (x4-x3)*(z2-z3)) +
1279 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))*
1280 ((x4-x3)*(y2-y3) - (x2-x3)*(y4-y3))))/2.0;
1283 area[index] = xarea ;
1289 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Areas on it !"));
1294 Area->setValue(MED_FULL_INTERLACE,area);
1299 FIELD<double>* MESH::getLength(const SUPPORT * Support) const throw (MEDEXCEPTION)
1301 const char * LOC = "MESH::getLength(SUPPORT*) : ";
1304 // Support must be on 1D elements
1306 // Make sure that the MESH class is the same as the MESH class attribut
1307 // in the class Support
1308 MESH* myMesh = Support->getMesh();
1310 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1312 int dim_space = getSpaceDimension();
1313 medEntityMesh support_entity = Support->getEntity();
1314 bool onAll = Support->isOnAllElements();
1316 int nb_type, length_values;
1317 const medGeometryElement* types;
1319 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1320 const int* global_connectivity;
1324 // nb_type = myMesh->getNumberOfTypes(support_entity);
1325 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1326 // types = getTypes(support_entity);
1330 nb_type = Support->getNumberOfTypes();
1331 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1332 types = Support->getTypes();
1336 FIELD<double>* Length;
1338 Length = new FIELD<double>::FIELD(Support,1);
1339 // double *length = new double [length_values];
1340 Length->setValueType(MED_REEL64);
1342 //const double *length = Length->getValue(MED_FULL_INTERLACE);
1343 MEDARRAY<double> * length = Length->getvalue();
1345 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1348 for (int i=0;i<nb_type;i++)
1350 medGeometryElement type = types[i] ;
1355 nb_entity_type = getNumberOfElements(support_entity,type);
1356 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1360 //throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1362 nb_entity_type = Support->getNumberOfElements(type);
1364 const int * supp_number = Support->getNumber(type);
1365 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1366 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1367 int* global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1369 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1370 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1371 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1374 global_connectivity = global_connectivity_tmp ;
1379 case MED_SEG2 : case MED_SEG3 :
1381 for (int edge=0;edge<nb_entity_type;edge++)
1383 int edge_index = (type%100)*edge;
1385 int N1 = global_connectivity[edge_index]-1;
1386 int N2 = global_connectivity[edge_index+1]-1;
1388 double x1 = coord[dim_space*N1];
1389 double x2 = coord[dim_space*N2];
1391 double y1 = coord[dim_space*N1+1];
1392 double y2 = coord[dim_space*N2+1];
1394 double z1, z2 ; z1 = 0.0 ; z2 = 0.0 ; if (dim_space==3) {
1395 z1 = coord[dim_space*N1+2]; z2 = coord[dim_space*N2+2];}
1397 xlength = sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2) +
1398 (z1 - z2)*(z1 - z2));
1400 length->setIJ(index,1,xlength) ;
1406 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Lengths on it !"));
1414 FIELD<double>* MESH::getNormal(const SUPPORT * Support) const throw (MEDEXCEPTION)
1416 const char * LOC = "MESH::getNormal(SUPPORT*) : ";
1419 // Support must be on 2D or 1D elements
1421 // Make sure that the MESH class is the same as the MESH class attribut
1422 // in the class Support
1423 MESH* myMesh = Support->getMesh();
1425 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : pointeur problem !"));
1427 int dim_space = getSpaceDimension();
1428 medEntityMesh support_entity = Support->getEntity();
1429 bool onAll = Support->isOnAllElements();
1431 int nb_type, length_values;
1432 const medGeometryElement* types;
1434 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1435 const int* global_connectivity;
1439 // nb_type = myMesh->getNumberOfTypes(support_entity);
1440 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1441 // types = getTypes(support_entity);
1445 nb_type = Support->getNumberOfTypes();
1446 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1447 types = Support->getTypes();
1452 FIELD<double>* Normal = new FIELD<double>::FIELD(Support,dim_space);
1453 Normal->setName("NORMAL");
1454 Normal->setDescription("cells or faces normal");
1455 for (int k=1;k<=dim_space;k++) {
1456 Normal->setComponentName(k,"normal");
1457 Normal->setComponentDescription(k,"desc-comp");
1458 Normal->setMEDComponentUnit(k,"unit");
1461 Normal->setValueType(MED_REEL64);
1463 Normal->setIterationNumber(MED_NOPDT);
1464 Normal->setOrderNumber(MED_NONOR);
1465 Normal->setTime(0.0);
1467 double * normal = new double [dim_space*length_values];
1468 //double *normal = Normal->getValue(MED_FULL_INTERLACE);
1470 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1473 for (int i=0;i<nb_type;i++)
1475 medGeometryElement type = types[i] ;
1477 // Make sure we trying to get Normals on
1478 // TRIA3 or TRIA6 or QUAD4 QUAD8 (2D) cells on a 3D mesh
1479 // or on SEG2 or SEG3 (1D) cells on a 2D mesh
1481 if ( (((type==MED_TRIA3) || (type==MED_TRIA6) ||
1482 (type==MED_QUAD4) || (type==MED_QUAD8)) &&
1483 (dim_space != 3)) || (((type==MED_SEG2) || (type==MED_SEG3)) &&
1485 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh : dimension problem !"));
1487 double xnormal1, xnormal2, xnormal3 ;
1491 nb_entity_type = getNumberOfElements(support_entity,type);
1492 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1496 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all for instance !"));
1497 nb_entity_type = Support->getNumberOfElements(type);
1499 const int * supp_number = Support->getNumber(type);
1500 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1501 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1502 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1504 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1505 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1506 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1510 global_connectivity = global_connectivity_tmp ;
1516 case MED_TRIA3 : case MED_TRIA6 :
1518 for (int tria=0;tria<nb_entity_type;tria++)
1520 int tria_index = (type%100)*tria;
1522 int N1 = global_connectivity[tria_index]-1;
1523 int N2 = global_connectivity[tria_index+1]-1;
1524 int N3 = global_connectivity[tria_index+2]-1;
1527 double x1 = coord[dim_space*N1];
1528 double x2 = coord[dim_space*N2];
1529 double x3 = coord[dim_space*N3];
1531 double y1 = coord[dim_space*N1+1];
1532 double y2 = coord[dim_space*N2+1];
1533 double y3 = coord[dim_space*N3+1];
1535 double z1 = coord[dim_space*N1+2];
1536 double z2 = coord[dim_space*N2+2];
1537 double z3 = coord[dim_space*N3+2];
1539 xnormal1 = ((y2-y1)*(z3-z1) - (y3-y1)*(z2-z1))/2.0;
1540 xnormal2 = ((x3-x1)*(z2-z1) - (x2-x1)*(z3-z1))/2.0;
1541 xnormal3 = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2.0;
1543 normal[3*index] = xnormal1 ;
1544 normal[3*index+1] = xnormal2 ;
1545 normal[3*index+2] = xnormal3 ;
1551 case MED_QUAD4 : case MED_QUAD8 :
1553 for (int quad=0;quad<nb_entity_type;quad++)
1555 int quad_index = (type%100)*quad;
1557 int N1 = global_connectivity[quad_index]-1;
1558 int N2 = global_connectivity[quad_index+1]-1;
1559 int N3 = global_connectivity[quad_index+2]-1;
1560 int N4 = global_connectivity[quad_index+3]-1;
1563 double x1 = coord[dim_space*N1];
1564 double x2 = coord[dim_space*N2];
1565 double x3 = coord[dim_space*N3];
1566 double x4 = coord[dim_space*N4];
1568 double y1 = coord[dim_space*N1+1];
1569 double y2 = coord[dim_space*N2+1];
1570 double y3 = coord[dim_space*N3+1];
1571 double y4 = coord[dim_space*N4+1];
1573 double z1 = coord[dim_space*N1+2];
1574 double z2 = coord[dim_space*N2+2];
1575 double z3 = coord[dim_space*N3+2];
1576 double z4 = coord[dim_space*N4+2];
1578 xnormal1 = (y2-y1)*(z4-z1) - (y4-y1)*(z2-z1);
1579 xnormal2 = (x4-x1)*(z2-z1) - (x2-x1)*(z4-z1);
1580 xnormal3 = (x2-x1)*(y4-y1) - (x4-x1)*(y2-y1);
1581 xarea = sqrt(xnormal1*xnormal1 + xnormal2*xnormal2 +
1584 xnormal1 = xnormal1/xarea;
1585 xnormal2 = xnormal2/xarea;
1586 xnormal3 = xnormal3/xarea;
1588 xarea = (sqrt(((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1))*
1589 ((y2-y1)*(z4-z1) - (y4-y1)*(z2-z1)) +
1590 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1))*
1591 ((x4-x1)*(z2-z1) - (x2-x1)*(z4-z1)) +
1592 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))*
1593 ((x2-x1)*(y4-y1) - (x4-x1)*(y2-y1))) +
1594 sqrt(((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2))*
1595 ((y3-y2)*(z4-z3) - (y4-y3)*(z3-z2)) +
1596 ((x4-x3)*(z2-z2) - (x3-x2)*(z4-z3))*
1597 ((x4-x3)*(z4-z2) - (x3-x2)*(z4-z3)) +
1598 ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))*
1599 ((x3-x2)*(y4-y3) - (x4-x3)*(y3-y2))))/2.0;
1601 xnormal1 = xnormal1*xarea;
1602 xnormal2 = xnormal2*xarea;
1603 xnormal3 = xnormal3*xarea;
1605 normal[3*index] = xnormal1 ;
1606 normal[3*index+1] = xnormal2 ;
1607 normal[3*index+2] = xnormal3 ;
1613 case MED_SEG2 : case MED_SEG3 :
1615 for (int edge=0;edge<nb_entity_type;edge++)
1617 int edge_index = (type%100)*edge;
1619 int N1 = global_connectivity[edge_index]-1;
1620 int N2 = global_connectivity[edge_index+1]-1;
1622 double x1 = coord[dim_space*N1];
1623 double x2 = coord[dim_space*N2];
1625 double y1 = coord[dim_space*N1+1];
1626 double y2 = coord[dim_space*N2+1];
1628 xnormal1 = -(y2-y1);
1631 normal[2*index] = xnormal1 ;
1632 normal[2*index+1] = xnormal2 ;
1639 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get Normals on it !"));
1644 Normal->setValue(MED_FULL_INTERLACE,normal);
1652 FIELD<double>* MESH::getBarycenter(const SUPPORT * Support) const throw (MEDEXCEPTION)
1654 const char * LOC = "MESH::getBarycenter(SUPPORT*) : ";
1657 // Make sure that the MESH class is the same as the MESH class attribut
1658 // in the class Support
1659 MESH* myMesh = Support->getMesh();
1661 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"no compatibility between *this and SUPPORT::_mesh !"));
1663 int dim_space = getSpaceDimension();
1664 medEntityMesh support_entity = Support->getEntity();
1665 bool onAll = Support->isOnAllElements();
1667 int nb_type, length_values;
1668 const medGeometryElement* types;
1670 // !!!! WARNING : use of nodal global numbering in the mesh !!!!
1671 const int* global_connectivity;
1675 // nb_type = myMesh->getNumberOfTypes(support_entity);
1676 // length_values = getNumberOfElements(support_entity,MED_ALL_ELEMENTS);
1677 // types = getTypes(support_entity);
1681 nb_type = Support->getNumberOfTypes();
1682 length_values = Support->getNumberOfElements(MED_ALL_ELEMENTS);
1683 types = Support->getTypes();
1687 FIELD<double>* Barycenter;
1689 Barycenter = new FIELD<double>::FIELD(Support,dim_space);
1690 Barycenter->setName("BARYCENTER");
1691 Barycenter->setDescription("cells or faces barycenter");
1693 for (int k=0;k<dim_space;k++) {
1695 Barycenter->setComponentName(kp1,myMesh->getCoordinatesNames()[k]);
1696 Barycenter->setComponentDescription(kp1,"desc-comp");
1697 Barycenter->setMEDComponentUnit(kp1,myMesh->getCoordinatesUnits()[k]);
1700 Barycenter->setValueType(MED_REEL64);
1702 Barycenter->setIterationNumber(0);
1703 Barycenter->setOrderNumber(0);
1704 Barycenter->setTime(0.0);
1706 double *barycenter = new double [dim_space*length_values];
1707 // double *barycenter = Barycenter->getValue(MED_FULL_INTERLACE);
1709 const double * coord = getCoordinates(MED_FULL_INTERLACE);
1712 for (int i=0;i<nb_type;i++)
1714 medGeometryElement type = types[i] ;
1715 double xbarycenter1, xbarycenter2, xbarycenter3;
1719 nb_entity_type = getNumberOfElements(support_entity,type);
1720 global_connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,type);
1724 // throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on all !"));
1725 nb_entity_type = Support->getNumberOfElements(type);
1727 const int * supp_number = Support->getNumber(type);
1728 const int * connectivity = getConnectivity(MED_FULL_INTERLACE,MED_NODAL,support_entity,MED_ALL_ELEMENTS);
1729 const int * connectivityIndex = getConnectivityIndex(MED_NODAL,support_entity);
1730 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
1732 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
1733 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
1734 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
1737 global_connectivity = global_connectivity_tmp;
1742 case MED_TETRA4 : case MED_TETRA10 :
1744 for (int tetra =0;tetra<nb_entity_type;tetra++)
1746 int tetra_index = (type%100)*tetra;
1748 int N1 = global_connectivity[tetra_index]-1;
1749 int N2 = global_connectivity[tetra_index+1]-1;
1750 int N3 = global_connectivity[tetra_index+2]-1;
1751 int N4 = global_connectivity[tetra_index+3]-1;
1753 double x1 = coord[dim_space*N1];
1754 double x2 = coord[dim_space*N2];
1755 double x3 = coord[dim_space*N3];
1756 double x4 = coord[dim_space*N4];
1758 double y1 = coord[dim_space*N1+1];
1759 double y2 = coord[dim_space*N2+1];
1760 double y3 = coord[dim_space*N3+1];
1761 double y4 = coord[dim_space*N4+1];
1763 double z1 = coord[dim_space*N1+2];
1764 double z2 = coord[dim_space*N2+2];
1765 double z3 = coord[dim_space*N3+2];
1766 double z4 = coord[dim_space*N4+2];
1768 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1769 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1770 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
1771 barycenter[3*index] = xbarycenter1 ;
1772 barycenter[3*index+1] = xbarycenter2 ;
1773 barycenter[3*index+2] = xbarycenter3 ;
1778 case MED_PYRA5 : case MED_PYRA13 :
1780 for (int pyra=0;pyra<nb_entity_type;pyra++)
1782 int pyra_index = (type%100)*pyra;
1784 int N1 = global_connectivity[pyra_index]-1;
1785 int N2 = global_connectivity[pyra_index+1]-1;
1786 int N3 = global_connectivity[pyra_index+2]-1;
1787 int N4 = global_connectivity[pyra_index+3]-1;
1788 int N5 = global_connectivity[pyra_index+4]-1;
1790 double x1 = coord[dim_space*N1];
1791 double x2 = coord[dim_space*N2];
1792 double x3 = coord[dim_space*N3];
1793 double x4 = coord[dim_space*N4];
1794 double x5 = coord[dim_space*N5];
1796 double y1 = coord[dim_space*N1+1];
1797 double y2 = coord[dim_space*N2+1];
1798 double y3 = coord[dim_space*N3+1];
1799 double y4 = coord[dim_space*N4+1];
1800 double y5 = coord[dim_space*N5+1];
1802 double z1 = coord[dim_space*N1+2];
1803 double z2 = coord[dim_space*N2+2];
1804 double z3 = coord[dim_space*N3+2];
1805 double z4 = coord[dim_space*N4+2];
1806 double z5 = coord[dim_space*N5+2];
1808 xbarycenter1 = (x1 + x2 + x3 + x4 + x5)/5.0;
1809 xbarycenter2 = (y1 + y2 + y3 + y4 + y5)/5.0;
1810 xbarycenter3 = (z1 + z2 + z3 + z4 + z5)/5.0;
1811 barycenter[3*index] = xbarycenter1 ;
1812 barycenter[3*index+1] = xbarycenter2 ;
1813 barycenter[3*index+2] = xbarycenter3 ;
1818 case MED_PENTA6 : case MED_PENTA15 :
1820 for (int penta=0;penta<nb_entity_type;penta++)
1822 int penta_index = (type%100)*penta;
1824 int N1 = global_connectivity[penta_index]-1;
1825 int N2 = global_connectivity[penta_index+1]-1;
1826 int N3 = global_connectivity[penta_index+2]-1;
1827 int N4 = global_connectivity[penta_index+3]-1;
1828 int N5 = global_connectivity[penta_index+4]-1;
1829 int N6 = global_connectivity[penta_index+5]-1;
1831 double x1 = coord[dim_space*N1];
1832 double x2 = coord[dim_space*N2];
1833 double x3 = coord[dim_space*N3];
1834 double x4 = coord[dim_space*N4];
1835 double x5 = coord[dim_space*N5];
1836 double x6 = coord[dim_space*N6];
1838 double y1 = coord[dim_space*N1+1];
1839 double y2 = coord[dim_space*N2+1];
1840 double y3 = coord[dim_space*N3+1];
1841 double y4 = coord[dim_space*N4+1];
1842 double y5 = coord[dim_space*N5+1];
1843 double y6 = coord[dim_space*N6+1];
1845 double z1 = coord[dim_space*N1+2];
1846 double z2 = coord[dim_space*N2+2];
1847 double z3 = coord[dim_space*N3+2];
1848 double z4 = coord[dim_space*N4+2];
1849 double z5 = coord[dim_space*N5+2];
1850 double z6 = coord[dim_space*N6+2];
1852 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6)/6.0;
1853 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6)/6.0;
1854 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6)/6.0;
1855 barycenter[3*index] = xbarycenter1 ;
1856 barycenter[3*index+1] = xbarycenter2 ;
1857 barycenter[3*index+2] = xbarycenter3 ;
1862 case MED_HEXA8 : case MED_HEXA20 :
1864 for (int hexa=0;hexa<nb_entity_type;hexa++)
1866 int hexa_index = (type%100)*hexa;
1868 int N1 = global_connectivity[hexa_index]-1;
1869 int N2 = global_connectivity[hexa_index+1]-1;
1870 int N3 = global_connectivity[hexa_index+2]-1;
1871 int N4 = global_connectivity[hexa_index+3]-1;
1872 int N5 = global_connectivity[hexa_index+4]-1;
1873 int N6 = global_connectivity[hexa_index+5]-1;
1874 int N7 = global_connectivity[hexa_index+6]-1;
1875 int N8 = global_connectivity[hexa_index+7]-1;
1877 double x1 = coord[dim_space*N1];
1878 double x2 = coord[dim_space*N2];
1879 double x3 = coord[dim_space*N3];
1880 double x4 = coord[dim_space*N4];
1881 double x5 = coord[dim_space*N5];
1882 double x6 = coord[dim_space*N6];
1883 double x7 = coord[dim_space*N7];
1884 double x8 = coord[dim_space*N8];
1886 double y1 = coord[dim_space*N1+1];
1887 double y2 = coord[dim_space*N2+1];
1888 double y3 = coord[dim_space*N3+1];
1889 double y4 = coord[dim_space*N4+1];
1890 double y5 = coord[dim_space*N5+1];
1891 double y6 = coord[dim_space*N6+1];
1892 double y7 = coord[dim_space*N7+1];
1893 double y8 = coord[dim_space*N8+1];
1895 double z1 = coord[dim_space*N1+2];
1896 double z2 = coord[dim_space*N2+2];
1897 double z3 = coord[dim_space*N3+2];
1898 double z4 = coord[dim_space*N4+2];
1899 double z5 = coord[dim_space*N5+2];
1900 double z6 = coord[dim_space*N6+2];
1901 double z7 = coord[dim_space*N7+2];
1902 double z8 = coord[dim_space*N8+2];
1904 xbarycenter1 = (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)/8.0;
1905 xbarycenter2 = (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8)/8.0;
1906 xbarycenter3 = (z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8)/8.0;
1908 barycenter[3*index] = xbarycenter1 ;
1909 barycenter[3*index+1] = xbarycenter2 ;
1910 barycenter[3*index+2] = xbarycenter3 ;
1916 case MED_TRIA3 : case MED_TRIA6 :
1918 for (int tria=0;tria<nb_entity_type;tria++)
1920 int tria_index = (type%100)*tria;
1922 int N1 = global_connectivity[tria_index]-1;
1923 int N2 = global_connectivity[tria_index+1]-1;
1924 int N3 = global_connectivity[tria_index+2]-1;
1926 double x1 = coord[dim_space*N1];
1927 double x2 = coord[dim_space*N2];
1928 double x3 = coord[dim_space*N3];
1930 double y1 = coord[dim_space*N1+1];
1931 double y2 = coord[dim_space*N2+1];
1932 double y3 = coord[dim_space*N3+1];
1934 xbarycenter1 = (x1 + x2 + x3)/3.0;
1935 xbarycenter2 = (y1 + y2 + y3)/3.0;
1939 barycenter[2*index] = xbarycenter1 ;
1940 barycenter[2*index+1] = xbarycenter2 ;
1945 coord[dim_space*N1+2];
1947 coord[dim_space*N2+2];
1949 coord[dim_space*N3+2];
1951 xbarycenter3 = (z1 + z2 + z3)/3.0;
1953 barycenter[3*index] = xbarycenter1 ;
1954 barycenter[3*index+1] = xbarycenter2 ;
1955 barycenter[3*index+2] = xbarycenter3 ;
1962 case MED_QUAD4 : case MED_QUAD8 :
1964 for (int quad=0;quad<nb_entity_type;quad++)
1966 int quad_index = (type%100)*quad;
1968 int N1 = global_connectivity[quad_index]-1;
1969 int N2 = global_connectivity[quad_index+1]-1;
1970 int N3 = global_connectivity[quad_index+2]-1;
1971 int N4 = global_connectivity[quad_index+3]-1;
1973 double x1 = coord[dim_space*N1];
1974 double x2 = coord[dim_space*N2];
1975 double x3 = coord[dim_space*N3];
1976 double x4 = coord[dim_space*N4];
1978 double y1 = coord[dim_space*N1+1];
1979 double y2 = coord[dim_space*N2+1];
1980 double y3 = coord[dim_space*N3+1];
1981 double y4 = coord[dim_space*N4+1];
1983 xbarycenter1 = (x1 + x2 + x3 + x4)/4.0;
1984 xbarycenter2 = (y1 + y2 + y3 + y4)/4.0;
1988 barycenter[2*index] = xbarycenter1 ;
1989 barycenter[2*index+1] = xbarycenter2 ;
1993 double z1 = coord[dim_space*N1+2];
1994 double z2 = coord[dim_space*N2+2];
1995 double z3 = coord[dim_space*N3+2];
1996 double z4 = coord[dim_space*N4+2];
1998 xbarycenter3 = (z1 + z2 + z3 + z4)/4.0;
2000 barycenter[3*index] = xbarycenter1 ;
2001 barycenter[3*index+1] = xbarycenter2 ;
2002 barycenter[3*index+2] = xbarycenter3 ;
2008 case MED_SEG2 : case MED_SEG3 :
2010 for (int edge=0;edge<nb_entity_type;edge++)
2012 int edge_index = (type%100)*edge;
2014 int N1 = global_connectivity[edge_index]-1;
2015 int N2 = global_connectivity[edge_index+1]-1;
2017 double x1 = coord[dim_space*N1];
2018 double x2 = coord[dim_space*N2];
2020 double y1 = coord[dim_space*N1+1];
2021 double y2 = coord[dim_space*N2+1];
2023 xbarycenter1 = (x1 + x2)/2.0;
2024 xbarycenter2 = (y1 + y2)/2.0;
2028 barycenter[2*index] = xbarycenter1 ;
2029 barycenter[2*index+1] = xbarycenter2 ;
2033 double z1 = coord[dim_space*N1+2];
2034 double z2 = coord[dim_space*N2+2];
2036 xbarycenter3 = (z1 + z2)/2.0;
2038 barycenter[3*index] = xbarycenter1 ;
2039 barycenter[3*index+1] = xbarycenter2 ;
2040 barycenter[3*index+2] = xbarycenter3 ;
2047 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad Support to get a barycenter on it (in fact unknown type) !"));
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();
2084 void MESH::read(int index)
2086 const char * LOC = "MESH::read(int index=0) : ";
2089 if (_drivers[index]) {
2090 _drivers[index]->open();
2091 _drivers[index]->read();
2092 _drivers[index]->close();
2095 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2096 << "The index given is invalid, index must be between 0 and |"
2101 ((GRID *) this)->fillMeshAfterRead();
2105 //=======================================================================
2106 //function : getSkin
2108 //=======================================================================
2110 SUPPORT * MESH::getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION)
2112 const char * LOC = "MESH::getSkin : " ;
2115 if (this != Support3D->getMesh())
2116 throw MEDEXCEPTION(STRING(LOC) << "no compatibility between *this and SUPPORT::_mesh !");
2117 if (_meshDimension != 3 || Support3D->getEntity() != MED_CELL)
2118 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Defined on 3D cells only"));
2120 // well, all rigth !
2121 SUPPORT * mySupport = new SUPPORT(this,"Skin",MED_FACE);
2122 mySupport->setAll(false);
2124 list<int> myElementsList ;
2127 calculateConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL);
2128 if (Support3D->isOnAllElements())
2130 int * myConnectivityValue = const_cast <int*> (getReverseConnectivity(MED_DESCENDING)) ;
2131 int nbFaces = getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS);
2132 for (i=0, j=1 ; j<=nbFaces; ++j, i += 2)
2134 int cellNb1 = myConnectivityValue [i];
2135 int cellNb2 = myConnectivityValue [i+1];
2136 //MESSAGE( " FACE # " << j << " -- Cells: " << cellNb1 << ", " << cellNb2 );
2137 if ((cellNb1 == 0 || cellNb2 == 0) && (cellNb1 + cellNb2 > 0))
2139 myElementsList.push_back( j ) ;
2146 map<int,int> FaceNbEncounterNb;
2148 int * myConnectivityValue = const_cast <int *>
2149 (getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING,
2150 MED_CELL, MED_ALL_ELEMENTS));
2151 int * myConnectivityIndex = const_cast <int *> (getConnectivityIndex(MED_DESCENDING, MED_CELL));
2152 int * myCellNbs = const_cast <int *> (Support3D->getnumber()->getValue());
2153 int nbCells = Support3D->getnumber()->getLength();
2154 for (i=0; i<nbCells; ++i)
2156 int cellNb = myCellNbs [ i ];
2157 int faceFirst = myConnectivityIndex[ cellNb-1 ];
2158 int faceLast = myConnectivityIndex[ cellNb ];
2159 for (j = faceFirst; j < faceLast; ++j)
2161 int faceNb = abs( myConnectivityValue [ j-1 ] );
2162 //MESSAGE( "Cell # " << i << " -- Face: " << faceNb);
2163 if (FaceNbEncounterNb.find( faceNb ) == FaceNbEncounterNb.end())
2164 FaceNbEncounterNb[ faceNb ] = 1;
2166 FaceNbEncounterNb[ faceNb ] ++;
2169 map<int,int>::iterator FaceNbEncounterNbItr;
2170 for (FaceNbEncounterNbItr = FaceNbEncounterNb.begin();
2171 FaceNbEncounterNbItr != FaceNbEncounterNb.end();
2172 FaceNbEncounterNbItr ++)
2173 if ((*FaceNbEncounterNbItr).second == 1)
2175 myElementsList.push_back( (*FaceNbEncounterNbItr).first) ;
2179 // Well, we must know how many geometric type we have found
2180 int * myListArray = new int[size] ;
2182 list<int>::iterator myElementsListIt ;
2183 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2184 myListArray[id]=(*myElementsListIt) ;
2188 int numberOfGeometricType ;
2189 medGeometryElement* geometricType ;
2190 int * numberOfGaussPoint ;
2191 int * geometricTypeNumber ;
2192 int * numberOfEntities ;
2193 // MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY() ;
2194 int * mySkyLineArrayIndex ;
2196 int numberOfType = getNumberOfTypes(MED_FACE) ;
2197 if (numberOfType == 1) { // wonderfull : it's easy !
2198 numberOfGeometricType = 1 ;
2199 geometricType = new medGeometryElement[1] ;
2200 const medGeometryElement * allType = getTypes(MED_FACE);
2201 geometricType[0] = allType[0] ;
2202 numberOfGaussPoint = new int[1] ;
2203 numberOfGaussPoint[0] = 1 ;
2204 geometricTypeNumber = new int[1] ; // not use, but initialized to nothing
2205 geometricTypeNumber[0] = 0 ;
2206 numberOfEntities = new int[1] ;
2207 numberOfEntities[0] = size ;
2208 mySkyLineArrayIndex = new int[2] ;
2209 mySkyLineArrayIndex[0]=1 ;
2210 mySkyLineArrayIndex[1]=1+size ;
2213 map<medGeometryElement,int> theType ;
2214 for (myElementsListIt=myElementsList.begin();myElementsListIt!=myElementsList.end();myElementsListIt++) {
2215 medGeometryElement myType = getElementType(MED_FACE,*myElementsListIt) ;
2216 if (theType.find(myType) != theType.end() )
2217 theType[myType]+=1 ;
2221 numberOfGeometricType = theType.size() ;
2222 geometricType = new medGeometryElement[numberOfGeometricType] ;
2223 const medGeometryElement * allType = getTypes(MED_FACE);
2224 numberOfGaussPoint = new int[numberOfGeometricType] ;
2225 geometricTypeNumber = new int[numberOfGeometricType] ; // not use, but initialized to nothing
2226 numberOfEntities = new int[numberOfGeometricType] ;
2227 mySkyLineArrayIndex = new int[numberOfGeometricType+1] ;
2229 mySkyLineArrayIndex[0]=1 ;
2230 map<medGeometryElement,int>::iterator theTypeIt ;
2231 for (theTypeIt=theType.begin();theTypeIt!=theType.end();theTypeIt++) {
2232 geometricType[index] = (*theTypeIt).first ;
2233 numberOfGaussPoint[index] = 1 ;
2234 geometricTypeNumber[index] = 0 ;
2235 numberOfEntities[index] = (*theTypeIt).second ;
2236 mySkyLineArrayIndex[index+1]=mySkyLineArrayIndex[index]+numberOfEntities[index] ;
2240 // mySkyLineArray->setMEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2241 MEDSKYLINEARRAY * mySkyLineArray = new MEDSKYLINEARRAY(numberOfGeometricType,size,mySkyLineArrayIndex,myListArray) ;
2243 mySupport->setNumberOfGeometricType(numberOfGeometricType) ;
2244 mySupport->setGeometricType(geometricType) ;
2245 mySupport->setNumberOfGaussPoint(numberOfGaussPoint) ;
2246 // mySupport->setGeometricTypeNumber(geometricTypeNumber) ;
2247 mySupport->setNumberOfElements(numberOfEntities) ;
2248 mySupport->setTotalNumberOfElements(size) ;
2249 mySupport->setNumber(mySkyLineArray) ;
2251 delete[] numberOfEntities;
2252 delete[] geometricTypeNumber;
2253 delete[] numberOfGaussPoint;
2254 delete[] geometricType;
2255 delete[] mySkyLineArrayIndex;
2256 delete[] myListArray;
2257 // delete mySkyLineArray;