src/GHS3DPlugin/GHS3DPlugin_Optimizer.cxx

   1 // Copyright (C) 2004-2016  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // 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, int studyId, SMESH_Gen* gen)
  43   : SMESH_3D_Algo(hypId, studyId, 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, meshDS->NbNodes() );
 170     int TypTab[] = { GmfSca };
 171     theMGInput->GmfSetKwd( sfile, GmfSolAtVertices, 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, 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                              tria->GetNode(0)->GetID(),
 192                              tria->GetNode(1)->GetID(),
 193                              tria->GetNode(2)->GetID(),
 194                              tag );
 195     }
 196
 197     // write all tetra
 198
 199     theMGInput->GmfSetKwd( mfile, GmfTetrahedra, 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                              tet->GetNode(0)->GetID(),
 206                              tet->GetNode(2)->GetID(),
 207                              tet->GetNode(1)->GetID(),
 208                              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     int nbNodesOld = meshDS->NbNodes();
 238     int   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 // 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
 304     // avoid "No mesh elements assigned to a sub-shape" error
 305     theHelper->GetMesh()->GetSubMesh( theHelper->GetSubShape() )->SetIsAlwaysComputed( true );
 306
 307     return true;
 308   }
 309
 310   void getNodeByGhsId( SMESH_Mesh& mesh, std::vector <const SMDS_MeshNode*> & nodeByGhsId )
 311   {
 312     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 313     const int nbNodes = meshDS->NbNodes();
 314     nodeByGhsId.resize( nbNodes + 1 );
 315     SMDS_NodeIteratorPtr nodeIt = meshDS->nodesIterator();
 316     while ( nodeIt->more() )
 317     {
 318       const SMDS_MeshNode* node = nodeIt->next();
 319       if ( node->GetID() <= nbNodes )
 320         nodeByGhsId[ node->GetID() ] = node;
 321 #ifdef _DEBUG_
 322       else
 323         throw SALOME_Exception(LOCALIZED ("bad ID -- not compacted mesh"));
 324 #endif
 325     }
 326   }
 327
 328   void removeFile(const std::string& name)
 329   {
 330     SMESH_File( name ).remove();
 331   }
 332 }
 333
 334 //================================================================================
 335 /*!
 336  * \brief Optimize an existing tetrahedral mesh
 337  */
 338 //================================================================================
 339
 340 bool GHS3DPlugin_Optimizer::Compute(SMESH_Mesh&         theMesh,
 341                                     SMESH_MesherHelper* theHelper)
 342 {
 343   if ( theMesh.NbTetras() == 0 )
 344     return error( COMPERR_BAD_INPUT_MESH, "No tetrahedra" );
 345   if ( theMesh.NbTetras( ORDER_QUADRATIC ) > 0 )
 346     return error( COMPERR_BAD_INPUT_MESH, "Quadratic mesh can't be optimized" );
 347   if ( theMesh.NbTriangles() == 0 )
 348     return error( COMPERR_BAD_INPUT_MESH, "2D mesh must exist around tetrahedra" );
 349
 350   std::string aGenericName    = GHS3DPlugin_Hypothesis::GetFileName(_hyp);
 351   std::string aLogFileName    = aGenericName + ".log";  // log
 352   std::string aGMFFileName    = aGenericName + ".mesh"; // input GMF mesh file
 353   std::string aSolFileName    = aGenericName + ".sol";  // input size map file
 354   std::string aResultFileName = aGenericName + "_Opt.mesh";  // out GMF mesh file
 355   std::string aResSolFileName = aGenericName + "_Opt.sol";   // out size map file
 356
 357   MG_Tetra_API mgTetra( _computeCanceled, _progress );
 358
 359   bool Ok = writeGMFFile( &mgTetra, theHelper, aGMFFileName, aSolFileName );
 360
 361   // -----------------
 362   // run MG-Tetra mesher
 363   // -----------------
 364
 365   bool logInStandardOutput = _hyp ? _hyp->GetStandardOutputLog() : false;
 366   bool removeLogOnSuccess  = _hyp ? _hyp->GetRemoveLogOnSuccess() : true;
 367   bool keepFiles           = _hyp ? _hyp->GetKeepFiles() : false;
 368
 369   std::string cmd = GHS3DPlugin_OptimizerHypothesis::CommandToRun( _hyp );
 370
 371   if ( mgTetra.IsExecutable() )
 372   {
 373     cmd += " --in " + aGMFFileName;
 374     cmd += " --out " + aResultFileName;
 375   }
 376   if ( !logInStandardOutput )
 377   {
 378     mgTetra.SetLogFile( aLogFileName.c_str() );
 379     cmd += " 1>" + aLogFileName;  // dump into file
 380   }
 381   std::cout << std::endl;
 382   std::cout << "MG-Tetra execution..." << std::endl;
 383   std::cout << cmd << std::endl;
 384
 385   _computeCanceled = false;
 386
 387   std::string errStr;
 388   Ok = mgTetra.Compute( cmd, errStr ); // run
 389
 390   if ( logInStandardOutput && mgTetra.IsLibrary() )
 391     std::cout << std::endl << mgTetra.GetLog() << std::endl;
 392   if ( Ok )
 393     std::cout << std::endl << "End of MG-Tetra execution !" << std::endl;
 394
 395   // --------------
 396   // read a result
 397   // --------------
 398   Ok = Ok && readGMFFile( &mgTetra, theHelper, aResultFileName );
 399
 400   // ---------------------
 401   // remove working files
 402   // ---------------------
 403   if ( Ok )
 404   {
 405     if ( removeLogOnSuccess )
 406       removeFile( aLogFileName );
 407   }
 408   if ( mgTetra.HasLog() )
 409   {
 410     if( _computeCanceled )
 411       error( "interruption initiated by user" );
 412     else
 413     {
 414       // get problem description from the log file
 415       std::vector <const SMDS_MeshNode*> nodeByGhsId;
 416       getNodeByGhsId( theMesh, nodeByGhsId );
 417       _Ghs2smdsConvertor conv( nodeByGhsId, SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( theMesh )));
 418       error( GHS3DPlugin_GHS3D::getErrorDescription( logInStandardOutput ? 0 : aLogFileName.c_str(),
 419                                                      mgTetra.GetLog(), conv, Ok ));
 420     }
 421   }
 422   else {
 423     // the log file is empty
 424     removeFile( aLogFileName );
 425     INFOS( "MG-Tetra Error, " << errStr);
 426     error(COMPERR_ALGO_FAILED, errStr);
 427   }
 428
 429   if ( !keepFiles )
 430   {
 431     if ( !Ok && _computeCanceled )
 432       removeFile( aLogFileName );
 433     removeFile( aGMFFileName );
 434     removeFile( aSolFileName );
 435     removeFile( aResultFileName );
 436     removeFile( aResSolFileName );
 437   }
 438   return Ok;
 439 }
 440
 441
 442 // double GHS3DPlugin_Optimizer::GetProgress() const
 443 // {
 444 // }