]> SALOME platform Git repositories - plugins/ghs3dplugin.git/blob - src/GHS3DPlugin/GHS3DPlugin_Optimizer.cxx
Salome HOME
ce54b2162d4c397ec69496f7bcbc0802b8b32173
[plugins/ghs3dplugin.git] / 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, 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, 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 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 int 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   std::cout << std::endl;
431   std::cout << "MG-Tetra execution..." << std::endl;
432   std::cout << cmd << std::endl;
433
434   _computeCanceled = false;
435
436   std::string errStr;
437   Ok = mgTetra.Compute( cmd, errStr ); // run
438
439   if ( logInStandardOutput && mgTetra.IsLibrary() )
440     std::cout << std::endl << mgTetra.GetLog() << std::endl;
441   if ( Ok )
442     std::cout << std::endl << "End of MG-Tetra execution !" << std::endl;
443
444   // --------------
445   // read a result
446   // --------------
447   Ok = Ok && readGMFFile( &mgTetra, theHelper, aResultFileName );
448
449   // ---------------------
450   // remove working files
451   // ---------------------
452   if ( mgTetra.HasLog() )
453   {
454     if( _computeCanceled )
455       error( "interruption initiated by user" );
456     else
457     {
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 ));
464     }
465   }
466   else {
467     // the log file is empty
468     removeFile( aLogFileName );
469     INFOS( "MG-Tetra Error, " << errStr);
470     error(COMPERR_ALGO_FAILED, errStr);
471   }
472
473   if ( Ok && removeLogOnSuccess )
474   {
475     removeFile( aLogFileName );
476   }
477   if ( !keepFiles )
478   {
479     if ( !Ok && _computeCanceled )
480       removeFile( aLogFileName );
481     removeFile( aGMFFileName );
482     removeFile( aSolFileName );
483     removeFile( aResultFileName );
484     removeFile( aResSolFileName );
485   }
486   return Ok;
487 }
488
489
490 // double GHS3DPlugin_Optimizer::GetProgress() const
491 // {
492 // }