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