Salome HOME
Porting HexoticPLUGIN SALOME module to the CMake build system: initial version.
[plugins/hexoticplugin.git] / src / HexoticPlugin / HexoticPlugin_Hexotic.cxx
1 // Copyright (C) 2007-2013  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 // ---
21 // File   : HexoticPlugin_Hexotic.cxx
22 // Author : Lioka RAZAFINDRAZAKA (CEA)
23 // ---
24 //
25 #include "HexoticPlugin_Hexotic.hxx"
26 #include "HexoticPlugin_Hypothesis.hxx"
27
28 #include "utilities.h"
29
30 #ifndef WIN32
31 #include <sys/sysinfo.h>
32 #else
33 #include <errno.h>
34 #endif
35
36 #ifdef _DEBUG_
37 #define DUMP(txt) \
38 //  cout << txt
39 #else
40 #define DUMP(txt)
41 #endif
42
43 #include <SMESHDS_Mesh.hxx>
44 #include <SMESHDS_GroupBase.hxx>
45 #include <SMESH_ComputeError.hxx>
46 #include <SMESH_File.hxx>
47 #include <SMESH_Gen.hxx>
48 #include <SMESH_HypoFilter.hxx>
49 #include <SMESH_MesherHelper.hxx>
50 #include <SMESH_subMesh.hxx>
51
52 #include <list>
53 #include <cstdlib>
54 #include <iostream>
55
56 #include <Standard_ProgramError.hxx>
57
58 #include <BRep_Tool.hxx>
59 #include <BRepBndLib.hxx>
60 #include <BRepClass3d_SolidClassifier.hxx>
61 #include <BRepBuilderAPI_MakeVertex.hxx>
62 #include <BRepExtrema_DistShapeShape.hxx>
63 #include <OSD_File.hxx>
64 #include <Precision.hxx>
65 #include <TopExp_Explorer.hxx>
66 #include <TopTools_MapOfShape.hxx>
67 #include <TopoDS.hxx>
68 #include <gp_Pnt.hxx>
69 #include <gp_Lin.hxx>
70 #include <GCPnts_UniformAbscissa.hxx>
71 #include <IntCurvesFace_Intersector.hxx>
72 #include <BRepMesh_IncrementalMesh.hxx>
73 #include <Poly_Triangulation.hxx>
74 #include <gp_Ax3.hxx>
75 #include <GeomAdaptor_Curve.hxx>
76
77 #include <Basics_Utils.hxx>
78 #include "GEOMImpl_Types.hxx"
79 #include "GEOM_wrap.hxx"
80
81 static void removeFile( const TCollection_AsciiString& fileName )
82 {
83   try {
84     OSD_File( fileName ).Remove();
85   }
86   catch ( Standard_ProgramError ) {
87     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
88   }
89 }
90
91 //=============================================================================
92 /*!
93  *  
94  */
95 //=============================================================================
96
97 HexoticPlugin_Hexotic::HexoticPlugin_Hexotic(int hypId, int studyId, SMESH_Gen* gen)
98   : SMESH_3D_Algo(hypId, studyId, gen)
99 {
100   MESSAGE("HexoticPlugin_Hexotic::HexoticPlugin_Hexotic");
101   _name = "Hexotic_3D";
102   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
103 //   _onlyUnaryInput = false;
104   _requireShape = false;
105   _iShape=0;
106   _nbShape=0;
107   _hexoticFilesKept=false;
108   _compatibleHypothesis.push_back("Hexotic_Parameters");
109 #ifdef WITH_BLSURFPLUGIN
110   _blsurfHypo = NULL;
111 #endif
112   _compute_canceled = false;
113   
114   // Copy of what is done in BLSURFPLugin TODO : share the code
115   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
116   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
117   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
118   
119   myStudy = NULL;
120   myStudy = aStudyMgr->GetStudyByID(_studyId);
121   if (myStudy)
122     MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
123 }
124
125 //=============================================================================
126 /*!
127  *  
128  */
129 //=============================================================================
130
131 HexoticPlugin_Hexotic::~HexoticPlugin_Hexotic()
132 {
133   MESSAGE("HexoticPlugin_Hexotic::~HexoticPlugin_Hexotic");
134 }
135
136
137 #ifdef WITH_BLSURFPLUGIN
138 bool HexoticPlugin_Hexotic::CheckBLSURFHypothesis( SMESH_Mesh&         aMesh,
139                                                    const TopoDS_Shape& aShape )
140 {
141   // MESSAGE("HexoticPlugin_Hexotic::CheckBLSURFHypothesis");
142   _blsurfHypo = NULL;
143
144   std::list<const SMESHDS_Hypothesis*>::const_iterator itl;
145   const SMESHDS_Hypothesis* theHyp;
146
147   // If a BLSURF hypothesis is applied, get it
148   SMESH_HypoFilter blsurfFilter;
149   blsurfFilter.Init( blsurfFilter.HasName( "BLSURF_Parameters" ));
150   std::list<const SMESHDS_Hypothesis *> appliedHyps;
151   aMesh.GetHypotheses( aShape, blsurfFilter, appliedHyps, false );
152
153   if ( appliedHyps.size() > 0 ) {
154     itl = appliedHyps.begin();
155     theHyp = (*itl); // use only the first hypothesis
156     std::string hypName = theHyp->GetName();
157     if (hypName == "BLSURF_Parameters") {
158       _blsurfHypo = static_cast<const BLSURFPlugin_Hypothesis*> (theHyp);
159       ASSERT(_blsurfHypo);
160       return true;
161     }
162   }
163   return false;
164 }
165 #endif
166
167 //=============================================================================
168 /*!
169  *  
170  */
171 //=============================================================================
172
173 bool HexoticPlugin_Hexotic::CheckHypothesis( SMESH_Mesh&                          aMesh,
174                                              const TopoDS_Shape&                  aShape,
175                                              SMESH_Hypothesis::Hypothesis_Status& aStatus )
176 {
177   // MESSAGE("HexoticPlugin_Hexotic::CheckHypothesis");
178   _hypothesis = NULL;
179
180   std::list<const SMESHDS_Hypothesis*>::const_iterator itl;
181   const SMESHDS_Hypothesis* theHyp;
182
183   const std::list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, false);
184   int nbHyp = hyps.size();
185   if (!nbHyp) {
186     aStatus = SMESH_Hypothesis::HYP_OK;
187     return true;  // can work with no hypothesis
188   }
189
190   itl = hyps.begin();
191   theHyp = (*itl); // use only the first hypothesis
192
193   std::string hypName = theHyp->GetName();
194   if (hypName == "Hexotic_Parameters") {
195     _hypothesis = static_cast<const HexoticPlugin_Hypothesis*> (theHyp);
196     ASSERT(_hypothesis);
197     aStatus = SMESH_Hypothesis::HYP_OK;
198   }
199   else
200     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
201   
202 #ifdef WITH_BLSURFPLUGIN
203   CheckBLSURFHypothesis(aMesh, aShape);
204 #endif
205   
206   return aStatus == SMESH_Hypothesis::HYP_OK;
207 }
208
209 //=======================================================================
210 //function : findShape
211 //purpose  :
212 //=======================================================================
213
214 static TopoDS_Shape findShape(SMDS_MeshNode**     t_Node,
215                               TopoDS_Shape        aShape,
216                               const TopoDS_Shape* t_Shape,
217                               double**            t_Box,
218                               const int           nShape)
219 {
220   double pntCoor[3];
221   int iShape, nbNode = 8;
222
223   for ( int i=0; i<3; i++ ) {
224     pntCoor[i] = 0;
225     for ( int j=0; j<nbNode; j++ ) {
226       if ( i == 0) pntCoor[i] += t_Node[j]->X();
227       if ( i == 1) pntCoor[i] += t_Node[j]->Y();
228       if ( i == 2) pntCoor[i] += t_Node[j]->Z();
229     }
230     pntCoor[i] /= nbNode;
231   }
232   gp_Pnt aPnt(pntCoor[0], pntCoor[1], pntCoor[2]);
233
234   if ( aShape.IsNull() ) aShape = t_Shape[0];
235   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
236   if ( !(SC.State() == TopAbs_IN) ) {
237     aShape.Nullify();
238     for (iShape = 0; iShape < nShape && aShape.IsNull(); iShape++) {
239       if ( !( pntCoor[0] < t_Box[iShape][0] || t_Box[iShape][1] < pntCoor[0] ||
240               pntCoor[1] < t_Box[iShape][2] || t_Box[iShape][3] < pntCoor[1] ||
241               pntCoor[2] < t_Box[iShape][4] || t_Box[iShape][5] < pntCoor[2]) ) {
242         BRepClass3d_SolidClassifier SC (t_Shape[iShape], aPnt, Precision::Confusion());
243         if (SC.State() == TopAbs_IN)
244           aShape = t_Shape[iShape];
245       }
246     }
247   }
248   return aShape;
249 }
250
251 //=======================================================================
252 //function : findEdge
253 //purpose  :
254 //=======================================================================
255
256 static int findEdge(const SMDS_MeshNode* aNode,
257                     const SMESHDS_Mesh*  aMesh,
258                     const int            nEdge,
259                     const TopoDS_Shape*  t_Edge) {
260
261   TopoDS_Shape aPntShape, foundEdge;
262   TopoDS_Vertex aVertex;
263   gp_Pnt aPnt( aNode->X(), aNode->Y(), aNode->Z() );
264
265   int foundInd, ind;
266   double nearest = RealLast(), *t_Dist;
267   double epsilon = Precision::Confusion();
268
269   t_Dist = new double[ nEdge ];
270   aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape();
271   aVertex   = TopoDS::Vertex( aPntShape );
272
273   for ( ind=0; ind < nEdge; ind++ ) {
274     BRepExtrema_DistShapeShape aDistance ( aVertex, t_Edge[ind] );
275     t_Dist[ind] = aDistance.Value();
276     if ( t_Dist[ind] < nearest ) {
277       nearest   = t_Dist[ind];
278       foundEdge = t_Edge[ind];
279       foundInd  = ind;
280       if ( nearest < epsilon )
281         ind = nEdge;
282     }
283   }
284
285   delete [] t_Dist;
286   return aMesh->ShapeToIndex( foundEdge );
287 }
288
289 //=======================================================================
290 //function : getNbShape
291 //purpose  :
292 //=======================================================================
293
294 static int getNbShape(std::string aFile, std::string aString, int defaultValue=0) {
295   int number = defaultValue;
296   std::string aLine;
297   std::ifstream file(aFile.c_str());
298   while ( !file.eof() ) {
299     getline( file, aLine);
300     if ( aLine == aString ) {
301       getline( file, aLine);
302       std::istringstream stringFlux( aLine );
303       stringFlux >> number;
304       number = ( number + defaultValue + std::abs(number - defaultValue) ) / 2;
305       break;
306     }
307   }
308   file.close();
309   return number;
310 }
311
312 //=======================================================================
313 //function : countShape
314 //purpose  :
315 //=======================================================================
316
317 template < class Mesh, class Shape >
318 static int countShape( Mesh* mesh, Shape shape ) {
319   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
320   TopTools_MapOfShape mapShape;
321   int nbShape = 0;
322   for ( ; expShape.More(); expShape.Next() ) {
323     if (mapShape.Add(expShape.Current())) {
324       nbShape++;
325     }
326   }
327   return nbShape;
328 }
329
330 //=======================================================================
331 //function : getShape
332 //purpose  :
333 //=======================================================================
334
335 template < class Mesh, class Shape, class Tab >
336 void getShape(Mesh* mesh, Shape shape, Tab *t_Shape) {
337   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
338   TopTools_MapOfShape mapShape;
339   for ( int i=0; expShape.More(); expShape.Next() ) {
340     if (mapShape.Add(expShape.Current())) {
341       t_Shape[i] = expShape.Current();
342       i++;
343     }
344   }
345   return;
346 }
347
348 //=======================================================================
349 //function : printWarning
350 //purpose  :
351 //=======================================================================
352
353 static void printWarning(const int nbExpected, std::string aString, const int nbFound) {
354   cout << std::endl;
355   cout << "WARNING : " << nbExpected << " " << aString << " expected, Hexotic has found " << nbFound << std::endl;
356   cout << "=======" << std::endl;
357   cout << std::endl;
358   return;
359 }
360
361 //=======================================================================
362 //function : removeHexoticFiles
363 //purpose  :
364 //=======================================================================
365
366 static void removeHexoticFiles(TCollection_AsciiString file_In, TCollection_AsciiString file_Out) {
367   removeFile( file_In );
368   removeFile( file_Out );
369 }
370
371 //=======================================================================
372 //function : readResult
373 //purpose  : Read GMF file in case of a mesh with geometry
374 //=======================================================================
375
376 static bool readResult(std::string         theFile,
377                        HexoticPlugin_Hexotic*  theAlgo,
378                        SMESHDS_Mesh*       theMesh,
379                        const int           nbShape,
380                        const TopoDS_Shape* tabShape,
381                        double**            tabBox)
382 {
383   // ---------------------------------
384   // Optimisation of the plugin ...
385   // Retrieve the correspondance edge --> shape
386   // (which is very costly) only when user
387   // has defined at least one group of edges
388   // which should be rare for a 3d mesh !
389   // ---------------------------------
390   
391   bool retrieve_edges = false;
392   const std::set<SMESHDS_GroupBase*>& aGroups = theMesh->GetGroups();
393   set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
394   for (; GrIt != aGroups.end(); GrIt++)
395     {
396       SMESHDS_GroupBase* aGrp = *GrIt;
397       if ( !aGrp )
398         continue;
399       if ( aGrp->GetType() == SMDSAbs_Edge )
400         {
401           retrieve_edges = true;
402           break;
403         }
404     }
405   
406   // ---------------------------------
407   // Read generated elements and nodes
408   // ---------------------------------
409
410   TopoDS_Shape aShape;
411   TopoDS_Vertex aVertex;
412   std::string token;
413   int EndOfFile = 0, nbElem = 0, nField = 9, nbRef = 0;
414   int aHexoticNodeID = 0, shapeID, hexoticShapeID;
415   const int IdShapeRef = 2;
416   int *tabID, *tabRef, *nodeAssigne;
417   bool *tabDummy, hasDummy = false;
418   double epsilon = Precision::Confusion();
419   std::map <std::string,int> mapField;
420   SMDS_MeshNode** HexoticNode;
421   TopoDS_Shape *tabCorner, *tabEdge;
422
423   const int nbDomains = countShape( theMesh, TopAbs_SHELL );
424   const int holeID = -1;
425
426   // tabID    = new int[nbShape];
427   tabID    = new int[nbDomains];
428   tabRef   = new int[nField];
429   tabDummy = new bool[nField];
430
431   for (int i=0; i<nbDomains; i++)
432     tabID[i] = 0;
433   if ( nbDomains == 1 )
434     tabID[0] = theMesh->ShapeToIndex( tabShape[0] );
435
436   mapField["MeshVersionFormatted"] = 0; tabRef[0] = 0; tabDummy[0] = false;
437   mapField["Dimension"]            = 1; tabRef[1] = 0; tabDummy[1] = false;
438   mapField["Vertices"]             = 2; tabRef[2] = 3; tabDummy[2] = true;
439   mapField["Corners"]              = 3; tabRef[3] = 1; tabDummy[3] = false;
440   mapField["Edges"]                = 4; tabRef[4] = 2; tabDummy[4] = true;
441   mapField["Ridges"]               = 5; tabRef[5] = 1; tabDummy[5] = false;
442   mapField["Quadrilaterals"]       = 6; tabRef[6] = 4; tabDummy[6] = true;
443   mapField["Hexahedra"]            = 7; tabRef[7] = 8; tabDummy[7] = true;
444   mapField["End"]                  = 8; tabRef[8] = 0; tabDummy[0] = false;
445
446   SMDS_NodeIteratorPtr itOnHexoticInputNode = theMesh->nodesIterator();
447   while ( itOnHexoticInputNode->more() )
448     theMesh->RemoveNode( itOnHexoticInputNode->next() );
449
450   int nbVertices   = getNbShape(theFile, "Vertices");
451   int nbCorners    = getNbShape(theFile, "Corners", countShape( theMesh, TopAbs_VERTEX ));
452   int nbShapeEdge  = countShape( theMesh, TopAbs_EDGE );
453
454   tabCorner   = new TopoDS_Shape[ nbCorners ];
455   tabEdge     = new TopoDS_Shape[ nbShapeEdge ];
456   nodeAssigne = new int[ nbVertices + 1 ];
457   HexoticNode = new SMDS_MeshNode*[ nbVertices + 1 ];
458
459   getShape(theMesh, TopAbs_VERTEX, tabCorner);
460   getShape(theMesh, TopAbs_EDGE,   tabEdge);
461
462   MESSAGE("Read " << theFile << " file");
463   std::ifstream fileRes(theFile.c_str());
464   ASSERT(fileRes);
465
466   while ( EndOfFile == 0  ) {
467     int dummy;
468     fileRes >> token;
469
470     if (mapField.count(token)) {
471       nField   = mapField[token];
472       nbRef    = tabRef[nField];
473       hasDummy = tabDummy[nField];
474     }
475     else {
476       nField = -1;
477       nbRef = 0;
478     }
479
480     nbElem = 0;
481     if ( nField < (mapField.size() - 1) && nField >= 0 )
482       fileRes >> nbElem;
483
484     switch (nField) {
485       case 0: { // "MeshVersionFormatted"
486         MESSAGE(token << " " << nbElem);
487         break;
488       }
489       case 1: { // "Dimension"
490         MESSAGE("Mesh dimension " << nbElem << "D");
491         break;
492       }
493       case 2: { // "Vertices"
494         MESSAGE("Read " << nbElem << " " << token);
495         int aHexoticID;
496         double *coord;
497         SMDS_MeshNode * aHexoticNode;
498
499         coord = new double[nbRef];
500         for ( int iElem = 0; iElem < nbElem; iElem++ ) {
501           if(theAlgo->computeCanceled())
502             {
503               return false;
504             }
505           aHexoticID = iElem + 1;
506           for ( int iCoord = 0; iCoord < 3; iCoord++ )
507             fileRes >> coord[ iCoord ];
508           fileRes >> dummy;
509           aHexoticNode = theMesh->AddNode(coord[0], coord[1], coord[2]);
510           HexoticNode[ aHexoticID ] = aHexoticNode;
511           nodeAssigne[ aHexoticID ] = 0;
512         }
513         delete [] coord;
514         break;
515       }
516       case 3: // "Corners"
517       case 4: // "Edges"
518       case 5: // "Ridges"
519       case 6: // "Quadrilaterals"
520       case 7: { // "Hexahedra"
521         MESSAGE("Read " << nbElem << " " << token);
522         SMDS_MeshNode** node;
523         int nodeDim, *nodeID;
524         SMDS_MeshElement * aHexoticElement = 0;
525
526         node   = new SMDS_MeshNode*[ nbRef ];
527         nodeID = new int[ nbRef ];
528         for ( int iElem = 0; iElem < nbElem; iElem++ ) {
529           if(theAlgo->computeCanceled())
530             {
531               return false;
532             }
533           for ( int iRef = 0; iRef < nbRef; iRef++ ) {
534             fileRes >> aHexoticNodeID;                          // read nbRef aHexoticNodeID
535             node[ iRef ]   = HexoticNode[ aHexoticNodeID ];
536             nodeID[ iRef ] = aHexoticNodeID;
537           }
538           if ( hasDummy )
539             fileRes >> dummy;
540           switch (nField) {
541             case 3: { // "Corners"
542               nodeDim = 1;
543               gp_Pnt HexoticPnt ( node[0]->X(), node[0]->Y(), node[0]->Z() );
544               for ( int i=0; i<nbElem; i++ ) {
545                 aVertex = TopoDS::Vertex( tabCorner[i] );
546                 gp_Pnt aPnt = BRep_Tool::Pnt( aVertex );
547                 if ( aPnt.Distance( HexoticPnt ) < epsilon )
548                   break;
549               }
550               break;
551             }
552             case 4: { // "Edges"
553               nodeDim = 2;
554               aHexoticElement = theMesh->AddEdge( node[0], node[1] );
555               int iNode = 1;
556               if ( nodeAssigne[ nodeID[0] ] == 0 || nodeAssigne[ nodeID[0] ] == 2 )
557                 iNode = 0;
558               if(retrieve_edges)
559                 shapeID = findEdge( node[iNode], theMesh, nbShapeEdge, tabEdge );
560               else
561                 shapeID = 0;
562               break;
563             }
564             case 5: { // "Ridges"
565               break;
566             }
567             case 6: { // "Quadrilaterals"
568               nodeDim = 3;
569               aHexoticElement = theMesh->AddFace( node[0], node[1], node[2], node[3] );
570               shapeID = dummy;
571               break;
572             }
573             case 7: { // "Hexahedra"
574               nodeDim = 4;
575               if ( nbDomains > 1 ) {
576                 hexoticShapeID = dummy - IdShapeRef;
577                 if ( tabID[ hexoticShapeID ] == 0 ) {
578                   aShape = findShape(node, aShape, tabShape, tabBox, nbShape);
579                   shapeID = aShape.IsNull() ? holeID : theMesh->ShapeToIndex( aShape );
580                   tabID[ hexoticShapeID ] = shapeID;
581                 }
582                 else
583                   shapeID = tabID[ hexoticShapeID ];
584                 if ( iElem == (nbElem - 1) ) {
585                   int shapeAssociated = 0;
586                   for ( int i=0; i<nbDomains; i++ ) {
587                     if (tabID[i] > 0 )
588                       shapeAssociated += 1;
589                   }
590                   if ( shapeAssociated != nbShape )
591                     printWarning(nbShape, "domains", shapeAssociated);
592                 }
593               }
594               else {
595                 shapeID = tabID[0];
596               }
597               if ( shapeID != holeID )
598                 aHexoticElement = theMesh->AddVolume( node[0], node[3], node[2], node[1], node[4], node[7], node[6], node[5] );
599               break;
600             }
601           } // switch (nField)
602
603           if ( token != "Ridges" && ( shapeID > 0 || token == "Corners")) {
604             for ( int i=0; i<nbRef; i++ ) {
605               if ( nodeAssigne[ nodeID[i] ] == 0 ) {
606                 if      ( token == "Corners" )        theMesh->SetNodeOnVertex( node[0], aVertex );
607                 else if ( token == "Edges" )          theMesh->SetNodeOnEdge( node[i], shapeID );
608                 else if ( token == "Quadrilaterals" ) theMesh->SetNodeOnFace( node[i], shapeID );
609                 else if ( token == "Hexahedra" )      theMesh->SetNodeInVolume( node[i], shapeID );
610                 nodeAssigne[ nodeID[i] ] = nodeDim;
611               }
612             }
613             if ( token != "Corners" && aHexoticElement )
614               theMesh->SetMeshElementOnShape( aHexoticElement, shapeID );
615           }
616         }
617         delete [] node;
618         delete [] nodeID;
619         break;
620       }
621       case 8: { // "End"
622         EndOfFile = 1;
623         MESSAGE("End of " << theFile << " file");
624         break;
625       }
626       default: {
627         MESSAGE("Unknown Token: " << token);
628       }
629     }
630   }
631   cout << std::endl;
632
633   // remove nodes in holes
634   if ( nbDomains > 1 )
635   {
636     SMESHDS_SubMesh* subMesh;
637     for ( int i = 1; i <= nbVertices; ++i )
638       if ( HexoticNode[i]->NbInverseElements() == 0 )
639       {
640         subMesh =  HexoticNode[i]->getshapeId() > 0 ? theMesh->MeshElements(HexoticNode[i]->getshapeId() ) : 0;
641         theMesh->RemoveFreeNode( HexoticNode[i], subMesh, /*fromGroups=*/false );
642       }
643   }
644   delete [] tabID;
645   delete [] tabRef;
646   delete [] tabDummy;
647   delete [] tabCorner;
648   delete [] tabEdge;
649   delete [] nodeAssigne;
650   delete [] HexoticNode;
651   return true;
652 }
653
654
655 //=======================================================================
656 //function : readResult
657 //purpose  : Read GMF file in case of a mesh w/o geometry
658 //=======================================================================
659
660 static bool readResult(std::string theFile,
661                        HexoticPlugin_Hexotic*  theAlgo,
662                        SMESH_MesherHelper* theHelper)
663 {
664   SMESHDS_Mesh* theMesh = theHelper->GetMeshDS();
665
666   // ---------------------------------
667   // Read generated elements and nodes
668   // ---------------------------------
669
670   std::string token;
671   const int nbField = 9;
672   int nField, EndOfFile = 0, nbElem = 0, nbRef = 0;
673   int aHexoticNodeID = 0, shapeID;
674   int tabRef[nbField], *nodeAssigne;
675   bool tabDummy[nbField], hasDummy = false;
676   std::map <std::string,int> mapField;
677   SMDS_MeshNode** HexoticNode;
678
679   mapField["MeshVersionFormatted"] = 0; tabRef[0] = 0; tabDummy[0] = false;
680   mapField["Dimension"]            = 1; tabRef[1] = 0; tabDummy[1] = false;
681   mapField["Vertices"]             = 2; tabRef[2] = 3; tabDummy[2] = true;
682   mapField["Corners"]              = 3; tabRef[3] = 1; tabDummy[3] = false;
683   mapField["Edges"]                = 4; tabRef[4] = 2; tabDummy[4] = true;
684   mapField["Ridges"]               = 5; tabRef[5] = 1; tabDummy[5] = false;
685   mapField["Quadrilaterals"]       = 6; tabRef[6] = 4; tabDummy[6] = true;
686   mapField["Hexahedra"]            = 7; tabRef[7] = 8; tabDummy[7] = true;
687   mapField["End"]                  = 8; tabRef[8] = 0; tabDummy[8] = false;
688
689   theHelper->GetMesh()->Clear();
690
691   int nbVertices = getNbShape(theFile, "Vertices");
692   HexoticNode = new SMDS_MeshNode*[ nbVertices + 1 ];
693   nodeAssigne = new int[ nbVertices + 1 ];
694
695   MESSAGE("Read " << theFile << " file");
696   std::ifstream fileRes(theFile.c_str());
697   ASSERT(fileRes);
698
699   while ( !EndOfFile  )
700   {
701     int dummy;
702     fileRes >> token;
703
704     if (mapField.count(token)) {
705       nField   = mapField[token];
706       nbRef    = tabRef[nField];
707       hasDummy = tabDummy[nField];
708     }
709     else {
710       nField = -1;
711       nbRef = 0;
712     }
713
714     nbElem = 0;
715     if ( nField < (mapField.size() - 1) && nField >= 0 )
716       fileRes >> nbElem;
717
718     switch (nField) {
719     case 0: { // "MeshVersionFormatted"
720       MESSAGE(token << " " << nbElem);
721       break;
722     }
723     case 1: { // "Dimension"
724       MESSAGE("Mesh dimension " << nbElem << "D");
725       break;
726     }
727     case 2: { // "Vertices"
728       MESSAGE("Read " << nbElem << " " << token);
729       int aHexoticID;
730       double coord[3];
731       SMDS_MeshNode * aHexoticNode;
732
733       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
734         if(theAlgo->computeCanceled())
735           {
736             return false;
737           }
738         aHexoticID = iElem + 1;
739         for ( int iCoord = 0; iCoord < 3; iCoord++ )
740           fileRes >> coord[ iCoord ];
741         fileRes >> dummy;
742         aHexoticNode = theMesh->AddNode(coord[0], coord[1], coord[2]);
743         HexoticNode[ aHexoticID ] = aHexoticNode;
744         nodeAssigne[ aHexoticID ] = 0;
745       }
746       break;
747     }
748     case 3: // "Corners"
749     case 4: // "Edges"
750     case 5: // "Ridges"
751     case 6: // "Quadrilaterals"
752     case 7: { // "Hexahedra"
753       MESSAGE("Read " << nbElem << " " << token);
754       std::vector< SMDS_MeshNode* > node( nbRef );
755       std::vector< int >          nodeID( nbRef );
756
757       for ( int iElem = 0; iElem < nbElem; iElem++ )
758       {
759         if(theAlgo->computeCanceled())
760           {
761             return false;
762           }
763         for ( int iRef = 0; iRef < nbRef; iRef++ )
764         {
765           fileRes >> aHexoticNodeID;                          // read nbRef aHexoticNodeID
766           node  [ iRef ] = HexoticNode[ aHexoticNodeID ];
767           nodeID[ iRef ] = aHexoticNodeID;
768         }
769         if ( hasDummy )
770           fileRes >> dummy;
771         switch (nField)
772         {
773         case 4: // "Edges"
774           theHelper->AddEdge( node[0], node[1] ); break;
775         case 6:  // "Quadrilaterals"
776           theMesh->AddFace( node[0], node[1], node[2], node[3] ); break;
777         case 7: // "Hexahedra"
778           theHelper->AddVolume( node[0], node[3], node[2], node[1],
779                                 node[4], node[7], node[6], node[5] ); break;
780         default: continue;
781         }
782         if ( nField == 6 )
783           for ( int iRef = 0; iRef < nbRef; iRef++ )
784             nodeAssigne[ nodeID[ iRef ]] = 1;
785       }
786       break;
787     }
788     case 8: { // "End"
789       EndOfFile = 1;
790       MESSAGE("End of " << theFile << " file");
791       break;
792     }
793     default: {
794       MESSAGE("Unknown Token: " << token);
795     }
796     }
797   }
798   cout << std::endl;
799
800   shapeID = theHelper->GetSubShapeID();
801   for ( int i = 0; i < nbVertices; ++i )
802     if ( !nodeAssigne[ i+1 ])
803       theMesh->SetNodeInVolume( HexoticNode[ i+1 ], shapeID );
804
805   delete [] HexoticNode;
806   delete [] nodeAssigne;
807   return true;
808 }
809
810 //=============================================================================
811 /*!
812  * Pass parameters to Hexotic
813  */
814 //=============================================================================
815
816 void HexoticPlugin_Hexotic::SetParameters(const HexoticPlugin_Hypothesis* hyp) {
817
818   MESSAGE("HexoticPlugin_Hexotic::SetParameters");
819   if (hyp) {
820     _hexesMinLevel = hyp->GetHexesMinLevel();
821     _hexesMaxLevel = hyp->GetHexesMaxLevel();
822     _hexesMinSize = hyp->GetMinSize();
823     _hexesMaxSize = hyp->GetMaxSize();
824     _hexoticIgnoreRidges = hyp->GetHexoticIgnoreRidges();
825     _hexoticInvalidElements = hyp->GetHexoticInvalidElements();
826     _hexoticSharpAngleThreshold = hyp->GetHexoticSharpAngleThreshold();
827     _hexoticNbProc = hyp->GetHexoticNbProc();
828     _hexoticWorkingDirectory = hyp->GetHexoticWorkingDirectory();
829     _hexoticVerbosity = hyp->GetHexoticVerbosity();
830     _hexoticMaxMemory = hyp->GetHexoticMaxMemory();
831     _hexoticSdMode = hyp->GetHexoticSdMode();
832     _sizeMaps = hyp->GetSizeMaps();
833   }
834   else {
835     cout << std::endl;
836     cout << "WARNING : The Hexotic default parameters are taken into account" << std::endl;
837     cout << "=======" << std::endl;
838     _hexesMinLevel = hyp->GetDefaultHexesMinLevel();
839     _hexesMaxLevel = hyp->GetDefaultHexesMaxLevel();
840     _hexesMinSize = hyp->GetDefaultMinSize();
841     _hexesMaxSize = hyp->GetDefaultMaxSize();
842     _hexoticIgnoreRidges = hyp->GetDefaultHexoticIgnoreRidges();
843     _hexoticInvalidElements = hyp->GetDefaultHexoticInvalidElements();
844     _hexoticSharpAngleThreshold = hyp->GetDefaultHexoticSharpAngleThreshold();
845     _hexoticNbProc = hyp->GetDefaultHexoticNbProc();
846     _hexoticWorkingDirectory = hyp->GetDefaultHexoticWorkingDirectory();
847     _hexoticVerbosity = hyp->GetDefaultHexoticVerbosity();
848     _hexoticMaxMemory = hyp->GetDefaultHexoticMaxMemory();
849     _hexoticSdMode = hyp->GetDefaultHexoticSdMode();
850     _sizeMaps = hyp->GetDefaultHexoticSizeMaps();
851   }
852 }
853
854 //=======================================================================
855 //function : getTmpDir
856 //purpose  :
857 //=======================================================================
858
859 static TCollection_AsciiString getTmpDir()
860 {
861   TCollection_AsciiString aTmpDir;
862
863   char *Tmp_dir = getenv("SALOME_TMP_DIR");
864 #ifdef WIN32
865   if(Tmp_dir == NULL) {
866     Tmp_dir = getenv("TEMP");
867     if( Tmp_dir== NULL )
868       Tmp_dir = getenv("TMP");
869 #endif
870
871   if(Tmp_dir != NULL) {
872     aTmpDir = Tmp_dir;
873 #ifdef WIN32
874     if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
875 #else
876     if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
877 #endif
878   }
879   else {
880 #ifdef WIN32
881     aTmpDir = TCollection_AsciiString("C:\\");
882 #else
883     aTmpDir = TCollection_AsciiString("/tmp/");
884 #endif
885   }
886   return aTmpDir;
887 }
888
889 //=======================================================================
890 //function : getSuffix
891 //purpose  : Returns a suffix that will be unique for the current process
892 //=======================================================================
893
894 static TCollection_AsciiString getSuffix()
895 {
896   TCollection_AsciiString aSuffix = "";
897   aSuffix += "_";
898   aSuffix += getenv("USER");
899   aSuffix += "_";
900   aSuffix += Kernel_Utils::GetHostname().c_str();
901   aSuffix += "_";
902   aSuffix += getpid();
903
904   return aSuffix;
905 }
906
907 //================================================================================
908 /*!
909  * \brief Returns a command to run Hexotic mesher
910  */
911 //================================================================================
912
913 std::string HexoticPlugin_Hexotic::getHexoticCommand(const TCollection_AsciiString& Hexotic_In,
914                                                      const TCollection_AsciiString& Hexotic_Out,
915                                                      const TCollection_AsciiString& Hexotic_SizeMap_Prefix) const
916 {
917   cout << std::endl;
918   cout << "Hexotic execution..." << std::endl;
919   cout << _name << " parameters :" << std::endl;
920   cout << "    " << _name << " Verbosity = " << _hexoticVerbosity << std::endl;
921   cout << "    " << _name << " Max Memory = " << _hexoticMaxMemory << std::endl;
922   cout << "    " << _name << " Segments Min Level = " << _hexesMinLevel << std::endl;
923   cout << "    " << _name << " Segments Max Level = " << _hexesMaxLevel << std::endl;
924   cout << "    " << _name << " Segments Min Size = " << _hexesMinSize << std::endl;
925   cout << "    " << _name << " Segments Max Size = " << _hexesMaxSize << std::endl;
926   cout << "    " << "Hexotic can ignore ridges : " << (_hexoticIgnoreRidges ? "yes":"no") << std::endl;
927   cout << "    " << "Hexotic authorize invalide elements : " << ( _hexoticInvalidElements ? "yes":"no") << std::endl;
928   cout << "    " << _name << " Sharp angle threshold = " << _hexoticSharpAngleThreshold << " degrees" << std::endl;
929   cout << "    " << _name << " Number of threads = " << _hexoticNbProc << std::endl;
930   cout << "    " << _name << " Working directory = \"" << _hexoticWorkingDirectory << "\"" << std::endl;
931   cout << "    " << _name << " Sub. Dom mode = " << _hexoticSdMode << std::endl;
932
933   TCollection_AsciiString run_Hexotic( "mg-hexa.exe" );
934
935   TCollection_AsciiString minl = " --min_level ", maxl = " --max_level ", angle = " --ridge_angle ";
936   TCollection_AsciiString mins = " --min_size ", maxs = " --max_size ";
937   TCollection_AsciiString in   = " --in ",   out  = " --out ";
938   TCollection_AsciiString sizeMap = " --read_sizemap ";
939   TCollection_AsciiString ignoreRidges = " --compute_ridges no ", invalideElements = " --allow_invalid_elements yes ";
940   TCollection_AsciiString subdom = " --components ";
941   TCollection_AsciiString proc = " --max_number_of_threads ";
942   TCollection_AsciiString verb = " --verbose ";
943   TCollection_AsciiString maxmem = " --max_memory ";
944
945   TCollection_AsciiString minLevel, maxLevel, minSize, maxSize, sharpAngle, mode, nbproc, verbosity, maxMemory;
946   minLevel = _hexesMinLevel;
947   maxLevel = _hexesMaxLevel;
948   minSize = _hexesMinSize;
949   maxSize = _hexesMaxSize;
950   sharpAngle = _hexoticSharpAngleThreshold;
951   // Mode translation for mg-tetra 1.1
952   switch ( _hexoticSdMode )
953   {
954     case 1:
955       mode = "outside_skin_only";
956       break;
957     case 2:
958       mode = "outside_components";
959       break;
960     case 3:
961       mode = "all";
962       break;
963     case 4:
964       mode = "all --manifold_geometry no";
965       break;
966   }
967   nbproc = _hexoticNbProc;
968   verbosity = _hexoticVerbosity;
969   maxMemory = _hexoticMaxMemory;
970
971   if (_hexoticIgnoreRidges)
972     run_Hexotic +=  ignoreRidges;
973
974   if (_hexoticInvalidElements)
975     run_Hexotic +=  invalideElements;
976
977   if (_hexesMinSize > 0)
978     run_Hexotic +=  mins + minSize;
979
980   if (_hexesMaxSize > 0)
981     run_Hexotic +=  maxs + maxSize;
982
983   if (_hexesMinLevel > 0)
984     run_Hexotic +=  minl + minLevel;
985
986   if (_hexesMaxLevel > 0)
987     run_Hexotic +=  maxl + maxLevel;
988
989   if (_hexoticSharpAngleThreshold > 0)
990     run_Hexotic +=  angle + sharpAngle;
991   
992   if (_sizeMaps.begin() != _sizeMaps.end())
993     run_Hexotic += sizeMap + Hexotic_SizeMap_Prefix;
994
995   run_Hexotic += in + Hexotic_In + out + Hexotic_Out;
996   run_Hexotic += subdom + mode;
997   run_Hexotic += proc + nbproc;
998   run_Hexotic += verb + verbosity;
999   run_Hexotic += maxmem + maxMemory;
1000
1001   return run_Hexotic.ToCString();
1002 }
1003
1004 // TODO : this is a duplication of some code found in BLSURFPlugin_BLSURF find a proper
1005 // way to share it
1006 TopoDS_Shape HexoticPlugin_Hexotic::entryToShape(std::string entry)
1007 {
1008   MESSAGE("HexoticPlugin_Hexotic::entryToShape "<<entry );
1009   GEOM::GEOM_Object_var aGeomObj;
1010   TopoDS_Shape S = TopoDS_Shape();
1011   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
1012   if (!aSObj->_is_nil()) {
1013     CORBA::Object_var obj = aSObj->GetObject();
1014     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
1015     aSObj->UnRegister();
1016   }
1017   if ( !aGeomObj->_is_nil() )
1018     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
1019   return S;
1020 }
1021
1022 //================================================================================
1023 /*!
1024  * \brief Produces a .mesh file with the size maps informations to give to Hexotic
1025  */
1026 //================================================================================
1027 std::vector<std::string> HexoticPlugin_Hexotic::writeSizeMapFile( std::string sizeMapPrefix )
1028 {
1029   HexoticPlugin_Hypothesis::THexoticSizeMaps::iterator it ;
1030   
1031   // The GMF driver will be used to write the size map file
1032   DriverGMF_Write aWriter;
1033   aWriter.SetSizeMapPrefix( sizeMapPrefix ); 
1034   
1035   std::vector<Control_Pnt> points;
1036   // Iterate on the size maps
1037   for (it=_sizeMaps.begin(); it!=_sizeMaps.end(); it++)
1038   {
1039     // Step 1 : Get the GEOM object entry and the size 
1040     // from the _sizeMaps infos
1041     std::string anEntry = it->first;
1042     double aLocalSize = it->second;
1043     TopoDS_Shape aShape = entryToShape( anEntry );
1044     
1045     // Step 2 : Create the points
1046     createControlPoints( aShape, aLocalSize, points );
1047   }
1048   // Write the .mesh size map file
1049   aWriter.PerformSizeMap(points); 
1050   return aWriter.GetSizeMapFiles();
1051 }
1052
1053 //================================================================================
1054 /*!
1055  * \brief Fills a vector of points from which a size map input file can be written
1056  */
1057 //================================================================================
1058 void HexoticPlugin_Hexotic::createControlPoints( const TopoDS_Shape& aShape, 
1059                                                  const double& theSize, 
1060                                                  std::vector<Control_Pnt>& thePoints )
1061
1062   if ( aShape.ShapeType() == TopAbs_VERTEX )
1063   {
1064     gp_Pnt aPnt = BRep_Tool::Pnt( TopoDS::Vertex(aShape) );
1065     Control_Pnt aControl_Pnt( aPnt, theSize );
1066     thePoints.push_back( aControl_Pnt );
1067   }
1068   if ( aShape.ShapeType() == TopAbs_EDGE )
1069   {
1070     createPointsSampleFromEdge( aShape, theSize, thePoints );  
1071   }
1072   else if ( aShape.ShapeType() == TopAbs_WIRE )
1073   {
1074     TopExp_Explorer Ex;
1075     for (Ex.Init(aShape,TopAbs_EDGE); Ex.More(); Ex.Next()) 
1076     {
1077       createPointsSampleFromEdge( Ex.Current(), theSize, thePoints );
1078     } 
1079   }
1080   else if ( aShape.ShapeType() ==  TopAbs_FACE )
1081   {
1082     createPointsSampleFromFace( aShape, theSize, thePoints ); 
1083   }
1084   else if ( aShape.ShapeType() ==  TopAbs_SOLID )
1085   {
1086     createPointsSampleFromSolid( aShape, theSize, thePoints ); 
1087   }
1088   else if ( aShape.ShapeType() == TopAbs_COMPOUND )
1089   {
1090     TopoDS_Iterator it( aShape );
1091     for(; it.More(); it.Next())
1092     {
1093       createControlPoints( it.Value(), theSize, thePoints );
1094     }
1095   }
1096 }
1097
1098 //================================================================================
1099 /*!
1100  * \brief Fills a vector of points with point samples approximately 
1101  * \brief spaced with a given size
1102  */
1103 //================================================================================
1104 void HexoticPlugin_Hexotic::createPointsSampleFromEdge( const TopoDS_Shape& aShape, 
1105                                                         const double& theSize, 
1106                                                         std::vector<Control_Pnt>& thePoints )
1107 {
1108   double step = theSize;
1109   double first, last;  
1110   Handle( Geom_Curve ) aCurve = BRep_Tool::Curve( TopoDS::Edge( aShape ), first, last );
1111   GeomAdaptor_Curve C ( aCurve );
1112   GCPnts_UniformAbscissa DiscretisationAlgo(C, step , first, last, Precision::Confusion());
1113   int nbPoints = DiscretisationAlgo.NbPoints();
1114   
1115   for ( int i = 1; i <= nbPoints; i++ )
1116   {
1117     double param = DiscretisationAlgo.Parameter( i );
1118     Control_Pnt aPnt;
1119     aCurve->D0( param, aPnt );
1120     aPnt.SetSize(theSize);
1121     thePoints.push_back( aPnt );
1122   }  
1123 }
1124
1125 //================================================================================
1126 /*!
1127  * \brief Fills a vector of points with point samples approximately 
1128  * \brief spaced with a given size
1129  */
1130 //================================================================================
1131 void HexoticPlugin_Hexotic::createPointsSampleFromFace( const TopoDS_Shape& aShape, 
1132                                                         const double& theSize, 
1133                                                         std::vector<Control_Pnt>& thePoints )
1134 {
1135   BRepMesh_IncrementalMesh M(aShape, 0.01, Standard_True);
1136   TopLoc_Location aLocation;
1137   TopoDS_Face aFace = TopoDS::Face(aShape);
1138
1139   // Triangulate the face
1140   Handle(Poly_Triangulation) aTri = BRep_Tool::Triangulation (aFace, aLocation);
1141   
1142   // Get the transformation associated to the face location
1143   gp_Trsf aTrsf = aLocation.Transformation();
1144   
1145   // Get triangles
1146   int nbTriangles = aTri->NbTriangles();
1147   Poly_Array1OfTriangle triangles(1,nbTriangles);
1148   triangles=aTri->Triangles();
1149   
1150   // GetNodes
1151   int nbNodes = aTri->NbNodes();
1152   TColgp_Array1OfPnt nodes(1,nbNodes);
1153   nodes = aTri->Nodes();
1154
1155   // Iterate on triangles and subdivide them
1156   for(int i=1; i<=nbTriangles; i++)
1157   {
1158      Poly_Triangle aTriangle = triangles.Value(i);
1159      gp_Pnt p1 = nodes.Value(aTriangle.Value(1));
1160      gp_Pnt p2 = nodes.Value(aTriangle.Value(2));
1161      gp_Pnt p3 = nodes.Value(aTriangle.Value(3));
1162      
1163      p1.Transform(aTrsf);
1164      p2.Transform(aTrsf);
1165      p3.Transform(aTrsf);
1166      
1167      subdivideTriangle(p1, p2, p3, theSize, thePoints);  
1168   }
1169 }
1170
1171 //================================================================================
1172 /*!
1173  * \brief Fills a vector of points with point samples approximately 
1174  * \brief spaced with a given size
1175  */
1176 //================================================================================
1177 void HexoticPlugin_Hexotic::createPointsSampleFromSolid( const TopoDS_Shape& aShape, 
1178                                                          const double& theSize, 
1179                                                          std::vector<Control_Pnt>& thePoints )
1180 {
1181   // Compute the bounding box
1182   double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
1183   Bnd_Box B;               
1184   BRepBndLib::Add(aShape, B);
1185   B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
1186   
1187   // Create the points
1188   double step = theSize;
1189   
1190   for ( double x=Xmin; x-Xmax<Precision::Confusion(); x=x+step )
1191   {
1192     for ( double y=Ymin; y-Ymax<Precision::Confusion(); y=y+step )
1193     {
1194       // Step1 : generate the Zmin -> Zmax line
1195       gp_Pnt startPnt(x, y, Zmin);
1196       gp_Pnt endPnt(x, y, Zmax);
1197       gp_Vec aVec(startPnt, endPnt);
1198       gp_Lin aLine(startPnt, aVec);
1199       double endParam = Zmax - Zmin;
1200       
1201       // Step2 : for each face of aShape:
1202       std::set<double> intersections;
1203       std::set<double>::iterator it = intersections.begin();
1204       
1205       TopExp_Explorer Ex;
1206       for (Ex.Init(aShape,TopAbs_FACE); Ex.More(); Ex.Next()) 
1207       { 
1208         // check if there is an intersection
1209         IntCurvesFace_Intersector anIntersector(TopoDS::Face(Ex.Current()), Precision::Confusion());
1210         anIntersector.Perform(aLine, 0, endParam);
1211         
1212         // get the intersection's parameter and store it
1213         int nbPoints = anIntersector.NbPnt();
1214         for(int i = 0 ; i < nbPoints ; i++ )
1215         {
1216           it = intersections.insert( it, anIntersector.WParameter(i+1) );
1217         }
1218       }
1219       // Step3 : go through the line chunk by chunk 
1220       if ( intersections.begin() != intersections.end() )
1221       {
1222         std::set<double>::iterator intersectionsIterator=intersections.begin();
1223         double first = *intersectionsIterator;
1224         intersectionsIterator++;
1225         bool innerPoints = true; 
1226         for ( ; intersectionsIterator!=intersections.end() ; intersectionsIterator++ )
1227         {
1228           double second = *intersectionsIterator;
1229           if ( innerPoints )
1230           {
1231             // If the last chunk was outside of the shape or this is the first chunk
1232             // add the points in the range [first, second] to the points vector
1233             double localStep = (second -first) / ceil( (second - first) / step );
1234             for ( double z = Zmin + first; z < Zmin + second; z = z + localStep )
1235             {
1236               thePoints.push_back(Control_Pnt( x, y, z, theSize ));
1237             }
1238             thePoints.push_back(Control_Pnt( x, y, Zmin + second, theSize ));
1239           }
1240           first = second;
1241           innerPoints = !innerPoints;
1242         }
1243       }
1244     }
1245   }
1246 }
1247
1248 //================================================================================
1249 /*!
1250  * \brief Subdivides a triangle until it reaches a certain size (recursive function)
1251  */
1252 //================================================================================
1253 void HexoticPlugin_Hexotic::subdivideTriangle( const gp_Pnt& p1, 
1254                                                const gp_Pnt& p2, 
1255                                                const gp_Pnt& p3, 
1256                                                const double& theSize, 
1257                                                std::vector<Control_Pnt>& thePoints)
1258 {
1259   // Size threshold to stop subdividing
1260   // This value ensures that two control points are distant no more than 2*theSize
1261   // as shown below
1262   //
1263   // The greater distance D of the mass center M to each Edge is 1/3 * Median 
1264   // and Median < sqrt(3/4) * a  where a is the greater side (by using Apollonius' thorem). 
1265   // So D < 1/3 * sqrt(3/4) * a and if a < sqrt(3) * S then D < S/2
1266   // and the distance between two mass centers of two neighbouring triangles 
1267   // sharing an edge is < 2 * 1/2 * S = S
1268   // If the traingles share a Vertex and no Edge the distance of the mass centers 
1269   // to the Vertices is 2*D < S so the mass centers are distant of less than 2*S 
1270   
1271   double threshold = sqrt( 3 ) * theSize;
1272   
1273   if ( (p1.Distance(p2) > threshold ||
1274         p2.Distance(p3) > threshold ||
1275         p3.Distance(p1) > threshold))
1276   { 
1277     std::vector<gp_Pnt> midPoints = computePointsForSplitting(p1, p2, p3);
1278  
1279     subdivideTriangle( midPoints[0], midPoints[1], midPoints[2], theSize, thePoints );
1280     subdivideTriangle( midPoints[0], p2, midPoints[1], theSize, thePoints );
1281     subdivideTriangle( midPoints[2], midPoints[1], p3, theSize, thePoints );
1282     subdivideTriangle( p1, midPoints[0], midPoints[2], theSize, thePoints );
1283   }
1284   else
1285   {
1286     double x = (p1.X() + p2.X() + p3.X()) / 3 ;
1287     double y = (p1.Y() + p2.Y() + p3.Y()) / 3 ;
1288     double z = (p1.Z() + p2.Z() + p3.Z()) / 3 ;
1289     
1290     Control_Pnt massCenter( x ,y ,z, theSize );
1291     thePoints.push_back( massCenter );
1292   }
1293 }
1294
1295 //================================================================================
1296 /*!
1297  * \brief Returns the appropriate points for splitting a triangle
1298  * \brief the tangency points of the incircle are used in order to have mostly
1299  * \brief well-shaped sub-triangles
1300  */
1301 //================================================================================
1302 std::vector<gp_Pnt> HexoticPlugin_Hexotic::computePointsForSplitting( const gp_Pnt& p1, 
1303                                                                       const gp_Pnt& p2, 
1304                                                                       const gp_Pnt& p3 )
1305 {
1306   std::vector<gp_Pnt> midPoints;
1307   //Change coordinates
1308   gp_Trsf Trsf_1;            // Identity transformation
1309   gp_Ax3 reference_system(gp::Origin(), gp::DZ(), gp::DX());   // OXY
1310  
1311   gp_Vec Vx(p1, p3);
1312   gp_Vec Vaux(p1, p2);
1313   gp_Dir Dx(Vx);
1314   gp_Dir Daux(Vaux);
1315   gp_Dir Dz = Dx.Crossed(Daux);
1316   gp_Ax3 current_system(p1, Dz, Dx);
1317   
1318   Trsf_1.SetTransformation( reference_system, current_system );
1319   
1320   gp_Pnt A = p1.Transformed(Trsf_1);
1321   gp_Pnt B = p2.Transformed(Trsf_1);
1322   gp_Pnt C = p3.Transformed(Trsf_1);
1323   
1324   double a =  B.Distance(C) ;
1325   double b =  A.Distance(C) ;
1326   double c =  B.Distance(A) ;
1327   
1328   // Incenter coordinates
1329   // see http://mathworld.wolfram.com/Incenter.html
1330   double Xi = ( b*B.X() + c*C.X() ) / ( a + b + c );
1331   double Yi = ( b*B.Y() ) / ( a + b + c );
1332   gp_Pnt Center(Xi, Yi, 0);
1333   
1334   // Calculate the tangency points of the incircle
1335   gp_Pnt T1 = tangencyPoint( A, B, Center);
1336   gp_Pnt T2 = tangencyPoint( B, C, Center);
1337   gp_Pnt T3 = tangencyPoint( C, A, Center);
1338   
1339   gp_Pnt p1_2 = T1.Transformed(Trsf_1.Inverted());
1340   gp_Pnt p2_3 = T2.Transformed(Trsf_1.Inverted());
1341   gp_Pnt p3_1 = T3.Transformed(Trsf_1.Inverted());
1342
1343   midPoints.push_back(p1_2);
1344   midPoints.push_back(p2_3);
1345   midPoints.push_back(p3_1);
1346   
1347   return midPoints;
1348 }
1349
1350 //================================================================================
1351 /*!
1352  * \brief Computes the tangency points of the circle of center Center with
1353  * \brief the straight line (p1 p2)
1354  */
1355 //================================================================================
1356 gp_Pnt HexoticPlugin_Hexotic::tangencyPoint(const gp_Pnt& p1,
1357                                             const gp_Pnt& p2,
1358                                             const gp_Pnt& Center)
1359 {
1360   double Xt = 0;
1361   double Yt = 0;
1362   
1363   // The tangency point is the intersection of the straight line (p1 p2)
1364   // and the straight line (Center T) which is orthogonal to (p1 p2)
1365   if ( fabs(p1.X() - p2.X()) <= Precision::Confusion() )
1366   {
1367     Xt=p1.X();     // T is on (p1 p2)
1368     Yt=Center.Y(); // (Center T) is orthogonal to (p1 p2)
1369   }
1370   else if ( fabs(p1.Y() - p2.Y()) <= Precision::Confusion() )
1371   {
1372     Yt=p1.Y();     // T is on (p1 p2) 
1373     Xt=Center.X(); // (Center T) is orthogonal to (p1 p2)
1374   }
1375   else
1376   {
1377     // First straight line coefficients (equation y=a*x+b)
1378     double a = (p2.Y() - p1.Y()) / (p2.X() - p1.X())  ;
1379     double b = p1.Y() - a*p1.X();         // p1 is on this straight line
1380     
1381     // Second straight line coefficients (equation y=c*x+d)
1382     double c = -1 / a;                    // The 2 lines are orthogonal
1383     double d = Center.Y() - c*Center.X(); // Center is on this straight line
1384     
1385     Xt = (d - b) / (a - c);
1386     Yt = a*Xt + b;
1387   }
1388   
1389   return gp_Pnt( Xt, Yt, 0 );
1390 }
1391
1392 //=============================================================================
1393 /*!
1394  * Here we are going to use the Hexotic mesher
1395  */
1396 //=============================================================================
1397
1398 bool HexoticPlugin_Hexotic::Compute(SMESH_Mesh&          aMesh,
1399                                      const TopoDS_Shape& aShape)
1400 {
1401   _compute_canceled = false;
1402   bool Ok = true;
1403   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1404   TCollection_AsciiString hexahedraMessage;
1405
1406   if (_iShape == 0 && _nbShape == 0) {
1407     _nbShape = countShape( meshDS, TopAbs_SOLID );  // we count the number of shapes
1408   }
1409
1410   // to prevent from displaying error message after computing,
1411   // SetIsAlwaysComputed( true ) to empty sub-meshes
1412   vector< SMESH_subMesh* > subMeshesAlwaysComp;
1413   for ( int i = 0; i < _nbShape; ++i )
1414     if ( SMESH_subMesh* sm = aMesh.GetSubMeshContaining( aShape ))
1415     {
1416       SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true,
1417                                                                /*complexShapeFirst=*/false);
1418       while ( smIt->more() )
1419       {
1420         sm = smIt->next();
1421         if ( !sm->IsMeshComputed() )
1422         {
1423           sm->SetIsAlwaysComputed( true );
1424           subMeshesAlwaysComp.push_back( sm );
1425         }
1426       }
1427     }
1428
1429   _iShape++;
1430
1431   if (_iShape == _nbShape ) {
1432
1433     // create bounding box for each shape of the compound
1434
1435     int iShape = 0;
1436     TopoDS_Shape *tabShape;
1437     double **tabBox;
1438
1439     tabShape = new TopoDS_Shape[_nbShape];
1440     tabBox   = new double*[_nbShape];
1441     for (int i=0; i<_nbShape; i++)
1442       tabBox[i] = new double[6];
1443     double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
1444
1445     TopExp_Explorer expBox (meshDS->ShapeToMesh(), TopAbs_SOLID);
1446     for (; expBox.More(); expBox.Next()) {
1447       tabShape[iShape] = expBox.Current();
1448       Bnd_Box BoundingBox;
1449       BRepBndLib::Add(expBox.Current(), BoundingBox);
1450       BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
1451       tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
1452       tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
1453       tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
1454       iShape++;
1455     }
1456
1457     SetParameters(_hypothesis);
1458
1459 //     TCollection_AsciiString aTmpDir = getTmpDir();
1460     TCollection_AsciiString aTmpDir = TCollection_AsciiString(_hexoticWorkingDirectory.c_str());
1461 #ifdef WIN32
1462     if ( aTmpDir.Value(aTmpDir.Length()) != '\\' ) aTmpDir += '\\';
1463 #else
1464     if ( aTmpDir.Value(aTmpDir.Length()) != '/' ) aTmpDir += '/';
1465 #endif
1466     TCollection_AsciiString Hexotic_In(""), Hexotic_Out, Hexotic_SizeMap_Prefix;
1467     TCollection_AsciiString modeFile_In( "chmod 666 " ), modeFile_Out( "chmod 666 " );
1468     TCollection_AsciiString aLogFileName = aTmpDir + "Hexotic"+getSuffix()+".log";    // log
1469
1470     std::map <int,int> aSmdsToHexoticIdMap;
1471     std::map <int,const SMDS_MeshNode*> aHexoticIdToNodeMap;
1472
1473     Hexotic_Out = aTmpDir + "Hexotic"+getSuffix()+"_Out.mesh";
1474 #ifdef WITH_BLSURFPLUGIN
1475     bool defaultInputFile = true;
1476     if (_blsurfHypo && !_blsurfHypo->GetQuadAllowed()) {
1477       Hexotic_In = TCollection_AsciiString(_blsurfHypo->GetGMFFile().c_str());
1478       if (Hexotic_In != "")
1479         defaultInputFile = false;
1480     }
1481     if (defaultInputFile) {
1482 #endif
1483       Hexotic_In  = aTmpDir + "Hexotic"+getSuffix()+"_In.mesh";
1484       removeHexoticFiles(Hexotic_In, Hexotic_Out);
1485       cout << std::endl;
1486       cout << "Creating Hexotic input mesh file : " << Hexotic_In << std::endl;
1487       aMesh.ExportGMF(Hexotic_In.ToCString(), meshDS, true);
1488 #ifdef WITH_BLSURFPLUGIN
1489     }
1490     else {
1491       removeFile( Hexotic_Out );
1492     }
1493 #endif
1494     
1495     Hexotic_SizeMap_Prefix = aTmpDir + "Hexotic_SizeMap" + getSuffix();
1496     std::vector<std::string> sizeMapFiles = writeSizeMapFile( Hexotic_SizeMap_Prefix.ToCString() );
1497     
1498     std::string run_Hexotic = getHexoticCommand(Hexotic_In, Hexotic_Out, Hexotic_SizeMap_Prefix);
1499     run_Hexotic += std::string(" 1> ") + aLogFileName.ToCString();  // dump into file
1500
1501     cout << std::endl;
1502     cout << "Hexotic command : " << run_Hexotic << std::endl;
1503
1504 #ifndef WIN32    
1505     modeFile_In += Hexotic_In;
1506     system( modeFile_In.ToCString() );
1507 #endif
1508     aSmdsToHexoticIdMap.clear();
1509     aHexoticIdToNodeMap.clear();
1510
1511     MESSAGE("HexoticPlugin_Hexotic::Compute");
1512
1513     int status = system( run_Hexotic.data() );
1514
1515     // --------------
1516     // read a result
1517     // --------------
1518
1519     std::ifstream fileRes( Hexotic_Out.ToCString() );
1520     modeFile_Out += Hexotic_Out;
1521     system( modeFile_Out.ToCString() );
1522     if ( ! fileRes.fail() ) {
1523       Ok = readResult( Hexotic_Out.ToCString(),
1524                        this,
1525                        meshDS, _nbShape, tabShape, tabBox );
1526       if(Ok) {
1527 /*********************
1528 // TODO: Detect and remove elements in holes in case of sd mode = 4
1529       // Remove previous nodes and elements
1530       SMDS_ElemIteratorPtr itElement = meshDS->elementsIterator();
1531       SMDS_NodeIteratorPtr itNode = meshDS->nodesIterator();
1532     
1533       while ( itElement->more() )
1534         meshDS->RemoveElement( itElement->next() );
1535       while ( itNode->more() )
1536         meshDS->RemoveNode( itNode->next() );
1537   
1538       SMESH_ComputeErrorPtr myError = aMesh.GMFToMesh(Hexotic_Out.ToCString());
1539       if (myError)
1540 */
1541         hexahedraMessage = "success";
1542         #ifndef _DEBUG_
1543         removeFile(Hexotic_Out);
1544         removeFile(Hexotic_In);
1545         removeFile(aLogFileName);
1546         for( int i=0; i<sizeMapFiles.size(); i++)
1547         {
1548           removeFile( TCollection_AsciiString( sizeMapFiles[i].c_str() ) );
1549         }
1550         #endif
1551       }
1552       else {
1553         hexahedraMessage = "failed"; 
1554       }
1555     }
1556     else {
1557       hexahedraMessage = "failed";
1558       cout << "Problem with Hexotic output file " << Hexotic_Out.ToCString() << std::endl;
1559       Ok = false;
1560       // analyse log file
1561       SMESH_File logFile( aLogFileName.ToCString() );
1562       if ( !logFile.eof() )
1563       {
1564         char msgLic[] = " Dlim ";
1565         const char* fileBeg = logFile.getPos(), *fileEnd = fileBeg + logFile.size();
1566         if ( std::search( fileBeg, fileEnd, msgLic, msgLic+strlen(msgLic)) != fileEnd )
1567           error("Licence problems.");
1568       }
1569 #ifndef WIN32
1570       if ( status > 0 && WEXITSTATUS(status) == 127 )
1571         error("hexotic: command not found");
1572 #else
1573       int err = errno;
1574       if ( status == 0 && err == ENOENT ) {
1575         error("hexotic: command not found");
1576       }
1577 #endif
1578     }
1579     cout << "Hexahedra meshing " << hexahedraMessage << std::endl;
1580     cout << std::endl;
1581
1582     // restore "always computed" flag of sub-meshes (0022127)
1583     for  ( size_t iSM = 0; iSM < subMeshesAlwaysComp.size(); ++iSM )
1584       subMeshesAlwaysComp[ iSM ]->SetIsAlwaysComputed( false );
1585
1586     delete [] tabShape;
1587     for (int i=0; i<_nbShape; i++)
1588       delete [] tabBox[i];
1589     delete [] tabBox;
1590     _nbShape = 0;
1591     _iShape  = 0;
1592   }
1593   if(_compute_canceled)
1594     return error(SMESH_Comment("interruption initiated by user"));
1595   return Ok;
1596 }
1597
1598 //=============================================================================
1599 /*!
1600  * \brief Computes mesh without geometry
1601  *  \param aMesh - the mesh
1602  *  \param aHelper - helper that must be used for adding elements to \aaMesh
1603  *  \retval bool - is a success
1604  *
1605  * The method is called if ( !aMesh->HasShapeToMesh() )
1606  */
1607 //=============================================================================
1608
1609 bool HexoticPlugin_Hexotic::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1610 {
1611   _compute_canceled = false;
1612 /*
1613   SMESH_ComputeErrorPtr myError = SMESH_ComputeError::New();
1614 */
1615   bool Ok = true;
1616   TCollection_AsciiString hexahedraMessage;
1617
1618   SetParameters(_hypothesis);
1619
1620   TCollection_AsciiString aTmpDir = getTmpDir();
1621   TCollection_AsciiString Hexotic_In, Hexotic_Out, Hexotic_SizeMap_Prefix;
1622   TCollection_AsciiString modeFile_In( "chmod 666 " ), modeFile_Out( "chmod 666 " );
1623   TCollection_AsciiString aLogFileName = aTmpDir + "Hexotic"+getSuffix()+".log";    // log
1624
1625   std::map <int,int> aSmdsToHexoticIdMap;
1626   std::map <int,const SMDS_MeshNode*> aHexoticIdToNodeMap;
1627
1628   Hexotic_In  = aTmpDir + "Hexotic"+getSuffix()+"_In.mesh";
1629   Hexotic_Out = aTmpDir + "Hexotic"+getSuffix()+"_Out.mesh";
1630   Hexotic_SizeMap_Prefix = aTmpDir + "Hexotic_SizeMap";
1631  
1632   
1633   std::vector<std::string> sizeMapFiles = writeSizeMapFile( Hexotic_SizeMap_Prefix.ToCString() );
1634
1635   std::string run_Hexotic = getHexoticCommand(Hexotic_In, Hexotic_Out, Hexotic_SizeMap_Prefix);
1636   run_Hexotic += std::string(" 1> ") + aLogFileName.ToCString();  // dump into file
1637
1638   removeHexoticFiles(Hexotic_In, Hexotic_Out);
1639
1640   cout << std::endl;
1641   cout << "Creating Hexotic input mesh file : " << Hexotic_In << std::endl;
1642   aMesh.ExportGMF(Hexotic_In.ToCString(), aHelper->GetMeshDS());
1643 #ifndef WIN32    
1644   modeFile_In += Hexotic_In;
1645   system( modeFile_In.ToCString() );
1646 #endif
1647   aSmdsToHexoticIdMap.clear();
1648   aHexoticIdToNodeMap.clear();
1649
1650   MESSAGE("HexoticPlugin_Hexotic::Compute");
1651
1652   cout << std::endl;
1653   cout << "Hexotic command : " << run_Hexotic << std::endl;
1654   system( run_Hexotic.data() );
1655
1656   // --------------
1657   // read a result
1658   // --------------
1659
1660   std::ifstream fileRes( Hexotic_Out.ToCString() );
1661   modeFile_Out += Hexotic_Out;
1662   system( modeFile_Out.ToCString() );
1663   if ( ! fileRes.fail() ) {
1664     Ok = readResult( Hexotic_Out.ToCString(),
1665                      this,
1666                      aHelper );
1667     if(Ok)
1668 /*
1669     // Remove previous nodes and elements
1670     SMDS_ElemIteratorPtr itElement = aHelper->GetMeshDS()->elementsIterator();
1671     SMDS_NodeIteratorPtr itNode = aHelper->GetMeshDS()->nodesIterator();
1672     
1673     while ( itElement->more() )
1674       aHelper->GetMeshDS()->RemoveElement( itElement->next() );
1675     while ( itNode->more() )
1676       aHelper->GetMeshDS()->RemoveNode( itNode->next() );
1677
1678     // Import GMF mesh
1679     myError = aMesh.GMFToMesh(Hexotic_Out.ToCString());
1680     
1681     itElement = aHelper->GetMeshDS()->elementsIterator();
1682     itNode = aHelper->GetMeshDS()->nodesIterator();
1683
1684     // Assign nodes and elements to the pseudo shape
1685     while ( itNode->more() )
1686       aHelper->GetMeshDS()->SetNodeInVolume(itNode->next(), 1);
1687     while ( itElement->more() )
1688       aHelper->GetMeshDS()->SetMeshElementOnShape(itElement->next(), 1);
1689
1690     if(myError->IsOK())
1691 */
1692       hexahedraMessage = "success";
1693     else
1694       hexahedraMessage = "failed";
1695   }
1696   else {
1697 /*
1698     myError->myName = COMPERR_EXCEPTION;
1699 */
1700     hexahedraMessage = "failed";
1701     cout << "Problem with Hexotic output file " << Hexotic_Out << std::endl;
1702     // analyse log file
1703     SMESH_File logFile( aLogFileName.ToCString() );
1704     if ( !logFile.eof() )
1705     {
1706       char msgLic[] = " Dlim ";
1707       const char* fileBeg = logFile.getPos(), *fileEnd = fileBeg + logFile.size();
1708       if ( std::search( fileBeg, fileEnd, msgLic, msgLic+strlen(msgLic)) != fileEnd )
1709         return error("Licence problems.");
1710     }
1711     return error(SMESH_Comment("Problem with Hexotic output file ")<<Hexotic_Out);
1712   }
1713   cout << "Hexahedra meshing " << hexahedraMessage << std::endl;
1714   cout << std::endl;
1715
1716   if(_compute_canceled)
1717     return error(SMESH_Comment("interruption initiated by user"));
1718   removeFile(Hexotic_Out);
1719   removeFile(Hexotic_In);
1720   removeFile(aLogFileName);
1721   for( int i=0; i<sizeMapFiles.size(); i++)
1722   {
1723     removeFile( TCollection_AsciiString(sizeMapFiles[i].c_str()) );
1724   }
1725   return Ok;
1726 /*
1727   return myError->IsOK();
1728 */
1729 }
1730
1731 //=============================================================================
1732 /*!
1733  *
1734  */
1735 //=============================================================================
1736
1737 bool HexoticPlugin_Hexotic::Evaluate(SMESH_Mesh& aMesh,
1738                                      const TopoDS_Shape& aShape,
1739                                      MapShapeNbElems& aResMap)
1740 {
1741   std::vector<int> aResVec(SMDSEntity_Last);
1742   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1743   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1744   aResMap.insert(std::make_pair(sm,aResVec));
1745
1746   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1747   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Evaluation is not implemented",this));
1748
1749   return true;
1750 }
1751
1752 void HexoticPlugin_Hexotic::CancelCompute()
1753 {
1754   _compute_canceled = true;
1755 #ifdef WNT
1756 #else
1757   TCollection_AsciiString aTmpDir = getTmpDir();
1758   TCollection_AsciiString Hexotic_In = aTmpDir + "Hexotic_In.mesh";
1759   TCollection_AsciiString cmd = TCollection_AsciiString("ps ux | grep ") + Hexotic_In;
1760   cmd += TCollection_AsciiString(" | grep -v grep | awk '{print $2}' | xargs kill -9 > /dev/null 2>&1");
1761   system( cmd.ToCString() );
1762 #endif
1763 }