src/NETGENPlugin/NETGENPlugin_NETGEN_1D2D3D_SA.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   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 // NETGENPlugin : C++ implementation
  21 // File   : NETGENPlugin_NETGEN_2D3D_SA.cxx
  22 // Author : Cesar Conopoima (OCC)
  23 // Date   : 01/11/2023
  24 // Project   : SALOME
  25 //=============================================================================
  26 //
  27
  28 #include "NETGENPlugin_DriverParam.hxx"
  29 #include "NETGENPlugin_Hypothesis.hxx"
  30 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
  31
  32 #include "NETGENPlugin_NETGEN_2D.hxx"
  33 #include "NETGENPlugin_NETGEN_1D2D3D_SA.hxx"
  34
  35 #include <SMESHDS_Mesh.hxx>
  36 #include <SMESH_ControlsDef.hxx>
  37 #include <SMESH_Gen.hxx>
  38 #include <SMESH_Mesh.hxx>
  39 #include <SMESH_MesherHelper.hxx>
  40 #include <SMESH_DriverShape.hxx>
  41 #include <SMESH_DriverMesh.hxx>
  42 #include <utilities.h>
  43
  44 #include <list>
  45
  46 namespace nglib {
  47 #include <nglib.h>
  48 }
  49 #include <meshing.hpp>
  50
  51 using namespace nglib;
  52
  53 //=============================================================================
  54 /*!
  55  *
  56  */
  57 //=============================================================================
  58
  59 NETGENPlugin_NETGEN_1D2D3D_SA::NETGENPlugin_NETGEN_1D2D3D_SA()
  60   : NETGENPlugin_NETGEN_2D3D(0, new SMESH_Gen())
  61 {
  62   _name = "NETGEN_1D2D3D_SA";
  63 }
  64
  65 //=============================================================================
  66 /*!
  67  *
  68  */
  69 //=============================================================================
  70
  71 NETGENPlugin_NETGEN_1D2D3D_SA::~NETGENPlugin_NETGEN_1D2D3D_SA()
  72 {
  73 }
  74
  75 /**
  76  * @brief Check presence and content of orientation file. Implemented for completness and future reference.
  77  *
  78  * @param element_orientation_file Binary file containing the orientation of surface elemnts
  79  * @return true, false
  80  */
  81 bool NETGENPlugin_NETGEN_1D2D3D_SA::checkOrientationFile( const std::string element_orientation_file )
  82 {
  83   if(element_orientation_file.empty()){
  84       MESSAGE("No element orientation file");
  85       return true;
  86     } else {
  87       MESSAGE("Reading from elements from file: " << element_orientation_file);
  88       // By construction the orientation file written by Remote version has a zero written to mark no need of orientation in 2D meshing
  89       int nbElement;
  90       std::ifstream df(element_orientation_file, ios::binary|ios::in);
  91       df.read((char*)&nbElement, sizeof(int));
  92       df.close();
  93       return (nbElement == 0);
  94     }
  95 }
  96
  97 /**
  98  * @brief fill plugin hypothesis from the netgen_params structure
  99  *
 100  * @param aParams the structure
 101  */
 102 void NETGENPlugin_NETGEN_1D2D3D_SA::fillHyp(netgen_params aParams)
 103 {
 104  if( aParams.myType == hypoType::Hypo )
 105   {
 106     NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, GetGen());
 107
 108     hypParameters->SetMaxSize(aParams.maxh);
 109     hypParameters->SetMinSize(aParams.minh);
 110     hypParameters->SetNbSegPerEdge(aParams.segmentsperedge);
 111     hypParameters->SetGrowthRate(aParams.grading);
 112     hypParameters->SetNbSegPerRadius(aParams.curvaturesafety);
 113     hypParameters->SetSecondOrder(aParams.secondorder);
 114     hypParameters->SetQuadAllowed(aParams.quad);
 115     hypParameters->SetOptimize(aParams.optimize);
 116     hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness);
 117     hypParameters->SetSurfaceCurvature(aParams.uselocalh);
 118     hypParameters->SetFuseEdges(aParams.merge_solids);
 119     hypParameters->SetChordalErrorEnabled(aParams.chordalError);
 120     if(aParams.optimize){
 121       hypParameters->SetNbSurfOptSteps(aParams.optsteps2d);
 122       hypParameters->SetNbVolOptSteps(aParams.optsteps3d);
 123     }
 124     hypParameters->SetElemSizeWeight(aParams.elsizeweight);
 125     hypParameters->SetWorstElemMeasure(aParams.opterrpow);
 126     hypParameters->SetUseDelauney(aParams.delaunay);
 127     hypParameters->SetCheckOverlapping(aParams.checkoverlap);
 128     hypParameters->SetCheckChartBoundary(aParams.checkchartboundary);
 129     hypParameters->SetMeshSizeFile(aParams.meshsizefilename);
 130
 131     _hypothesis = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters);
 132   }
 133   else if ( aParams.myType == hypoType::Simple2D )
 134   {
 135     NETGENPlugin_SimpleHypothesis_2D * hypParameters = new NETGENPlugin_SimpleHypothesis_2D(0, GetGen());
 136
 137     // mandatory to fill in this branch case!
 138     // Number of segments   (int)
 139     // localLenght          (double)
 140     // maxElement area      (double)
 141     // GetAllowQuadrangles  (bool)
 142     hypParameters->SetNumberOfSegments( aParams.numberOfSegments );
 143     if ( !aParams.numberOfSegments )
 144       hypParameters->SetLocalLength( aParams.localLength );
 145     hypParameters->SetMaxElementArea( aParams.maxElementArea );
 146     hypParameters->SetAllowQuadrangles( aParams.allowQuadrangles );
 147     _hypothesis = dynamic_cast< const NETGENPlugin_SimpleHypothesis_2D *> (hypParameters);
 148   }
 149   else if ( aParams.myType == hypoType::Simple3D )
 150   {
 151     NETGENPlugin_SimpleHypothesis_3D * hypParameters = new NETGENPlugin_SimpleHypothesis_3D(0, GetGen());
 152
 153     // mandatory to fill in this branch case!
 154     // Number of segments   (int)
 155     // localLenght          (double)
 156     // maxElement area      (double)
 157     // maxElement volume    (double)
 158     // GetAllowQuadrangles  (bool)
 159     hypParameters->SetNumberOfSegments( aParams.numberOfSegments );
 160     if ( !aParams.numberOfSegments )
 161       hypParameters->SetLocalLength( aParams.localLength );
 162     hypParameters->SetMaxElementArea( aParams.maxElementArea );
 163     hypParameters->SetMaxElementVolume( aParams.maxElementVol );
 164     hypParameters->SetAllowQuadrangles( aParams.allowQuadrangles );
 165
 166     _hypothesis = dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D *> (hypParameters);
 167   }
 168 }
 169
 170 /**
 171  * @brief Write a binary file containing information on the elements/nodes
 172  *        created by the netgen mesher
 173  *
 174  * @param nodeVec mapping between the mesh id and the netgen structure id
 175  * @param ngLib Wrapper on netgen library
 176  * @param new_element_file Name of the output file
 177  * @param NumOfPremeshedNodes Number of nodes in the netgen structure
 178  * @return true if there are some error
 179  */
 180 bool NETGENPlugin_NETGEN_1D2D3D_SA::FillNewElementFile( std::vector< const SMDS_MeshNode* > &nodeVec, NETGENPlugin_NetgenLibWrapper &ngLib,
 181                                                         std::string new_element_file, const NETGENPlugin_Mesher::DIM dim )
 182 {
 183
 184   // Particularities: As submeshing is not supported nodeVect is empty and NumberOfPremeshedNodes is also zero
 185   Ng_Mesh* NetgenMesh  = ngLib.ngMesh();
 186   int NetgenNodes       = Ng_GetNP(NetgenMesh);
 187   int NetgenSeg2D       = Ng_GetNSeg_2D( NetgenMesh );
 188   int NetgenFaces       = Ng_GetNSE(NetgenMesh);
 189   int NetgenVols        = Ng_GetNE(NetgenMesh);
 190   int segmentId;
 191   bool isOK = ( NetgenNodes > 0 && NetgenSeg2D > 0 );
 192   if ( isOK && !new_element_file.empty() )
 193   {
 194     MESSAGE("Writting new elements")
 195     std::ofstream df(new_element_file, ios::out|ios::binary);
 196
 197     double NetgenPoint[3];
 198     int    NetgenSegment[2];
 199     int    NetgenSurface[8];
 200     int    NetgenVolumens[10];
 201
 202     // Writing nodevec (correspondence netgen numbering mesh numbering)
 203     // Number of nodes
 204     const int NumOfPremeshedNodes = nodeVec.size();
 205     df.write((char*) &NumOfPremeshedNodes, sizeof(int));
 206     df.write((char*) &NetgenNodes, sizeof(int));
 207
 208     for (int nodeIndex = 1 ; nodeIndex <= NumOfPremeshedNodes; ++nodeIndex )
 209     {
 210       //Id of the point
 211       int id = nodeVec.at(nodeIndex)->GetID();
 212       df.write((char*) &id, sizeof(int));
 213     }
 214
 215     // Writing all new points
 216     for (int nodeIndex = NumOfPremeshedNodes + 1; nodeIndex <= NetgenNodes; ++nodeIndex )
 217     {
 218       Ng_GetPoint( NetgenMesh, nodeIndex, NetgenPoint );
 219       // Coordinates of the point
 220       df.write((char *) &NetgenPoint, sizeof(double)*3);
 221     }
 222
 223     if ( dim >= NETGENPlugin_Mesher::D1 )
 224     {
 225       // create segments at boundaries.
 226       df.write((char*) &NetgenSeg2D, sizeof(int));
 227       for ( int elemIndex = 1; elemIndex <= NetgenSeg2D; ++elemIndex )
 228       {
 229         Ng_GetSegment_2D( NetgenMesh, elemIndex, NetgenSegment, &segmentId );
 230         df.write((char*) &NetgenSegment, sizeof(int) * 2 );
 231       }
 232     }
 233     if ( dim >= NETGENPlugin_Mesher::D2 )
 234     {
 235       // create surface elements.
 236       df.write((char*) &NetgenFaces, sizeof(int));
 237       for ( int elemIndex = 1; elemIndex <= NetgenFaces; ++elemIndex )
 238       {
 239         nglib::Ng_Surface_Element_Type elemType = Ng_GetSurfaceElement( NetgenMesh, elemIndex, NetgenSurface );
 240         switch (elemType)
 241         {
 242           case nglib::NG_TRIG:
 243           {
 244             df.write((char*) &NetgenSurface, sizeof(int) * 3 );
 245             break;
 246           }
 247           case nglib::NG_QUAD:
 248           {
 249             df.write((char*) &NetgenSurface, sizeof(int) * 4 );
 250             break;
 251           }
 252           case nglib::NG_TRIG6:
 253           {
 254             df.write((char*) &NetgenSurface, sizeof(int) * 6 );
 255             break;
 256           }
 257           case nglib::NG_QUAD8:
 258           {
 259             df.write((char*) &NetgenSurface, sizeof(int) * 8 );
 260             break;
 261           }
 262           default:
 263           { break; }
 264         }
 265       }
 266     }
 267     if ( dim >= NETGENPlugin_Mesher::D3 )
 268     {
 269       // create volume elements.
 270       df.write((char*) &NetgenVols, sizeof(int));
 271       for ( int elemIndex = 1; elemIndex <= NetgenVols; ++elemIndex )
 272       {
 273         nglib::Ng_Volume_Element_Type elemType =  Ng_GetVolumeElement( NetgenMesh, elemIndex, NetgenVolumens );
 274         switch (elemType)
 275         {
 276           case nglib::NG_TET:
 277           {
 278             df.write((char*) &NetgenVolumens, sizeof(int) * 4 );
 279             break;
 280           }
 281           case nglib::NG_PYRAMID:
 282           {
 283             df.write((char*) &NetgenVolumens, sizeof(int) * 5 );
 284             break;
 285           }
 286           case nglib::NG_PRISM:
 287           {
 288             df.write((char*) &NetgenVolumens, sizeof(int) * 6 );
 289             break;
 290           }
 291           case nglib::NG_TET10:
 292           {
 293             df.write((char*) &NetgenVolumens, sizeof(int) * 10 );
 294             break;
 295           }
 296           default:
 297           { break; }
 298         }
 299       }
 300     }
 301     df.close();
 302   }
 303   return false;
 304 }
 305
 306 /**
 307  * @brief Compute the mesh based on the
 308  *
 309  * @param aMesh the read Mesh
 310  * @param aShape the loaded shape
 311  * @param new_element_file output file containing info the elements created by the mesher
 312  * @param output_mesh whether or not write the created elements into the mesh
 313  * @param dim the dimension to be meshed.
 314  * @return negation of mesh fail: true, false
 315  */
 316 bool NETGENPlugin_NETGEN_1D2D3D_SA::Compute(SMESH_Mesh& aMesh, TopoDS_Shape &aShape, std::string new_element_file, bool output_mesh, NETGENPlugin_Mesher::DIM dim )
 317 {
 318   //
 319   netgen::multithread.terminate = 0;
 320   NETGENPlugin_Mesher mesher(&aMesh, aShape, /*is3D = */ false );
 321   mesher.SetParameters(dynamic_cast<const NETGENPlugin_Hypothesis*>(_hypothesis));
 322   if ( dim == NETGENPlugin_Mesher::D3 )
 323     mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_3D*>(_hypothesis));
 324   else
 325     mesher.SetParameters(dynamic_cast<const NETGENPlugin_SimpleHypothesis_2D*>(_hypothesis));
 326   NETGENPlugin_NetgenLibWrapper ngLib;
 327   vector< const SMDS_MeshNode* > nodeVec;
 328   bool err = mesher.Compute( ngLib, nodeVec, output_mesh, dim );
 329   FillNewElementFile( nodeVec, ngLib, new_element_file, dim );
 330   return err;
 331 }
 332
 333 /**
 334  * @brief Running the mesher on the given files
 335  *
 336  * @param input_mesh_file Mesh file (containing 2D elements)
 337  * @param shape_file Shape file (BREP or STEP format)
 338  * @param hypo_file Ascii file containing the netgen parameters
 339  * @param element_orientation_file Binary file containing the orientation of surface elemnts
 340  * @param new_element_file output file containing info the elements created by the mesher
 341  * @param output_mesh_file output mesh file (if empty it will not be created)
 342  * @return negation of mesh fail: true, false
 343  */
 344 int NETGENPlugin_NETGEN_1D2D3D_SA::run(const std::string input_mesh_file,
 345                                       const std::string shape_file,
 346                                       const std::string hypo_file,
 347                                       const std::string element_orientation_file,
 348                                       const std::string new_element_file,
 349                                       const std::string output_mesh_file,
 350                                       const NETGENPlugin_Mesher::DIM dim )
 351   {
 352
 353   // _element_orientation_file = element_orientation_file;
 354
 355   std::unique_ptr<SMESH_Mesh> myMesh(_gen->CreateMesh(false));
 356
 357   // Importing mesh
 358   SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh);
 359   // Importing shape
 360   TopoDS_Shape myShape;
 361   SMESH_DriverShape::importShape(shape_file, myShape);
 362   // Importing hypothesis
 363   netgen_params myParams;
 364   importNetgenParams(hypo_file, myParams);
 365   fillHyp(myParams);
 366   int ret = 1;
 367
 368   if ( checkOrientationFile(element_orientation_file) )
 369   {
 370     ret = Compute( *myMesh, myShape, new_element_file, !output_mesh_file.empty(), dim );
 371     if(ret){
 372       std::cerr << "Meshing failed" << std::endl;
 373       return ret;
 374     }
 375
 376     if(!output_mesh_file.empty()){
 377       std::string meshName = "MESH";
 378       SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName);
 379     }
 380   }
 381   else
 382     std::cerr << "For NETGENPlugin_NETGEN_1D2D3D_SA, orientation file should be market with 0 or be empty!" << std::endl;
 383
 384
 385   return ret;
 386 }