1 // Copyright (C) 2004-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : GHS3DPlugin_Optimizer.cxx
20 // Created : Mon Oct 31 20:05:28 2016
22 #include "GHS3DPlugin_Optimizer.hxx"
24 #include "GHS3DPlugin_GHS3D.hxx"
25 #include "GHS3DPlugin_OptimizerHypothesis.hxx"
26 #include "MG_Tetra_API.hxx"
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>
34 #include <TopoDS_Shape.hxx>
36 //================================================================================
40 //================================================================================
42 GHS3DPlugin_Optimizer::GHS3DPlugin_Optimizer(int hypId, SMESH_Gen* gen)
43 : SMESH_3D_Algo(hypId, gen)
46 _compatibleHypothesis.push_back( GHS3DPlugin_OptimizerHypothesis::GetHypType());
47 _requireShape = false; // work without shape only
50 //================================================================================
52 * \brief Get a hypothesis
54 //================================================================================
56 bool GHS3DPlugin_Optimizer::CheckHypothesis(SMESH_Mesh& aMesh,
57 const TopoDS_Shape& aShape,
58 Hypothesis_Status& aStatus)
62 const std::list <const SMESHDS_Hypothesis * >& hyps = GetUsedHypothesis(aMesh, aShape);
64 _hyp = dynamic_cast< const GHS3DPlugin_OptimizerHypothesis* >( hyps.front() );
70 //================================================================================
72 * \brief Terminate a process
74 //================================================================================
76 void GHS3DPlugin_Optimizer::CancelCompute()
78 _computeCanceled = true;
81 //================================================================================
83 * \brief Evaluate nb of elements
85 //================================================================================
87 bool GHS3DPlugin_Optimizer::Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape,
88 MapShapeNbElems& aResMap)
92 //================================================================================
96 //================================================================================
98 bool GHS3DPlugin_Optimizer::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
105 //================================================================================
107 * \brief Compute average size of tetrahedra araound a node
109 //================================================================================
111 double getSizeAtNode( const SMDS_MeshNode* node )
117 vt.SetExternalNormal();
119 SMDS_ElemIteratorPtr volIt = node->GetInverseElementIterator(SMDSAbs_Volume);
120 while ( volIt->more() )
122 const SMDS_MeshElement* vol = volIt->next();
123 if ( vol->GetGeomType() != SMDSGeom_TETRA )
127 int dN = vol->IsQuadratic() ? 2 : 1;
129 for ( int iF = 0, nbF = vt.NbFaces(); iF < nbF; ++iF )
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
135 volSize += SMESH_TNodeXYZ( nodes[ iN ]).Distance( nodes[ iN + dN ]);
139 size += volSize / nbLinks;
142 return nbTet > 0 ? size/nbTet : 0.;
145 //================================================================================
147 * \brief Write a mesh file and a solution file
149 //================================================================================
151 bool writeGMFFile( MG_Tetra_API* theMGInput,
152 SMESH_MesherHelper* theHelper,
153 const std::string& theMeshFileName,
154 const std::string& theSolFileName )
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 )
163 // write all nodes and volume size at them
165 SMESHDS_Mesh* meshDS = theHelper->GetMeshDS();
166 if ( meshDS->NbNodes() != meshDS->MaxNodeID() )
167 meshDS->CompactMesh();
169 theMGInput->GmfSetKwd( mfile, GmfVertices, meshDS->NbNodes() );
170 int TypTab[] = { GmfSca };
171 theMGInput->GmfSetKwd( sfile, GmfSolAtVertices, meshDS->NbNodes(), 1, TypTab);
173 SMDS_NodeIteratorPtr nodeIt = theHelper->GetMeshDS()->nodesIterator();
174 while ( nodeIt->more() )
176 const SMDS_MeshNode* node = nodeIt->next();
177 theMGInput->GmfSetLin( mfile, GmfVertices, node->X(), node->Y(), node->Z(), tag );
179 double size = getSizeAtNode( node );
180 theMGInput->GmfSetLin( sfile, GmfSolAtVertices, &size );
183 // write all triangles
185 theMGInput->GmfSetKwd( mfile, GmfTriangles, meshDS->GetMeshInfo().NbTriangles() );
186 SMDS_ElemIteratorPtr triaIt = meshDS->elementGeomIterator( SMDSGeom_TRIANGLE );
187 while ( triaIt->more() )
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(),
199 theMGInput->GmfSetKwd( mfile, GmfTetrahedra, meshDS->GetMeshInfo().NbTetras() );
200 SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
201 while ( tetIt->more() )
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(),
212 theMGInput->GmfCloseMesh( mfile );
213 theMGInput->GmfCloseMesh( sfile );
218 //================================================================================
220 * \brief Read optimized tetrahedra
222 //================================================================================
224 bool readGMFFile( MG_Tetra_API* theMGOutput,
225 SMESH_MesherHelper* theHelper,
226 const std::string& theMeshFileName )
229 int inFile = theMGOutput->GmfOpenMesh( theMeshFileName.c_str(), GmfRead, &ver, &dim);
233 SMESHDS_Mesh* meshDS = theHelper->GetMeshDS();
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;
246 theMGOutput->GmfGotoKwd( inFile, GmfVertices );
248 if ( nbNodes == nbNodesOld && nbTet == nbTetOld )
251 SMDS_NodeIteratorPtr nodeIt = theHelper->GetMeshDS()->nodesIterator();
252 while ( nodeIt->more() )
254 const SMDS_MeshNode* node = nodeIt->next();
255 theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
256 meshDS->MoveNode( node, x,y,z );
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 )
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 );
274 else if ( nbNodes >= nbNodesOld ) // tetra added/removed
277 for ( int iN = 1; iN <= nbNodes; ++iN )
279 theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
280 const SMDS_MeshNode* node = meshDS->FindNode( iN );
282 node = meshDS->AddNode( x,y,z );
284 meshDS->MoveNode( node, x,y,z );
288 SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
289 while ( tetIt->more() )
290 meshDS->RemoveFreeElement( tetIt->next(), /*sm=*/0 );
293 theMGOutput->GmfGotoKwd( inFile, GmfTetrahedra );
294 for ( int i = 0; i < nbTet; ++i )
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 ));
303 else if ( nbNodes < nbNodesOld ) // nodes and tetra removed
306 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
307 while ( nIt->more() )
308 nIt->next()->setIsMarked( false );
310 // remove tetrahedra and mark nodes used by them
311 SMDS_ElemIteratorPtr tetIt = meshDS->elementGeomIterator( SMDSGeom_TETRA );
312 while ( tetIt->more() )
314 const SMDS_MeshElement* tet = tetIt->next();
315 SMDS_ElemIteratorPtr nIter = tet->nodesIterator();
316 while ( nIter->more() )
317 nIter->next()->setIsMarked( true );
319 meshDS->RemoveFreeElement( tet, /*sm=*/0 );
323 for ( int iN = 1; iN <= nbNodes; ++iN )
325 theMGOutput->GmfGetLin( inFile, GmfVertices, &x, &y, &z, &tag );
326 const SMDS_MeshNode* node = meshDS->FindNode( iN );
328 node = meshDS->AddNode( x,y,z );
330 meshDS->MoveNode( node, x,y,z );
334 theMGOutput->GmfGotoKwd( inFile, GmfTetrahedra );
335 for ( int i = 0; i < nbTet; ++i )
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 ));
344 // remove free marked nodes
345 for ( nIt = meshDS->nodesIterator(); nIt->more(); )
347 const SMDS_MeshNode* node = nIt->next();
348 if ( node->NbInverseElements() == 0 && node->isMarked() )
349 meshDS->RemoveFreeNode( node, 0 );
353 // avoid "No mesh elements assigned to a sub-shape" error
354 theHelper->GetMesh()->GetSubMesh( theHelper->GetSubShape() )->SetIsAlwaysComputed( true );
359 void getNodeByGhsId( SMESH_Mesh& mesh, std::vector <const SMDS_MeshNode*> & nodeByGhsId )
361 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
362 const int nbNodes = meshDS->NbNodes();
363 nodeByGhsId.resize( nbNodes + 1 );
364 SMDS_NodeIteratorPtr nodeIt = meshDS->nodesIterator();
365 while ( nodeIt->more() )
367 const SMDS_MeshNode* node = nodeIt->next();
368 if ( node->GetID() <= nbNodes )
369 nodeByGhsId[ node->GetID() ] = node;
372 throw SALOME_Exception(LOCALIZED ("bad ID -- not compacted mesh"));
377 void removeFile(const std::string& name)
379 SMESH_File( name ).remove();
383 //================================================================================
385 * \brief Optimize an existing tetrahedral mesh
387 //================================================================================
389 bool GHS3DPlugin_Optimizer::Compute(SMESH_Mesh& theMesh,
390 SMESH_MesherHelper* theHelper)
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" );
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
406 MG_Tetra_API mgTetra( _computeCanceled, _progress );
408 bool Ok = writeGMFFile( &mgTetra, theHelper, aGMFFileName, aSolFileName );
411 // run MG-Tetra mesher
414 bool logInStandardOutput = _hyp ? _hyp->GetStandardOutputLog() : false;
415 bool removeLogOnSuccess = _hyp ? _hyp->GetRemoveLogOnSuccess() : true;
416 bool keepFiles = _hyp ? _hyp->GetKeepFiles() : false;
418 std::string cmd = GHS3DPlugin_OptimizerHypothesis::CommandToRun( _hyp );
420 if ( mgTetra.IsExecutable() )
422 cmd += " --in " + aGMFFileName;
423 cmd += " --out " + aResultFileName;
425 if ( !logInStandardOutput )
427 mgTetra.SetLogFile( aLogFileName.c_str() );
428 cmd += " 1>" + aLogFileName; // dump into file
430 std::cout << std::endl;
431 std::cout << "MG-Tetra execution..." << std::endl;
432 std::cout << cmd << std::endl;
434 _computeCanceled = false;
437 Ok = mgTetra.Compute( cmd, errStr ); // run
439 if ( logInStandardOutput && mgTetra.IsLibrary() )
440 std::cout << std::endl << mgTetra.GetLog() << std::endl;
442 std::cout << std::endl << "End of MG-Tetra execution !" << std::endl;
447 Ok = Ok && readGMFFile( &mgTetra, theHelper, aResultFileName );
449 // ---------------------
450 // remove working files
451 // ---------------------
452 if ( mgTetra.HasLog() )
454 if( _computeCanceled )
455 error( "interruption initiated by user" );
458 // get problem description from the log file
459 std::vector <const SMDS_MeshNode*> nodeByGhsId;
460 getNodeByGhsId( theMesh, nodeByGhsId );
461 _Ghs2smdsConvertor conv( nodeByGhsId, SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( theMesh )));
462 error( GHS3DPlugin_GHS3D::getErrorDescription( logInStandardOutput ? 0 : aLogFileName.c_str(),
463 mgTetra.GetLog(), conv, Ok ));
467 // the log file is empty
468 removeFile( aLogFileName );
469 INFOS( "MG-Tetra Error, " << errStr);
470 error(COMPERR_ALGO_FAILED, errStr);
473 if ( Ok && removeLogOnSuccess )
475 removeFile( aLogFileName );
479 if ( !Ok && _computeCanceled )
480 removeFile( aLogFileName );
481 removeFile( aGMFFileName );
482 removeFile( aSolFileName );
483 removeFile( aResultFileName );
484 removeFile( aResSolFileName );
490 // double GHS3DPlugin_Optimizer::GetProgress() const