]> SALOME platform Git repositories - plugins/ghs3dplugin.git/blob - src/GHS3DPlugin_GHS3D.cxx
Salome HOME
0019296: EDF 681 SMESH - Pre-evaluation of the number of elements before mesh
[plugins/ghs3dplugin.git] / src / GHS3DPlugin_GHS3D.cxx
1 //  Copyright (C) 2004-2008  CEA/DEN, EDF R&D
2 //
3 //  This library is free software; you can redistribute it and/or
4 //  modify it under the terms of the GNU Lesser General Public
5 //  License as published by the Free Software Foundation; either
6 //  version 2.1 of the License.
7 //
8 //  This library is distributed in the hope that it will be useful,
9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 //  Lesser General Public License for more details.
12 //
13 //  You should have received a copy of the GNU Lesser General Public
14 //  License along with this library; if not, write to the Free Software
15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 //=============================================================================
20 // File      : GHS3DPlugin_GHS3D.cxx
21 // Created   : 
22 // Author    : Edward AGAPOV, modified by Lioka RAZAFINDRAZAKA (CEA) 09/02/2007
23 // Project   : SALOME
24 // $Header$
25 //=============================================================================
26 //
27 using namespace std;
28
29 #include "GHS3DPlugin_GHS3D.hxx"
30 #include "GHS3DPlugin_Hypothesis.hxx"
31
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
34 #include "SMESH_Comment.hxx"
35 #include "SMESH_MesherHelper.hxx"
36 #include "SMESH_MeshEditor.hxx"
37
38 #include "SMDS_MeshElement.hxx"
39 #include "SMDS_MeshNode.hxx"
40 #include "SMDS_FaceOfNodes.hxx"
41 #include "SMDS_VolumeOfNodes.hxx"
42
43 #include <BRepAdaptor_Surface.hxx>
44 #include <BRepBndLib.hxx>
45 #include <BRepClass3d_SolidClassifier.hxx>
46 #include <BRepTools.hxx>
47 #include <BRep_Tool.hxx>
48 #include <Bnd_Box.hxx>
49 #include <GeomAPI_ProjectPointOnSurf.hxx>
50 #include <OSD_File.hxx>
51 #include <Precision.hxx>
52 #include <Quantity_Parameter.hxx>
53 #include <Standard_ErrorHandler.hxx>
54 #include <Standard_Failure.hxx>
55 #include <TopExp.hxx>
56 #include <TopExp_Explorer.hxx>
57 #include <TopTools_IndexedMapOfShape.hxx>
58 #include <TopTools_ListIteratorOfListOfShape.hxx>
59 #include <TopoDS.hxx>
60 //#include <BRepClass_FaceClassifier.hxx>
61 #include <TopTools_MapOfShape.hxx>
62 #include <BRepGProp.hxx>
63 #include <GProp_GProps.hxx>
64
65 #include "utilities.h"
66
67 #ifndef WIN32
68 #include <sys/sysinfo.h>
69 #endif
70
71 //#include <Standard_Stream.hxx>
72
73
74 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
75
76 #ifdef _DEBUG_
77 #define DUMP(txt) \
78 //  cout << txt
79 #else
80 #define DUMP(txt)
81 #endif
82
83 extern "C"
84 {
85 #ifndef WNT
86 #include <unistd.h>
87 #include <sys/mman.h>
88 #endif
89 #include <sys/stat.h>
90 #include <fcntl.h>
91 }
92
93 #define HOLE_ID -1
94
95 //=============================================================================
96 /*!
97  *  
98  */
99 //=============================================================================
100
101 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
102   : SMESH_3D_Algo(hypId, studyId, gen)
103 {
104   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
105   _name = "GHS3D_3D";
106   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
107   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
108   _iShape=0;
109   _nbShape=0;
110   _compatibleHypothesis.push_back("GHS3D_Parameters");
111   _requireShape = false; // can work without shape
112 }
113
114 //=============================================================================
115 /*!
116  *  
117  */
118 //=============================================================================
119
120 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
121 {
122   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
123 }
124
125 //=============================================================================
126 /*!
127  *  
128  */
129 //=============================================================================
130
131 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
132                                           const TopoDS_Shape& aShape,
133                                           Hypothesis_Status&  aStatus )
134 {
135   aStatus = SMESH_Hypothesis::HYP_OK;
136
137   // there is only one compatible Hypothesis so far
138   _hyp = 0;
139   _keepFiles = false;
140   const list <const SMESHDS_Hypothesis * >& hyps = GetUsedHypothesis(aMesh, aShape);
141   if ( !hyps.empty() )
142     _hyp = static_cast<const GHS3DPlugin_Hypothesis*> ( hyps.front() );
143   if ( _hyp )
144     _keepFiles = _hyp->GetKeepFiles();
145
146   return true;
147 }
148
149 //=======================================================================
150 //function : findShape
151 //purpose  : 
152 //=======================================================================
153
154 static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
155                               TopoDS_Shape        aShape,
156                               const TopoDS_Shape  shape[],
157                               double**            box,
158                               const int           nShape,
159                               TopAbs_State *      state = 0)
160 {
161   gp_XYZ aPnt(0,0,0);
162   int j, iShape, nbNode = 4;
163
164   for ( j=0; j<nbNode; j++ ) {
165     gp_XYZ p ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
166     if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) {
167       aPnt = p;
168       break;
169     }
170     aPnt += p;
171   }
172
173   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
174   if (state) *state = SC.State();
175   if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
176     for (iShape = 0; iShape < nShape; iShape++) {
177       aShape = shape[iShape];
178       if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
179               aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
180               aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
181         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
182         if (state) *state = SC.State();
183         if (SC.State() == TopAbs_IN)
184           break;
185       }
186     }
187   }
188   return aShape;
189 }
190
191 //=======================================================================
192 //function : readMapIntLine
193 //purpose  : 
194 //=======================================================================
195
196 static char* readMapIntLine(char* ptr, int tab[]) {
197   long int intVal;
198   cout << endl;
199
200   for ( int i=0; i<17; i++ ) {
201     intVal = strtol(ptr, &ptr, 10);
202     if ( i < 3 )
203       tab[i] = intVal;
204   }
205   return ptr;
206 }
207
208 //=======================================================================
209 //function : countShape
210 //purpose  :
211 //=======================================================================
212
213 template < class Mesh, class Shape >
214 static int countShape( Mesh* mesh, Shape shape ) {
215   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
216   int nbShape = 0;
217   for ( ; expShape.More(); expShape.Next() ) {
218       nbShape++;
219   }
220   return nbShape;
221 }
222
223 //=======================================================================
224 //function : writeFaces
225 //purpose  : 
226 //=======================================================================
227
228 static bool writeFaces (ofstream &            theFile,
229                         SMESHDS_Mesh *        theMesh,
230                         const map <int,int> & theSmdsToGhs3dIdMap)
231 {
232   // record structure:
233   //
234   // NB_ELEMS DUMMY_INT
235   // Loop from 1 to NB_ELEMS
236   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
237
238   int nbShape = countShape( theMesh, TopAbs_FACE );
239
240   int *tabID;             tabID    = new int[nbShape];
241   TopoDS_Shape *tabShape; tabShape = new TopoDS_Shape[nbShape];
242   TopoDS_Shape aShape;
243   SMESHDS_SubMesh* theSubMesh;
244   const SMDS_MeshElement* aFace;
245   const char* space    = "  ";
246   const int   dummyint = 0;
247   map<int,int>::const_iterator itOnMap;
248   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
249   int shapeID, nbNodes, aSmdsID;
250   bool idFound;
251
252   cout << "    " << theMesh->NbFaces() << " shapes of 2D dimension" << endl;
253   cout << endl;
254
255   theFile << space << theMesh->NbFaces() << space << dummyint << endl;
256
257   TopExp_Explorer expface( theMesh->ShapeToMesh(), TopAbs_FACE );
258   for ( int i = 0; expface.More(); expface.Next(), i++ ) {
259     tabID[i] = 0;
260     aShape   = expface.Current();
261     shapeID  = theMesh->ShapeToIndex( aShape );
262     idFound  = false;
263     for ( int j=0; j<=i; j++) {
264       if ( shapeID == tabID[j] ) {
265         idFound = true;
266         break;
267       }
268     }
269     if ( ! idFound ) {
270       tabID[i]    = shapeID;
271       tabShape[i] = aShape;
272     }
273   }
274   for ( int i =0; i < nbShape; i++ ) {
275     if ( tabID[i] != 0 ) {
276       aShape      = tabShape[i];
277       shapeID     = tabID[i];
278       theSubMesh  = theMesh->MeshElements( aShape );
279       if ( !theSubMesh ) continue;
280       itOnSubMesh = theSubMesh->GetElements();
281       while ( itOnSubMesh->more() ) {
282         aFace   = itOnSubMesh->next();
283         nbNodes = aFace->NbNodes();
284
285         theFile << space << nbNodes;
286
287         itOnSubFace = aFace->nodesIterator();
288         while ( itOnSubFace->more() ) {
289           // find GHS3D ID
290           aSmdsID = itOnSubFace->next()->GetID();
291           itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
292           ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
293
294           theFile << space << (*itOnMap).second;
295         }
296
297         // (NB_NODES + 1) times: DUMMY_INT
298         for ( int j=0; j<=nbNodes; j++)
299           theFile << space << dummyint;
300
301         theFile << endl;
302       }
303     }
304   }
305
306   delete [] tabID;
307   delete [] tabShape;
308
309   return true;
310 }
311
312 //=======================================================================
313 //function : writeFaces
314 //purpose  : Write Faces in case if generate 3D mesh w/o geometry
315 //=======================================================================
316
317 static bool writeFaces (ofstream &            theFile,
318                         SMESHDS_Mesh *        theMesh,
319                         vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
320 {
321   // record structure:
322   //
323   // NB_ELEMS DUMMY_INT
324   // Loop from 1 to NB_ELEMS
325   //   NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
326
327
328   int nbFaces = 0;
329   list< const SMDS_MeshElement* > faces;
330   list< const SMDS_MeshElement* >::iterator f;
331   map< const SMDS_MeshNode*,int >::iterator it;
332   SMDS_ElemIteratorPtr nodeIt;
333   const SMDS_MeshElement* elem;
334   int nbNodes;
335
336   const char* space    = "  ";
337   const int   dummyint = 0;
338
339   //get all faces from mesh
340   SMDS_FaceIteratorPtr eIt = theMesh->facesIterator();
341   while ( eIt->more() ) {
342     const SMDS_MeshElement* elem = eIt->next();
343     if ( !elem )
344       return false;
345     faces.push_back( elem );
346     nbFaces++;
347   }
348
349   if ( nbFaces == 0 )
350     return false;
351   
352   cout << "The initial 2D mesh contains " << nbFaces << " faces and ";
353
354   // NB_ELEMS DUMMY_INT
355   theFile << space << nbFaces << space << dummyint << endl;
356
357   // Loop from 1 to NB_ELEMS
358
359   map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
360   f = faces.begin();
361   for ( ; f != faces.end(); ++f )
362   {
363     // NB_NODES PER FACE
364     elem = *f;
365     nbNodes = elem->NbNodes();
366     theFile << space << nbNodes;
367
368     // NODE_NB_1 NODE_NB_2 ...
369     nodeIt = elem->nodesIterator();
370     while ( nodeIt->more() )
371     {
372       // find GHS3D ID
373       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
374       int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
375       it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
376       theFile << space << it->second;
377     }
378
379     // (NB_NODES + 1) times: DUMMY_INT
380     for ( int i=0; i<=nbNodes; i++)
381       theFile << space << dummyint;
382
383     theFile << endl;
384   }
385
386   // put nodes to theNodeByGhs3dId vector
387   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
388   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
389   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
390   {
391     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
392   }
393
394   return true;
395 }
396
397 //=======================================================================
398 //function : writePoints
399 //purpose  : 
400 //=======================================================================
401
402 static bool writePoints (ofstream &                       theFile,
403                          SMESHDS_Mesh *                   theMesh,
404                          map <int,int> &                  theSmdsToGhs3dIdMap,
405                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap)
406 {
407   // record structure:
408   //
409   // NB_NODES
410   // Loop from 1 to NB_NODES
411   //   X Y Z DUMMY_INT
412
413   int nbNodes = theMesh->NbNodes();
414   if ( nbNodes == 0 )
415     return false;
416
417   const char* space    = "  ";
418   const int   dummyint = 0;
419
420   int aGhs3dID = 1;
421   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
422   const SMDS_MeshNode* node;
423
424   // NB_NODES
425   theFile << space << nbNodes << endl;
426   cout << endl;
427   cout << "The initial 2D mesh contains :" << endl;
428   cout << "    " << nbNodes << " vertices" << endl;
429
430   // Loop from 1 to NB_NODES
431
432   while ( it->more() )
433   {
434     node = it->next();
435     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
436     theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
437     aGhs3dID++;
438
439     // X Y Z DUMMY_INT
440     theFile
441       << space << node->X()
442       << space << node->Y()
443       << space << node->Z()
444       << space << dummyint;
445
446     theFile << endl;
447   }
448
449   return true;
450 }
451
452 //=======================================================================
453 //function : writePoints
454 //purpose  : 
455 //=======================================================================
456
457 static bool writePoints (ofstream &                            theFile,
458                          SMESHDS_Mesh *                        theMesh,
459                          const vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
460 {
461   // record structure:
462   //
463   // NB_NODES
464   // Loop from 1 to NB_NODES
465   //   X Y Z DUMMY_INT
466
467   //int nbNodes = theMesh->NbNodes();
468   int nbNodes = theNodeByGhs3dId.size();
469   if ( nbNodes == 0 )
470     return false;
471
472   const char* space    = "  ";
473   const int   dummyint = 0;
474
475   const SMDS_MeshNode* node;
476
477   // NB_NODES
478   theFile << space << nbNodes << endl;
479   cout << nbNodes << " nodes" << endl;
480
481   // Loop from 1 to NB_NODES
482
483   vector<const SMDS_MeshNode*>::const_iterator nodeIt = theNodeByGhs3dId.begin();
484   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
485   for ( ; nodeIt != after; ++nodeIt )
486   {
487     node = *nodeIt;
488
489     // X Y Z DUMMY_INT
490     theFile
491       << space << node->X()
492       << space << node->Y()
493       << space << node->Z()
494       << space << dummyint;
495
496     theFile << endl;
497   }
498
499   return true;
500 }
501
502 //=======================================================================
503 //function : findShapeID
504 //purpose  : find the solid corresponding to GHS3D sub-domain following
505 //           the technique proposed in GHS3D manual in chapter
506 //           "B.4 Subdomain (sub-region) assignment"
507 //=======================================================================
508
509 static int findShapeID(SMESH_Mesh&          mesh,
510                        const SMDS_MeshNode* node1,
511                        const SMDS_MeshNode* node2,
512                        const SMDS_MeshNode* node3,
513                        const bool           toMeshHoles)
514 {
515   const int invalidID = 0;
516   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
517
518   // face the nodes belong to
519   const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
520   if ( !face )
521     return invalidID;
522
523   // geom face the face assigned to
524   SMESH_MeshEditor editor(&mesh);
525   int geomFaceID = editor.FindShape( face );
526   if ( !geomFaceID )
527     return invalidID;
528   TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
529   if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
530     return invalidID;
531   TopoDS_Face geomFace = TopoDS::Face( shape );
532
533   // solids bounded by geom face
534   TopTools_IndexedMapOfShape solids, shells;
535   TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
536   for ( ; ansIt.More(); ansIt.Next() ) {
537     switch ( ansIt.Value().ShapeType() ) {
538     case TopAbs_SOLID:
539       solids.Add( ansIt.Value() ); break;
540     case TopAbs_SHELL:
541       shells.Add( ansIt.Value() ); break;
542     default:;
543     }
544   }
545   // analyse found solids
546   if ( solids.Extent() == 0 || shells.Extent() == 0)
547     return invalidID;
548
549   const TopoDS_Solid& solid1 = TopoDS::Solid( solids(1) );
550   if ( solids.Extent() == 1 )
551   {
552     if ( toMeshHoles )
553       return meshDS->ShapeToIndex( solid1 );
554
555     // - Are we at a hole boundary face?
556     if ( shells(1).IsSame( BRepTools::OuterShell( solid1 )) )
557     { // - No, but maybe a hole is bound by two shapes? Does shells(1) touches another shell?
558       bool touch = false;
559       TopExp_Explorer eExp( shells(1), TopAbs_EDGE );
560       // check if any edge of shells(1) belongs to another shell
561       for ( ; eExp.More() && !touch; eExp.Next() ) {
562         ansIt = mesh.GetAncestors( eExp.Current() );
563         for ( ; ansIt.More() && !touch; ansIt.Next() ) {
564           if ( ansIt.Value().ShapeType() == TopAbs_SHELL )
565             touch = ( !ansIt.Value().IsSame( shells(1) ));
566         }
567       }
568       if (!touch)
569         return meshDS->ShapeToIndex( solid1 );
570     }
571   }
572   // find orientation of geom face within the first solid
573   TopExp_Explorer fExp( solid1, TopAbs_FACE );
574   for ( ; fExp.More(); fExp.Next() )
575     if ( geomFace.IsSame( fExp.Current() )) {
576       geomFace = TopoDS::Face( fExp.Current() );
577       break;
578     }
579   if ( !fExp.More() )
580     return invalidID; // face not found
581
582   // normale to triangle
583   gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
584   gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
585   gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
586   gp_Vec vec12( node1Pnt, node2Pnt );
587   gp_Vec vec13( node1Pnt, node3Pnt );
588   gp_Vec meshNormal = vec12 ^ vec13;
589   if ( meshNormal.SquareMagnitude() < DBL_MIN )
590     return invalidID;
591
592   // find UV of node1 on geomFace
593   SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
594   const SMDS_MeshNode* nNotOnSeamEdge = 0;
595   if ( helper.IsSeamShape( node1->GetPosition()->GetShapeId() ))
596     if ( helper.IsSeamShape( node2->GetPosition()->GetShapeId() ))
597       nNotOnSeamEdge = node3;
598     else
599       nNotOnSeamEdge = node2;
600   gp_XY uv = helper.GetNodeUV( geomFace, node1, nNotOnSeamEdge );
601   // check that uv is correct
602   BRepAdaptor_Surface surface( geomFace );
603   double tol = BRep_Tool::Tolerance( geomFace );
604   if ( node1Pnt.Distance( surface.Value( uv.X(), uv.Y() )) > 2 * tol ) {
605     GeomAPI_ProjectPointOnSurf projector( node1Pnt, surface.Surface().Surface(), 2 * tol );
606     if ( !projector.IsDone() || projector.NbPoints() < 1 || projector.LowerDistance() > 2 * tol)
607       return invalidID;
608     Quantity_Parameter U,V;
609     projector.LowerDistanceParameters(U,V);
610     if ( node1Pnt.Distance( surface.Value( U, V )) > 2 * tol )
611       return invalidID;
612     uv.SetCoord( U,V );
613   }
614   // normale to geomFace at UV
615   gp_Vec du, dv;
616   surface.D1( uv.X(), uv.Y(), node1Pnt, du, dv );
617   gp_Vec geomNormal = du ^ dv;
618   if ( geomNormal.SquareMagnitude() < DBL_MIN )
619     return findShapeID( mesh, node2, node3, node1, toMeshHoles );
620   if ( geomFace.Orientation() == TopAbs_REVERSED )
621     geomNormal.Reverse();
622
623   // compare normals
624   bool isReverse = ( meshNormal * geomNormal ) < 0;
625   if ( !isReverse )
626     return meshDS->ShapeToIndex( solid1 );
627
628   if ( solids.Extent() == 1 )
629     return HOLE_ID; // we are inside a hole
630   else
631     return meshDS->ShapeToIndex( solids(2) );
632 }
633                        
634 //=======================================================================
635 //function : readResultFile
636 //purpose  : 
637 //=======================================================================
638
639 static bool readResultFile(const int                       fileOpen,
640                            SMESH_Mesh&                     theMesh,
641                            TopoDS_Shape                    tabShape[],
642                            double**                        tabBox,
643                            const int                       nbShape,
644                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
645                            bool                            toMeshHoles)
646 {
647   struct stat status;
648   size_t      length;
649
650   char *ptr, *mapPtr;
651   char *tetraPtr;
652   char *shapePtr;
653
654   SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
655
656   int fileStat;
657   int nbElems, nbNodes, nbInputNodes;
658   int nodeId/*, triangleId*/;
659   int nbTriangle;
660   int ID, shapeID, ghs3dShapeID;
661   int IdShapeRef = 1;
662   int compoundID =
663     nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
664
665   int *tab, *tabID, *nodeID, *nodeAssigne;
666   double *coord;
667   const SMDS_MeshNode **node;
668
669   tab    = new int[3];
670   //tabID  = new int[nbShape];
671   nodeID = new int[4];
672   coord  = new double[3];
673   node   = new const SMDS_MeshNode*[4];
674
675   TopoDS_Shape aSolid;
676   SMDS_MeshNode * aNewNode;
677   map <int,const SMDS_MeshNode*>::iterator itOnNode;
678   SMDS_MeshElement* aTet;
679 #ifdef _DEBUG_
680   set<int> shapeIDs;
681 #endif
682
683   // Read the file state
684   fileStat = fstat(fileOpen, &status);
685   length   = status.st_size;
686
687   // Mapping the result file into memory
688   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
689   mapPtr = ptr;
690
691   ptr      = readMapIntLine(ptr, tab);
692   tetraPtr = ptr;
693
694   nbElems      = tab[0];
695   nbNodes      = tab[1];
696   nbInputNodes = tab[2];
697
698   nodeAssigne = new int[ nbNodes+1 ];
699
700   if (nbShape > 0)
701     aSolid = tabShape[0];
702
703   // Reading the nodeId
704   for (int i=0; i < 4*nbElems; i++)
705     nodeId = strtol(ptr, &ptr, 10);
706
707   // Reading the nodeCoor and update the nodeMap
708   for (int iNode=1; iNode <= nbNodes; iNode++) {
709     for (int iCoor=0; iCoor < 3; iCoor++)
710       coord[ iCoor ] = strtod(ptr, &ptr);
711     nodeAssigne[ iNode ] = 1;
712     if ( iNode > nbInputNodes ) {
713       nodeAssigne[ iNode ] = 0;
714       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
715       theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode ));
716     }
717   }
718
719   // Reading the number of triangles which corresponds to the number of sub-domains
720   nbTriangle = strtol(ptr, &ptr, 10);
721
722   tabID = new int[nbTriangle];
723   for (int i=0; i < nbTriangle; i++) {
724     tabID[i] = 0;
725     // find the solid corresponding to GHS3D sub-domain following
726     // the technique proposed in GHS3D manual in chapter
727     // "B.4 Subdomain (sub-region) assignment"
728     int nodeId1 = strtol(ptr, &ptr, 10);
729     int nodeId2 = strtol(ptr, &ptr, 10);
730     int nodeId3 = strtol(ptr, &ptr, 10);
731     if ( nbTriangle > 1 ) {
732       const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
733       const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
734       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
735       try {
736         OCC_CATCH_SIGNALS;
737         tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
738         // -- 0020330: Pb with ghs3d as a submesh
739         // check that found shape is to be meshed
740         if ( tabID[i] > 0 ) {
741           const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] );
742           bool isToBeMeshed = false;
743           for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS )
744             isToBeMeshed = foundShape.IsSame( tabShape[ iS ]);
745           if ( !isToBeMeshed )
746             tabID[i] = HOLE_ID;
747         }
748         // END -- 0020330: Pb with ghs3d as a submesh
749 #ifdef _DEBUG_
750         cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << endl;
751 #endif
752       } catch ( Standard_Failure ) {
753       } catch (...) {}
754     }
755   }
756
757   shapePtr = ptr;
758
759   if ( nbTriangle <= nbShape ) // no holes
760     toMeshHoles = true; // not avoid creating tetras in holes
761
762   // Associating the tetrahedrons to the shapes
763   shapeID = compoundID;
764   for (int iElem = 0; iElem < nbElems; iElem++) {
765     for (int iNode = 0; iNode < 4; iNode++) {
766       ID = strtol(tetraPtr, &tetraPtr, 10);
767       itOnNode = theGhs3dIdToNodeMap.find(ID);
768       node[ iNode ] = itOnNode->second;
769       nodeID[ iNode ] = ID;
770     }
771     // We always run GHS3D with "to mesh holes'==TRUE but we must not create
772     // tetras within holes depending on hypo option,
773     // so we first check if aTet is inside a hole and then create it 
774     //aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
775     if ( nbTriangle > 1 ) {
776       shapeID = HOLE_ID; // negative shapeID means not to create tetras if !toMeshHoles
777       ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
778       if ( tabID[ ghs3dShapeID ] == 0 ) {
779         TopAbs_State state;
780         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
781         if ( toMeshHoles || state == TopAbs_IN )
782           shapeID = theMeshDS->ShapeToIndex( aSolid );
783         tabID[ ghs3dShapeID ] = shapeID;
784       }
785       else
786         shapeID = tabID[ ghs3dShapeID ];
787     }
788     else if ( nbShape > 1 ) {
789       // Case where nbTriangle == 1 while nbShape == 2 encountered
790       // with compound of 2 boxes and "To mesh holes"==False,
791       // so there are no subdomains specified for each tetrahedron.
792       // Try to guess a solid by a node already bound to shape
793       shapeID = 0;
794       for ( int i=0; i<4 && shapeID==0; i++ ) {
795         if ( nodeAssigne[ nodeID[i] ] == 1 &&
796              node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
797              node[i]->GetPosition()->GetShapeId() > 1 )
798         {
799           shapeID = node[i]->GetPosition()->GetShapeId();
800         }
801       }
802       if ( shapeID==0 ) {
803         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
804         shapeID = theMeshDS->ShapeToIndex( aSolid );
805       }
806     }
807     // set new nodes and tetrahedron onto the shape
808     for ( int i=0; i<4; i++ ) {
809       if ( nodeAssigne[ nodeID[i] ] == 0 ) {
810         if ( shapeID != HOLE_ID )
811           theMeshDS->SetNodeInVolume( node[i], shapeID );
812         nodeAssigne[ nodeID[i] ] = shapeID;
813       }
814     }
815     if ( toMeshHoles || shapeID != HOLE_ID ) {
816       aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
817       theMeshDS->SetMeshElementOnShape( aTet, shapeID );
818     }
819 #ifdef _DEBUG_
820     shapeIDs.insert( shapeID );
821 #endif
822   }
823   // Remove nodes of tetras inside holes if !toMeshHoles
824   if ( !toMeshHoles ) {
825     itOnNode = theGhs3dIdToNodeMap.find( nbInputNodes );
826     for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
827       ID = itOnNode->first;
828       if ( nodeAssigne[ ID ] == HOLE_ID )
829         theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
830     }
831   }
832
833   if ( nbElems )
834     cout << nbElems << " tetrahedrons have been associated to " << nbShape << " shapes" << endl;
835   munmap(mapPtr, length);
836   close(fileOpen);
837
838   delete [] tab;
839   delete [] tabID;
840   delete [] nodeID;
841   delete [] coord;
842   delete [] node;
843   delete [] nodeAssigne;
844
845 #ifdef _DEBUG_
846   if ( shapeIDs.size() != nbShape ) {
847     cout << "Only " << shapeIDs.size() << " solids of " << nbShape << " found" << endl;
848     for (int i=0; i<nbShape; i++) {
849       shapeID = theMeshDS->ShapeToIndex( tabShape[i] );
850       if ( shapeIDs.find( shapeID ) == shapeIDs.end() )
851         cout << "  Solid #" << shapeID << " not found" << endl;
852     }
853   }
854 #endif
855
856   return true;
857 }
858
859 //=======================================================================
860 //function : readResultFile
861 //purpose  : 
862 //=======================================================================
863
864 static bool readResultFile(const int                      fileOpen,
865                            SMESHDS_Mesh*                  theMeshDS,
866                            TopoDS_Shape                   aSolid,
867                            vector <const SMDS_MeshNode*>& theNodeByGhs3dId) {
868
869   struct stat  status;
870   size_t       length;
871
872   char *ptr, *mapPtr;
873   char *tetraPtr;
874   char *shapePtr;
875
876   int fileStat;
877   int nbElems, nbNodes, nbInputNodes;
878   int nodeId, triangleId;
879   int nbTriangle;
880   int ID, shapeID;
881
882   int *tab;
883   double *coord;
884   const SMDS_MeshNode **node;
885
886   tab   = new int[3];
887   coord = new double[3];
888   node  = new const SMDS_MeshNode*[4];
889
890   SMDS_MeshNode * aNewNode;
891   map <int,const SMDS_MeshNode*>::iterator IdNode;
892   SMDS_MeshElement* aTet;
893
894   // Read the file state
895   fileStat = fstat(fileOpen, &status);
896   length   = status.st_size;
897
898   // Mapping the result file into memory
899   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
900   mapPtr = ptr;
901
902   ptr      = readMapIntLine(ptr, tab);
903   tetraPtr = ptr;
904
905   nbElems      = tab[0];
906   nbNodes      = tab[1];
907   nbInputNodes = tab[2];
908
909   theNodeByGhs3dId.resize( nbNodes );
910
911   // Reading the nodeId
912   for (int i=0; i < 4*nbElems; i++)
913     nodeId = strtol(ptr, &ptr, 10);
914
915   // Reading the nodeCoor and update the nodeMap
916   shapeID = theMeshDS->ShapeToIndex( aSolid );
917   for (int iNode=0; iNode < nbNodes; iNode++) {
918     for (int iCoor=0; iCoor < 3; iCoor++)
919       coord[ iCoor ] = strtod(ptr, &ptr);
920     if ((iNode+1) > nbInputNodes) {
921       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
922       theMeshDS->SetNodeInVolume( aNewNode, shapeID );
923       theNodeByGhs3dId[ iNode ] = aNewNode;
924     }
925   }
926
927   // Reading the triangles
928   nbTriangle = strtol(ptr, &ptr, 10);
929
930   for (int i=0; i < 3*nbTriangle; i++)
931     triangleId = strtol(ptr, &ptr, 10);
932
933   shapePtr = ptr;
934
935   // Associating the tetrahedrons to the shapes
936   for (int iElem = 0; iElem < nbElems; iElem++) {
937     for (int iNode = 0; iNode < 4; iNode++) {
938       ID = strtol(tetraPtr, &tetraPtr, 10);
939       node[ iNode ] = theNodeByGhs3dId[ ID-1 ];
940     }
941     aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
942     shapeID = theMeshDS->ShapeToIndex( aSolid );
943     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
944   }
945   if ( nbElems )
946     cout << nbElems << " tetrahedrons have been associated to " << nbTriangle << " shapes" << endl;
947   munmap(mapPtr, length);
948   close(fileOpen);
949
950   delete [] tab;
951   delete [] coord;
952   delete [] node;
953
954   return true;
955 }
956
957 //=============================================================================
958 /*!
959  *Here we are going to use the GHS3D mesher
960  */
961 //=============================================================================
962
963 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
964                                 const TopoDS_Shape& theShape)
965 {
966   bool Ok(false);
967   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
968
969   // we count the number of shapes
970   // _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with ghs3d as a submesh
971   _nbShape = 0;
972   TopExp_Explorer expBox ( theShape, TopAbs_SOLID );
973   for ( ; expBox.More(); expBox.Next() )
974     _nbShape++;
975
976   // create bounding box for every shape inside the compound
977
978   int iShape = 0;
979   TopoDS_Shape* tabShape;
980   double**      tabBox;
981   tabShape = new TopoDS_Shape[_nbShape];
982   tabBox   = new double*[_nbShape];
983   for (int i=0; i<_nbShape; i++)
984     tabBox[i] = new double[6];
985   Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
986
987   // TopExp_Explorer expBox (meshDS->ShapeToMesh(), TopAbs_SOLID); -- 0020330:...ghs3d as a submesh
988   for (expBox.ReInit(); expBox.More(); expBox.Next()) {
989     tabShape[iShape] = expBox.Current();
990     Bnd_Box BoundingBox;
991     BRepBndLib::Add(expBox.Current(), BoundingBox);
992     BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
993     tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
994     tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
995     tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
996     iShape++;
997   }
998
999   // a unique working file name
1000   // to avoid access to the same files by eg different users
1001   TCollection_AsciiString aGenericName
1002     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
1003
1004   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
1005   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
1006   aFacesFileName  = aGenericName + ".faces";  // in faces
1007   aPointsFileName = aGenericName + ".points"; // in points
1008   aResultFileName = aGenericName + ".noboite";// out points and volumes
1009   aBadResFileName = aGenericName + ".boite";  // out bad result
1010   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
1011   aLogFileName    = aGenericName + ".log";    // log
1012
1013   // -----------------
1014   // make input files
1015   // -----------------
1016
1017   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
1018   ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
1019
1020   Ok =
1021 #ifdef WIN32
1022     aFacesFile->is_open() && aPointsFile->is_open();
1023 #else
1024     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
1025 #endif
1026   if (!Ok) {
1027     INFOS( "Can't write into " << aFacesFileName);
1028     return error(SMESH_Comment("Can't write into ") << aFacesFileName);
1029   }
1030   map <int,int> aSmdsToGhs3dIdMap;
1031   map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
1032
1033   Ok = writePoints( aPointsFile, meshDS, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap ) &&
1034        writeFaces ( aFacesFile,  meshDS, aSmdsToGhs3dIdMap );
1035
1036   aFacesFile.close();
1037   aPointsFile.close();
1038
1039   if ( ! Ok ) {
1040     if ( !_keepFiles ) {
1041       OSD_File( aFacesFileName ).Remove();
1042       OSD_File( aPointsFileName ).Remove();
1043     }
1044     return error(COMPERR_BAD_INPUT_MESH);
1045   }
1046   OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
1047
1048   // -----------------
1049   // run ghs3d mesher
1050   // -----------------
1051
1052   TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
1053   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
1054   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
1055
1056   cout << endl;
1057   cout << "Ghs3d execution..." << endl;
1058   cout << cmd << endl;
1059
1060   system( cmd.ToCString() ); // run
1061
1062   cout << endl;
1063   cout << "End of Ghs3d execution !" << endl;
1064
1065   // --------------
1066   // read a result
1067   // --------------
1068
1069   // Mapping the result file
1070
1071   int fileOpen;
1072   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1073   if ( fileOpen < 0 ) {
1074     cout << endl;
1075     cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << endl;
1076     cout << "Log: " << aLogFileName << endl;
1077     Ok = false;
1078   }
1079   else {
1080     bool toMeshHoles =
1081       _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
1082     Ok = readResultFile( fileOpen, theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
1083                          toMeshHoles );
1084   }
1085
1086   // ---------------------
1087   // remove working files
1088   // ---------------------
1089
1090   if ( Ok )
1091   {
1092     if ( !_keepFiles )
1093       OSD_File( aLogFileName ).Remove();
1094   }
1095   else if ( OSD_File( aLogFileName ).Size() > 0 )
1096   {
1097     // get problem description from the log file
1098     _Ghs2smdsConvertor conv( aGhs3dIdToNodeMap );
1099     storeErrorDescription( aLogFileName, conv );
1100   }
1101   else
1102   {
1103     // the log file is empty
1104     OSD_File( aLogFileName ).Remove();
1105     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
1106     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
1107   }
1108
1109   if ( !_keepFiles ) {
1110     OSD_File( aFacesFileName ).Remove();
1111     OSD_File( aPointsFileName ).Remove();
1112     OSD_File( aResultFileName ).Remove();
1113     OSD_File( aBadResFileName ).Remove();
1114     OSD_File( aBbResFileName ).Remove();
1115   }
1116   cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
1117   if ( !Ok )
1118     cout << "not ";
1119   cout << "treated !" << endl;
1120   cout << endl;
1121
1122   _nbShape = 0;    // re-initializing _nbShape for the next Compute() method call
1123   delete [] tabShape;
1124   delete [] tabBox;
1125
1126   return Ok;
1127 }
1128
1129 //=============================================================================
1130 /*!
1131  *Here we are going to use the GHS3D mesher w/o geometry
1132  */
1133 //=============================================================================
1134 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
1135                                 SMESH_MesherHelper* aHelper)
1136 {
1137   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
1138
1139   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
1140   TopoDS_Shape theShape = aHelper->GetSubShape();
1141
1142   // a unique working file name
1143   // to avoid access to the same files by eg different users
1144   TCollection_AsciiString aGenericName
1145     = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str();
1146
1147   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
1148   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
1149   aFacesFileName  = aGenericName + ".faces";  // in faces
1150   aPointsFileName = aGenericName + ".points"; // in points
1151   aResultFileName = aGenericName + ".noboite";// out points and volumes
1152   aBadResFileName = aGenericName + ".boite";  // out bad result
1153   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
1154   aLogFileName    = aGenericName + ".log";    // log
1155
1156   // -----------------
1157   // make input files
1158   // -----------------
1159
1160   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
1161   ofstream aPointsFile  ( aPointsFileName.ToCString()  , ios::out);
1162   bool Ok =
1163 #ifdef WIN32
1164     aFacesFile->is_open() && aPointsFile->is_open();
1165 #else
1166     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
1167 #endif
1168
1169   if (!Ok)
1170     return error( SMESH_Comment("Can't write into ") << aPointsFileName);
1171
1172   vector <const SMDS_MeshNode*> aNodeByGhs3dId;
1173
1174   Ok = (writeFaces ( aFacesFile, meshDS, aNodeByGhs3dId ) &&
1175         writePoints( aPointsFile, meshDS, aNodeByGhs3dId));
1176   
1177   aFacesFile.close();
1178   aPointsFile.close();
1179   
1180   if ( ! Ok ) {
1181     if ( !_keepFiles ) {
1182       OSD_File( aFacesFileName ).Remove();
1183       OSD_File( aPointsFileName ).Remove();
1184     }
1185     return error(COMPERR_BAD_INPUT_MESH);
1186   }
1187   OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
1188
1189   // -----------------
1190   // run ghs3d mesher
1191   // -----------------
1192
1193   TCollection_AsciiString cmd =
1194     (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str();
1195   cmd += TCollection_AsciiString(" -f ") + aGenericName;  // file to read
1196   cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
1197   
1198   system( cmd.ToCString() ); // run
1199
1200   // --------------
1201   // read a result
1202   // --------------
1203   int fileOpen;
1204   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1205   if ( fileOpen < 0 ) {
1206     cout << endl;
1207     cout << "Error when opening the " << aResultFileName.ToCString() << " file" << endl;
1208     cout << "Log: " << aLogFileName << endl;
1209     cout << endl;
1210     Ok = false;
1211   }
1212   else {
1213     Ok = readResultFile( fileOpen, meshDS, theShape ,aNodeByGhs3dId );
1214   }
1215   
1216   // ---------------------
1217   // remove working files
1218   // ---------------------
1219
1220   if ( Ok )
1221   {
1222     if ( !_keepFiles )
1223       OSD_File( aLogFileName ).Remove();
1224   }
1225   else if ( OSD_File( aLogFileName ).Size() > 0 )
1226   {
1227     // get problem description from the log file
1228     _Ghs2smdsConvertor conv( aNodeByGhs3dId );
1229     storeErrorDescription( aLogFileName, conv );
1230   }
1231   else {
1232     // the log file is empty
1233     OSD_File( aLogFileName ).Remove();
1234     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
1235     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
1236   }
1237
1238   if ( !_keepFiles )
1239   {
1240     OSD_File( aFacesFileName ).Remove();
1241     OSD_File( aPointsFileName ).Remove();
1242     OSD_File( aResultFileName ).Remove();
1243     OSD_File( aBadResFileName ).Remove();
1244     OSD_File( aBbResFileName ).Remove();
1245   }
1246   
1247   return Ok;
1248 }
1249
1250 //================================================================================
1251 /*!
1252  * \brief Provide human readable text by error code reported by ghs3d
1253  */
1254 //================================================================================
1255
1256 static string translateError(const int errNum)
1257 {
1258   switch ( errNum ) {
1259   case 0:
1260     return "The surface mesh includes a face of type other than edge, "
1261       "triangle or quadrilateral. This face type is not supported.";
1262   case 1:
1263     return "Not enough memory for the face table.";
1264   case 2:
1265     return "Not enough memory.";
1266   case 3:
1267     return "Not enough memory.";
1268   case 4:
1269     return "Face is ignored.";
1270   case 5:
1271     return "End of file. Some data are missing in the file.";
1272   case 6:
1273     return "Read error on the file. There are wrong data in the file.";
1274   case 7:
1275     return "the metric file is inadequate (dimension other than 3).";
1276   case 8:
1277     return "the metric file is inadequate (values not per vertices).";
1278   case 9:
1279     return "the metric file contains more than one field.";
1280   case 10:
1281     return "the number of values in the \".bb\" (metric file) is incompatible with the expected"
1282       "value of number of mesh vertices in the \".noboite\" file.";
1283   case 12:
1284     return "Too many sub-domains.";
1285   case 13:
1286     return "the number of vertices is negative or null.";
1287   case 14:
1288     return "the number of faces is negative or null.";
1289   case 15:
1290     return "A face has a null vertex.";
1291   case 22:
1292     return "incompatible data.";
1293   case 131:
1294     return "the number of vertices is negative or null.";
1295   case 132:
1296     return "the number of vertices is negative or null (in the \".mesh\" file).";
1297   case 133:
1298     return "the number of faces is negative or null.";
1299   case 1000:
1300     return "A face appears more than once in the input surface mesh.";
1301   case 1001:
1302     return "An edge appears more than once in the input surface mesh.";
1303   case 1002:
1304     return "A face has a vertex negative or null.";
1305   case 1003:
1306     return "NOT ENOUGH MEMORY.";
1307   case 2000:
1308     return "Not enough available memory.";
1309   case 2002:
1310     return "Some initial points cannot be inserted. The surface mesh is probably very bad "
1311       "in terms of quality or the input list of points is wrong.";
1312   case 2003:
1313     return "Some vertices are too close to one another or coincident.";
1314   case 2004:
1315     return "Some vertices are too close to one another or coincident.";
1316   case 2012:
1317     return "A vertex cannot be inserted.";
1318   case 2014:
1319     return "There are at least two points considered as coincident.";
1320   case 2103:
1321     return "Some vertices are too close to one another or coincident.";
1322   case 3000:
1323     return "The surface mesh regeneration step has failed.";
1324   case 3009:
1325     return "Constrained edge cannot be enforced.";
1326   case 3019:
1327     return "Constrained face cannot be enforced.";
1328   case 3029:
1329     return "Missing faces.";
1330   case 3100:
1331     return "No guess to start the definition of the connected component(s).";
1332   case 3101:
1333     return "The surface mesh includes at least one hole. The domain is not well defined.";
1334   case 3102:
1335     return "Impossible to define a component.";
1336   case 3103:
1337     return "The surface edge intersects another surface edge.";
1338   case 3104:
1339     return "The surface edge intersects the surface face.";
1340   case 3105:
1341     return "One boundary point lies within a surface face.";
1342   case 3106:
1343     return "One surface edge intersects a surface face.";
1344   case 3107:
1345     return "One boundary point lies within a surface edge.";
1346   case 3108:
1347     return "Insufficient memory ressources detected due to a bad quality surface mesh leading "
1348       "to too many swaps.";
1349   case 3109:
1350     return "Edge is unique (i.e., bounds a hole in the surface).";
1351   case 3122:
1352     return "Presumably, the surface mesh is not compatible with the domain being processed.";
1353   case 3123:
1354     return "Too many components, too many sub-domain.";
1355   case 3209:
1356     return "The surface mesh includes at least one hole. "
1357       "Therefore there is no domain properly defined.";
1358   case 3300:
1359     return "Statistics.";
1360   case 3400:
1361     return "Statistics.";
1362   case 3500:
1363     return "Warning, it is dramatically tedious to enforce the boundary items.";
1364   case 4000:
1365     return "Not enough memory at this time, nevertheless, the program continues. "
1366       "The expected mesh will be correct but not really as large as required.";
1367   case 4002:
1368     return "see above error code, resulting quality may be poor.";
1369   case 4003:
1370     return "Not enough memory at this time, nevertheless, the program continues (warning).";
1371   case 8000:
1372     return "Unknown face type.";
1373   case 8005:
1374   case 8006:
1375     return "End of file. Some data are missing in the file.";
1376   case 9000:
1377     return "A too small volume element is detected.";
1378   case 9001:
1379     return "There exists at least a null or negative volume element.";
1380   case 9002:
1381     return "There exist null or negative volume elements.";
1382   case 9003:
1383     return "A too small volume element is detected. A face is considered being degenerated.";
1384   case 9100:
1385     return "Some element is suspected to be very bad shaped or wrong.";
1386   case 9102:
1387     return "A too bad quality face is detected. This face is considered degenerated.";
1388   case 9112:
1389     return "A too bad quality face is detected. This face is degenerated.";
1390   case 9122:
1391     return "Presumably, the surface mesh is not compatible with the domain being processed.";
1392   case 9999:
1393     return "Abnormal error occured, contact hotline.";
1394   case 23600:
1395     return "Not enough memory for the face table.";
1396   case 23601:
1397     return "The algorithm cannot run further. "
1398       "The surface mesh is probably very bad in terms of quality.";
1399   case 23602:
1400     return "Bad vertex number.";
1401   }
1402   return "";
1403 }
1404
1405 //================================================================================
1406 /*!
1407  * \brief Retrieve from a string given number of integers
1408  */
1409 //================================================================================
1410
1411 static char* getIds( char* ptr, int nbIds, vector<int>& ids )
1412 {
1413   ids.clear();
1414   ids.reserve( nbIds );
1415   while ( nbIds )
1416   {
1417     while ( !isdigit( *ptr )) ++ptr;
1418     if ( ptr[-1] == '-' ) --ptr;
1419     ids.push_back( strtol( ptr, &ptr, 10 ));
1420     --nbIds;
1421   }
1422   return ptr;
1423 }
1424
1425 //================================================================================
1426 /*!
1427  * \brief Retrieve problem description form a log file
1428  *  \retval bool - always false
1429  */
1430 //================================================================================
1431
1432 bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& logFile,
1433                                               const _Ghs2smdsConvertor &     toSmdsConvertor )
1434 {
1435   // open file
1436 #ifdef WNT
1437   int file = ::_open (logFile.ToCString(), _O_RDONLY|_O_BINARY);
1438 #else
1439   int file = ::open (logFile.ToCString(), O_RDONLY);
1440 #endif
1441   if ( file < 0 )
1442     return error( SMESH_Comment("See ") << logFile << " for problem description");
1443
1444   // get file size
1445 //   struct stat status;
1446 //   fstat(file, &status);
1447 //   size_t length = status.st_size;
1448   off_t length = lseek( file, 0, SEEK_END);
1449   lseek( file, 0, SEEK_SET);
1450
1451   // read file
1452   vector< char > buf( length );
1453   int nBytesRead = ::read (file, & buf[0], length);
1454   ::close (file);
1455   char* ptr = & buf[0];
1456   char* bufEnd = ptr + nBytesRead;
1457
1458   SMESH_Comment errDescription;
1459
1460   enum { NODE = 1, EDGE, TRIA, VOL, ID = 1 };
1461
1462   // look for errors "ERR #"
1463
1464   set<string> foundErrorStr; // to avoid reporting same error several times
1465   set<int>    elemErrorNums; // not to report different types of errors with bad elements
1466   while ( ++ptr < bufEnd )
1467   {
1468     if ( strncmp( ptr, "ERR ", 4 ) != 0 )
1469       continue;
1470
1471     list<const SMDS_MeshElement*> badElems;
1472     vector<int> nodeIds;
1473
1474     ptr += 4;
1475     char* errBeg = ptr;
1476     int   errNum = strtol(ptr, &ptr, 10);
1477     switch ( errNum ) { // we treat errors enumerated in [SALOME platform 0019316] issue
1478     case 0015:
1479       // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex.
1480       ptr = getIds(ptr, NODE, nodeIds);
1481       ptr = getIds(ptr, TRIA, nodeIds);
1482       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1483       break;
1484     case 1000: // ERR  1000 :  1 3 2
1485       // Face (f 1, f 2, f 3) appears more than once in the input surface mesh.
1486       ptr = getIds(ptr, TRIA, nodeIds);
1487       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1488       break;
1489     case 1001:
1490       // Edge (e1, e2) appears more than once in the input surface mesh
1491       ptr = getIds(ptr, EDGE, nodeIds);
1492       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1493       break;
1494     case 1002:
1495       // Face (f 1, f 2, f 3) has a vertex negative or null
1496       ptr = getIds(ptr, TRIA, nodeIds);
1497       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1498       break;
1499     case 2004:
1500       // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
1501       ptr = getIds(ptr, NODE, nodeIds);
1502       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1503       ptr = getIds(ptr, NODE, nodeIds);
1504       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1505       break;
1506     case 2012:
1507       // Vertex v1 cannot be inserted (warning).
1508       ptr = getIds(ptr, NODE, nodeIds);
1509       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1510       break;
1511     case 2014:
1512       // There are at least two points whose distance is dist, i.e., considered as coincident
1513     case 2103: // ERR  2103 :  16 WITH  3
1514       // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
1515       ptr = getIds(ptr, NODE, nodeIds);
1516       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1517       ptr = getIds(ptr, NODE, nodeIds);
1518       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1519       break;
1520     case 3009:
1521       // Constrained edge (e1, e2) cannot be enforced (warning).
1522       ptr = getIds(ptr, EDGE, nodeIds);
1523       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1524       break;
1525     case 3019:
1526       // Constrained face (f 1, f 2, f 3) cannot be enforced
1527       ptr = getIds(ptr, TRIA, nodeIds);
1528       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1529       break;
1530     case 3103: // ERR  3103 :  1 2 WITH  7 3
1531       // The surface edge (e1, e2) intersects another surface edge (e3, e4)
1532       ptr = getIds(ptr, EDGE, nodeIds);
1533       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1534       ptr = getIds(ptr, EDGE, nodeIds);
1535       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1536       break;
1537     case 3104: // ERR  3104 :  9 10 WITH  1 2 3
1538       // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3)
1539       ptr = getIds(ptr, EDGE, nodeIds);
1540       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1541       ptr = getIds(ptr, TRIA, nodeIds);
1542       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1543       break;
1544     case 3105: // ERR  3105 :  8 IN  2 3 5
1545       // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3)
1546       ptr = getIds(ptr, NODE, nodeIds);
1547       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1548       ptr = getIds(ptr, TRIA, nodeIds);
1549       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1550       break;
1551     case 3106:
1552       // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3)
1553       ptr = getIds(ptr, EDGE, nodeIds);
1554       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1555       ptr = getIds(ptr, TRIA, nodeIds);
1556       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1557       break;
1558     case 3107: // ERR  3107 :  2 IN  4 1
1559       // One boundary point (say p1) lies within a surface edge (e1, e2) (stop).
1560       ptr = getIds(ptr, NODE, nodeIds);
1561       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1562       ptr = getIds(ptr, EDGE, nodeIds);
1563       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1564       break;
1565     case 3109: // ERR  3109 :  EDGE  5 6 UNIQUE
1566       // Edge (e1, e2) is unique (i.e., bounds a hole in the surface)
1567       ptr = getIds(ptr, EDGE, nodeIds);
1568       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1569       break;
1570     case 9000: // ERR  9000 
1571       //  ELEMENT  261 WITH VERTICES :  7 396 -8 242 
1572       //  VOLUME   : -1.11325045E+11 W.R.T. EPSILON   0.
1573       // A too small volume element is detected. Are reported the index of the element,
1574       // its four vertex indices, its volume and the tolerance threshold value
1575       ptr = getIds(ptr, ID, nodeIds);
1576       ptr = getIds(ptr, VOL, nodeIds);
1577       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1578       // even if all nodes found, volume it most probably invisible,
1579       // add its faces to demenstrate it anyhow
1580       {
1581         vector<int> faceNodes( nodeIds.begin(), --nodeIds.end() ); // 012
1582         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1583         faceNodes[2] = nodeIds[3]; // 013
1584         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1585         faceNodes[1] = nodeIds[2]; // 023
1586         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1587         faceNodes[0] = nodeIds[1]; // 123
1588         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
1589       }
1590       break;
1591     case 9001: // ERR  9001
1592       //  %% NUMBER OF NEGATIVE VOLUME TETS  :  1 
1593       //  %% THE LARGEST NEGATIVE TET        :   1.75376581E+11 
1594       //  %%  NUMBER OF NULL VOLUME TETS     :  0
1595       // There exists at least a null or negative volume element
1596       break;
1597     case 9002:
1598       // There exist n null or negative volume elements
1599       break;
1600     case 9003:
1601       // A too small volume element is detected
1602       break;
1603     case 9102:
1604       // A too bad quality face is detected. This face is considered degenerated,
1605       // its index, its three vertex indices together with its quality value are reported
1606       break; // same as next
1607     case 9112: // ERR  9112 
1608       //  FACE   2 WITH VERTICES :  4 2 5 
1609       //  SMALL INRADIUS :   0.
1610       // A too bad quality face is detected. This face is degenerated,
1611       // its index, its three vertex indices together with its inradius are reported
1612       ptr = getIds(ptr, ID, nodeIds);
1613       ptr = getIds(ptr, TRIA, nodeIds);
1614       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
1615       // add triangle edges as it most probably has zero area and hence invisible
1616       {
1617         vector<int> edgeNodes(2);
1618         edgeNodes[0] = nodeIds[0]; edgeNodes[1] = nodeIds[1]; // 0-1
1619         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1620         edgeNodes[1] = nodeIds[2]; // 0-2
1621         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1622         edgeNodes[0] = nodeIds[1]; // 1-2
1623         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
1624       }
1625       break;
1626     }
1627
1628     bool isNewError = foundErrorStr.insert( string( errBeg, ptr )).second;
1629     if ( !isNewError )
1630       continue; // not to report same error several times
1631
1632 //     const SMDS_MeshElement* nullElem = 0;
1633 //     bool allElemsOk = ( find( badElems.begin(), badElems.end(), nullElem) == badElems.end());
1634
1635 //     if ( allElemsOk && !badElems.empty() && !elemErrorNums.empty() ) {
1636 //       bool oneMoreErrorType = elemErrorNums.insert( errNum ).second;
1637 //       if ( oneMoreErrorType )
1638 //         continue; // not to report different types of errors with bad elements
1639 //     }
1640
1641     // store bad elements
1642     //if ( allElemsOk ) {
1643       list<const SMDS_MeshElement*>::iterator elem = badElems.begin();
1644       for ( ; elem != badElems.end(); ++elem )
1645         addBadInputElement( *elem );
1646       //}
1647
1648     // make error text
1649     string text = translateError( errNum );
1650     if ( errDescription.find( text ) == text.npos ) {
1651       if ( !errDescription.empty() )
1652         errDescription << "\n";
1653       errDescription << text;
1654     }
1655
1656   } // end while
1657
1658   if ( errDescription.empty() ) { // no errors found
1659     char msg[] = "connection to server failed";
1660     if ( search( &buf[0], bufEnd, msg, msg + strlen(msg)) != bufEnd )
1661       errDescription << "Licence problems.";
1662   }
1663
1664   if ( errDescription.empty() )
1665     errDescription << "See " << logFile << " for problem description";
1666   else
1667     errDescription << "\nSee " << logFile << " for more information";
1668
1669   return error( errDescription );
1670 }
1671
1672 //================================================================================
1673 /*!
1674  * \brief Creates _Ghs2smdsConvertor
1675  */
1676 //================================================================================
1677
1678 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const map <int,const SMDS_MeshNode*> & ghs2NodeMap)
1679   :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 )
1680 {
1681 }
1682
1683 //================================================================================
1684 /*!
1685  * \brief Creates _Ghs2smdsConvertor
1686  */
1687 //================================================================================
1688
1689 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector <const SMDS_MeshNode*> &  nodeByGhsId)
1690   : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId )
1691 {
1692 }
1693
1694 //================================================================================
1695 /*!
1696  * \brief Return SMDS element by ids of GHS3D nodes
1697  */
1698 //================================================================================
1699
1700 const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector<int>& ghsNodes) const
1701 {
1702   size_t nbNodes = ghsNodes.size();
1703   vector<const SMDS_MeshNode*> nodes( nbNodes, 0 );
1704   for ( size_t i = 0; i < nbNodes; ++i ) {
1705     int ghsNode = ghsNodes[ i ];
1706     if ( _ghs2NodeMap ) {
1707       map <int,const SMDS_MeshNode*>::const_iterator in = _ghs2NodeMap->find( ghsNode);
1708       if ( in == _ghs2NodeMap->end() )
1709         return 0;
1710       nodes[ i ] = in->second;
1711     }
1712     else {
1713       if ( ghsNode < 1 || ghsNode > _nodeByGhsId->size() )
1714         return 0;
1715       nodes[ i ] = (*_nodeByGhsId)[ ghsNode-1 ];
1716     }
1717   }
1718   if ( nbNodes == 1 )
1719     return nodes[0];
1720
1721   if ( nbNodes == 2 ) {
1722     const SMDS_MeshElement* edge= SMDS_Mesh::FindEdge( nodes[0], nodes[1] );
1723     if ( !edge )
1724       edge = new SMDS_MeshEdge( nodes[0], nodes[1] );
1725     return edge;
1726   }
1727   if ( nbNodes == 3 ) {
1728     const SMDS_MeshElement* face = SMDS_Mesh::FindFace( nodes );
1729     if ( !face )
1730       face = new SMDS_FaceOfNodes( nodes[0], nodes[1], nodes[2] );
1731     return face;
1732   }
1733   if ( nbNodes == 4 )
1734     return new SMDS_VolumeOfNodes( nodes[0], nodes[1], nodes[2], nodes[3] );
1735
1736   return 0;
1737 }
1738
1739
1740 //=============================================================================
1741 /*!
1742  *  
1743  */
1744 //=============================================================================
1745 bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
1746                                  const TopoDS_Shape& aShape,
1747                                  MapShapeNbElems& aResMap)
1748 {
1749   int nbtri = 0, nbqua = 0;
1750   double fullArea = 0.0;
1751   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
1752     TopoDS_Face F = TopoDS::Face( exp.Current() );
1753     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1754     MapShapeNbElemsItr anIt = aResMap.find(sm);
1755     std::vector<int> aVec = (*anIt).second;
1756     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1757     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1758     GProp_GProps G;
1759     BRepGProp::SurfaceProperties(F,G);
1760     double anArea = G.Mass();
1761     fullArea += anArea;
1762   }
1763
1764   // collect info from edges
1765   int nb0d_e = 0, nb1d_e = 0;
1766   bool IsQuadratic = false;
1767   bool IsFirst = true;
1768   TopTools_MapOfShape tmpMap;
1769   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
1770     TopoDS_Edge E = TopoDS::Edge(exp.Current());
1771     if( tmpMap.Contains(E) )
1772       continue;
1773     tmpMap.Add(E);
1774     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
1775     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
1776     std::vector<int> aVec = (*anIt).second;
1777     nb0d_e += aVec[SMDSEntity_Node];
1778     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1779     if(IsFirst) {
1780       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1781       IsFirst = false;
1782     }
1783   }
1784   tmpMap.Clear();
1785
1786   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
1787
1788   GProp_GProps G;
1789   BRepGProp::VolumeProperties(aShape,G);
1790   double aVolume = G.Mass();
1791   double tetrVol = 0.1179*ELen*ELen*ELen;
1792   double CoeffQuality = 0.9;
1793   int nbVols = (int)aVolume/tetrVol/CoeffQuality;
1794   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
1795   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
1796   std::vector<int> aVec(SMDSEntity_Last);
1797   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1798   if( IsQuadratic ) {
1799     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
1800     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
1801     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
1802   }
1803   else {
1804     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
1805     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
1806     aVec[SMDSEntity_Pyramid] = nbqua;
1807   }
1808   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
1809   aResMap.insert(std::make_pair(sm,aVec));
1810
1811   return true;
1812 }