src/GHS3DPlugin/GHS3DPlugin_Optimizer.cxx

   1 // Copyright (C) 2004-2024  CEA, EDF
   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 // File      : GHS3DPlugin_Optimizer.cxx
  20 // Created   : Mon Oct 31 20:05:28 2016
  21
  22 #include "GHS3DPlugin_Optimizer.hxx"
  23
  24 #include "GHS3DPlugin_GHS3D.hxx"
  25 #include "GHS3DPlugin_OptimizerHypothesis.hxx"
  26 #include "MG_Tetra_API.hxx"
  27
  28 #include <SMDS_VolumeTool.hxx>
  29 #include <SMESHDS_Mesh.hxx>
  30 #include <SMESH_File.hxx>
  31 #include <SMESH_Mesh.hxx>
  32 #include <SMESH_MesherHelper.hxx>
  33
  34 #include <TopoDS_Shape.hxx>
  35
  36 //================================================================================
  37 /*!
  38  * \brief Constructor
  39  */
  40 //================================================================================
  41
  42 GHS3DPlugin_Optimizer::GHS3DPlugin_Optimizer(int hypId, SMESH_Gen* gen)
  43   : SMESH_3D_Algo(hypId, gen)
  44 {
  45   _name = Name();
  46   _compatibleHypothesis.push_back( GHS3DPlugin_OptimizerHypothesis::GetHypType());
  47   _requireShape = false; // work without shape only
  48 }
  49
  50 //================================================================================
  51 /*!
  52  * \brief Get a hypothesis
  53  */
  54 //================================================================================
  55
  56 bool GHS3DPlugin_Optimizer::CheckHypothesis(SMESH_Mesh&         aMesh,
  57                                             const TopoDS_Shape& aShape,
  58                                             Hypothesis_Status&  aStatus)
  59 {
  60   _hyp = NULL;
  61
  62   const std::list <const SMESHDS_Hypothesis * >& hyps = GetUsedHypothesis(aMesh, aShape);
  63   if ( !hyps.empty() )
  64     _hyp = dynamic_cast< const GHS3DPlugin_OptimizerHypothesis* >( hyps.front() );
  65
  66   aStatus = HYP_OK;
  67   return true;
  68 }
  69
  70 //================================================================================
  71 /*!
  72  * \brief Terminate a process
  73  */
  74 //================================================================================
  75
  76 void GHS3DPlugin_Optimizer::CancelCompute()
  77 {
  78   _computeCanceled = true;
  79 }
  80
  81 //================================================================================
  82 /*!
  83  * \brief Evaluate nb of elements
  84  */
  85 //================================================================================
  86
  87 bool GHS3DPlugin_Optimizer::Evaluate(SMESH_Mesh& /*aMesh*/, const TopoDS_Shape& /*aShape*/,
  88                                      MapShapeNbElems& /*aResMap*/)
  89 {
  90   return false;
  91 }
  92 //================================================================================
  93 /*!
  94  * \brief Does nothing
  95  */
  96 //================================================================================
  97
  98 bool GHS3DPlugin_Optimizer::Compute(SMESH_Mesh & /*aMesh*/, const TopoDS_Shape & /*aShape*/)
  99 {
 100   return false;
 101 }
 102
 103 namespace
 104 {
 105   //================================================================================
 106   /*!
 107    * \brief Compute average size of tetrahedra araound a node
 108    */
 109   //================================================================================
 110
 111   double getSizeAtNode( const SMDS_MeshNode* node )
 112   {
 113     double size = 0;
 114     int nbTet = 0;
 115
 116     SMDS_VolumeTool vt;
 117     vt.SetExternalNormal();
 118
 119     SMDS_ElemIteratorPtr volIt = node->GetInverseElementIterator(SMDSAbs_Volume);
 120     while ( volIt->more() )
 121     {
 122       const SMDS_MeshElement* vol = volIt->next();
 123       if ( vol->GetGeomType() != SMDSGeom_TETRA )
 124         continue;
 125       double volSize = 0.;
 126       int    nbLinks = 0;
 127       int dN = vol->IsQuadratic() ? 2 : 1;
 128       vt.Set( vol );
 129       for ( int iF = 0, nbF = vt.NbFaces(); iF < nbF; ++iF )
 130       {
 131         const SMDS_MeshNode** nodes = vt. GetFaceNodes( iF );
 132         for ( int iN = 0, nbN = vt.NbFaceNodes( iF ); iN < nbN; iN += dN )
 133           if ( nodes[ iN ] < nodes[ iN + dN ]) // use a link once
 134           {
 135             volSize += SMESH_TNodeXYZ( nodes[ iN ]).Distance( nodes[ iN + dN ]);
 136             ++nbLinks;
 137           }
 138       }
 139       size += volSize / nbLinks;
 140       ++nbTet;
 141     }
 142     return nbTet > 0 ? size/nbTet : 0.;
 143   }
 144
 145   //================================================================================
 146   /*!
 147    * \brief Write a mesh file and a solution file
 148    */
 149   //================================================================================
 150
 151   bool writeGMFFile( MG_Tetra_API*       theMGInput,
 152                      SMESH_MesherHelper* theHelper,
 153                      const std::string&  theMeshFileName,
 154                      const std::string&  theSolFileName )
 155   {
 156     const int tag = 0;
 157
 158     int mfile = theMGInput->GmfOpenMesh( theMeshFileName.c_str(), GmfWrite, GMFVERSION,GMFDIMENSION );
 159     int sfile = theMGInput->GmfOpenMesh( theSolFileName.c_str(), GmfWrite, GMFVERSION,GMFDIMENSION );
 160     if ( !mfile || !sfile )
 161       return false;
 162
 163     // write all nodes and volume size at them
 164
 165     SMESHDS_Mesh* meshDS = theHelper->GetMeshDS();
 166     if ( meshDS->NbNodes() != meshDS->MaxNodeID() )
 167       meshDS->CompactMesh();
 168
 169     theMGInput->GmfSetKwd( mfile, GmfVertices, (int) meshDS->NbNodes() );
 170     int TypTab[] = { GmfSca };
 171     theMGInput->GmfSetKwd( sfile, GmfSolAtVertices, (int) meshDS->NbNodes(), 1, TypTab);
 172
 173     SMDS_NodeIteratorPtr nodeIt = theHelper->GetMeshDS()->nodesIterator();
 174     while ( nodeIt->more() )
 175     {
 176       const SMDS_MeshNode* node = nodeIt->next();
 177       theMGInput->GmfSetLin( mfile, GmfVertices, node->X(), node->Y(), node->Z(), tag );
 178
 179       double size = getSizeAtNode( node );
 180       theMGInput->GmfSetLin( sfile, GmfSolAtVertices, &size );
 181     }
 182
 183     // write all triangles
 184
 185     theMGInput->GmfSetKwd( mfile, GmfTriangles, (int) meshDS->GetMeshInfo().NbTriangles() );
 186     SMDS_ElemIteratorPtr triaIt = meshDS->elementGeomIterator( SMDSGeom_TRIANGLE );
 187     while ( triaIt->more() )
 188     {
 189       const SMDS_MeshElement* tria = triaIt->next();
 190       theMGInput->GmfSetLin( mfile, GmfTriangles,
 191                              static_cast<int>(tria->GetNode(0)->GetID()),
 192                              static_cast<int>(tria->GetNode(1)->GetID()),
 193                              static_cast<int>(tria->GetNode(2)->GetID()),
 194                              tag );
 195     }
 196
 197     // write all tetra
 198
 199     theMGInput->GmfSetKwd( mfile, GmfTetrahedra, (int) meshDS->GetMeshInfo().NbTetras() );
 200     SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
 201     while ( tetIt->more() )
 202     {
 203       const SMDS_MeshElement* tet = tetIt->next();
 204       theMGInput->GmfSetLin( mfile, GmfTetrahedra,
 205                              static_cast<int>(tet->GetNode(0)->GetID()),
 206                              static_cast<int>(tet->GetNode(2)->GetID()),
 207                              static_cast<int>(tet->GetNode(1)->GetID()),
 208                              static_cast<int>(tet->GetNode(3)->GetID()),
 209                              tag );
 210     }
 211
 212     theMGInput->GmfCloseMesh( mfile );
 213     theMGInput->GmfCloseMesh( sfile );
 214
 215     return true;
 216   }
 217
 218   //================================================================================
 219   /*!
 220    * \brief Read optimized tetrahedra
 221    */
 222   //================================================================================
 223
 224   bool readGMFFile( MG_Tetra_API*       theMGOutput,
 225                     SMESH_MesherHelper* theHelper,
 226                     const std::string&  theMeshFileName )
 227   {
 228     int ver, dim, tag;
 229     int inFile = theMGOutput->GmfOpenMesh( theMeshFileName.c_str(), GmfRead, &ver, &dim);
 230     if ( !inFile )
 231       return false;
 232
 233     SMESHDS_Mesh* meshDS = theHelper->GetMeshDS();
 234
 235     int nbNodes = theMGOutput->GmfStatKwd( inFile, GmfVertices );
 236     int   nbTet = theMGOutput->GmfStatKwd( inFile, GmfTetrahedra );
 237     smIdType nbNodesOld = meshDS->NbNodes();
 238     smIdType nbTetOld = meshDS->GetMeshInfo().NbTetras();
 239     std::cout << "Optimization input:  "
 240               << nbNodesOld << " nodes, \t" << nbTetOld << " tetra" << std::endl;
 241     std::cout << "Optimization output: "
 242               << nbNodes << " nodes, \t" << nbTet << " tetra" << std::endl;
 243
 244     double x,y,z;
 245     int i1, i2, i3, i4;
 246     theMGOutput->GmfGotoKwd( inFile, GmfVertices );
 247
 248     if ( nbNodes == nbNodesOld && nbTet == nbTetOld )
 249     {
 250       // move nodes
 251       SMDS_NodeIteratorPtr nodeIt = theHelper->GetMeshDS()->nodesIterator();
 252       while ( nodeIt->more() )
 253       {
 254         const SMDS_MeshNode* node = nodeIt->next();
 255         theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
 256         meshDS->MoveNode( node, x,y,z );
 257       }
 258
 259       // update tetra
 260       const SMDS_MeshNode* nodes[ 4 ];
 261       theMGOutput->GmfGotoKwd( inFile, GmfTetrahedra );
 262       SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
 263       for ( int i = 0; i < nbTet; ++i )
 264       {
 265         const SMDS_MeshElement* tet = tetIt->next();
 266         theMGOutput->GmfGetLin( inFile, GmfTetrahedra, &i1, &i2, &i3, &i4, &tag);
 267         nodes[ 0 ] = meshDS->FindNode( i1 );
 268         nodes[ 1 ] = meshDS->FindNode( i3 );
 269         nodes[ 2 ] = meshDS->FindNode( i2 );
 270         nodes[ 3 ] = meshDS->FindNode( i4 );
 271         meshDS->ChangeElementNodes( tet, &nodes[0], 4 );
 272       }
 273     }
 274     else if ( nbNodes >= nbNodesOld ) // tetra added/removed
 275     {
 276       // move or add nodes
 277       for ( int iN = 1; iN <= nbNodes; ++iN )
 278       {
 279         theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
 280         const SMDS_MeshNode* node = meshDS->FindNode( iN );
 281         if ( !node )
 282           node = meshDS->AddNode( x,y,z );
 283         else
 284           meshDS->MoveNode( node, x,y,z );
 285       }
 286
 287       // remove tetrahedra
 288       SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
 289       while ( tetIt->more() )
 290         meshDS->RemoveFreeElement( tetIt->next(), /*sm=*/0 );
 291
 292       // add tetrahedra
 293       theMGOutput->GmfGotoKwd( inFile, GmfTetrahedra );
 294       for ( int i = 0; i < nbTet; ++i )
 295       {
 296         theMGOutput->GmfGetLin( inFile, GmfTetrahedra, &i1, &i2, &i3, &i4, &tag);
 297         meshDS->AddVolume( meshDS->FindNode( i1 ),
 298                            meshDS->FindNode( i3 ),
 299                            meshDS->FindNode( i2 ),
 300                            meshDS->FindNode( i4 ));
 301       }
 302     }
 303     else if ( nbNodes < nbNodesOld ) // nodes and tetra removed
 304     {
 305       // unmark nodes
 306       SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
 307       while ( nIt->more() )
 308         nIt->next()->setIsMarked( false );
 309
 310       // remove tetrahedra and mark nodes used by them
 311       SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
 312       while ( tetIt->more() )
 313       {
 314         const SMDS_MeshElement* tet = tetIt->next();
 315         SMDS_ElemIteratorPtr nIter = tet->nodesIterator();
 316         while ( nIter->more() )
 317           nIter->next()->setIsMarked( true );
 318
 319         meshDS->RemoveFreeElement( tet, /*sm=*/0 );
 320       }
 321
 322       // move or add nodes
 323       for ( int iN = 1; iN <= nbNodes; ++iN )
 324       {
 325         theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
 326         const SMDS_MeshNode* node = meshDS->FindNode( iN );
 327         if ( !node )
 328           node = meshDS->AddNode( x,y,z );
 329         else
 330           meshDS->MoveNode( node, x,y,z );
 331       }
 332
 333       // add tetrahedra
 334       theMGOutput->GmfGotoKwd( inFile, GmfTetrahedra );
 335       for ( int i = 0; i < nbTet; ++i )
 336       {
 337         theMGOutput->GmfGetLin( inFile, GmfTetrahedra, &i1, &i2, &i3, &i4, &tag);
 338         meshDS->AddVolume( meshDS->FindNode( i1 ),
 339                            meshDS->FindNode( i3 ),
 340                            meshDS->FindNode( i2 ),
 341                            meshDS->FindNode( i4 ));
 342       }
 343
 344       // remove free marked nodes
 345       for ( nIt = meshDS->nodesIterator(); nIt->more(); )
 346       {
 347         const SMDS_MeshNode* node = nIt->next();
 348         if ( node->NbInverseElements() == 0 && node->isMarked() )
 349           meshDS->RemoveFreeNode( node, 0 );
 350       }
 351     }
 352
 353     // avoid "No mesh elements assigned to a sub-shape" error
 354     theHelper->GetMesh()->GetSubMesh( theHelper->GetSubShape() )->SetIsAlwaysComputed( true );
 355
 356     return true;
 357   }
 358
 359   void getNodeByGhsId( SMESH_Mesh& mesh, std::vector <const SMDS_MeshNode*> & nodeByGhsId )
 360   {
 361     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 362     const smIdType nbNodes = meshDS->NbNodes();
 363     nodeByGhsId.resize( nbNodes + 1 );
 364     SMDS_NodeIteratorPtr nodeIt = meshDS->nodesIterator();
 365     while ( nodeIt->more() )
 366     {
 367       const SMDS_MeshNode* node = nodeIt->next();
 368       if ( node->GetID() <= nbNodes )
 369         nodeByGhsId[ node->GetID() ] = node;
 370 #ifdef _DEBUG_
 371       else
 372         throw SALOME_Exception(LOCALIZED ("bad ID -- not compacted mesh"));
 373 #endif
 374     }
 375   }
 376
 377   void removeFile(const std::string& name)
 378   {
 379     SMESH_File( name ).remove();
 380   }
 381 }
 382
 383 //================================================================================
 384 /*!
 385  * \brief Optimize an existing tetrahedral mesh
 386  */
 387 //================================================================================
 388
 389 bool GHS3DPlugin_Optimizer::Compute(SMESH_Mesh&         theMesh,
 390                                     SMESH_MesherHelper* theHelper)
 391 {
 392   if ( theMesh.NbTetras() == 0 )
 393     return error( COMPERR_BAD_INPUT_MESH, "No tetrahedra" );
 394   if ( theMesh.NbTetras( ORDER_QUADRATIC ) > 0 )
 395     return error( COMPERR_BAD_INPUT_MESH, "Quadratic mesh can't be optimized" );
 396   if ( theMesh.NbTriangles() == 0 )
 397     return error( COMPERR_BAD_INPUT_MESH, "2D mesh must exist around tetrahedra" );
 398
 399   std::string aGenericName    = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
 400   std::string aLogFileName    = aGenericName + ".log";  // log
 401   std::string aGMFFileName    = aGenericName + ".mesh"; // input GMF mesh file
 402   std::string aSolFileName    = aGenericName + ".sol";  // input size map file
 403   std::string aResultFileName = aGenericName + "_Opt.mesh";  // out GMF mesh file
 404   std::string aResSolFileName = aGenericName + "_Opt.sol";   // out size map file
 405
 406   MG_Tetra_API mgTetra( _computeCanceled, _progress );
 407
 408   bool Ok = writeGMFFile( &mgTetra, theHelper, aGMFFileName, aSolFileName );
 409
 410   // -----------------
 411   // run MG-Tetra mesher
 412   // -----------------
 413
 414   bool logInStandardOutput = _hyp ? _hyp->GetStandardOutputLog() : false;
 415   bool removeLogOnSuccess  = _hyp ? _hyp->GetRemoveLogOnSuccess() : true;
 416   bool keepFiles           = _hyp ? _hyp->GetKeepFiles() : false;
 417
 418   std::string cmd = GHS3DPlugin_OptimizerHypothesis::CommandToRun( _hyp );
 419
 420   if ( mgTetra.IsExecutable() )
 421   {
 422     cmd += " --in " + aGMFFileName;
 423     cmd += " --out " + aResultFileName;
 424   }
 425   if ( !logInStandardOutput )
 426   {
 427     mgTetra.SetLogFile( aLogFileName.c_str() );
 428     cmd += " 1>" + aLogFileName;  // dump into file
 429   }
 430
 431   BRIEF_INFOS("")
 432   BRIEF_INFOS("MG-Tetra execution...")
 433   BRIEF_INFOS(cmd)
 434
 435   _computeCanceled = false;
 436
 437   std::string errStr;
 438   Ok = mgTetra.Compute( cmd, errStr ); // run
 439
 440   if ( logInStandardOutput && mgTetra.IsLibrary() ) {
 441     BRIEF_INFOS("");
 442     BRIEF_INFOS(mgTetra.GetLog());
 443     BRIEF_INFOS("")
 444   }
 445   if ( Ok ) {
 446     BRIEF_INFOS("");
 447     BRIEF_INFOS("End of MG-Tetra execution !");
 448     BRIEF_INFOS("")
 449   }
 450
 451
 452   // --------------
 453   // read a result
 454   // --------------
 455   Ok = Ok && readGMFFile( &mgTetra, theHelper, aResultFileName );
 456
 457   // ---------------------
 458   // remove working files
 459   // ---------------------
 460   if ( mgTetra.HasLog() )
 461   {
 462     if( _computeCanceled )
 463       error( "interruption initiated by user" );
 464     else
 465     {
 466       // get problem description from the log file
 467       std::vector <const SMDS_MeshNode*> nodeByGhsId;
 468       getNodeByGhsId( theMesh, nodeByGhsId );
 469       _Ghs2smdsConvertor conv( nodeByGhsId, SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( theMesh )));
 470       error( GHS3DPlugin_GHS3D::getErrorDescription( logInStandardOutput ? 0 : aLogFileName.c_str(),
 471                                                      mgTetra.GetLog(), conv, Ok ));
 472     }
 473   }
 474   else {
 475     // the log file is empty
 476     removeFile( aLogFileName );
 477     INFOS( "MG-Tetra Error, " << errStr);
 478     error(COMPERR_ALGO_FAILED, errStr);
 479   }
 480
 481   if ( Ok && removeLogOnSuccess )
 482   {
 483     removeFile( aLogFileName );
 484   }
 485   if ( !keepFiles )
 486   {
 487     if ( !Ok && _computeCanceled )
 488       removeFile( aLogFileName );
 489     removeFile( aGMFFileName );
 490     removeFile( aSolFileName );
 491     removeFile( aResultFileName );
 492     removeFile( aResSolFileName );
 493   }
 494   return Ok;
 495 }
 496
 497
 498 // double GHS3DPlugin_Optimizer::GetProgress() const
 499 // {
 500 // }