src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx

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