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