1 // File : MEDMEM_Grid.hxx
2 // Created : Wed Dec 18 08:35:26 2002
3 // Descr : class containing structured mesh data
5 // Author : Edward AGAPOV (eap)
6 // Project : SALOME Pro
8 // Copyright : Open CASCADE
12 #include "MEDMEM_Grid.hxx"
13 #include "MEDMEM_CellModel.hxx"
14 #include "MEDMEM_SkyLineArray.hxx"
17 using namespace MEDMEM;
18 using namespace MED_EN;
20 //=======================================================================
22 //purpose : empty constructor
23 //=======================================================================
27 MESSAGE("A GRID CREATED");
30 //=======================================================================
32 ////purpose : array constructor
33 ////=======================================================================
34 GRID::GRID(const std::vector<std::vector<double> >& xyz_array,const std::vector<std::string>& coord_name,
35 const std::vector<std::string>& coord_unit, const med_grid_type type) : _gridType(type)
37 _spaceDimension = xyz_array.size();
39 // compute & set _numberOfNodes
41 for(int i=0;i!=xyz_array.size();++i)
42 NumberOfNodes*=xyz_array[i].size();
43 _numberOfNodes = NumberOfNodes ;
45 // create a non allocated COORDINATE
46 _coordinate = new COORDINATE(_spaceDimension, &coord_name[0], &coord_unit[0]);
47 string coordinateSystem = "UNDEFINED";
48 if( _gridType == MED_CARTESIAN)
49 coordinateSystem = "CARTESIAN";
50 else if ( _gridType == MED_POLAR)
51 coordinateSystem = "CYLINDRICAL";
52 _coordinate->setCoordinatesSystem(coordinateSystem);
55 if (_spaceDimension>=1)
57 _iArrayLength=xyz_array[0].size();
58 _iArray=new double[_iArrayLength];
59 std::copy(xyz_array[0].begin(),xyz_array[0].end(),_iArray);
61 if (_spaceDimension>=2)
63 _jArrayLength=xyz_array[1].size();
64 _jArray=new double[_jArrayLength];
65 std::copy(xyz_array[1].begin(),xyz_array[1].end(),_jArray);
67 if (_spaceDimension>=3)
69 _kArrayLength=xyz_array[2].size();
70 _kArray=new double[_kArrayLength];
71 std::copy(xyz_array[2].begin(),xyz_array[2].end(),_kArray);
74 _is_coordinates_filled = false;
75 _is_connectivity_filled = false;
79 //=======================================================================
81 //purpose : empty constructor
82 //=======================================================================
84 GRID::GRID(const med_grid_type type)
88 MESSAGE("A TYPED GRID CREATED");
91 //=======================================================================
93 //purpose : copying constructor
94 //=======================================================================
96 GRID::GRID(const GRID& otherGrid) {
100 //=======================================================================
102 //purpose : destructor
103 //=======================================================================
106 MESSAGE("GRID::~GRID() : Destroying the Grid");
107 if ( _iArray != (double* ) NULL) delete [] _iArray;
108 if ( _jArray != (double* ) NULL) delete [] _jArray;
109 if ( _kArray != (double* ) NULL) delete [] _kArray;
112 //=======================================================================
115 //=======================================================================
119 _gridType = MED_CARTESIAN;
121 _iArray = _jArray = _kArray = (double* ) NULL;
122 _iArrayLength = _jArrayLength = _kArrayLength = 0;
124 _is_coordinates_filled = false;
125 _is_connectivity_filled = false;
132 //=======================================================================
133 //function: operator=
134 //purpose : operator=
135 //=======================================================================
137 GRID & GRID::operator=(const GRID & otherGrid)
141 MESH::operator=(otherGrid);
145 //=======================================================================
147 //purpose : Create a GRID object using a MESH driver of type
148 // (MED_DRIVER, ....) associated with file <fileName>.
149 //=======================================================================
151 GRID::GRID(driverTypes driverType, const string & fileName,
152 const string & driverName) : MESH(driverType, fileName, driverName)
154 const char * LOC ="GRID::GRID(driverTypes , const string & , const string &) : ";
161 // MESH(driverType,fileName,driverName);
163 // current = addDriver(driverType,fileName,driverName);
165 // switch(_drivers[current]->getAccessMode() ) {
167 // MESSAGE(LOC << "driverType must have a MED_RDWR accessMode");
168 // rmDriver(current);
172 // _drivers[current]->open();
173 // _drivers[current]->read();
174 // _drivers[current]->close();
176 // // fill some fields of MESH
177 // fillMeshAfterRead();
185 return the GRID Geometric type, without computing all connectivity
187 const medGeometryElement * GRID::getTypes(medEntityMesh entity) const
189 static const medGeometryElement _gridGeometry[4]={MED_HEXA8,MED_QUAD4,MED_SEG2,MED_POINT1};
195 else if(entity==MED_FACE && _spaceDimension>2 )
197 else if(entity==MED_EDGE && _spaceDimension>1 )
199 else if(entity==MED_NODE && _spaceDimension>0)
202 throw MEDEXCEPTION(LOCALIZED("CONNECTIVITY::getGeometricTypes : Entity not defined !"));
203 return &_gridGeometry[i];
206 MED_EN::medGeometryElement * GRID::getTypesWithPoly(MED_EN::medEntityMesh Entity) const
208 int size=getNumberOfTypesWithPoly(Entity);
209 MED_EN::medGeometryElement *ret=new MED_EN::medGeometryElement[size];
210 memcpy(ret,getTypes(Entity),size*sizeof(MED_EN::medGeometryElement));
214 //=======================================================================
215 //function : fillMeshAfterRead
217 //=======================================================================
219 void GRID::fillMeshAfterRead()
221 // fill not only MESH (:-)
222 _is_coordinates_filled = false;
223 _is_connectivity_filled = false;
225 if (_gridType == MED_BODY_FITTED)
227 _is_coordinates_filled = true;
229 // nb of nodes in each direction is not known, set anything
230 // in order to be able to work anyhow
231 // INFOS("GRID::fillMeshAfterRead(): This stub must be removed");
232 // switch (_spaceDimension) {
234 // _iArrayLength = _numberOfNodes;
235 // _jArrayLength = 0;
236 // _kArrayLength = 0;
239 // _iArrayLength = 2;
240 // _jArrayLength = _numberOfNodes / _iArrayLength;
241 // _kArrayLength = 0;
244 // _iArrayLength = 2;
245 // _jArrayLength = _numberOfNodes / _iArrayLength / 2;
246 // _kArrayLength = _numberOfNodes / _iArrayLength / _jArrayLength;
248 //cout << "ARRAY LENGTHS: " << _iArrayLength << " " << _jArrayLength
249 // << " " << _kArrayLength << endl;
252 // int NbNodes = _iArrayLength;
253 // if (_jArrayLength)
254 // NbNodes *= _jArrayLength;
255 // if (_kArrayLength)
256 // NbNodes *= _kArrayLength;
257 // MESH::_numberOfNodes = NbNodes;
260 MESH::_meshDimension = MESH::_spaceDimension;
263 //=======================================================================
264 //function : fillCoordinates
266 //=======================================================================
268 void GRID::fillCoordinates() const
270 if (_is_coordinates_filled)
273 const char * LOC ="GRID::fillCoordinates()";
276 // if coordonate has not been allocated, perform shalow copy, transfer ownership of matrix
277 if(_coordinate->getSpaceDimension()*_coordinate->getNumberOfNodes() == 0)
278 _coordinate->setCoordinates(new MEDARRAY<double>(_spaceDimension,_numberOfNodes,MED_FULL_INTERLACE),true);
280 double* myCoord = const_cast <double *> ( _coordinate->getCoordinates(MED_FULL_INTERLACE) );
282 bool hasJ = _jArrayLength, hasK = _kArrayLength;
283 int J = hasJ ? _jArrayLength : 1;
284 int K = hasK ? _kArrayLength : 1;
285 //int nb, !! UNUSED VARIABLE !!
287 for (k=0; k < K; ++k) {
288 for (j=0; j < J; ++j) {
289 for (i=0; i < _iArrayLength; ++i) {
291 * myCoord = _iArray[ i ];
296 * myCoord = _jArray[ j ];
301 * myCoord = _kArray[ k ];
309 (const_cast <GRID *> (this))->_is_coordinates_filled = true;
313 //=======================================================================
314 //function : makeConnectivity
316 //=======================================================================
318 CONNECTIVITY * GRID::makeConnectivity (medEntityMesh Entity,
319 medGeometryElement Geometry,
323 const int * NodeNumbers)
326 CONNECTIVITY * Connectivity = new CONNECTIVITY(Entity) ;
327 Connectivity->_numberOfNodes = nbMeshNodes ;
328 Connectivity->_entityDimension = Geometry/100 ;
330 int numberOfGeometricType = 1;
331 Connectivity->_numberOfTypes = numberOfGeometricType;
333 Connectivity->_count = new int [numberOfGeometricType + 1] ;
334 Connectivity->_count[0] = 1;
335 Connectivity->_count[1] = 1 + NbEntities;
337 Connectivity->_type = new CELLMODEL [numberOfGeometricType];
338 Connectivity->_type[0] = CELLMODEL( Geometry ) ;
340 Connectivity->_geometricTypes = new medGeometryElement [ numberOfGeometricType ];
341 Connectivity->_geometricTypes[0] = Geometry;
343 // Connectivity->_nodal = new MEDSKYLINEARRAY() ;
344 int * skyLineArrayIndex = new int [NbEntities + 1];
345 int i, j, nbEntityNodes = Connectivity->_type[0].getNumberOfNodes();
346 for (i=0, j=1; i <= NbEntities; ++i, j += nbEntityNodes)
347 skyLineArrayIndex [ i ] = j;
349 // Connectivity->_nodal->setMEDSKYLINEARRAY (NbEntities, NbNodes,
350 //skyLineArrayIndex, NodeNumbers);
352 Connectivity->_nodal = new MEDSKYLINEARRAY (NbEntities, NbNodes,
353 skyLineArrayIndex, NodeNumbers);
355 delete [] skyLineArrayIndex;
357 // test connectivity right away
358 // for (med_int j = 0; j < numberOfGeometricType; j++)
360 // int node_number = Connectivity->_type[j].getNumberOfNodes();
361 // for (med_int k = Connectivity->_count[j]; k < Connectivity->_count[j+1]; k++)
362 // for (med_int local_node_number = 1 ; local_node_number < node_number+1; local_node_number++)
364 // cout << "MEDSKYLINEARRAY::getIJ(" << k << ", " << local_node_number << endl;
365 // med_int global_node = Connectivity->_nodal->getIJ(k,local_node_number) ;
366 // cout << "...= " << global_node << endl;
373 //=======================================================================
374 //function : fillConnectivity
375 //purpose : fill _coordinates and _connectivity of MESH if not yet done
376 //=======================================================================
378 void GRID::fillConnectivity() const
380 if (_is_connectivity_filled)
382 MESSAGE("GRID::fillConnectivity(): Already filled");
386 const char * LOC = "GRID::fillConnectivity() ";
389 int nbCells, nbFaces, nbEdges;
390 int nbCNodes, nbFNodes, nbENodes, nbMeshNodes;
391 int indexC, indexF, indexE;
392 int * nodeCNumbers, * nodeFNumbers, * nodeENumbers;
393 // about descending connectivity
394 int nbSub, nbRevSub, indexSub, indexRevSub, * subNumbers, * subRevNumbers;
396 bool hasFaces = _kArrayLength, hasEdges = _jArrayLength;
398 int iLenMin1 = _iArrayLength-1, jLenMin1 = _jArrayLength-1;
399 int kLenMin1 = _kArrayLength-1, ijLenMin1 = iLenMin1 * jLenMin1;
400 if (!hasEdges) jLenMin1 = 1;
401 if (!hasFaces) kLenMin1 = 1;
403 // nb of cells and of their connectivity nodes
405 nbCells = iLenMin1 * jLenMin1 * kLenMin1;
406 nbMeshNodes = _iArrayLength * (_jArrayLength ? _jArrayLength : 1) * (_kArrayLength ? _kArrayLength : 1);
407 nbCNodes = nbCells * 2 * (hasEdges ? 2 : 1) * (hasFaces ? 2 : 1);
408 nodeCNumbers = new int [ nbCNodes ];
410 // nb of faces and of their connectivity nodes
413 nbFaces = _iArrayLength * jLenMin1 * kLenMin1;
414 nbFaces += _jArrayLength * kLenMin1 * iLenMin1;
415 nbFaces += _kArrayLength * iLenMin1 * jLenMin1;
416 nbFNodes = nbFaces * 4;
417 nodeFNumbers = new int [ nbFNodes ];
419 nbFaces = nbFNodes = 0;
421 // nb of edges and of their connectivity nodes
424 if (_kArrayLength) { // 3d grid
425 nbEdges = iLenMin1 * _jArrayLength * _kArrayLength;
426 nbEdges += jLenMin1 * _kArrayLength * _iArrayLength;
427 nbEdges += kLenMin1 * _iArrayLength * _jArrayLength;
429 else if (_jArrayLength) { // 2d
430 nbEdges = iLenMin1 * _jArrayLength;
431 nbEdges += jLenMin1 * _iArrayLength;
433 nbENodes = nbEdges * 2;
434 nodeENumbers = new int [ nbENodes ];
436 nbEdges = nbENodes = 0;
438 // nb of descending connectivity Elements
443 nbRevSub = nbFaces * 2;
448 nbRevSub = nbEdges * 2;
450 subNumbers = new int [ nbSub ];
451 subRevNumbers = new int [ nbRevSub ];
452 for (int r=0; r<nbRevSub; ++r)
453 subRevNumbers[ r ] = 0;
455 int nSubI = 1, nSubJ, nSubK; // subelement numbers
458 nSubJ = getFaceNumber(2, 0, 0, 0);
459 nSubK = getFaceNumber(3, 0, 0, 0);
462 nSubJ = getEdgeNumber(2, 0, 0, 0);
465 // fill cell node numbers and descending element numbers
468 indexC = indexF = indexE = indexSub = indexRevSub = -1;
469 int iNode = 0, iCell = 0;
470 int ijLen = _iArrayLength * _jArrayLength;
471 int i, j, k, n1, n2, n3 ,n4;
473 int I = _iArrayLength;
474 int J = _jArrayLength ? _jArrayLength : 2; // pass at least once
475 int K = _kArrayLength ? _kArrayLength : 2;
476 // pass by all but last nodes in all directions
477 for (k = 1; k < K; ++k ) {
478 for (j = 1; j < J; ++j ) {
479 for (i = 1; i < I; ++i ) {
483 n1 = ++iNode; // iNode
485 nodeCNumbers [ ++indexC ] = n1;
486 nodeCNumbers [ ++indexC ] = n2;
488 if (hasEdges) { // at least 2d
491 nodeCNumbers [ ++indexC ] = n3;
492 nodeCNumbers [ ++indexC ] = n4;
494 if (hasFaces) { // 3d
495 nodeCNumbers [ ++indexC ] = n1 + ijLen;
496 nodeCNumbers [ ++indexC ] = n2 + ijLen;
497 nodeCNumbers [ ++indexC ] = n3 + ijLen;
498 nodeCNumbers [ ++indexC ] = n4 + ijLen;
505 subNumbers [ ++indexSub ] = n1;
506 subRevNumbers[ n3 ] = iCell;
507 subNumbers [ ++indexSub ] = n2;
508 subRevNumbers[ n4 ] = -iCell;
513 subNumbers [ ++indexSub ] = n1;
514 subRevNumbers[ n3 ] = iCell;
515 subNumbers [ ++indexSub ] = n2;
516 subRevNumbers[ n4 ] = -iCell;
521 subNumbers [ ++indexSub ] = n1;
522 subRevNumbers[ n3 ] = iCell;
523 subNumbers [ ++indexSub ] = n2;
524 subRevNumbers[ n4 ] = -iCell;
533 subNumbers [ ++indexSub ] = n1;
534 subRevNumbers[ n3 ] = iCell;
535 subNumbers [ ++indexSub ] = n2;
536 subRevNumbers[ n4 ] = -iCell;
541 subNumbers [ ++indexSub ] = n1;
542 subRevNumbers[ n3 ] = iCell;
543 subNumbers [ ++indexSub ] = n2;
544 subRevNumbers[ n4 ] = -iCell;
546 ++nSubI; ++nSubJ; ++nSubK;
548 ++iNode; // skip the last node in a row
554 iNode += I; // skip the whole last row
557 // fill face node numbers
559 int ax, AX = hasFaces ? 3 : 0;
560 for (ax = 1; ax <= AX; ++ax) {
568 case 1: --J; --K; break;
569 case 2: --K; --I; break;
572 for (k = 1; k <= K; ++k ) {
573 for (j = 1; j <= J; ++j ) {
574 for (i = 1; i <= I; ++i ) {
579 case 1: // nodes for faces normal to i direction
584 case 2: // nodes for faces normal to j direction
589 default: // nodes for faces normal to k direction
594 nodeFNumbers [ ++indexF ] = n1;
595 nodeFNumbers [ ++indexF ] = n2;
596 nodeFNumbers [ ++indexF ] = n3;
597 nodeFNumbers [ ++indexF ] = n4;
599 if (ax != 1) ++iNode;
601 if (ax != 2) iNode += _iArrayLength;
604 if (nbFNodes != indexF+1) {
605 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Wrong nbFNodes : " \
606 << nbFNodes << " indexF : " << indexF ));
609 // fill edge node numbers
611 AX = hasEdges ? _spaceDimension : 0;
612 for (ax = 1; ax <= AX; ++ax) {
626 for (k = 1; k <= K; ++k ) {
627 for (j = 1; j <= J; ++j ) {
628 for (i = 1; i <= I; ++i ) {
633 case 1: // nodes for edges going along i direction
636 case 2: // nodes for edges going along j direction
637 n2 = n1 + _iArrayLength;
639 default: // nodes for edges going along k direction
642 nodeENumbers [ ++indexE ] = n1;
643 nodeENumbers [ ++indexE ] = n2;
645 if (ax == 1) ++iNode;
647 if (ax == 2) iNode += _iArrayLength;
650 if (nbENodes != indexE+1) {
651 throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Wrong nbFNodes : " \
652 << nbENodes << " indexE : " << indexE ));
656 CONNECTIVITY * CellCNCT, * FaceCNCT, * EdgeCNCT;
658 // make connectivity for CELL
660 medGeometryElement aCellGeometry;
661 if (_kArrayLength) aCellGeometry = MED_HEXA8;
662 else if (_jArrayLength) aCellGeometry = MED_QUAD4;
663 else aCellGeometry = MED_SEG2;
666 CellCNCT = makeConnectivity (MED_CELL, aCellGeometry, nbCells, nbCNodes, nbMeshNodes, nodeCNumbers);
668 delete [] nodeCNumbers;
672 // CellCNCT->_descending = new MEDSKYLINEARRAY() ;
673 int * skyLineArrayIndex = new int [nbCells + 1];
674 int nbEntitySub = CellCNCT->_type[0].getNumberOfConstituents(1);
675 for (i=0, j=1; i <= nbCells; ++i, j += nbEntitySub)
676 skyLineArrayIndex [ i ] = j;
677 // CellCNCT->_descending->setMEDSKYLINEARRAY (nbCells, nbSub,
678 // skyLineArrayIndex, subNumbers);
679 CellCNCT->_descending = new MEDSKYLINEARRAY (nbCells, nbSub,
682 delete [] skyLineArrayIndex;
684 delete [] subNumbers;
686 // reverse descending
688 // CellCNCT->_reverseDescendingConnectivity = new MEDSKYLINEARRAY() ;
690 int * skyLineArrayIndex = new int [nbSub + 1];
691 for (i=0, j=1; i <= nbSub; ++i, j += 2)
692 skyLineArrayIndex [ i ] = j;
693 // CellCNCT->_reverseDescendingConnectivity->setMEDSKYLINEARRAY
694 // (nbSub, nbRevSub, skyLineArrayIndex, subRevNumbers);
696 CellCNCT->_reverseDescendingConnectivity =
697 new MEDSKYLINEARRAY(nbSub, nbRevSub, skyLineArrayIndex, subRevNumbers);
699 delete [] skyLineArrayIndex;
701 delete [] subRevNumbers;
703 // make connectivity for FACE and/or EDGE
706 FaceCNCT = makeConnectivity (MED_FACE, MED_QUAD4, nbFaces, nbFNodes, nbMeshNodes, nodeFNumbers);
708 delete [] nodeFNumbers;
710 CellCNCT->_constituent = FaceCNCT;
713 EdgeCNCT = makeConnectivity (MED_EDGE, MED_SEG2, nbEdges, nbENodes, nbMeshNodes, nodeENumbers);
715 delete [] nodeENumbers;
718 FaceCNCT->_constituent = EdgeCNCT;
720 CellCNCT->_constituent = EdgeCNCT;
723 MESH::_connectivity = CellCNCT;
725 (const_cast <GRID *> (this))->_is_connectivity_filled = true;
730 //=======================================================================
731 //function : getArrayLength
732 //purpose : return array length. Axis = [1,2,3] meaning [i,j,k],
733 //=======================================================================
735 int GRID::getArrayLength( const int Axis ) const throw (MEDEXCEPTION)
738 case 1: return _iArrayLength;
739 case 2: return _jArrayLength;
740 case 3: return _kArrayLength;
742 throw MED_EXCEPTION ( LOCALIZED( STRING("GRID::getArrayLength ( ") << Axis << ")"));
747 //=======================================================================
748 //function : getArrayValue
749 //purpose : return i-th array component. Axis = [1,2,3] meaning [i,j,k],
750 // exception if Axis out of [1-3] range
751 // exception if i is out of range 0 <= i < getArrayLength(Axis);
752 //=======================================================================
754 const double GRID::getArrayValue (const int Axis, const int i) const throw (MEDEXCEPTION)
756 if (i < 0 || i >= getArrayLength(Axis))
758 ( LOCALIZED(STRING("GRID::getArrayValue ( ") << Axis << ", " << i << ")"));
761 case 1: return _iArray[ i ];
762 case 2: return _jArray[ i ];
763 default: return _kArray[ i ];
768 //=======================================================================
769 //function : getEdgeNumber
771 //=======================================================================
773 int GRID::getEdgeNumber(const int Axis, const int i, const int j, const int k)
774 const throw (MEDEXCEPTION)
776 const char * LOC = "GRID::getEdgeNumber(Axis, i,j,k) :";
780 int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
781 int maxAxis = Len[ K ] ? 3 : 2;
783 if (Axis <= 0 || Axis > maxAxis)
784 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis out of range: " << Axis));
787 int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
793 if (Axis > 1) { // add all edges in i direction
795 Nb += Len[ I ]*Len[ J ]*Len[ K ];
798 if (Axis > 2) { // add all edges in j direction
800 Nb += Len[ I ]*Len[ J ]*Len[ K ];
808 //=======================================================================
809 //function : getFaceNumber
810 //purpose : return a NODE, EDGE, FACE, CELL number by its position in the grid.
811 // Axis [1,2,3] means one of directions: along i, j or k
812 // For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
813 // * a FACE which is normal to direction along given Axis;
814 // * an EDGE going along given Axis.
815 // Exception for Axis out of range
816 //=======================================================================
818 int GRID::getFaceNumber(const int Axis, const int i, const int j, const int k)
819 const throw (MEDEXCEPTION)
821 const char * LOC = "GRID::getFaceNumber(Axis, i,j,k) :";
825 // if (Axis <= 0 || Axis > 3)
826 if (Axis < 0 || Axis > 3)
827 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Axis = " << Axis));
829 int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
832 int Nb = 1 + i + j*Len[ I ] + k*Len[ J ]*Len[ K ];
835 if (Axis > 1) { // add all faces in i direction
837 Nb += Len[ I ]*Len[ J ]*Len[ K ];
840 if (Axis > 2) { // add all faces in j direction
842 Nb += Len[ I ]*Len[ J ]*Len[ K ];
850 //=======================================================================
851 //function : getNodePosition
853 //=======================================================================
855 void GRID::getNodePosition(const int Number, int& i, int& j, int& k) const
858 const char * LOC = "GRID::getNodePosition(Number, i,j,k) :";
862 if (Number <= 0 || Number > _numberOfNodes)
863 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
865 int Len[] = {_iArrayLength, _jArrayLength, _kArrayLength }, I=0, J=1;
866 // , K=2; !! UNUSED VARIABLE !!
868 int ijLen = Len[I] * Len[J]; // nb in a full k layer
869 int kLen = (Number - 1) % ijLen; // nb in the non full k layer
873 k = (Number - 1) / ijLen;
875 ////cout <<" NODE POS: " << Number << " - " << i << ", " << j << ", " << k << endl;
881 //=======================================================================
882 //function : getCellPosition
884 //=======================================================================
886 void GRID::getCellPosition(const int Number, int& i, int& j, int& k) const
889 const char * LOC = "GRID::getCellPosition(Number, i,j,k) :";
893 int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2;
894 // , K=3; !! UNUSED VARIABLE !!
896 // if (Number <= 0 || Number > getCellNumber(Len[I]-1, Len[J]-1, Len[K]-1))
897 // throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
899 int ijLen = Len[I] * Len[J]; // nb in a full k layer
900 int kLen = (Number - 1) % ijLen; // nb in the non full k layer
904 k = (Number - 1) / ijLen;
909 //=======================================================================
910 //function : getEdgePosition
912 //=======================================================================
914 void GRID::getEdgePosition(const int Number, int& Axis, int& i, int& j, int& k)
915 const throw (MEDEXCEPTION)
917 const char * LOC = "GRID::getEdgePosition(Number, i,j,k) :";
922 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no edges in the grid: "));
925 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
927 int Len[4] = {0,_iArrayLength, _jArrayLength, _kArrayLength }, I=1, J=2, K=3;
930 int maxAxis = _kArrayLength ? 3 : 2;
932 for (Axis = 1; Axis <= maxAxis; ++Axis)
934 Len[Axis]--; // one less edge in Axis direction
936 // max edge number in Axis direction
937 int maxNb = getEdgeNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
954 int ijLen = Len[I] * Len[J]; // nb in a full k layer
955 int kLen = (theNb - 1) % ijLen; // nb in the non full k layer
958 k = (theNb - 1) / ijLen;
966 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
969 //=======================================================================
970 //function : getFacePosition
971 //purpose : return position (i,j,k) of an entity #Number
972 // Axis [1,2,3] means one of directions: along i, j or k
973 // For Cell contituents (FACE or EDGE), Axis selects one of those having same (i,j,k):
974 // * a FACE which is normal to direction along given Axis;
975 // * an EDGE going along given Axis.
976 // Exception for Number out of range
977 //=======================================================================
979 void GRID::getFacePosition(const int Number, int& Axis, int& i, int& j, int& k)
980 const throw (MEDEXCEPTION)
982 const char * LOC = "GRID::getFacePosition(Number, i,j,k) :";
986 if (_kArrayLength == 0) {
987 getCellPosition(Number, i, j, k);
993 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "no faces in the grid: "));
996 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
998 int Len[4] = {0,_iArrayLength-1, _jArrayLength-1, _kArrayLength-1 }, I=1, J=2, K=3;
1001 for (Axis = 1; Axis <= 3; ++Axis)
1005 // max face number in Axis direction
1006 int maxNb = getFaceNumber (Axis, Len[I]-1, Len[J]-1, Len[K]-1);
1023 int ijLen = Len[I] * Len[J]; // nb in a full k layer
1024 int kLen = (theNb - 1) % ijLen; // nb in the non full k layer
1027 k = (theNb - 1) / ijLen;
1035 throw MED_EXCEPTION ( LOCALIZED(STRING(LOC) << "Number is out of range: " << Number));
1038 //=======================================================================
1039 //function : writeUnstructured
1040 //purpose : write a Grid as an Unstructured mesh
1041 //=======================================================================
1043 void GRID::writeUnstructured(int index, const string & driverName)
1045 const char * LOC = "GRID::writeUnstructured(int index=0, const string & driverName) : ";
1048 if ( _drivers[index] ) {
1053 _drivers[index]->open();
1054 if (driverName != "") _drivers[index]->setMeshName(driverName);
1055 _drivers[index]->write();
1056 _drivers[index]->close();
1061 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1062 << "The index given is invalid, index must be between 0 and |"
1069 void GRID::read(int index)
1071 const char * LOC = "GRID::read(int index=0) : ";
1074 if (_drivers[index]) {
1075 _drivers[index]->open();
1076 _drivers[index]->read();
1077 _drivers[index]->close();
1080 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1081 << "The index given is invalid, index must be between 0 and |"
1086 fillMeshAfterRead();