src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

   1 // Copyright (C) 2004-2016  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, or (at your option) any later version.
   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      : GHS3DPlugin_GHS3D.cxx
  22 // Created   :
  23 // Author    : Edward AGAPOV, modified by Lioka RAZAFINDRAZAKA (CEA) 09/02/2007
  24 // Project   : SALOME
  25 //=============================================================================
  26 //
  27 #include "GHS3DPlugin_GHS3D.hxx"
  28 #include "GHS3DPlugin_Hypothesis.hxx"
  29
  30 #include <SMDS_FaceOfNodes.hxx>
  31 #include <SMDS_MeshElement.hxx>
  32 #include <SMDS_MeshNode.hxx>
  33 #include <SMDS_VolumeOfNodes.hxx>
  34 #include <SMESHDS_Group.hxx>
  35 #include <SMESH_Comment.hxx>
  36 #include <SMESH_Group.hxx>
  37 #include <SMESH_HypoFilter.hxx>
  38 #include <SMESH_Mesh.hxx>
  39 #include <SMESH_MeshAlgos.hxx>
  40 #include <SMESH_MeshEditor.hxx>
  41 #include <SMESH_MesherHelper.hxx>
  42 #include <SMESH_OctreeNode.hxx>
  43 #include <SMESH_subMeshEventListener.hxx>
  44 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  45 #include <StdMeshers_ViscousLayers.hxx>
  46
  47 #include <BRepAdaptor_Surface.hxx>
  48 #include <BRepBndLib.hxx>
  49 #include <BRepBuilderAPI_MakeVertex.hxx>
  50 #include <BRepClass3d.hxx>
  51 #include <BRepClass3d_SolidClassifier.hxx>
  52 #include <BRepExtrema_DistShapeShape.hxx>
  53 #include <BRepGProp.hxx>
  54 #include <BRepTools.hxx>
  55 #include <BRep_Tool.hxx>
  56 #include <Bnd_Box.hxx>
  57 #include <GProp_GProps.hxx>
  58 #include <GeomAPI_ProjectPointOnSurf.hxx>
  59 #include <OSD_File.hxx>
  60 #include <Precision.hxx>
  61 #include <Standard_ErrorHandler.hxx>
  62 #include <Standard_Failure.hxx>
  63 #include <Standard_ProgramError.hxx>
  64 #include <TopExp.hxx>
  65 #include <TopExp_Explorer.hxx>
  66 #include <TopTools_IndexedMapOfShape.hxx>
  67 #include <TopTools_ListIteratorOfListOfShape.hxx>
  68 #include <TopTools_MapOfShape.hxx>
  69 #include <TopoDS.hxx>
  70 #include <TopoDS_Shell.hxx>
  71 #include <TopoDS_Solid.hxx>
  72
  73 #include <Basics_Utils.hxx>
  74 #include <utilities.h>
  75
  76 #ifdef WIN32
  77 #include <io.h>
  78 #else
  79 #include <sys/sysinfo.h>
  80 #endif
  81 #include <algorithm>
  82 #include <errno.h>
  83
  84 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
  85
  86 extern "C"
  87 {
  88 #ifndef WIN32
  89 #include <unistd.h>
  90 #include <sys/mman.h>
  91 #endif
  92 #include <sys/stat.h>
  93 #include <fcntl.h>
  94 }
  95
  96 using namespace std;
  97
  98 #define HOLE_ID -1
  99
 100 // flags returning state of enforced entities, returned from writeGMFFile
 101 enum InvalidEnforcedFlags { FLAG_BAD_ENF_VERT = 1,
 102                             FLAG_BAD_ENF_NODE = 2,
 103                             FLAG_BAD_ENF_EDGE = 4,
 104                             FLAG_BAD_ENF_TRIA = 8
 105 };
 106 static std::string flagsToErrorStr( int anInvalidEnforcedFlags )
 107 {
 108   std::string str;
 109   if ( anInvalidEnforcedFlags != 0 )
 110   {
 111     if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_VERT )
 112       str = "There are enforced vertices incorrectly defined.\n";
 113     if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_NODE )
 114       str += "There are enforced nodes incorrectly defined.\n";
 115     if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_EDGE )
 116       str += "There are enforced edge incorrectly defined.\n";
 117     if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_TRIA )
 118       str += "There are enforced triangles incorrectly defined.\n";
 119   }
 120   return str;
 121 }
 122
 123 typedef const list<const SMDS_MeshFace*> TTriaList;
 124
 125 static const char theDomainGroupNamePrefix[] = "Domain_";
 126
 127 static void removeFile( const TCollection_AsciiString& fileName )
 128 {
 129   try {
 130     OSD_File( fileName ).Remove();
 131   }
 132   catch ( Standard_ProgramError ) {
 133     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
 134   }
 135 }
 136
 137 //=============================================================================
 138 /*!
 139  *
 140  */
 141 //=============================================================================
 142
 143 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
 144   : SMESH_3D_Algo(hypId, studyId, gen)
 145 {
 146   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
 147   _name = Name();
 148   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
 149   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
 150   _iShape=0;
 151   _nbShape=0;
 152   _compatibleHypothesis.push_back( GHS3DPlugin_Hypothesis::GetHypType());
 153   _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
 154   _requireShape = false; // can work without shape_studyId
 155
 156   smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 157   CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 158   SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 159
 160   MESSAGE("studyid = " << _studyId);
 161
 162   myStudy = NULL;
 163   myStudy = aStudyMgr->GetStudyByID(_studyId);
 164   if (myStudy)
 165     MESSAGE("myStudy->StudyId() = " << myStudy->StudyId());
 166
 167   _compute_canceled = false;
 168 }
 169
 170 //=============================================================================
 171 /*!
 172  *
 173  */
 174 //=============================================================================
 175
 176 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
 177 {
 178   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
 179 }
 180
 181 //=============================================================================
 182 /*!
 183  *
 184  */
 185 //=============================================================================
 186
 187 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&         aMesh,
 188                                           const TopoDS_Shape& aShape,
 189                                           Hypothesis_Status&  aStatus )
 190 {
 191   aStatus = SMESH_Hypothesis::HYP_OK;
 192
 193   _hyp = 0;
 194   _viscousLayersHyp = 0;
 195   _keepFiles = false;
 196   _removeLogOnSuccess = true;
 197   _logInStandardOutput = false;
 198
 199   const list <const SMESHDS_Hypothesis * >& hyps =
 200     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 201   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 202   for ( ; h != hyps.end(); ++h )
 203   {
 204     if ( !_hyp )
 205       _hyp = dynamic_cast< const GHS3DPlugin_Hypothesis*> ( *h );
 206     if ( !_viscousLayersHyp )
 207       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 208   }
 209   if ( _hyp )
 210   {
 211     _keepFiles = _hyp->GetKeepFiles();
 212     _removeLogOnSuccess = _hyp->GetRemoveLogOnSuccess();
 213     _logInStandardOutput = _hyp->GetStandardOutputLog();
 214   }
 215
 216   if ( _viscousLayersHyp )
 217     error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus ));
 218
 219   return aStatus == HYP_OK;
 220 }
 221
 222
 223 //=======================================================================
 224 //function : entryToShape
 225 //purpose  :
 226 //=======================================================================
 227
 228 TopoDS_Shape GHS3DPlugin_GHS3D::entryToShape(std::string entry)
 229 {
 230   MESSAGE("GHS3DPlugin_GHS3D::entryToShape "<<entry );
 231   GEOM::GEOM_Object_var aGeomObj;
 232   TopoDS_Shape S = TopoDS_Shape();
 233   SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 234   if (!aSObj->_is_nil() ) {
 235     CORBA::Object_var obj = aSObj->GetObject();
 236     aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 237     aSObj->UnRegister();
 238   }
 239   if ( !aGeomObj->_is_nil() )
 240     S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 241   return S;
 242 }
 243
 244 //=======================================================================
 245 //function : findShape
 246 //purpose  :
 247 //=======================================================================
 248
 249 // static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
 250 //                               TopoDS_Shape        aShape,
 251 //                               const TopoDS_Shape  shape[],
 252 //                               double**            box,
 253 //                               const int           nShape,
 254 //                               TopAbs_State *      state = 0)
 255 // {
 256 //   gp_XYZ aPnt(0,0,0);
 257 //   int j, iShape, nbNode = 4;
 258
 259 //   for ( j=0; j<nbNode; j++ ) {
 260 //     gp_XYZ p ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
 261 //     if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) {
 262 //       aPnt = p;
 263 //       break;
 264 //     }
 265 //     aPnt += p / nbNode;
 266 //   }
 267
 268 //   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 269 //   if (state) *state = SC.State();
 270 //   if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
 271 //     for (iShape = 0; iShape < nShape; iShape++) {
 272 //       aShape = shape[iShape];
 273 //       if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
 274 //               aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
 275 //               aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
 276 //         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 277 //         if (state) *state = SC.State();
 278 //         if (SC.State() == TopAbs_IN)
 279 //           break;
 280 //       }
 281 //     }
 282 //   }
 283 //   return aShape;
 284 // }
 285
 286 //================================================================================
 287 /*!
 288  * \brief returns id of a solid if a triangle defined by the nodes is a temporary face on a
 289  * side facet of pyramid and defines sub-domian outside the pyramid; else returns HOLE_ID
 290  */
 291 //================================================================================
 292
 293 static int checkTmpFace(const SMDS_MeshNode* node1,
 294                         const SMDS_MeshNode* node2,
 295                         const SMDS_MeshNode* node3)
 296 {
 297   // find a pyramid sharing the 3 nodes
 298   SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume);
 299   while ( vIt1->more() )
 300   {
 301     const SMDS_MeshElement* pyram = vIt1->next();
 302     if ( pyram->NbCornerNodes() != 5 ) continue;
 303     int i2, i3;
 304     if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 &&
 305          (i3 = pyram->GetNodeIndex( node3 )) >= 0 )
 306     {
 307       // Triangle defines sub-domian inside the pyramid if it's
 308       // normal points out of the pyram
 309
 310       // make i2 and i3 hold indices of base nodes of the pyram while
 311       // keeping the nodes order in the triangle
 312       const int iApex = 4;
 313       if ( i2 == iApex )
 314         i2 = i3, i3 = pyram->GetNodeIndex( node1 );
 315       else if ( i3 == iApex )
 316         i3 = i2, i2 = pyram->GetNodeIndex( node1 );
 317
 318       int i3base = (i2+1) % 4; // next index after i2 within the pyramid base
 319       bool isDomainInPyramid = ( i3base != i3 );
 320       return isDomainInPyramid ? HOLE_ID : pyram->getshapeId();
 321     }
 322   }
 323   return HOLE_ID;
 324 }
 325
 326 //=======================================================================
 327 //function : findShapeID
 328 //purpose  : find the solid corresponding to MG-Tetra sub-domain following
 329 //           the technique proposed in MG-Tetra manual (available within
 330 //           MG-Tetra installation) in chapter "B.4 Subdomain (sub-region) assignment".
 331 //           In brief: normal of the triangle defined by the given nodes
 332 //           points out of the domain it is associated to
 333 //=======================================================================
 334
 335 static int findShapeID(SMESH_Mesh&          mesh,
 336                        const SMDS_MeshNode* node1,
 337                        const SMDS_MeshNode* node2,
 338                        const SMDS_MeshNode* node3,
 339                        const bool           toMeshHoles)
 340 {
 341   const int invalidID = 0;
 342   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 343
 344   // face the nodes belong to
 345   vector<const SMDS_MeshNode *> nodes(3);
 346   nodes[0] = node1;
 347   nodes[1] = node2;
 348   nodes[2] = node3;
 349   const SMDS_MeshElement * face = meshDS->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/true);
 350   if ( !face )
 351     return checkTmpFace(node1, node2, node3);
 352 #ifdef _DEBUG_
 353   std::cout << "bnd face " << face->GetID() << " - ";
 354 #endif
 355   // geom face the face assigned to
 356   SMESH_MeshEditor editor(&mesh);
 357   int geomFaceID = editor.FindShape( face );
 358   if ( !geomFaceID )
 359     return checkTmpFace(node1, node2, node3);
 360   TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
 361   if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
 362     return invalidID;
 363   TopoDS_Face geomFace = TopoDS::Face( shape );
 364
 365   // solids bounded by geom face
 366   TopTools_IndexedMapOfShape solids, shells;
 367   TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
 368   for ( ; ansIt.More(); ansIt.Next() ) {
 369     switch ( ansIt.Value().ShapeType() ) {
 370     case TopAbs_SOLID:
 371       solids.Add( ansIt.Value() ); break;
 372     case TopAbs_SHELL:
 373       shells.Add( ansIt.Value() ); break;
 374     default:;
 375     }
 376   }
 377   // analyse found solids
 378   if ( solids.Extent() == 0 || shells.Extent() == 0)
 379     return invalidID;
 380
 381   const TopoDS_Solid& solid1 = TopoDS::Solid( solids(1) );
 382   if ( solids.Extent() == 1 )
 383   {
 384     if ( toMeshHoles )
 385       return meshDS->ShapeToIndex( solid1 );
 386
 387     // - Are we at a hole boundary face?
 388     if ( shells(1).IsSame( BRepClass3d::OuterShell( solid1 )) )
 389     { // - No, but maybe a hole is bound by two shapes? Does shells(1) touch another shell?
 390       bool touch = false;
 391       TopExp_Explorer eExp( shells(1), TopAbs_EDGE );
 392       // check if any edge of shells(1) belongs to another shell
 393       for ( ; eExp.More() && !touch; eExp.Next() ) {
 394         ansIt = mesh.GetAncestors( eExp.Current() );
 395         for ( ; ansIt.More() && !touch; ansIt.Next() ) {
 396           if ( ansIt.Value().ShapeType() == TopAbs_SHELL )
 397             touch = ( !ansIt.Value().IsSame( shells(1) ));
 398         }
 399       }
 400       if (!touch)
 401         return meshDS->ShapeToIndex( solid1 );
 402     }
 403   }
 404   // find orientation of geom face within the first solid
 405   TopExp_Explorer fExp( solid1, TopAbs_FACE );
 406   for ( ; fExp.More(); fExp.Next() )
 407     if ( geomFace.IsSame( fExp.Current() )) {
 408       geomFace = TopoDS::Face( fExp.Current() );
 409       break;
 410     }
 411   if ( !fExp.More() )
 412     return invalidID; // face not found
 413
 414   // normale to triangle
 415   gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
 416   gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
 417   gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
 418   gp_Vec vec12( node1Pnt, node2Pnt );
 419   gp_Vec vec13( node1Pnt, node3Pnt );
 420   gp_Vec meshNormal = vec12 ^ vec13;
 421   if ( meshNormal.SquareMagnitude() < DBL_MIN )
 422     return invalidID;
 423
 424   // get normale to geomFace at any node
 425   bool geomNormalOK = false;
 426   gp_Vec geomNormal;
 427   SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
 428   for ( int i = 0; !geomNormalOK && i < 3; ++i )
 429   {
 430     // find UV of i-th node on geomFace
 431     const SMDS_MeshNode* nNotOnSeamEdge = 0;
 432     if ( helper.IsSeamShape( nodes[i]->getshapeId() )) {
 433       if ( helper.IsSeamShape( nodes[(i+1)%3]->getshapeId() ))
 434         nNotOnSeamEdge = nodes[(i+2)%3];
 435       else
 436         nNotOnSeamEdge = nodes[(i+1)%3];
 437     }
 438     bool uvOK = true;
 439     gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK );
 440     // check that uv is correct
 441     if (uvOK) {
 442       double tol = 1e-6;
 443       TopoDS_Shape nodeShape = helper.GetSubShapeByNode( nodes[i], meshDS );
 444       if ( !nodeShape.IsNull() )
 445         switch ( nodeShape.ShapeType() )
 446         {
 447         case TopAbs_FACE:   tol = BRep_Tool::Tolerance( TopoDS::Face( nodeShape )); break;
 448         case TopAbs_EDGE:   tol = BRep_Tool::Tolerance( TopoDS::Edge( nodeShape )); break;
 449         case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
 450         default:;
 451         }
 452       gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
 453       BRepAdaptor_Surface surface( geomFace );
 454       uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
 455       if ( uvOK ) {
 456         // normale to geomFace at UV
 457         gp_Vec du, dv;
 458         surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
 459         geomNormal = du ^ dv;
 460         if ( geomFace.Orientation() == TopAbs_REVERSED )
 461           geomNormal.Reverse();
 462         geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
 463       }
 464     }
 465   }
 466   if ( !geomNormalOK)
 467     return invalidID;
 468
 469   // compare normals
 470   bool isReverse = ( meshNormal * geomNormal ) < 0;
 471   if ( !isReverse )
 472     return meshDS->ShapeToIndex( solid1 );
 473
 474   if ( solids.Extent() == 1 )
 475     return HOLE_ID; // we are inside a hole
 476   else
 477     return meshDS->ShapeToIndex( solids(2) );
 478 }
 479
 480 //=======================================================================
 481 //function : addElemInMeshGroup
 482 //purpose  : Update or create groups in mesh
 483 //=======================================================================
 484
 485 static void addElemInMeshGroup(SMESH_Mesh*             theMesh,
 486                                const SMDS_MeshElement* anElem,
 487                                std::string&            groupName,
 488                                std::set<std::string>&  groupsToRemove)
 489 {
 490   if ( !anElem ) return; // issue 0021776
 491
 492   bool groupDone = false;
 493   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
 494   while (grIt->more()) {
 495     SMESH_Group * group = grIt->next();
 496     if ( !group ) continue;
 497     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
 498     if ( !groupDS ) continue;
 499     if ( groupDS->GetType()==anElem->GetType() &&groupName.compare(group->GetName())==0) {
 500       SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( groupDS );
 501       aGroupDS->SMDSGroup().Add(anElem);
 502       groupDone = true;
 503 //       MESSAGE("Successfully added enforced element to existing group " << groupName);
 504       break;
 505     }
 506   }
 507
 508   if (!groupDone)
 509   {
 510     int groupId;
 511     SMESH_Group* aGroup = theMesh->AddGroup(anElem->GetType(), groupName.c_str(), groupId);
 512     aGroup->SetName( groupName.c_str() );
 513     SMESHDS_Group* aGroupDS = static_cast<SMESHDS_Group*>( aGroup->GetGroupDS() );
 514     aGroupDS->SMDSGroup().Add(anElem);
 515 //     MESSAGE("Successfully created enforced vertex group " << groupName);
 516     groupDone = true;
 517   }
 518   if (!groupDone)
 519     throw SALOME_Exception(LOCALIZED("A given element was not added to a group"));
 520 }
 521
 522
 523 //=======================================================================
 524 //function : updateMeshGroups
 525 //purpose  : Update or create groups in mesh
 526 //=======================================================================
 527
 528 static void updateMeshGroups(SMESH_Mesh* theMesh, std::set<std::string> groupsToRemove)
 529 {
 530   SMESH_Mesh::GroupIteratorPtr grIt = theMesh->GetGroups();
 531   while (grIt->more()) {
 532     SMESH_Group * group = grIt->next();
 533     if ( !group ) continue;
 534     SMESHDS_GroupBase* groupDS = group->GetGroupDS();
 535     if ( !groupDS ) continue;
 536     std::string currentGroupName = (string)group->GetName();
 537     if (groupDS->IsEmpty() && groupsToRemove.find(currentGroupName) != groupsToRemove.end()) {
 538       // Previous group created by enforced elements
 539       MESSAGE("Delete previous group created by removed enforced elements: " << group->GetName())
 540       theMesh->RemoveGroup(groupDS->GetID());
 541     }
 542   }
 543 }
 544
 545 //=======================================================================
 546 //function : removeEmptyGroupsOfDomains
 547 //purpose  : remove empty groups named "Domain_nb" created due to
 548 //           "To make groups of domains" option.
 549 //=======================================================================
 550
 551 static void removeEmptyGroupsOfDomains(SMESH_Mesh* mesh,
 552                                        bool        notEmptyAsWell = false)
 553 {
 554   const char* refName = theDomainGroupNamePrefix;
 555   const size_t refLen = strlen( theDomainGroupNamePrefix );
 556
 557   std::list<int> groupIDs = mesh->GetGroupIds();
 558   std::list<int>::const_iterator id = groupIDs.begin();
 559   for ( ; id != groupIDs.end(); ++id )
 560   {
 561     SMESH_Group* group = mesh->GetGroup( *id );
 562     if ( !group || ( !group->GetGroupDS()->IsEmpty() && !notEmptyAsWell ))
 563       continue;
 564     const char* name = group->GetName();
 565     char* end;
 566     // check the name
 567     if ( strncmp( name, refName, refLen ) == 0 && // starts from refName;
 568          isdigit( *( name + refLen )) &&          // refName is followed by a digit;
 569          strtol( name + refLen, &end, 10) >= 0 && // there are only digits ...
 570          *end == '\0')                            // ... till a string end.
 571     {
 572       mesh->RemoveGroup( *id );
 573     }
 574   }
 575 }
 576
 577 //================================================================================
 578 /*!
 579  * \brief Create the groups corresponding to domains
 580  */
 581 //================================================================================
 582
 583 static void makeDomainGroups( std::vector< std::vector< const SMDS_MeshElement* > >& elemsOfDomain,
 584                               SMESH_MesherHelper*                                    theHelper)
 585 {
 586   // int nbDomains = 0;
 587   // for ( size_t i = 0; i < elemsOfDomain.size(); ++i )
 588   //   nbDomains += ( elemsOfDomain[i].size() > 0 );
 589
 590   // if ( nbDomains > 1 )
 591   for ( size_t iDomain = 0; iDomain < elemsOfDomain.size(); ++iDomain )
 592   {
 593     std::vector< const SMDS_MeshElement* > & elems = elemsOfDomain[ iDomain ];
 594     if ( elems.empty() ) continue;
 595
 596     // find existing groups
 597     std::vector< SMESH_Group* > groupOfType( SMDSAbs_NbElementTypes, (SMESH_Group*)NULL );
 598     const std::string domainName = ( SMESH_Comment( theDomainGroupNamePrefix ) << iDomain );
 599     SMESH_Mesh::GroupIteratorPtr groupIt = theHelper->GetMesh()->GetGroups();
 600     while ( groupIt->more() )
 601     {
 602       SMESH_Group* group = groupIt->next();
 603       if ( domainName == group->GetName() &&
 604            dynamic_cast< SMESHDS_Group* >( group->GetGroupDS()) )
 605         groupOfType[ group->GetGroupDS()->GetType() ] = group;
 606     }
 607     // create and fill the groups
 608     size_t iElem = 0;
 609     int groupID;
 610     do
 611     {
 612       SMESH_Group* group = groupOfType[ elems[ iElem ]->GetType() ];
 613       if ( !group )
 614         group = theHelper->GetMesh()->AddGroup( elems[ iElem ]->GetType(),
 615                                                 domainName.c_str(), groupID );
 616       SMDS_MeshGroup& groupDS =
 617         static_cast< SMESHDS_Group* >( group->GetGroupDS() )->SMDSGroup();
 618
 619       while ( iElem < elems.size() && groupDS.Add( elems[iElem] ))
 620         ++iElem;
 621
 622     } while ( iElem < elems.size() );
 623   }
 624 }
 625
 626 //=======================================================================
 627 //function : readGMFFile
 628 //purpose  : read GMF file w/o geometry associated to mesh
 629 //=======================================================================
 630
 631 static bool readGMFFile(const char*                     theFile,
 632                         GHS3DPlugin_GHS3D*              theAlgo,
 633                         SMESH_MesherHelper*             theHelper,
 634                         std::vector <const SMDS_MeshNode*> &    theNodeByGhs3dId,
 635                         std::vector <const SMDS_MeshElement*> & theFaceByGhs3dId,
 636                         map<const SMDS_MeshNode*,int> & theNodeToGhs3dIdMap,
 637                         std::vector<std::string> &      aNodeGroupByGhs3dId,
 638                         std::vector<std::string> &      anEdgeGroupByGhs3dId,
 639                         std::vector<std::string> &      aFaceGroupByGhs3dId,
 640                         std::set<std::string> &         groupsToRemove,
 641                         bool                            toMakeGroupsOfDomains=false,
 642                         bool                            toMeshHoles=true)
 643 {
 644   std::string tmpStr;
 645   SMESHDS_Mesh* theMeshDS = theHelper->GetMeshDS();
 646   const bool hasGeom = ( theHelper->GetMesh()->HasShapeToMesh() );
 647
 648   int nbInitialNodes = theNodeByGhs3dId.size();
 649   int nbMeshNodes = theMeshDS->NbNodes();
 650
 651   const bool isQuadMesh =
 652     theHelper->GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
 653     theHelper->GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
 654     theHelper->GetMesh()->NbVolumes( ORDER_QUADRATIC );
 655
 656 #ifdef _DEBUG_
 657   std::cout << "theNodeByGhs3dId.size(): " << nbInitialNodes << std::endl;
 658   std::cout << "theHelper->GetMesh()->NbNodes(): " << nbMeshNodes << std::endl;
 659   std::cout << "isQuadMesh: " << isQuadMesh << std::endl;
 660 #endif
 661
 662   // ---------------------------------
 663   // Read generated elements and nodes
 664   // ---------------------------------
 665
 666   int nbElem = 0, nbRef = 0;
 667   int aGMFNodeID = 0;
 668   const SMDS_MeshNode** GMFNode;
 669 #ifdef _DEBUG_
 670   std::map<int, std::set<int> > subdomainId2tetraId;
 671 #endif
 672   std::map <GmfKwdCod,int> tabRef;
 673   const bool force3d = !hasGeom;
 674   const int  noID    = 0;
 675
 676   tabRef[GmfVertices]       = 3; // for new nodes and enforced nodes
 677   tabRef[GmfCorners]        = 1;
 678   tabRef[GmfEdges]          = 2; // for enforced edges
 679   tabRef[GmfRidges]         = 1;
 680   tabRef[GmfTriangles]      = 3; // for enforced faces
 681   tabRef[GmfQuadrilaterals] = 4;
 682   tabRef[GmfTetrahedra]     = 4; // for new tetras
 683   tabRef[GmfHexahedra]      = 8;
 684
 685   int ver, dim;
 686   MESSAGE("Read " << theFile << " file");
 687   int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim);
 688   if (!InpMsh)
 689     return false;
 690   MESSAGE("Done ");
 691
 692   // Read ids of domains
 693   vector< int > solidIDByDomain;
 694   if ( hasGeom )
 695   {
 696     int solid1; // id used in case of 1 domain or some reading failure
 697     if ( theHelper->GetSubShape().ShapeType() == TopAbs_SOLID )
 698       solid1 = theHelper->GetSubShapeID();
 699     else
 700       solid1 = theMeshDS->ShapeToIndex
 701         ( TopExp_Explorer( theHelper->GetSubShape(), TopAbs_SOLID ).Current() );
 702
 703     int nbDomains = GmfStatKwd( InpMsh, GmfSubDomainFromGeom );
 704     if ( nbDomains > 1 )
 705     {
 706       solidIDByDomain.resize( nbDomains+1, theHelper->GetSubShapeID() );
 707       int faceNbNodes, faceIndex, orientation, domainNb;
 708       GmfGotoKwd( InpMsh, GmfSubDomainFromGeom );
 709       for ( int i = 0; i < nbDomains; ++i )
 710       {
 711         faceIndex = 0;
 712         GmfGetLin( InpMsh, GmfSubDomainFromGeom,
 713                    &faceNbNodes, &faceIndex, &orientation, &domainNb);
 714         solidIDByDomain[ domainNb ] = 1;
 715         if ( 0 < faceIndex && faceIndex-1 < (int)theFaceByGhs3dId.size() )
 716         {
 717           const SMDS_MeshElement* face = theFaceByGhs3dId[ faceIndex-1 ];
 718           const SMDS_MeshNode* nn[3] = { face->GetNode(0),
 719                                          face->GetNode(1),
 720                                          face->GetNode(2) };
 721           if ( orientation < 0 )
 722             std::swap( nn[1], nn[2] );
 723           solidIDByDomain[ domainNb ] =
 724             findShapeID( *theHelper->GetMesh(), nn[0], nn[1], nn[2], toMeshHoles );
 725           if ( solidIDByDomain[ domainNb ] > 0 )
 726           {
 727 #ifdef _DEBUG_
 728             std::cout << "solid " << solidIDByDomain[ domainNb ] << std::endl;
 729 #endif
 730             const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( solidIDByDomain[ domainNb ] );
 731             if ( ! theHelper->IsSubShape( foundShape, theHelper->GetSubShape() ))
 732               solidIDByDomain[ domainNb ] = HOLE_ID;
 733           }
 734         }
 735       }
 736     }
 737     if ( solidIDByDomain.size() < 2 )
 738       solidIDByDomain.resize( 2, solid1 );
 739   }
 740
 741   // Issue 0020682. Avoid creating nodes and tetras at place where
 742   // volumic elements already exist
 743   SMESH_ElementSearcher* elemSearcher = 0;
 744   std::vector< const SMDS_MeshElement* > foundVolumes;
 745   if ( !hasGeom && theHelper->GetMesh()->NbVolumes() > 0 )
 746     elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *theMeshDS );
 747   auto_ptr< SMESH_ElementSearcher > elemSearcherDeleter( elemSearcher );
 748
 749   // IMP 0022172: [CEA 790] create the groups corresponding to domains
 750   std::vector< std::vector< const SMDS_MeshElement* > > elemsOfDomain;
 751
 752   int nbVertices = GmfStatKwd(InpMsh, GmfVertices) - nbInitialNodes;
 753   GMFNode = new const SMDS_MeshNode*[ nbVertices + 1 ];
 754
 755   std::map <GmfKwdCod,int>::const_iterator it = tabRef.begin();
 756   for ( ; it != tabRef.end() ; ++it)
 757   {
 758     if(theAlgo->computeCanceled()) {
 759       GmfCloseMesh(InpMsh);
 760       delete [] GMFNode;
 761       return false;
 762     }
 763     int dummy, solidID;
 764     GmfKwdCod token = it->first;
 765     nbRef           = it->second;
 766
 767     nbElem = GmfStatKwd(InpMsh, token);
 768     if (nbElem > 0) {
 769       GmfGotoKwd(InpMsh, token);
 770       std::cout << "Read " << nbElem;
 771     }
 772     else
 773       continue;
 774
 775     std::vector<int> id (nbElem*tabRef[token]); // node ids
 776     std::vector<int> domainID( nbElem ); // domain
 777
 778     if (token == GmfVertices) {
 779       (nbElem <= 1) ? tmpStr = " vertex" : tmpStr = " vertices";
 780 //       std::cout << nbInitialNodes << " from input mesh " << std::endl;
 781
 782       // Remove orphan nodes from previous enforced mesh which was cleared
 783 //       if ( nbElem < nbMeshNodes ) {
 784 //         const SMDS_MeshNode* node;
 785 //         SMDS_NodeIteratorPtr nodeIt = theMeshDS->nodesIterator();
 786 //         while ( nodeIt->more() )
 787 //         {
 788 //           node = nodeIt->next();
 789 //           if (theNodeToGhs3dIdMap.find(node) != theNodeToGhs3dIdMap.end())
 790 //             theMeshDS->RemoveNode(node);
 791 //         }
 792 //       }
 793
 794
 795       int aGMFID;
 796
 797       float VerTab_f[3];
 798       double x, y, z;
 799       const SMDS_MeshNode * aGMFNode;
 800
 801       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 802         if(theAlgo->computeCanceled()) {
 803           GmfCloseMesh(InpMsh);
 804           delete [] GMFNode;
 805           return false;
 806         }
 807         if (ver == GmfFloat) {
 808           GmfGetLin(InpMsh, token, &VerTab_f[0], &VerTab_f[1], &VerTab_f[2], &dummy);
 809           x = VerTab_f[0];
 810           y = VerTab_f[1];
 811           z = VerTab_f[2];
 812         }
 813         else {
 814           GmfGetLin(InpMsh, token, &x, &y, &z, &dummy);
 815         }
 816         if (iElem >= nbInitialNodes) {
 817           if ( elemSearcher &&
 818                 elemSearcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_Volume, foundVolumes))
 819             aGMFNode = 0;
 820           else
 821             aGMFNode = theHelper->AddNode(x, y, z);
 822
 823           aGMFID = iElem -nbInitialNodes +1;
 824           GMFNode[ aGMFID ] = aGMFNode;
 825           if (aGMFID-1 < (int)aNodeGroupByGhs3dId.size() && !aNodeGroupByGhs3dId.at(aGMFID-1).empty())
 826             addElemInMeshGroup(theHelper->GetMesh(), aGMFNode, aNodeGroupByGhs3dId.at(aGMFID-1), groupsToRemove);
 827         }
 828       }
 829     }
 830     else if (token == GmfCorners && nbElem > 0) {
 831       (nbElem <= 1) ? tmpStr = " corner" : tmpStr = " corners";
 832       for ( int iElem = 0; iElem < nbElem; iElem++ )
 833         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 834     }
 835     else if (token == GmfRidges && nbElem > 0) {
 836       (nbElem <= 1) ? tmpStr = " ridge" : tmpStr = " ridges";
 837       for ( int iElem = 0; iElem < nbElem; iElem++ )
 838         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]);
 839     }
 840     else if (token == GmfEdges && nbElem > 0) {
 841       (nbElem <= 1) ? tmpStr = " edge" : tmpStr = " edges";
 842       for ( int iElem = 0; iElem < nbElem; iElem++ )
 843         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &domainID[iElem]);
 844     }
 845     else if (token == GmfTriangles && nbElem > 0) {
 846       (nbElem <= 1) ? tmpStr = " triangle" : tmpStr = " triangles";
 847       for ( int iElem = 0; iElem < nbElem; iElem++ )
 848         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &domainID[iElem]);
 849     }
 850     else if (token == GmfQuadrilaterals && nbElem > 0) {
 851       (nbElem <= 1) ? tmpStr = " Quadrilateral" : tmpStr = " Quadrilaterals";
 852       for ( int iElem = 0; iElem < nbElem; iElem++ )
 853         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]);
 854     }
 855     else if (token == GmfTetrahedra && nbElem > 0) {
 856       (nbElem <= 1) ? tmpStr = " Tetrahedron" : tmpStr = " Tetrahedra";
 857       for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 858         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]);
 859 #ifdef _DEBUG_
 860         subdomainId2tetraId[dummy].insert(iElem+1);
 861 //         MESSAGE("subdomainId2tetraId["<<dummy<<"].insert("<<iElem+1<<")");
 862 #endif
 863       }
 864     }
 865     else if (token == GmfHexahedra && nbElem > 0) {
 866       (nbElem <= 1) ? tmpStr = " Hexahedron" : tmpStr = " Hexahedra";
 867       for ( int iElem = 0; iElem < nbElem; iElem++ )
 868         GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3],
 869                   &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &domainID[iElem]);
 870     }
 871     std::cout << tmpStr << std::endl;
 872     std::cout << std::endl;
 873
 874     switch (token) {
 875     case GmfCorners:
 876     case GmfRidges:
 877     case GmfEdges:
 878     case GmfTriangles:
 879     case GmfQuadrilaterals:
 880     case GmfTetrahedra:
 881     case GmfHexahedra:
 882     {
 883       std::vector< const SMDS_MeshNode* > node( nbRef );
 884       std::vector< int >          nodeID( nbRef );
 885       std::vector< SMDS_MeshNode* > enfNode( nbRef );
 886       const SMDS_MeshElement* aCreatedElem;
 887
 888       for ( int iElem = 0; iElem < nbElem; iElem++ )
 889       {
 890         if(theAlgo->computeCanceled()) {
 891           GmfCloseMesh(InpMsh);
 892           delete [] GMFNode;
 893           return false;
 894         }
 895         // Check if elem is already in input mesh. If yes => skip
 896         bool fullyCreatedElement = false; // if at least one of the nodes was created
 897         for ( int iRef = 0; iRef < nbRef; iRef++ )
 898         {
 899           aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID
 900           if (aGMFNodeID <= nbInitialNodes) // input nodes
 901           {
 902             aGMFNodeID--;
 903             node[ iRef ] = theNodeByGhs3dId[aGMFNodeID];
 904           }
 905           else
 906           {
 907             fullyCreatedElement = true;
 908             aGMFNodeID -= nbInitialNodes;
 909             nodeID[ iRef ] = aGMFNodeID ;
 910             node  [ iRef ] = GMFNode[ aGMFNodeID ];
 911           }
 912         }
 913         aCreatedElem = 0;
 914         switch (token)
 915         {
 916         case GmfEdges:
 917           if (fullyCreatedElement) {
 918             aCreatedElem = theHelper->AddEdge( node[0], node[1], noID, force3d );
 919             if (anEdgeGroupByGhs3dId.size() && !anEdgeGroupByGhs3dId[iElem].empty())
 920               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, anEdgeGroupByGhs3dId[iElem], groupsToRemove);
 921           }
 922           break;
 923         case GmfTriangles:
 924           if (fullyCreatedElement) {
 925             aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], noID, force3d );
 926             if (aFaceGroupByGhs3dId.size() && !aFaceGroupByGhs3dId[iElem].empty())
 927               addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, aFaceGroupByGhs3dId[iElem], groupsToRemove);
 928           }
 929           break;
 930         case GmfQuadrilaterals:
 931           if (fullyCreatedElement) {
 932             aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], node[3], noID, force3d );
 933           }
 934           break;
 935         case GmfTetrahedra:
 936           if ( hasGeom )
 937           {
 938             solidID = solidIDByDomain[ domainID[iElem]];
 939             if ( solidID != HOLE_ID )
 940             {
 941               aCreatedElem = theHelper->AddVolume( node[1], node[0], node[2], node[3],
 942                                                    noID, force3d );
 943               theMeshDS->SetMeshElementOnShape( aCreatedElem, solidID );
 944               for ( int iN = 0; iN < 4; ++iN )
 945                 if ( node[iN]->getshapeId() < 1 )
 946                   theMeshDS->SetNodeInVolume( node[iN], solidID );
 947             }
 948           }
 949           else
 950           {
 951             if ( elemSearcher ) {
 952               // Issue 0020682. Avoid creating nodes and tetras at place where
 953               // volumic elements already exist
 954               if ( !node[1] || !node[0] || !node[2] || !node[3] )
 955                 continue;
 956               if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
 957                                                       SMESH_TNodeXYZ(node[1]) +
 958                                                       SMESH_TNodeXYZ(node[2]) +
 959                                                       SMESH_TNodeXYZ(node[3]) ) / 4.,
 960                                                      SMDSAbs_Volume, foundVolumes ))
 961                 break;
 962             }
 963             aCreatedElem = theHelper->AddVolume( node[1], node[0], node[2], node[3],
 964                                                  noID, force3d );
 965           }
 966           break;
 967         case GmfHexahedra:
 968           if ( hasGeom )
 969           {
 970             solidID = solidIDByDomain[ domainID[iElem]];
 971             if ( solidID != HOLE_ID )
 972             {
 973               aCreatedElem = theHelper->AddVolume( node[0], node[3], node[2], node[1],
 974                                                    node[4], node[7], node[6], node[5],
 975                                                    noID, force3d );
 976               theMeshDS->SetMeshElementOnShape( aCreatedElem, solidID );
 977               for ( int iN = 0; iN < 8; ++iN )
 978                 if ( node[iN]->getshapeId() < 1 )
 979                   theMeshDS->SetNodeInVolume( node[iN], solidID );
 980             }
 981           }
 982           else
 983           {
 984             if ( elemSearcher ) {
 985               // Issue 0020682. Avoid creating nodes and tetras at place where
 986               // volumic elements already exist
 987               if ( !node[1] || !node[0] || !node[2] || !node[3] || !node[4] || !node[5] || !node[6] || !node[7])
 988                 continue;
 989               if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) +
 990                                                       SMESH_TNodeXYZ(node[1]) +
 991                                                       SMESH_TNodeXYZ(node[2]) +
 992                                                       SMESH_TNodeXYZ(node[3]) +
 993                                                       SMESH_TNodeXYZ(node[4]) +
 994                                                       SMESH_TNodeXYZ(node[5]) +
 995                                                       SMESH_TNodeXYZ(node[6]) +
 996                                                       SMESH_TNodeXYZ(node[7])) / 8.,
 997                                                      SMDSAbs_Volume, foundVolumes ))
 998                 break;
 999             }
1000             aCreatedElem = theHelper->AddVolume( node[0], node[3], node[2], node[1],
1001                                                  node[4], node[7], node[6], node[5],
1002                                                  noID, force3d );
1003           }
1004           break;
1005         default: continue;
1006         } // switch (token)
1007
1008         if ( aCreatedElem && toMakeGroupsOfDomains )
1009         {
1010           if ( domainID[iElem] >= (int) elemsOfDomain.size() )
1011             elemsOfDomain.resize( domainID[iElem] + 1 );
1012           elemsOfDomain[ domainID[iElem] ].push_back( aCreatedElem );
1013         }
1014       } // loop on elements of one type
1015       break;
1016     } // case ...
1017     default:;
1018     } // switch (token)
1019   } // loop on tabRef
1020
1021   // remove nodes in holes
1022   if ( hasGeom )
1023   {
1024     for ( int i = 1; i <= nbVertices; ++i )
1025       if ( GMFNode[i]->NbInverseElements() == 0 )
1026         theMeshDS->RemoveFreeNode( GMFNode[i], /*sm=*/0, /*fromGroups=*/false );
1027   }
1028
1029   GmfCloseMesh(InpMsh);
1030   delete [] GMFNode;
1031
1032   // 0022172: [CEA 790] create the groups corresponding to domains
1033   if ( toMakeGroupsOfDomains )
1034     makeDomainGroups( elemsOfDomain, theHelper );
1035
1036 #ifdef _DEBUG_
1037   MESSAGE("Nb subdomains " << subdomainId2tetraId.size());
1038   std::map<int, std::set<int> >::const_iterator subdomainIt = subdomainId2tetraId.begin();
1039   TCollection_AsciiString aSubdomainFileName = theFile;
1040   aSubdomainFileName = aSubdomainFileName + ".subdomain";
1041   ofstream aSubdomainFile  ( aSubdomainFileName.ToCString()  , ios::out);
1042
1043   aSubdomainFile << "Nb subdomains " << subdomainId2tetraId.size() << std::endl;
1044   for(;subdomainIt != subdomainId2tetraId.end() ; ++subdomainIt) {
1045     int subdomainId = subdomainIt->first;
1046     std::set<int> tetraIds = subdomainIt->second;
1047     MESSAGE("Subdomain #"<<subdomainId<<": "<<tetraIds.size()<<" tetrahedrons");
1048     std::set<int>::const_iterator tetraIdsIt = tetraIds.begin();
1049     aSubdomainFile << subdomainId << std::endl;
1050     for(;tetraIdsIt != tetraIds.end() ; ++tetraIdsIt) {
1051       aSubdomainFile << (*tetraIdsIt) << " ";
1052     }
1053     aSubdomainFile << std::endl;
1054   }
1055   aSubdomainFile.close();
1056 #endif
1057
1058   return true;
1059 }
1060
1061
1062 static bool writeGMFFile(const char*                                     theMeshFileName,
1063                          const char*                                     theRequiredFileName,
1064                          const char*                                     theSolFileName,
1065                          const SMESH_ProxyMesh&                          theProxyMesh,
1066                          SMESH_MesherHelper&                             theHelper,
1067                          std::vector <const SMDS_MeshNode*> &            theNodeByGhs3dId,
1068                          std::vector <const SMDS_MeshElement*> &         theFaceByGhs3dId,
1069                          std::map<const SMDS_MeshNode*,int> &            aNodeToGhs3dIdMap,
1070                          std::vector<std::string> &                      aNodeGroupByGhs3dId,
1071                          std::vector<std::string> &                      anEdgeGroupByGhs3dId,
1072                          std::vector<std::string> &                      aFaceGroupByGhs3dId,
1073                          GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap & theEnforcedNodes,
1074                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges,
1075                          GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles,
1076                          std::map<std::vector<double>, std::string> &    enfVerticesWithGroup,
1077                          GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices,
1078                          int &                                           theInvalidEnforcedFlags)
1079 {
1080   MESSAGE("writeGMFFile w/o geometry");
1081   std::string tmpStr;
1082   int idx, idxRequired = 0, idxSol = 0;
1083   const int dummyint = 0;
1084   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt;
1085   std::vector<double> enfVertexSizes;
1086   const SMDS_MeshElement* elem;
1087   TIDSortedElemSet anElemSet, theKeptEnforcedEdges, theKeptEnforcedTriangles;
1088   SMDS_ElemIteratorPtr nodeIt;
1089   std::vector <const SMDS_MeshNode*> theEnforcedNodeByGhs3dId;
1090   map<const SMDS_MeshNode*,int> anEnforcedNodeToGhs3dIdMap, anExistingEnforcedNodeToGhs3dIdMap;
1091   std::vector< const SMDS_MeshElement* > foundElems;
1092   map<const SMDS_MeshNode*,TopAbs_State> aNodeToTopAbs_StateMap;
1093   int nbFoundElems;
1094   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::iterator elemIt;
1095   TIDSortedElemSet::iterator elemSetIt;
1096   bool isOK;
1097   SMESH_Mesh* theMesh = theHelper.GetMesh();
1098   const bool hasGeom = theMesh->HasShapeToMesh();
1099   SMESHUtils::Deleter< SMESH_ElementSearcher > pntCls
1100     ( SMESH_MeshAlgos::GetElementSearcher(*theMesh->GetMeshDS()));
1101
1102   int nbEnforcedVertices = theEnforcedVertices.size();
1103   theInvalidEnforcedFlags = 0;
1104
1105   // count faces
1106   int nbFaces = theProxyMesh.NbFaces();
1107   int nbNodes;
1108   theFaceByGhs3dId.reserve( nbFaces );
1109
1110   // groups management
1111   int usedEnforcedNodes = 0;
1112   std::string gn = "";
1113
1114   if ( nbFaces == 0 )
1115     return false;
1116
1117   idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1118   if (!idx)
1119     return false;
1120
1121   /* ========================== FACES ========================== */
1122   /* TRIANGLES ========================== */
1123   SMDS_ElemIteratorPtr eIt =
1124     hasGeom ? theProxyMesh.GetFaces( theHelper.GetSubShape()) : theProxyMesh.GetFaces();
1125   while ( eIt->more() )
1126   {
1127     elem = eIt->next();
1128     anElemSet.insert(elem);
1129     nodeIt = elem->nodesIterator();
1130     nbNodes = elem->NbCornerNodes();
1131     while ( nodeIt->more() && nbNodes--)
1132     {
1133       // find MG-Tetra ID
1134       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1135       int newId = aNodeToGhs3dIdMap.size() + 1; // MG-Tetra ids count from 1
1136       aNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1137     }
1138   }
1139
1140   /* EDGES ========================== */
1141
1142   // Iterate over the enforced edges
1143   for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) {
1144     elem = elemIt->first;
1145     isOK = true;
1146     nodeIt = elem->nodesIterator();
1147     nbNodes = 2;
1148     while ( nodeIt->more() && nbNodes-- ) {
1149       // find MG-Tetra ID
1150       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1151       // Test if point is inside shape to mesh
1152       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1153       TopAbs_State result = pntCls->GetPointState( myPoint );
1154       if ( result == TopAbs_OUT ) {
1155         isOK = false;
1156         theInvalidEnforcedFlags |= FLAG_BAD_ENF_EDGE;
1157         break;
1158       }
1159       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1160     }
1161     if (isOK) {
1162       nodeIt = elem->nodesIterator();
1163       nbNodes = 2;
1164       int newId = -1;
1165       while ( nodeIt->more() && nbNodes-- ) {
1166         // find MG-Tetra ID
1167         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1168         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1169         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1170 #ifdef _DEBUG_
1171         std::cout << "Node at "<<node->X()<<", "<<node->Y()<<", "<<node->Z()<<std::endl;
1172         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1173 #endif
1174         if (nbFoundElems ==0) {
1175           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1176             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // MG-Tetra ids count from 1
1177             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1178           }
1179         }
1180         else if (nbFoundElems ==1) {
1181           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1182           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1183           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1184         }
1185         else
1186           isOK = false;
1187 #ifdef _DEBUG_
1188         std::cout << "MG-Tetra node ID: "<<newId<<std::endl;
1189 #endif
1190       }
1191       if (isOK)
1192         theKeptEnforcedEdges.insert(elem);
1193       else
1194         theInvalidEnforcedFlags |= FLAG_BAD_ENF_EDGE;
1195     }
1196   }
1197
1198   /* ENFORCED TRIANGLES ========================== */
1199
1200   // Iterate over the enforced triangles
1201   for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) {
1202     elem = elemIt->first;
1203     isOK = true;
1204     nodeIt = elem->nodesIterator();
1205     nbNodes = 3;
1206     while ( nodeIt->more() && nbNodes--) {
1207       // find MG-Tetra ID
1208       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1209       // Test if point is inside shape to mesh
1210       gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1211       TopAbs_State result = pntCls->GetPointState( myPoint );
1212       if ( result == TopAbs_OUT ) {
1213         isOK = false;
1214         theInvalidEnforcedFlags |= FLAG_BAD_ENF_TRIA;
1215         break;
1216       }
1217       aNodeToTopAbs_StateMap.insert( make_pair( node, result ));
1218     }
1219     if (isOK) {
1220       nodeIt = elem->nodesIterator();
1221       nbNodes = 3;
1222       int newId = -1;
1223       while ( nodeIt->more() && nbNodes--) {
1224         // find MG-Tetra ID
1225         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1226         gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1227         nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1228 #ifdef _DEBUG_
1229         std::cout << "Nb nodes found : "<<nbFoundElems<<std::endl;
1230 #endif
1231         if (nbFoundElems ==0) {
1232           if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1233             newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // MG-Tetra ids count from 1
1234             anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1235           }
1236         }
1237         else if (nbFoundElems ==1) {
1238           const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1239           newId = (*aNodeToGhs3dIdMap.find(existingNode)).second;
1240           anExistingEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId ));
1241         }
1242         else
1243           isOK = false;
1244 #ifdef _DEBUG_
1245         std::cout << "MG-Tetra node ID: "<<newId<<std::endl;
1246 #endif
1247       }
1248       if (isOK)
1249         theKeptEnforcedTriangles.insert(elem);
1250       else
1251         theInvalidEnforcedFlags |= FLAG_BAD_ENF_TRIA;
1252     }
1253   }
1254
1255   // put nodes to theNodeByGhs3dId vector
1256 #ifdef _DEBUG_
1257   std::cout << "aNodeToGhs3dIdMap.size(): "<<aNodeToGhs3dIdMap.size()<<std::endl;
1258 #endif
1259   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
1260   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
1261   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
1262   {
1263 //     std::cout << "n2id->first: "<<n2id->first<<std::endl;
1264     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // MG-Tetra ids count from 1
1265   }
1266
1267   // put nodes to anEnforcedNodeToGhs3dIdMap vector
1268 #ifdef _DEBUG_
1269   std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<<anEnforcedNodeToGhs3dIdMap.size()<<std::endl;
1270 #endif
1271   theEnforcedNodeByGhs3dId.resize( anEnforcedNodeToGhs3dIdMap.size());
1272   n2id = anEnforcedNodeToGhs3dIdMap.begin();
1273   for ( ; n2id != anEnforcedNodeToGhs3dIdMap.end(); ++ n2id)
1274   {
1275     if (n2id->second > (int)aNodeToGhs3dIdMap.size()) {
1276       theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // MG-Tetra ids count from 1
1277     }
1278   }
1279
1280
1281   /* ========================== NODES ========================== */
1282   vector<const SMDS_MeshNode*> theOrderedNodes, theRequiredNodes;
1283   std::set< std::vector<double> > nodesCoords;
1284   vector<const SMDS_MeshNode*>::const_iterator ghs3dNodeIt = theNodeByGhs3dId.begin();
1285   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
1286
1287   (theNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1288   std::cout << theNodeByGhs3dId.size() << tmpStr << " from mesh ..." << std::endl;
1289   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1290   {
1291     const SMDS_MeshNode* node = *ghs3dNodeIt;
1292     std::vector<double> coords;
1293     coords.push_back(node->X());
1294     coords.push_back(node->Y());
1295     coords.push_back(node->Z());
1296     nodesCoords.insert(coords);
1297     theOrderedNodes.push_back(node);
1298   }
1299
1300   // Iterate over the enforced nodes given by enforced elements
1301   ghs3dNodeIt = theEnforcedNodeByGhs3dId.begin();
1302   after  = theEnforcedNodeByGhs3dId.end();
1303   (theEnforcedNodeByGhs3dId.size() <= 1) ? tmpStr = " node" : " nodes";
1304   std::cout << theEnforcedNodeByGhs3dId.size() << tmpStr << " from enforced elements ..." << std::endl;
1305   for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt )
1306   {
1307     const SMDS_MeshNode* node = *ghs3dNodeIt;
1308     std::vector<double> coords;
1309     coords.push_back(node->X());
1310     coords.push_back(node->Y());
1311     coords.push_back(node->Z());
1312 #ifdef _DEBUG_
1313     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1314 #endif
1315
1316     if (nodesCoords.find(coords) != nodesCoords.end()) {
1317       // node already exists in original mesh
1318 #ifdef _DEBUG_
1319       std::cout << " found" << std::endl;
1320 #endif
1321       continue;
1322     }
1323
1324     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1325       // node already exists in enforced vertices
1326 #ifdef _DEBUG_
1327       std::cout << " found" << std::endl;
1328 #endif
1329       continue;
1330     }
1331
1332 //     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1333 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1334 //     if (nbFoundElems ==0) {
1335 //       std::cout << " not found" << std::endl;
1336 //       if ((*aNodeToTopAbs_StateMap.find(node)).second == TopAbs_IN) {
1337 //         nodesCoords.insert(coords);
1338 //         theOrderedNodes.push_back(node);
1339 //       }
1340 //     }
1341 //     else {
1342 //       std::cout << " found in initial mesh" << std::endl;
1343 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1344 //       nodesCoords.insert(coords);
1345 //       theOrderedNodes.push_back(existingNode);
1346 //     }
1347
1348 #ifdef _DEBUG_
1349     std::cout << " not found" << std::endl;
1350 #endif
1351
1352     nodesCoords.insert(coords);
1353     theOrderedNodes.push_back(node);
1354 //     theRequiredNodes.push_back(node);
1355   }
1356
1357
1358   // Iterate over the enforced nodes
1359   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt;
1360   (theEnforcedNodes.size() <= 1) ? tmpStr = " node" : " nodes";
1361   std::cout << theEnforcedNodes.size() << tmpStr << " from enforced nodes ..." << std::endl;
1362   for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt)
1363   {
1364     const SMDS_MeshNode* node = enfNodeIt->first;
1365     std::vector<double> coords;
1366     coords.push_back(node->X());
1367     coords.push_back(node->Y());
1368     coords.push_back(node->Z());
1369 #ifdef _DEBUG_
1370     std::cout << "Node at " << node->X()<<", " <<node->Y()<<", " <<node->Z();
1371 #endif
1372
1373     // Test if point is inside shape to mesh
1374     gp_Pnt myPoint(node->X(),node->Y(),node->Z());
1375     TopAbs_State result = pntCls->GetPointState( myPoint );
1376     if ( result == TopAbs_OUT ) {
1377 #ifdef _DEBUG_
1378       std::cout << " out of volume" << std::endl;
1379 #endif
1380       theInvalidEnforcedFlags |= FLAG_BAD_ENF_NODE;
1381       continue;
1382     }
1383
1384     if (nodesCoords.find(coords) != nodesCoords.end()) {
1385 #ifdef _DEBUG_
1386       std::cout << " found in nodesCoords" << std::endl;
1387 #endif
1388 //       theRequiredNodes.push_back(node);
1389       continue;
1390     }
1391
1392     if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) {
1393 #ifdef _DEBUG_
1394       std::cout << " found in theEnforcedVertices" << std::endl;
1395 #endif
1396       continue;
1397     }
1398
1399 //     nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems);
1400 //     if (nbFoundElems ==0) {
1401 //       std::cout << " not found" << std::endl;
1402 //       if (result == TopAbs_IN) {
1403 //         nodesCoords.insert(coords);
1404 //         theRequiredNodes.push_back(node);
1405 //       }
1406 //     }
1407 //     else {
1408 //       std::cout << " found in initial mesh" << std::endl;
1409 //       const SMDS_MeshNode* existingNode = (SMDS_MeshNode*) foundElems.at(0);
1410 // //       nodesCoords.insert(coords);
1411 //       theRequiredNodes.push_back(existingNode);
1412 //     }
1413 //
1414 //
1415 //
1416 //     if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0)
1417 //       continue;
1418
1419 //     if ( result != TopAbs_IN )
1420 //       continue;
1421
1422 #ifdef _DEBUG_
1423     std::cout << " not found" << std::endl;
1424 #endif
1425     nodesCoords.insert(coords);
1426 //     theOrderedNodes.push_back(node);
1427     theRequiredNodes.push_back(node);
1428   }
1429   int requiredNodes = theRequiredNodes.size();
1430
1431   int solSize = 0;
1432   std::vector<std::vector<double> > ReqVerTab;
1433   if (nbEnforcedVertices) {
1434 //    ReqVerTab.clear();
1435     (nbEnforcedVertices <= 1) ? tmpStr = " node" : " nodes";
1436     std::cout << nbEnforcedVertices << tmpStr << " from enforced vertices ..." << std::endl;
1437     // Iterate over the enforced vertices
1438     for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
1439       double x = vertexIt->first[0];
1440       double y = vertexIt->first[1];
1441       double z = vertexIt->first[2];
1442       // Test if point is inside shape to mesh
1443       gp_Pnt myPoint(x,y,z);
1444       TopAbs_State result = pntCls->GetPointState( myPoint );
1445       if ( result == TopAbs_OUT )
1446       {
1447         std::cout << "Warning: enforced vertex at ( " << x << "," << y << "," << z << " ) is out of the meshed domain!!!" << std::endl;
1448         theInvalidEnforcedFlags |= FLAG_BAD_ENF_VERT;
1449         //continue;
1450       }
1451       std::vector<double> coords;
1452       coords.push_back(x);
1453       coords.push_back(y);
1454       coords.push_back(z);
1455       ReqVerTab.push_back(coords);
1456       enfVertexSizes.push_back(vertexIt->second);
1457       solSize++;
1458     }
1459   }
1460
1461
1462   // GmfVertices
1463   std::cout << "Begin writting required nodes in GmfVertices" << std::endl;
1464   std::cout << "Nb vertices: " << theOrderedNodes.size() << std::endl;
1465   GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()/*+solSize*/);
1466   for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt) {
1467     GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
1468   }
1469
1470   std::cout << "End writting required nodes in GmfVertices" << std::endl;
1471
1472   if (requiredNodes + solSize) {
1473     std::cout << "Begin writting in req and sol file" << std::endl;
1474     aNodeGroupByGhs3dId.resize( requiredNodes + solSize );
1475     idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1476     if (!idxRequired) {
1477       GmfCloseMesh(idx);
1478       return false;
1479     }
1480     idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1481     if (!idxSol) {
1482       GmfCloseMesh(idx);
1483       if (idxRequired)
1484         GmfCloseMesh(idxRequired);
1485       return false;
1486     }
1487     int TypTab[] = {GmfSca};
1488     double ValTab[] = {0.0};
1489     GmfSetKwd(idxRequired, GmfVertices, requiredNodes + solSize);
1490     GmfSetKwd(idxSol, GmfSolAtVertices, requiredNodes + solSize, 1, TypTab);
1491 //     int usedEnforcedNodes = 0;
1492 //     std::string gn = "";
1493     for (ghs3dNodeIt = theRequiredNodes.begin();ghs3dNodeIt != theRequiredNodes.end();++ghs3dNodeIt) {
1494       GmfSetLin(idxRequired, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint);
1495       GmfSetLin(idxSol, GmfSolAtVertices, ValTab);
1496       if (theEnforcedNodes.find((*ghs3dNodeIt)) != theEnforcedNodes.end())
1497         gn = theEnforcedNodes.find((*ghs3dNodeIt))->second;
1498       aNodeGroupByGhs3dId[usedEnforcedNodes] = gn;
1499       usedEnforcedNodes++;
1500     }
1501
1502     for (int i=0;i<solSize;i++) {
1503       std::cout << ReqVerTab[i][0] <<" "<< ReqVerTab[i][1] << " "<< ReqVerTab[i][2] << std::endl;
1504 #ifdef _DEBUG_
1505       std::cout << "enfVertexSizes.at("<<i<<"): " << enfVertexSizes.at(i) << std::endl;
1506 #endif
1507       double solTab[] = {enfVertexSizes.at(i)};
1508       GmfSetLin(idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint);
1509       GmfSetLin(idxSol, GmfSolAtVertices, solTab);
1510       aNodeGroupByGhs3dId[usedEnforcedNodes] = enfVerticesWithGroup.find(ReqVerTab[i])->second;
1511 #ifdef _DEBUG_
1512       std::cout << "aNodeGroupByGhs3dId["<<usedEnforcedNodes<<"] = \""<<aNodeGroupByGhs3dId[usedEnforcedNodes]<<"\""<<std::endl;
1513 #endif
1514       usedEnforcedNodes++;
1515     }
1516     std::cout << "End writting in req and sol file" << std::endl;
1517   }
1518
1519   int nedge[2], ntri[3];
1520
1521   // GmfEdges
1522   int usedEnforcedEdges = 0;
1523   if (theKeptEnforcedEdges.size()) {
1524     anEdgeGroupByGhs3dId.resize( theKeptEnforcedEdges.size() );
1525 //    idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION);
1526 //    if (!idxRequired)
1527 //      return false;
1528     GmfSetKwd(idx, GmfEdges, theKeptEnforcedEdges.size());
1529 //    GmfSetKwd(idxRequired, GmfEdges, theKeptEnforcedEdges.size());
1530     for(elemSetIt = theKeptEnforcedEdges.begin() ; elemSetIt != theKeptEnforcedEdges.end() ; ++elemSetIt) {
1531       elem = (*elemSetIt);
1532       nodeIt = elem->nodesIterator();
1533       int index=0;
1534       while ( nodeIt->more() ) {
1535         // find MG-Tetra ID
1536         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1537         map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
1538         if (it == anEnforcedNodeToGhs3dIdMap.end()) {
1539           it = anExistingEnforcedNodeToGhs3dIdMap.find(node);
1540           if (it == anEnforcedNodeToGhs3dIdMap.end())
1541             throw "Node not found";
1542         }
1543         nedge[index] = it->second;
1544         index++;
1545       }
1546       GmfSetLin(idx, GmfEdges, nedge[0], nedge[1], dummyint);
1547       anEdgeGroupByGhs3dId[usedEnforcedEdges] = theEnforcedEdges.find(elem)->second;
1548 //      GmfSetLin(idxRequired, GmfEdges, nedge[0], nedge[1], dummyint);
1549       usedEnforcedEdges++;
1550     }
1551 //    GmfCloseMesh(idxRequired);
1552   }
1553
1554
1555   if (usedEnforcedEdges) {
1556     GmfSetKwd(idx, GmfRequiredEdges, usedEnforcedEdges);
1557     for (int enfID=1;enfID<=usedEnforcedEdges;enfID++) {
1558       GmfSetLin(idx, GmfRequiredEdges, enfID);
1559     }
1560   }
1561
1562   // GmfTriangles
1563   int usedEnforcedTriangles = 0;
1564   if (anElemSet.size()+theKeptEnforcedTriangles.size()) {
1565     aFaceGroupByGhs3dId.resize( anElemSet.size()+theKeptEnforcedTriangles.size() );
1566     GmfSetKwd(idx, GmfTriangles, anElemSet.size()+theKeptEnforcedTriangles.size());
1567     int k=0;
1568     for(elemSetIt = anElemSet.begin() ; elemSetIt != anElemSet.end() ; ++elemSetIt,++k) {
1569       elem = (*elemSetIt);
1570       theFaceByGhs3dId.push_back( elem );
1571       nodeIt = elem->nodesIterator();
1572       int index=0;
1573       for ( int j = 0; j < 3; ++j ) {
1574         // find MG-Tetra ID
1575         const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1576         map< const SMDS_MeshNode*,int >::iterator it = aNodeToGhs3dIdMap.find(node);
1577         if (it == aNodeToGhs3dIdMap.end())
1578           throw "Node not found";
1579         ntri[index] = it->second;
1580         index++;
1581       }
1582       GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
1583       aFaceGroupByGhs3dId[k] = "";
1584     }
1585     if ( !theHelper.GetMesh()->HasShapeToMesh() )
1586       SMESHUtils::FreeVector( theFaceByGhs3dId );
1587     if (theKeptEnforcedTriangles.size()) {
1588       for(elemSetIt = theKeptEnforcedTriangles.begin() ; elemSetIt != theKeptEnforcedTriangles.end() ; ++elemSetIt,++k) {
1589         elem = (*elemSetIt);
1590         nodeIt = elem->nodesIterator();
1591         int index=0;
1592         for ( int j = 0; j < 3; ++j ) {
1593           // find MG-Tetra ID
1594           const SMDS_MeshNode* node = castToNode( nodeIt->next() );
1595           map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node);
1596           if (it == anEnforcedNodeToGhs3dIdMap.end()) {
1597             it = anExistingEnforcedNodeToGhs3dIdMap.find(node);
1598             if (it == anEnforcedNodeToGhs3dIdMap.end())
1599               throw "Node not found";
1600           }
1601           ntri[index] = it->second;
1602           index++;
1603         }
1604         GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint);
1605         aFaceGroupByGhs3dId[k] = theEnforcedTriangles.find(elem)->second;
1606         usedEnforcedTriangles++;
1607       }
1608     }
1609   }
1610
1611
1612   if (usedEnforcedTriangles) {
1613     GmfSetKwd(idx, GmfRequiredTriangles, usedEnforcedTriangles);
1614     for (int enfID=1;enfID<=usedEnforcedTriangles;enfID++)
1615       GmfSetLin(idx, GmfRequiredTriangles, anElemSet.size()+enfID);
1616   }
1617
1618   GmfCloseMesh(idx);
1619   if (idxRequired)
1620     GmfCloseMesh(idxRequired);
1621   if (idxSol)
1622     GmfCloseMesh(idxSol);
1623
1624   return true;
1625
1626 }
1627
1628 //=============================================================================
1629 /*!
1630  *Here we are going to use the MG-Tetra mesher with geometry
1631  */
1632 //=============================================================================
1633
1634 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
1635                                 const TopoDS_Shape& theShape)
1636 {
1637   bool Ok(false);
1638   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
1639
1640   // we count the number of shapes
1641   // _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with MG-Tetra as a submesh
1642   // _nbShape = 0;
1643   TopExp_Explorer expBox ( theShape, TopAbs_SOLID );
1644   // for ( ; expBox.More(); expBox.Next() )
1645   //   _nbShape++;
1646
1647   // create bounding box for every shape inside the compound
1648
1649   // int iShape = 0;
1650   // TopoDS_Shape* tabShape;
1651   // double**      tabBox;
1652   // tabShape = new TopoDS_Shape[_nbShape];
1653   // tabBox   = new double*[_nbShape];
1654   // for (int i=0; i<_nbShape; i++)
1655   //   tabBox[i] = new double[6];
1656   // Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
1657
1658   // for (expBox.ReInit(); expBox.More(); expBox.Next()) {
1659   //   tabShape[iShape] = expBox.Current();
1660   //   Bnd_Box BoundingBox;
1661   //   BRepBndLib::Add(expBox.Current(), BoundingBox);
1662   //   BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
1663   //   tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
1664   //   tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
1665   //   tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
1666   //   iShape++;
1667   // }
1668
1669   // a unique working file name
1670   // to avoid access to the same files by eg different users
1671   _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
1672   TCollection_AsciiString aGenericName((char*) _genericName.c_str() );
1673   TCollection_AsciiString aGenericNameRequired = aGenericName + "_required";
1674
1675   TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
1676   TCollection_AsciiString aResultFileName;
1677
1678   TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName, aResSolFileName;
1679 //#ifdef _DEBUG_
1680   aGMFFileName              = aGenericName + ".mesh"; // GMF mesh file
1681   aResultFileName           = aGenericName + "Vol.mesh"; // GMF mesh file
1682   aResSolFileName           = aGenericName + "Vol.sol"; // GMF mesh file
1683   aRequiredVerticesFileName = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file
1684   aSolFileName              = aGenericNameRequired + ".sol"; // GMF solution file
1685 //#else
1686 //  aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
1687 //  aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
1688 //  aRequiredVerticesFileName    = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file
1689 //  aSolFileName    = aGenericNameRequired + ".solb"; // GMF solution file
1690 //#endif
1691
1692   std::map <int,int> aNodeId2NodeIndexMap, aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap;
1693   //std::map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
1694   std::map <int, int> nodeID2nodeIndexMap;
1695   std::map<std::vector<double>, std::string> enfVerticesWithGroup;
1696   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap = GHS3DPlugin_Hypothesis::GetEnforcedVerticesCoordsSize(_hyp);
1697   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
1698   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
1699   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
1700 //   TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
1701   GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
1702
1703   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList enfVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
1704   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList::const_iterator enfVerIt = enfVertices.begin();
1705   std::vector<double> coords;
1706
1707   for ( ; enfVerIt != enfVertices.end() ; ++enfVerIt)
1708   {
1709     GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertex* enfVertex = (*enfVerIt);
1710 //     if (enfVertex->geomEntry.empty() && enfVertex->coords.size()) {
1711     if (enfVertex->coords.size()) {
1712       coordsSizeMap.insert(make_pair(enfVertex->coords,enfVertex->size));
1713       enfVerticesWithGroup.insert(make_pair(enfVertex->coords,enfVertex->groupName));
1714 //       MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<enfVertex->coords[0]<<","<<enfVertex->coords[1]<<","<<enfVertex->coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
1715     }
1716     else {
1717 //     if (!enfVertex->geomEntry.empty()) {
1718       TopoDS_Shape GeomShape = entryToShape(enfVertex->geomEntry);
1719 //       GeomType = GeomShape.ShapeType();
1720
1721 //       if (!enfVertex->isCompound) {
1722 // //       if (GeomType == TopAbs_VERTEX) {
1723 //         coords.clear();
1724 //         aPnt = BRep_Tool::Pnt(TopoDS::Vertex(GeomShape));
1725 //         coords.push_back(aPnt.X());
1726 //         coords.push_back(aPnt.Y());
1727 //         coords.push_back(aPnt.Z());
1728 //         if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
1729 //           coordsSizeMap.insert(make_pair(coords,enfVertex->size));
1730 //           enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
1731 //         }
1732 //       }
1733 //
1734 //       // Group Management
1735 //       else {
1736 //       if (GeomType == TopAbs_COMPOUND){
1737         for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
1738           coords.clear();
1739           if (it.Value().ShapeType() == TopAbs_VERTEX){
1740             gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value()));
1741             coords.push_back(aPnt.X());
1742             coords.push_back(aPnt.Y());
1743             coords.push_back(aPnt.Z());
1744             if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
1745               coordsSizeMap.insert(make_pair(coords,enfVertex->size));
1746               enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
1747 //               MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<coords[0]<<","<<coords[1]<<","<<coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
1748             }
1749           }
1750         }
1751 //       }
1752     }
1753   }
1754   int nbEnforcedVertices = coordsSizeMap.size();
1755   int nbEnforcedNodes = enforcedNodes.size();
1756
1757   std::string tmpStr;
1758   (nbEnforcedNodes <= 1) ? tmpStr = "node" : "nodes";
1759   std::cout << nbEnforcedNodes << " enforced " << tmpStr << " from hypo" << std::endl;
1760   (nbEnforcedVertices <= 1) ? tmpStr = "vertex" : "vertices";
1761   std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl;
1762
1763   SMESH_MesherHelper helper( theMesh );
1764   helper.SetSubShape( theShape );
1765
1766   std::vector <const SMDS_MeshNode*> aNodeByGhs3dId, anEnforcedNodeByGhs3dId;
1767   std::vector <const SMDS_MeshElement*> aFaceByGhs3dId;
1768   std::map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
1769   std::vector<std::string> aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId;
1770
1771   // proxyMesh must live till readGMFFile() as a proxy face can be used by
1772   // MG-Tetra for domain indication
1773   //{
1774     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
1775
1776     // make prisms on quadrangles
1777     if ( theMesh.NbQuadrangles() > 0 )
1778     {
1779       vector<SMESH_ProxyMesh::Ptr> components;
1780       for (expBox.ReInit(); expBox.More(); expBox.Next())
1781       {
1782         if ( _viscousLayersHyp )
1783         {
1784           proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() );
1785           if ( !proxyMesh )
1786             return false;
1787         }
1788         StdMeshers_QuadToTriaAdaptor* q2t = new StdMeshers_QuadToTriaAdaptor;
1789         q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() );
1790         components.push_back( SMESH_ProxyMesh::Ptr( q2t ));
1791       }
1792       proxyMesh.reset( new SMESH_ProxyMesh( components ));
1793     }
1794     // build viscous layers
1795     else if ( _viscousLayersHyp )
1796     {
1797       proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape );
1798       if ( !proxyMesh )
1799         return false;
1800     }
1801
1802     // Ok = (writePoints( aPointsFile, helper,
1803     //                    aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, aGhs3dIdToNodeMap,
1804     //                    nodeIDToSizeMap,
1805     //                    coordsSizeMap, enforcedNodes, enforcedEdges, enforcedTriangles)
1806     //       &&
1807     //       writeFaces ( aFacesFile, *proxyMesh, theShape,
1808     //                    aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap,
1809     //                    enforcedEdges, enforcedTriangles ));
1810     int anInvalidEnforcedFlags = 0;
1811     Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
1812                       *proxyMesh, helper,
1813                       aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap,
1814                       aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
1815                       enforcedNodes, enforcedEdges, enforcedTriangles, /*enforcedQuadrangles,*/
1816                       enfVerticesWithGroup, coordsSizeMap, anInvalidEnforcedFlags);
1817     //}
1818
1819   // Write aSmdsToGhs3dIdMap to temp file
1820   TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
1821   aSmdsToGhs3dIdMapFileName = aGenericName + ".ids";  // ids relation
1822   ofstream aIdsFile  ( aSmdsToGhs3dIdMapFileName.ToCString()  , ios::out);
1823   Ok = aIdsFile.rdbuf()->is_open();
1824   if (!Ok) {
1825     INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName);
1826     return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName);
1827   }
1828   INFOS( "Writing ids relation into " << aSmdsToGhs3dIdMapFileName);
1829   aIdsFile << "Smds MG-Tetra" << std::endl;
1830   map <int,int>::const_iterator myit;
1831   for (myit=aSmdsToGhs3dIdMap.begin() ; myit != aSmdsToGhs3dIdMap.end() ; ++myit) {
1832     aIdsFile << myit->first << " " << myit->second << std::endl;
1833   }
1834
1835   aIdsFile.close();
1836
1837   if ( ! Ok ) {
1838     if ( !_keepFiles ) {
1839       removeFile( aGMFFileName );
1840       removeFile( aRequiredVerticesFileName );
1841       removeFile( aSolFileName );
1842       removeFile( aSmdsToGhs3dIdMapFileName );
1843     }
1844     return error(COMPERR_BAD_INPUT_MESH);
1845   }
1846   removeFile( aResultFileName ); // needed for boundary recovery module usage
1847
1848   // -----------------
1849   // run MG-Tetra mesher
1850   // -----------------
1851
1852   TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
1853
1854   cmd += TCollection_AsciiString(" --in ") + aGMFFileName;
1855   if ( nbEnforcedVertices + nbEnforcedNodes)
1856     cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired;
1857   cmd += TCollection_AsciiString(" --out ") + aResultFileName;
1858   if ( !_logInStandardOutput )
1859     cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
1860
1861   std::cout << std::endl;
1862   std::cout << "MG-Tetra execution..." << std::endl;
1863   std::cout << cmd << std::endl;
1864
1865   _compute_canceled = false;
1866
1867   int err = system( cmd.ToCString() ); // run
1868
1869   std::string errStr;
1870   if ( err )
1871     errStr = SMESH_Comment("system(mg-tetra.exe ...) command failed with error: ")
1872       << strerror( errno );
1873   else
1874     std::cout << std::endl << "End of MG-Tetra execution !" << std::endl;
1875
1876   // --------------
1877   // read a result
1878   // --------------
1879
1880   // Mapping the result file
1881
1882   // int fileOpen;
1883   // fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1884   // if ( fileOpen < 0 ) {
1885   //   std::cout << std::endl;
1886   //   std::cout << "Can't open the " << aResultFileName.ToCString() << " MG-Tetra output file" << std::endl;
1887   //   std::cout << "Log: " << aLogFileName << std::endl;
1888   //   Ok = false;
1889   // }
1890   // else {
1891     GHS3DPlugin_Hypothesis::TSetStrings groupsToRemove = GHS3DPlugin_Hypothesis::GetGroupsToRemove(_hyp);
1892     bool toMeshHoles =
1893       _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
1894     const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp );
1895
1896     helper.IsQuadraticSubMesh( theShape );
1897     helper.SetElementsOnShape( false );
1898
1899 //     Ok = readResultFile( fileOpen,
1900 // #ifdef WIN32
1901 //                          aResultFileName.ToCString(),
1902 // #endif
1903 //                          this,
1904 //                          /*theMesh, */helper, tabShape, tabBox, _nbShape,
1905 //                          aGhs3dIdToNodeMap, aNodeId2NodeIndexMap,
1906 //                          toMeshHoles,
1907 //                          nbEnforcedVertices, nbEnforcedNodes,
1908 //                          enforcedEdges, enforcedTriangles,
1909 //                          toMakeGroupsOfDomains );
1910
1911     Ok = readGMFFile(aResultFileName.ToCString(),
1912                      this,
1913                      &helper, aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap,
1914                      aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
1915                      groupsToRemove, toMakeGroupsOfDomains, toMeshHoles);
1916
1917     removeEmptyGroupsOfDomains( helper.GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains );
1918     //}
1919
1920
1921
1922
1923   // ---------------------
1924   // remove working files
1925   // ---------------------
1926
1927   if ( Ok )
1928   {
1929     if ( anInvalidEnforcedFlags )
1930       error( COMPERR_WARNING, flagsToErrorStr( anInvalidEnforcedFlags ));
1931     if ( _removeLogOnSuccess )
1932       removeFile( aLogFileName );
1933     // if ( _hyp && _hyp->GetToMakeGroupsOfDomains() )
1934     //   error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since the mesh is on shape" );
1935   }
1936   else if ( OSD_File( aLogFileName ).Size() > 0 )
1937   {
1938     // get problem description from the log file
1939     _Ghs2smdsConvertor conv( aNodeByGhs3dId, proxyMesh );
1940     storeErrorDescription( aLogFileName, conv );
1941   }
1942   else
1943   {
1944     // the log file is empty
1945     removeFile( aLogFileName );
1946     INFOS( "MG-Tetra Error, " << errStr);
1947     error(COMPERR_ALGO_FAILED, errStr);
1948   }
1949
1950   if ( !_keepFiles ) {
1951     if (! Ok && _compute_canceled)
1952       removeFile( aLogFileName );
1953     removeFile( aGMFFileName );
1954     removeFile( aRequiredVerticesFileName );
1955     removeFile( aSolFileName );
1956     removeFile( aResSolFileName );
1957     removeFile( aResultFileName );
1958     removeFile( aSmdsToGhs3dIdMapFileName );
1959   }
1960   std::cout << "<" << aResultFileName.ToCString() << "> MG-Tetra output file ";
1961   if ( !Ok )
1962     std::cout << "not ";
1963   std::cout << "treated !" << std::endl;
1964   std::cout << std::endl;
1965
1966   // _nbShape = 0;    // re-initializing _nbShape for the next Compute() method call
1967   // delete [] tabShape;
1968   // delete [] tabBox;
1969
1970   return Ok;
1971 }
1972
1973 //=============================================================================
1974 /*!
1975  *Here we are going to use the MG-Tetra mesher w/o geometry
1976  */
1977 //=============================================================================
1978 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
1979                                 SMESH_MesherHelper* theHelper)
1980 {
1981   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
1982
1983   theHelper->IsQuadraticSubMesh( theHelper->GetSubShape() );
1984
1985   // a unique working file name
1986   // to avoid access to the same files by eg different users
1987   _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
1988   TCollection_AsciiString aGenericName((char*) _genericName.c_str() );
1989   TCollection_AsciiString aGenericNameRequired = aGenericName + "_required";
1990
1991   TCollection_AsciiString aLogFileName    = aGenericName + ".log";    // log
1992   TCollection_AsciiString aResultFileName;
1993   bool Ok;
1994
1995   TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName, aResSolFileName;
1996 //#ifdef _DEBUG_
1997   aGMFFileName              = aGenericName + ".mesh"; // GMF mesh file
1998   aResultFileName           = aGenericName + "Vol.mesh"; // GMF mesh file
1999   aResSolFileName           = aGenericName + "Vol.sol"; // GMF mesh file
2000   aRequiredVerticesFileName = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file
2001   aSolFileName              = aGenericNameRequired + ".sol"; // GMF solution file
2002 //#else
2003 //  aGMFFileName    = aGenericName + ".meshb"; // GMF mesh file
2004 //  aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file
2005 //  aRequiredVerticesFileName    = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file
2006 //  aSolFileName    = aGenericNameRequired + ".solb"; // GMF solution file
2007 //#endif
2008
2009   std::map <int, int> nodeID2nodeIndexMap;
2010   std::map<std::vector<double>, std::string> enfVerticesWithGroup;
2011   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap;
2012   TopoDS_Shape GeomShape;
2013 //   TopAbs_ShapeEnum GeomType;
2014   std::vector<double> coords;
2015   gp_Pnt aPnt;
2016   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertex* enfVertex;
2017
2018   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList enfVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
2019   GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList::const_iterator enfVerIt = enfVertices.begin();
2020
2021   for ( ; enfVerIt != enfVertices.end() ; ++enfVerIt)
2022   {
2023     enfVertex = (*enfVerIt);
2024 //     if (enfVertex->geomEntry.empty() && enfVertex->coords.size()) {
2025     if (enfVertex->coords.size()) {
2026       coordsSizeMap.insert(make_pair(enfVertex->coords,enfVertex->size));
2027       enfVerticesWithGroup.insert(make_pair(enfVertex->coords,enfVertex->groupName));
2028 //       MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<enfVertex->coords[0]<<","<<enfVertex->coords[1]<<","<<enfVertex->coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
2029     }
2030     else {
2031 //     if (!enfVertex->geomEntry.empty()) {
2032       GeomShape = entryToShape(enfVertex->geomEntry);
2033 //       GeomType = GeomShape.ShapeType();
2034
2035 //       if (!enfVertex->isCompound) {
2036 // //       if (GeomType == TopAbs_VERTEX) {
2037 //         coords.clear();
2038 //         aPnt = BRep_Tool::Pnt(TopoDS::Vertex(GeomShape));
2039 //         coords.push_back(aPnt.X());
2040 //         coords.push_back(aPnt.Y());
2041 //         coords.push_back(aPnt.Z());
2042 //         if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
2043 //           coordsSizeMap.insert(make_pair(coords,enfVertex->size));
2044 //           enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
2045 //         }
2046 //       }
2047 //
2048 //       // Group Management
2049 //       else {
2050 //       if (GeomType == TopAbs_COMPOUND){
2051         for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){
2052           coords.clear();
2053           if (it.Value().ShapeType() == TopAbs_VERTEX){
2054             aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value()));
2055             coords.push_back(aPnt.X());
2056             coords.push_back(aPnt.Y());
2057             coords.push_back(aPnt.Z());
2058             if (coordsSizeMap.find(coords) == coordsSizeMap.end()) {
2059               coordsSizeMap.insert(make_pair(coords,enfVertex->size));
2060               enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName));
2061 //               MESSAGE("enfVerticesWithGroup.insert(make_pair(("<<coords[0]<<","<<coords[1]<<","<<coords[2]<<"),\""<<enfVertex->groupName<<"\"))");
2062             }
2063           }
2064         }
2065 //       }
2066     }
2067   }
2068
2069 //   const SMDS_MeshNode* enfNode;
2070   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp);
2071 //   GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt = enforcedNodes.begin();
2072 //   for ( ; enfNodeIt != enforcedNodes.end() ; ++enfNodeIt)
2073 //   {
2074 //     enfNode = enfNodeIt->first;
2075 //     coords.clear();
2076 //     coords.push_back(enfNode->X());
2077 //     coords.push_back(enfNode->Y());
2078 //     coords.push_back(enfNode->Z());
2079 //     if (enfVerticesWithGro
2080 //       enfVerticesWithGroup.insert(make_pair(coords,enfNodeIt->second));
2081 //   }
2082
2083
2084   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp);
2085   GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp);
2086 //   TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp);
2087   GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp);
2088
2089   std::string tmpStr;
2090
2091   int nbEnforcedVertices = coordsSizeMap.size();
2092   int nbEnforcedNodes = enforcedNodes.size();
2093   (nbEnforcedNodes <= 1) ? tmpStr = "node" : tmpStr = "nodes";
2094   std::cout << nbEnforcedNodes << " enforced " << tmpStr << " from hypo" << std::endl;
2095   (nbEnforcedVertices <= 1) ? tmpStr = "vertex" : tmpStr = "vertices";
2096   std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl;
2097
2098   std::vector <const SMDS_MeshNode*> aNodeByGhs3dId, anEnforcedNodeByGhs3dId;
2099   std::vector <const SMDS_MeshElement*> aFaceByGhs3dId;
2100   std::map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
2101   std::vector<std::string> aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId;
2102
2103   // proxyMesh must live till readGMFFile() as a proxy face can be used by
2104   // MG-Tetra for domain indication
2105   //{
2106     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh ));
2107     if ( theMesh.NbQuadrangles() > 0 )
2108     {
2109       StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor;
2110       aQuad2Trias->Compute( theMesh );
2111       proxyMesh.reset( aQuad2Trias );
2112     }
2113
2114     int anInvalidEnforcedFlags = 0;
2115     Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(),
2116                       *proxyMesh, *theHelper,
2117                       aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap,
2118                       aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
2119                       enforcedNodes, enforcedEdges, enforcedTriangles,
2120                       enfVerticesWithGroup, coordsSizeMap, anInvalidEnforcedFlags);
2121     //}
2122
2123   // -----------------
2124   // run MG-Tetra mesher
2125   // -----------------
2126
2127   TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str());
2128
2129   cmd += TCollection_AsciiString(" --in ") + aGMFFileName;
2130   if ( nbEnforcedVertices + nbEnforcedNodes)
2131     cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired;
2132   cmd += TCollection_AsciiString(" --out ") + aResultFileName;
2133   if ( !_logInStandardOutput )
2134     cmd += TCollection_AsciiString(" 1>" ) + aLogFileName;  // dump into file
2135
2136   std::cout << std::endl;
2137   std::cout << "MG-Tetra execution..." << std::endl;
2138   std::cout << cmd << std::endl;
2139
2140   _compute_canceled = false;
2141
2142   int err = system( cmd.ToCString() ); // run
2143
2144   std::string errStr;
2145   if ( err )
2146     errStr = SMESH_Comment("system(mg-tetra.exe ...) command failed with error: ")
2147       << strerror( errno );
2148   else
2149     std::cout << std::endl << "End of MG-Tetra execution !" << std::endl;
2150
2151   // --------------
2152   // read a result
2153   // --------------
2154   GHS3DPlugin_Hypothesis::TSetStrings groupsToRemove = GHS3DPlugin_Hypothesis::GetGroupsToRemove(_hyp);
2155   const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp );
2156
2157   Ok = readGMFFile(aResultFileName.ToCString(),
2158                    this,
2159                    theHelper, aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap,
2160                    aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId,
2161                    groupsToRemove, toMakeGroupsOfDomains);
2162
2163   updateMeshGroups(theHelper->GetMesh(), groupsToRemove);
2164   removeEmptyGroupsOfDomains( theHelper->GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains );
2165
2166   if ( Ok ) {
2167     GHS3DPlugin_Hypothesis* that = (GHS3DPlugin_Hypothesis*)this->_hyp;
2168     if (that)
2169       that->ClearGroupsToRemove();
2170   }
2171   // ---------------------
2172   // remove working files
2173   // ---------------------
2174
2175   if ( Ok )
2176   {
2177     if ( anInvalidEnforcedFlags )
2178       error( COMPERR_WARNING, flagsToErrorStr( anInvalidEnforcedFlags ));
2179     if ( _removeLogOnSuccess )
2180       removeFile( aLogFileName );
2181
2182     //if ( !toMakeGroupsOfDomains && _hyp && _hyp->GetToMakeGroupsOfDomains() )
2183     //error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since 'toMeshHoles' is OFF." );
2184   }
2185   else if ( OSD_File( aLogFileName ).Size() > 0 )
2186   {
2187     // get problem description from the log file
2188     _Ghs2smdsConvertor conv( aNodeByGhs3dId, proxyMesh );
2189     storeErrorDescription( aLogFileName, conv );
2190   }
2191   else {
2192     // the log file is empty
2193     removeFile( aLogFileName );
2194     INFOS( "MG-Tetra Error, " << errStr);
2195     error(COMPERR_ALGO_FAILED, errStr);
2196   }
2197
2198   if ( !_keepFiles )
2199   {
2200     if (! Ok && _compute_canceled)
2201       removeFile( aLogFileName );
2202     removeFile( aGMFFileName );
2203     removeFile( aResultFileName );
2204     removeFile( aRequiredVerticesFileName );
2205     removeFile( aSolFileName );
2206     removeFile( aResSolFileName );
2207   }
2208   return Ok;
2209 }
2210
2211 void GHS3DPlugin_GHS3D::CancelCompute()
2212 {
2213   _compute_canceled = true;
2214 #ifdef WIN32
2215 #else
2216   std::string cmd = "ps xo pid,args | grep " + _genericName;
2217   //cmd += " | grep -e \"^ *[0-9]\\+ \\+" + GHS3DPlugin_Hypothesis::GetExeName() + "\"";
2218   cmd += " | awk '{print $1}' | xargs kill -9 > /dev/null 2>&1";
2219   system( cmd.c_str() );
2220 #endif
2221 }
2222
2223 //================================================================================
2224 /*!
2225  * \brief Provide human readable text by error code reported by MG-Tetra
2226  */
2227 //================================================================================
2228
2229 static const char* translateError(const int errNum)
2230 {
2231   switch ( errNum ) {
2232   case 0:
2233     return "The surface mesh includes a face of type other than edge, "
2234       "triangle or quadrilateral. This face type is not supported.";
2235   case 1:
2236     return "Not enough memory for the face table.";
2237   case 2:
2238     return "Not enough memory.";
2239   case 3:
2240     return "Not enough memory.";
2241   case 4:
2242     return "Face is ignored.";
2243   case 5:
2244     return "End of file. Some data are missing in the file.";
2245   case 6:
2246     return "Read error on the file. There are wrong data in the file.";
2247   case 7:
2248     return "the metric file is inadequate (dimension other than 3).";
2249   case 8:
2250     return "the metric file is inadequate (values not per vertices).";
2251   case 9:
2252     return "the metric file contains more than one field.";
2253   case 10:
2254     return "the number of values in the \".bb\" (metric file) is incompatible with the expected"
2255       "value of number of mesh vertices in the \".noboite\" file.";
2256   case 12:
2257     return "Too many sub-domains.";
2258   case 13:
2259     return "the number of vertices is negative or null.";
2260   case 14:
2261     return "the number of faces is negative or null.";
2262   case 15:
2263     return "A face has a null vertex.";
2264   case 22:
2265     return "incompatible data.";
2266   case 131:
2267     return "the number of vertices is negative or null.";
2268   case 132:
2269     return "the number of vertices is negative or null (in the \".mesh\" file).";
2270   case 133:
2271     return "the number of faces is negative or null.";
2272   case 1000:
2273     return "A face appears more than once in the input surface mesh.";
2274   case 1001:
2275     return "An edge appears more than once in the input surface mesh.";
2276   case 1002:
2277     return "A face has a vertex negative or null.";
2278   case 1003:
2279     return "NOT ENOUGH MEMORY.";
2280   case 2000:
2281     return "Not enough available memory.";
2282   case 2002:
2283     return "Some initial points cannot be inserted. The surface mesh is probably very bad "
2284       "in terms of quality or the input list of points is wrong.";
2285   case 2003:
2286     return "Some vertices are too close to one another or coincident.";
2287   case 2004:
2288     return "Some vertices are too close to one another or coincident.";
2289   case 2012:
2290     return "A vertex cannot be inserted.";
2291   case 2014:
2292     return "There are at least two points considered as coincident.";
2293   case 2103:
2294     return "Some vertices are too close to one another or coincident.";
2295   case 3000:
2296     return "The surface mesh regeneration step has failed.";
2297   case 3009:
2298     return "Constrained edge cannot be enforced.";
2299   case 3019:
2300     return "Constrained face cannot be enforced.";
2301   case 3029:
2302     return "Missing faces.";
2303   case 3100:
2304     return "No guess to start the definition of the connected component(s).";
2305   case 3101:
2306     return "The surface mesh includes at least one hole. The domain is not well defined.";
2307   case 3102:
2308     return "Impossible to define a component.";
2309   case 3103:
2310     return "The surface edge intersects another surface edge.";
2311   case 3104:
2312     return "The surface edge intersects the surface face.";
2313   case 3105:
2314     return "One boundary point lies within a surface face.";
2315   case 3106:
2316     return "One surface edge intersects a surface face.";
2317   case 3107:
2318     return "One boundary point lies within a surface edge.";
2319   case 3108:
2320     return "Insufficient memory ressources detected due to a bad quality surface mesh leading "
2321       "to too many swaps.";
2322   case 3109:
2323     return "Edge is unique (i.e., bounds a hole in the surface).";
2324   case 3122:
2325     return "Presumably, the surface mesh is not compatible with the domain being processed.";
2326   case 3123:
2327     return "Too many components, too many sub-domain.";
2328   case 3209:
2329     return "The surface mesh includes at least one hole. "
2330       "Therefore there is no domain properly defined.";
2331   case 3300:
2332     return "Statistics.";
2333   case 3400:
2334     return "Statistics.";
2335   case 3500:
2336     return "Warning, it is dramatically tedious to enforce the boundary items.";
2337   case 4000:
2338     return "Not enough memory at this time, nevertheless, the program continues. "
2339       "The expected mesh will be correct but not really as large as required.";
2340   case 4002:
2341     return "see above error code, resulting quality may be poor.";
2342   case 4003:
2343     return "Not enough memory at this time, nevertheless, the program continues (warning).";
2344   case 8000:
2345     return "Unknown face type.";
2346   case 8005:
2347   case 8006:
2348     return "End of file. Some data are missing in the file.";
2349   case 9000:
2350     return "A too small volume element is detected.";
2351   case 9001:
2352     return "There exists at least a null or negative volume element.";
2353   case 9002:
2354     return "There exist null or negative volume elements.";
2355   case 9003:
2356     return "A too small volume element is detected. A face is considered being degenerated.";
2357   case 9100:
2358     return "Some element is suspected to be very bad shaped or wrong.";
2359   case 9102:
2360     return "A too bad quality face is detected. This face is considered degenerated.";
2361   case 9112:
2362     return "A too bad quality face is detected. This face is degenerated.";
2363   case 9122:
2364     return "Presumably, the surface mesh is not compatible with the domain being processed.";
2365   case 9999:
2366     return "Abnormal error occured, contact hotline.";
2367   case 23600:
2368     return "Not enough memory for the face table.";
2369   case 23601:
2370     return "The algorithm cannot run further. "
2371       "The surface mesh is probably very bad in terms of quality.";
2372   case 23602:
2373     return "Bad vertex number.";
2374   case 1001200:
2375     return "Cannot close mesh file NomFil.";
2376   case 1002010:
2377     return "There are wrong data.";
2378   case 1002120:
2379     return "The number of faces is negative or null.";
2380   case 1002170:
2381     return "The number of vertices is negative or null in the '.sol' file.";
2382   case 1002190:
2383     return "The number of tetrahedra is negative or null.";
2384   case 1002210:
2385     return "The number of vertices is negative or null.";
2386   case 1002211:
2387     return "A face has a vertex negative or null.";
2388   case 1002270:
2389     return "The field is not a size in file NomFil.";
2390   case 1002280:
2391     return "A count is wrong in the enclosing box in the .boite.mesh input "
2392       "file (option '--read_boite').";
2393   case 1002290:
2394     return "A tetrahedron has a vertex with a negative number.";
2395   case 1002300:
2396     return "the 'MeshVersionFormatted' is not 1 or 2 in the '.mesh' file or the '.sol'.";
2397  case 1002370:
2398    return "The number of values in the '.sol' (metric file) is incompatible with "
2399      "the expected value of number of mesh vertices in the '.mesh' file.";
2400   case 1003000:
2401     return "Not enough memory.";
2402   case 1003020:
2403     return "Not enough memory for the face table.";
2404   case 1003050:
2405     return "Insufficient memory ressources detected due to a bad quality "
2406       "surface mesh leading to too many swaps.";
2407   case 1005010:
2408     return "The surface coordinates of a vertex are differing from the "
2409       "volume coordinates, probably due to a precision problem.";
2410   case 1005050:
2411     return "Invalid dimension. Dimension 3 expected.";
2412   case 1005100:
2413     return "A point has a tag 0. This point is probably outside the domain which has been meshed.";
2414   case 1005103:
2415     return "The vertices of an element are too close to one another or coincident.";
2416   case 1005104:
2417     return "There are at least two points whose distance is very small, and considered as coincident.";
2418   case 1005105:
2419     return "Two vertices are too close to one another or coincident.";
2420   case 1005106:
2421     return "A vertex cannot be inserted.";
2422   case 1005107:
2423     return "Two vertices are too close to one another or coincident. Note : When "
2424       "this error occurs during the overconstrained processing phase, this is only "
2425       "a warning which means that it is difficult to break some overconstrained facets.";
2426   case 1005110:
2427     return "Two surface edges are intersecting.";
2428   case 1005120:
2429     return "A surface edge intersects a surface face.";
2430   case 1005150:
2431     return "A boundary point lies within a surface face.";
2432   case 1005160:
2433     return "A boundary point lies within a surface edge.";
2434   case 1005200:
2435     return "A surface mesh appears more than once in the input surface mesh.";
2436   case 1005210:
2437     return "An edge appears more than once in the input surface mesh.";
2438   case 1005225:
2439     return "Surface with unvalid triangles.";
2440   case 1005270:
2441     return "The metric in the '.sol' file contains more than one field.";
2442   case 1005300:
2443     return "The surface mesh includes at least one hole. The domain is not well defined.";
2444   case 1005301:
2445     return "Presumably, the surface mesh is not compatible with the domain being processed (warning).";
2446   case 1005302:
2447     return "Probable faces overlapping somewher.";
2448   case 1005320:
2449     return "The quadratic version does not work with prescribed free edges.";
2450   case 1005321:
2451     return "The quadratic version does not work with a volume mesh.";
2452   case 1005370:
2453     return "The metric in the '.sol' file is inadequate (values not per vertices).";
2454   case 1005371:
2455     return "The number of vertices in the '.sol' is different from the one in the "
2456       "'.mesh' file for the required vertices (option '--required_vertices').";
2457   case 1005372:
2458     return "More than one type in file NomFil. The type must be equal to 1 in the '.sol'"
2459       "for the required vertices (option '--required_vertices').";
2460   case 1005515:
2461     return "Bad vertex number.";
2462   case 1005560:
2463     return "No guess to start the definition of the connected component(s).";
2464   case 1005602:
2465     return "Some initial points cannot be inserted.";
2466   case 1005620:
2467     return "A too bad quality face is detected. This face is considered degenerated.";
2468   case 1005621:
2469     return "A too bad quality face is detected. This face is degenerated.";
2470   case 1005622:
2471     return "The algorithm cannot run further.";
2472   case 1005690:
2473     return "A too small volume element is detected.";
2474   case 1005691:
2475     return "A tetrahedra is suspected to be very bad shaped or wrong.";
2476   case 1005692:
2477     return "There is at least a null or negative volume element. The resulting mesh"
2478       "may be inappropriate.";
2479   case 1005693:
2480     return "There are some null or negative volume element. The resulting mesh may"
2481       "be inappropriate.";
2482   case 1005820:
2483     return "An edge is unique (i.e., bounds a hole in the surface).";
2484   case 1007000:
2485     return "Abnormal or internal error.";
2486   case 1007010:
2487     return "Too many components with respect to too many sub-domain.";
2488   case 1007400:
2489     return "An internal error has been encountered or a signal has been received. "
2490       "Current mesh will not be saved.";
2491   case 1008491:
2492     return "Impossible to define a component.";
2493   case 1008410:
2494     return "There are some overconstrained edges.";
2495   case 1008420:
2496     return "There are some overconstrained facets.";
2497   case 1008422:
2498     return "Give the number of missing faces (information given when regeneration phase failed).";
2499   case 1008423:
2500     return "A constrained face cannot be enforced (information given when regeneration phase failed).";
2501   case 1008441:
2502     return "A constrained edge cannot be enforced.";
2503   case 1008460:
2504     return "It is dramatically tedious to enforce the boundary items.";
2505   case 1008480:
2506     return "The surface mesh regeneration step has failed. A .boite.mesh and .boite.map files are created.";
2507   case 1008490:
2508     return "Invalid resulting mesh.";
2509   case 1008495:
2510     return "P2 correction not successful.";
2511   case 1009000:
2512     return "Program has received an interruption or a termination signal sent by the "
2513       "user or the system administrator. Current mesh will not be saved.";
2514   }
2515   return "";
2516 }
2517
2518 //================================================================================
2519 /*!
2520  * \brief Retrieve from a string given number of integers
2521  */
2522 //================================================================================
2523
2524 static char* getIds( char* ptr, int nbIds, vector<int>& ids )
2525 {
2526   ids.clear();
2527   ids.reserve( nbIds );
2528   while ( nbIds )
2529   {
2530     while ( !isdigit( *ptr )) ++ptr;
2531     if ( ptr[-1] == '-' ) --ptr;
2532     ids.push_back( strtol( ptr, &ptr, 10 ));
2533     --nbIds;
2534   }
2535   return ptr;
2536 }
2537
2538 //================================================================================
2539 /*!
2540  * \brief Retrieve problem description form a log file
2541  *  \retval bool - always false
2542  */
2543 //================================================================================
2544
2545 bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& logFile,
2546                                               const _Ghs2smdsConvertor &     toSmdsConvertor )
2547 {
2548   if(_compute_canceled)
2549     return error(SMESH_Comment("interruption initiated by user"));
2550   // open file
2551 #ifdef WIN32
2552   int file = ::_open (logFile.ToCString(), _O_RDONLY|_O_BINARY);
2553 #else
2554   int file = ::open (logFile.ToCString(), O_RDONLY);
2555 #endif
2556   if ( file < 0 )
2557     return error( SMESH_Comment("See ") << logFile << " for problem description");
2558
2559   // get file size
2560   off_t length = lseek( file, 0, SEEK_END);
2561   lseek( file, 0, SEEK_SET);
2562
2563   // read file
2564   vector< char > buf( length );
2565   int nBytesRead = ::read (file, & buf[0], length);
2566   ::close (file);
2567   char* ptr = & buf[0];
2568   char* bufEnd = ptr + nBytesRead;
2569
2570   SMESH_Comment errDescription;
2571
2572   enum { NODE = 1, EDGE, TRIA, VOL, SKIP_ID = 1 };
2573
2574   // look for MeshGems version
2575   // Since "MG-TETRA -- MeshGems 1.1-3 (January, 2013)" error codes change.
2576   // To discriminate old codes from new ones we add 1000000 to the new codes.
2577   // This way value of the new codes is same as absolute value of codes printed
2578   // in the log after "MGMESSAGE" string.
2579   int versionAddition = 0;
2580   {
2581     char* verPtr = ptr;
2582     while ( ++verPtr < bufEnd )
2583     {
2584       if ( strncmp( verPtr, "MG-TETRA -- MeshGems ", 21 ) != 0 )
2585         continue;
2586       if ( strcmp( verPtr, "MG-TETRA -- MeshGems 1.1-3 " ) >= 0 )
2587         versionAddition = 1000000;
2588       ptr = verPtr;
2589       break;
2590     }
2591   }
2592
2593   // look for errors "ERR #"
2594
2595   set<string> foundErrorStr; // to avoid reporting same error several times
2596   set<int>    elemErrorNums; // not to report different types of errors with bad elements
2597   while ( ++ptr < bufEnd )
2598   {
2599     if ( strncmp( ptr, "ERR ", 4 ) != 0 )
2600       continue;
2601
2602     list<const SMDS_MeshElement*> badElems;
2603     vector<int> nodeIds;
2604
2605     ptr += 4;
2606     char* errBeg = ptr;
2607     int   errNum = strtol(ptr, &ptr, 10) + versionAddition;
2608     // we treat errors enumerated in [SALOME platform 0019316] issue
2609     // and all errors from a new (Release 1.1) MeshGems User Manual
2610     switch ( errNum ) {
2611     case 0015: // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex.
2612     case 1005620 : // a too bad quality face is detected. This face is considered degenerated.
2613       ptr = getIds(ptr, SKIP_ID, nodeIds);
2614       ptr = getIds(ptr, TRIA, nodeIds);
2615       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2616       break;
2617     case 1005621 : // a too bad quality face is detected. This face is degenerated.
2618       // hence the is degenerated it is invisible, add its edges in addition
2619       ptr = getIds(ptr, SKIP_ID, nodeIds);
2620       ptr = getIds(ptr, TRIA, nodeIds);
2621       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2622       {
2623         vector<int> edgeNodes( nodeIds.begin(), --nodeIds.end() ); // 01
2624         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
2625         edgeNodes[1] = nodeIds[2]; // 02
2626         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
2627         edgeNodes[0] = nodeIds[1]; // 12
2628       }
2629       break;
2630     case 1000: // Face (f 1, f 2, f 3) appears more than once in the input surface mesh.
2631       // ERR  1000 :  1 3 2
2632     case 1002: // Face (f 1, f 2, f 3) has a vertex negative or null
2633     case 3019: // Constrained face (f 1, f 2, f 3) cannot be enforced
2634     case 1002211: // a face has a vertex negative or null.
2635     case 1005200 : // a surface mesh appears more than once in the input surface mesh.
2636     case 1008423 : // a constrained face cannot be enforced (regeneration phase failed).
2637       ptr = getIds(ptr, TRIA, nodeIds);
2638       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2639       break;
2640     case 1001: // Edge (e1, e2) appears more than once in the input surface mesh
2641     case 3009: // Constrained edge (e1, e2) cannot be enforced (warning).
2642       // ERR  3109 :  EDGE  5 6 UNIQUE
2643     case 3109: // Edge (e1, e2) is unique (i.e., bounds a hole in the surface)
2644     case 1005210 : // an edge appears more than once in the input surface mesh.
2645     case 1005820 : // an edge is unique (i.e., bounds a hole in the surface).
2646     case 1008441 : // a constrained edge cannot be enforced.
2647       ptr = getIds(ptr, EDGE, nodeIds);
2648       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2649       break;
2650     case 2004: // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
2651     case 2014: // at least two points whose distance is dist, i.e., considered as coincident
2652     case 2103: // Vertex v1 and vertex v2 are too close to one another or coincident (warning).
2653       // ERR  2103 :  16 WITH  3
2654     case 1005105 : // two vertices are too close to one another or coincident.
2655     case 1005107: // Two vertices are too close to one another or coincident.
2656       ptr = getIds(ptr, NODE, nodeIds);
2657       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2658       ptr = getIds(ptr, NODE, nodeIds);
2659       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2660       break;
2661     case 2012: // Vertex v1 cannot be inserted (warning).
2662     case 1005106 : // a vertex cannot be inserted.
2663       ptr = getIds(ptr, NODE, nodeIds);
2664       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2665       break;
2666     case 3103: // The surface edge (e1, e2) intersects another surface edge (e3, e4)
2667     case 1005110 : // two surface edges are intersecting.
2668       // ERR  3103 :  1 2 WITH  7 3
2669       ptr = getIds(ptr, EDGE, nodeIds);
2670       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2671       ptr = getIds(ptr, EDGE, nodeIds);
2672       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2673       break;
2674     case 3104: // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3)
2675       // ERR  3104 :  9 10 WITH  1 2 3
2676     case 3106: // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3)
2677     case 1005120 : // a surface edge intersects a surface face.
2678       ptr = getIds(ptr, EDGE, nodeIds);
2679       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2680       ptr = getIds(ptr, TRIA, nodeIds);
2681       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2682       break;
2683     case 3105: // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3)
2684       // ERR  3105 :  8 IN  2 3 5
2685     case 1005150 : // a boundary point lies within a surface face.
2686       ptr = getIds(ptr, NODE, nodeIds);
2687       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2688       ptr = getIds(ptr, TRIA, nodeIds);
2689       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2690       break;
2691     case 3107: // One boundary point (say p1) lies within a surface edge (e1, e2) (stop).
2692       // ERR  3107 :  2 IN  4 1
2693     case 1005160 : // a boundary point lies within a surface edge.
2694       ptr = getIds(ptr, NODE, nodeIds);
2695       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2696       ptr = getIds(ptr, EDGE, nodeIds);
2697       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2698       break;
2699     case 9000: // ERR  9000
2700       //  ELEMENT  261 WITH VERTICES :  7 396 -8 242
2701       //  VOLUME   : -1.11325045E+11 W.R.T. EPSILON   0.
2702       // A too small volume element is detected. Are reported the index of the element,
2703       // its four vertex indices, its volume and the tolerance threshold value
2704       ptr = getIds(ptr, SKIP_ID, nodeIds);
2705       ptr = getIds(ptr, VOL, nodeIds);
2706       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2707       // even if all nodes found, volume it most probably invisible,
2708       // add its faces to demonstrate it anyhow
2709       {
2710         vector<int> faceNodes( nodeIds.begin(), --nodeIds.end() ); // 012
2711         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
2712         faceNodes[2] = nodeIds[3]; // 013
2713         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
2714         faceNodes[1] = nodeIds[2]; // 023
2715         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
2716         faceNodes[0] = nodeIds[1]; // 123
2717         badElems.push_back( toSmdsConvertor.getElement(faceNodes));
2718       }
2719       break;
2720     case 9001: // ERR  9001
2721       //  %% NUMBER OF NEGATIVE VOLUME TETS  :  1
2722       //  %% THE LARGEST NEGATIVE TET        :   1.75376581E+11
2723       //  %%  NUMBER OF NULL VOLUME TETS     :  0
2724       // There exists at least a null or negative volume element
2725       break;
2726     case 9002:
2727       // There exist n null or negative volume elements
2728       break;
2729     case 9003:
2730       // A too small volume element is detected
2731       break;
2732     case 9102:
2733       // A too bad quality face is detected. This face is considered degenerated,
2734       // its index, its three vertex indices together with its quality value are reported
2735       break; // same as next
2736     case 9112: // ERR  9112
2737       //  FACE   2 WITH VERTICES :  4 2 5
2738       //  SMALL INRADIUS :   0.
2739       // A too bad quality face is detected. This face is degenerated,
2740       // its index, its three vertex indices together with its inradius are reported
2741       ptr = getIds(ptr, SKIP_ID, nodeIds);
2742       ptr = getIds(ptr, TRIA, nodeIds);
2743       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2744       // add triangle edges as it most probably has zero area and hence invisible
2745       {
2746         vector<int> edgeNodes(2);
2747         edgeNodes[0] = nodeIds[0]; edgeNodes[1] = nodeIds[1]; // 0-1
2748         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
2749         edgeNodes[1] = nodeIds[2]; // 0-2
2750         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
2751         edgeNodes[0] = nodeIds[1]; // 1-2
2752         badElems.push_back( toSmdsConvertor.getElement(edgeNodes));
2753       }
2754       break;
2755     case 1005103 : // the vertices of an element are too close to one another or coincident.
2756       ptr = getIds(ptr, TRIA, nodeIds);
2757       if ( nodeIds.back() == 0 ) // index of the third vertex of the element (0 for an edge)
2758         nodeIds.resize( EDGE );
2759       badElems.push_back( toSmdsConvertor.getElement(nodeIds));
2760       break;
2761     }
2762
2763     bool isNewError = foundErrorStr.insert( string( errBeg, ptr )).second;
2764     if ( !isNewError )
2765       continue; // not to report same error several times
2766
2767 //     const SMDS_MeshElement* nullElem = 0;
2768 //     bool allElemsOk = ( find( badElems.begin(), badElems.end(), nullElem) == badElems.end());
2769
2770 //     if ( allElemsOk && !badElems.empty() && !elemErrorNums.empty() ) {
2771 //       bool oneMoreErrorType = elemErrorNums.insert( errNum ).second;
2772 //       if ( oneMoreErrorType )
2773 //         continue; // not to report different types of errors with bad elements
2774 //     }
2775
2776     // store bad elements
2777     //if ( allElemsOk ) {
2778       list<const SMDS_MeshElement*>::iterator elem = badElems.begin();
2779       for ( ; elem != badElems.end(); ++elem )
2780         addBadInputElement( *elem );
2781       //}
2782
2783     // make error text
2784     string text = translateError( errNum );
2785     if ( errDescription.find( text ) == text.npos ) {
2786       if ( !errDescription.empty() )
2787         errDescription << "\n";
2788       errDescription << text;
2789     }
2790
2791   } // end while
2792
2793   if ( errDescription.empty() ) { // no errors found
2794     char msgLic1[] = "connection to server failed";
2795     char msgLic2[] = " Dlim ";
2796     if ( search( &buf[0], bufEnd, msgLic1, msgLic1 + strlen(msgLic1)) != bufEnd ||
2797          search( &buf[0], bufEnd, msgLic2, msgLic2 + strlen(msgLic2)) != bufEnd )
2798       errDescription << "Licence problems.";
2799     else
2800     {
2801       char msg2[] = "SEGMENTATION FAULT";
2802       if ( search( &buf[0], bufEnd, msg2, msg2 + strlen(msg2)) != bufEnd )
2803         errDescription << "MG-Tetra: SEGMENTATION FAULT. ";
2804     }
2805   }
2806
2807   if ( errDescription.empty() )
2808     errDescription << "See " << logFile << " for problem description";
2809   else
2810     errDescription << "\nSee " << logFile << " for more information";
2811
2812   return error( errDescription );
2813 }
2814
2815 //================================================================================
2816 /*!
2817  * \brief Creates _Ghs2smdsConvertor
2818  */
2819 //================================================================================
2820
2821 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const map <int,const SMDS_MeshNode*> & ghs2NodeMap,
2822                                         SMESH_ProxyMesh::Ptr                   mesh)
2823   :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 ), _mesh( mesh )
2824 {
2825 }
2826
2827 //================================================================================
2828 /*!
2829  * \brief Creates _Ghs2smdsConvertor
2830  */
2831 //================================================================================
2832
2833 _Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector <const SMDS_MeshNode*> &  nodeByGhsId,
2834                                         SMESH_ProxyMesh::Ptr                   mesh)
2835   : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId ), _mesh( mesh )
2836 {
2837 }
2838
2839 //================================================================================
2840 /*!
2841  * \brief Return SMDS element by ids of MG-Tetra nodes
2842  */
2843 //================================================================================
2844
2845 const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector<int>& ghsNodes) const
2846 {
2847   size_t nbNodes = ghsNodes.size();
2848   vector<const SMDS_MeshNode*> nodes( nbNodes, 0 );
2849   for ( size_t i = 0; i < nbNodes; ++i ) {
2850     int ghsNode = ghsNodes[ i ];
2851     if ( _ghs2NodeMap ) {
2852       map <int,const SMDS_MeshNode*>::const_iterator in = _ghs2NodeMap->find( ghsNode);
2853       if ( in == _ghs2NodeMap->end() )
2854         return 0;
2855       nodes[ i ] = in->second;
2856     }
2857     else {
2858       if ( ghsNode < 1 || ghsNode > (int)_nodeByGhsId->size() )
2859         return 0;
2860       nodes[ i ] = (*_nodeByGhsId)[ ghsNode-1 ];
2861     }
2862   }
2863   if ( nbNodes == 1 )
2864     return nodes[0];
2865
2866   if ( nbNodes == 2 ) {
2867     const SMDS_MeshElement* edge= SMDS_Mesh::FindEdge( nodes[0], nodes[1] );
2868     if ( !edge || edge->GetID() < 1 || _mesh->IsTemporary( edge ))
2869       edge = new SMDS_LinearEdge( nodes[0], nodes[1] );
2870     return edge;
2871   }
2872   if ( nbNodes == 3 ) {
2873     const SMDS_MeshElement* face = SMDS_Mesh::FindFace( nodes );
2874     if ( !face || face->GetID() < 1 || _mesh->IsTemporary( face ))
2875       face = new SMDS_FaceOfNodes( nodes[0], nodes[1], nodes[2] );
2876     return face;
2877   }
2878   if ( nbNodes == 4 )
2879     return new SMDS_VolumeOfNodes( nodes[0], nodes[1], nodes[2], nodes[3] );
2880
2881   return 0;
2882 }
2883
2884
2885 //=============================================================================
2886 /*!
2887  *
2888  */
2889 //=============================================================================
2890 bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
2891                                  const TopoDS_Shape& aShape,
2892                                  MapShapeNbElems& aResMap)
2893 {
2894   int nbtri = 0, nbqua = 0;
2895   double fullArea = 0.0;
2896   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
2897     TopoDS_Face F = TopoDS::Face( exp.Current() );
2898     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
2899     MapShapeNbElemsItr anIt = aResMap.find(sm);
2900     if( anIt==aResMap.end() ) {
2901       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2902       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
2903                                             "Submesh can not be evaluated",this));
2904       return false;
2905     }
2906     std::vector<int> aVec = (*anIt).second;
2907     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
2908     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
2909     GProp_GProps G;
2910     BRepGProp::SurfaceProperties(F,G);
2911     double anArea = G.Mass();
2912     fullArea += anArea;
2913   }
2914
2915   // collect info from edges
2916   int nb0d_e = 0, nb1d_e = 0;
2917   bool IsQuadratic = false;
2918   bool IsFirst = true;
2919   TopTools_MapOfShape tmpMap;
2920   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
2921     TopoDS_Edge E = TopoDS::Edge(exp.Current());
2922     if( tmpMap.Contains(E) )
2923       continue;
2924     tmpMap.Add(E);
2925     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
2926     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
2927     std::vector<int> aVec = (*anIt).second;
2928     nb0d_e += aVec[SMDSEntity_Node];
2929     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
2930     if(IsFirst) {
2931       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
2932       IsFirst = false;
2933     }
2934   }
2935   tmpMap.Clear();
2936
2937   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
2938
2939   GProp_GProps G;
2940   BRepGProp::VolumeProperties(aShape,G);
2941   double aVolume = G.Mass();
2942   double tetrVol = 0.1179*ELen*ELen*ELen;
2943   double CoeffQuality = 0.9;
2944   int nbVols = int(aVolume/tetrVol/CoeffQuality);
2945   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
2946   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
2947   std::vector<int> aVec(SMDSEntity_Last);
2948   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
2949   if( IsQuadratic ) {
2950     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
2951     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
2952     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
2953   }
2954   else {
2955     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
2956     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
2957     aVec[SMDSEntity_Pyramid] = nbqua;
2958   }
2959   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
2960   aResMap.insert(std::make_pair(sm,aVec));
2961
2962   return true;
2963 }
2964
2965 bool GHS3DPlugin_GHS3D::importGMFMesh(const char* theGMFFileName, SMESH_Mesh& theMesh)
2966 {
2967   SMESH_MesherHelper* helper  = new SMESH_MesherHelper(theMesh );
2968   std::vector <const SMDS_MeshNode*> dummyNodeVector;
2969   std::vector <const SMDS_MeshElement*> aFaceByGhs3dId;
2970   std::map<const SMDS_MeshNode*,int> dummyNodeMap;
2971   std::map<std::vector<double>, std::string> dummyEnfVertGroup;
2972   std::vector<std::string> dummyElemGroup;
2973   std::set<std::string> dummyGroupsToRemove;
2974
2975   bool ok = readGMFFile(theGMFFileName,
2976                         this,
2977                         helper, dummyNodeVector, aFaceByGhs3dId, dummyNodeMap, dummyElemGroup, dummyElemGroup, dummyElemGroup, dummyGroupsToRemove);
2978   theMesh.GetMeshDS()->Modified();
2979   return ok;
2980 }
2981
2982 namespace
2983 {
2984   //================================================================================
2985   /*!
2986    * \brief Sub-mesh event listener setting enforced elements as soon as an enforced
2987    *        mesh is loaded
2988    */
2989   struct _EnforcedMeshRestorer : public SMESH_subMeshEventListener
2990   {
2991     _EnforcedMeshRestorer():
2992       SMESH_subMeshEventListener( /*isDeletable = */true, Name() )
2993     {}
2994
2995     //================================================================================
2996     /*!
2997      * \brief Returns an ID of listener
2998      */
2999     static const char* Name() { return "GHS3DPlugin_GHS3D::_EnforcedMeshRestorer"; }
3000
3001     //================================================================================
3002     /*!
3003      * \brief Treat events of the subMesh
3004      */
3005     void ProcessEvent(const int                       event,
3006                       const int                       eventType,
3007                       SMESH_subMesh*                  subMesh,
3008                       SMESH_subMeshEventListenerData* data,
3009                       const SMESH_Hypothesis*         hyp)
3010     {
3011       if ( SMESH_subMesh::SUBMESH_LOADED == event &&
3012            SMESH_subMesh::COMPUTE_EVENT  == eventType &&
3013            data &&
3014            !data->mySubMeshes.empty() )
3015       {
3016         // An enforced mesh (subMesh->_father) has been loaded from hdf file
3017         if ( GHS3DPlugin_Hypothesis* hyp = GetGHSHypothesis( data->mySubMeshes.front() ))
3018           hyp->RestoreEnfElemsByMeshes();
3019       }
3020     }
3021     //================================================================================
3022     /*!
3023      * \brief Returns GHS3DPlugin_Hypothesis used to compute a subMesh
3024      */
3025     static GHS3DPlugin_Hypothesis* GetGHSHypothesis( SMESH_subMesh* subMesh )
3026     {
3027       SMESH_HypoFilter ghsHypFilter
3028         ( SMESH_HypoFilter::HasName( GHS3DPlugin_Hypothesis::GetHypType() ));
3029       return (GHS3DPlugin_Hypothesis* )
3030         subMesh->GetFather()->GetHypothesis( subMesh->GetSubShape(),
3031                                              ghsHypFilter,
3032                                              /*visitAncestors=*/true);
3033     }
3034   };
3035
3036   //================================================================================
3037   /*!
3038    * \brief Sub-mesh event listener removing empty groups created due to "To make
3039    *        groups of domains".
3040    */
3041   struct _GroupsOfDomainsRemover : public SMESH_subMeshEventListener
3042   {
3043     _GroupsOfDomainsRemover():
3044       SMESH_subMeshEventListener( /*isDeletable = */true,
3045                                   "GHS3DPlugin_GHS3D::_GroupsOfDomainsRemover" ) {}
3046     /*!
3047      * \brief Treat events of the subMesh
3048      */
3049     void ProcessEvent(const int                       event,
3050                       const int                       eventType,
3051                       SMESH_subMesh*                  subMesh,
3052                       SMESH_subMeshEventListenerData* data,
3053                       const SMESH_Hypothesis*         hyp)
3054     {
3055       if (SMESH_subMesh::ALGO_EVENT == eventType &&
3056           !subMesh->GetAlgo() )
3057       {
3058         removeEmptyGroupsOfDomains( subMesh->GetFather(), /*notEmptyAsWell=*/true );
3059       }
3060     }
3061   };
3062 }
3063
3064 //================================================================================
3065 /*!
3066  * \brief Set an event listener to set enforced elements as soon as an enforced
3067  *        mesh is loaded
3068  */
3069 //================================================================================
3070
3071 void GHS3DPlugin_GHS3D::SubmeshRestored(SMESH_subMesh* subMesh)
3072 {
3073   if ( GHS3DPlugin_Hypothesis* hyp = _EnforcedMeshRestorer::GetGHSHypothesis( subMesh ))
3074   {
3075     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList enfMeshes = hyp->_GetEnforcedMeshes();
3076     GHS3DPlugin_Hypothesis::TGHS3DEnforcedMeshList::iterator it = enfMeshes.begin();
3077     for(;it != enfMeshes.end();++it) {
3078       GHS3DPlugin_Hypothesis::TGHS3DEnforcedMesh* enfMesh = *it;
3079       if ( SMESH_Mesh* mesh = GetMeshByPersistentID( enfMesh->persistID ))
3080       {
3081         SMESH_subMesh* smToListen = mesh->GetSubMesh( mesh->GetShapeToMesh() );
3082         // a listener set to smToListen will care of hypothesis stored in SMESH_EventListenerData
3083         subMesh->SetEventListener( new _EnforcedMeshRestorer(),
3084                                    SMESH_subMeshEventListenerData::MakeData( subMesh ),
3085                                    smToListen);
3086       }
3087     }
3088   }
3089 }
3090
3091 //================================================================================
3092 /*!
3093  * \brief Sets an event listener removing empty groups created due to "To make
3094  *        groups of domains".
3095  * \param subMesh - submesh where algo is set
3096  *
3097  * This method is called when a submesh gets HYP_OK algo_state.
3098  * After being set, event listener is notified on each event of a submesh.
3099  */
3100 //================================================================================
3101
3102 void GHS3DPlugin_GHS3D::SetEventListener(SMESH_subMesh* subMesh)
3103 {
3104   subMesh->SetEventListener( new _GroupsOfDomainsRemover(), 0, subMesh );
3105 }