src/GMSHPlugin/GMSHPlugin_GMSH_3D_SA.cxx

   1 // Copyright (C) 2012-2015  ALNEOS
   2 // Copyright (C) 2016-2023  EDF
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License, or (at your option) any later version.
   8 //
   9 // This library is distributed in the hope that it will be useful,
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.alneos.com/ or email : contact@alneos.fr
  19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 //
  21 #include "GMSHPlugin_GMSH_3D.hxx"
  22 #include "GMSHPlugin_GMSH_3D_SA.hxx"
  23 #include "GMSHPlugin_Hypothesis_2D.hxx"
  24 #include "GMSHPlugin_Mesher.hxx"
  25
  26 #include <SMESH_Gen.hxx>
  27 #include <SMESH_ControlsDef.hxx>
  28 #include <SMESHDS_Mesh.hxx>
  29 #include <utilities.h>
  30
  31 #include <list>
  32
  33 //=============================================================================
  34 /*!
  35  *
  36  */
  37 //=============================================================================
  38
  39 GMSHPlugin_GMSH_3D_SA::GMSHPlugin_GMSH_3D_SA()
  40   : GMSHPlugin_GMSH_3D(0, new SMESH_Gen() )
  41 {
  42   MESSAGE("GMSHPlugin_GMSH_3D_SA::GMSHPlugin_GMSH_3D_SA");
  43   _name = "GMSH_3D_SA";
  44 }
  45
  46 //=============================================================================
  47 /*!
  48  *
  49  */
  50 //=============================================================================
  51
  52 GMSHPlugin_GMSH_3D_SA::~GMSHPlugin_GMSH_3D_SA()
  53 {
  54   MESSAGE("GMSHPlugin_GMSH_3D_SA::~GMSHPlugin_GMSH_3D_SA");
  55 }
  56
  57 void GMSHPlugin_GMSH_3D_SA::importGMSHParameters( const std::string hypo_file )
  58 {
  59   GMSHPlugin_Hypothesis * hypParameters = new GMSHPlugin_Hypothesis(0, GetGen());
  60   std::ifstream myfile(hypo_file);
  61   std::string line;
  62
  63   std::getline(myfile, line);
  64   int hasParams = std::stoi(line);
  65   if ( hasParams )
  66   {
  67     std::getline(myfile, line);
  68     hypParameters->Set2DAlgo( (GMSHPlugin_Hypothesis::Algo2D)(std::stoi(line)) );
  69     std::getline(myfile, line);
  70     hypParameters->Set3DAlgo( (GMSHPlugin_Hypothesis::Algo3D)(std::stoi(line)) );
  71     std::getline(myfile, line);
  72     hypParameters->SetRecomb2DAlgo( (GMSHPlugin_Hypothesis::Recomb2DAlgo)(std::stoi(line)) );
  73     std::getline(myfile, line);
  74     hypParameters->SetRecombineAll( std::stoi(line) );
  75     std::getline(myfile, line);
  76     hypParameters->SetSubdivAlgo( (GMSHPlugin_Hypothesis::SubdivAlgo)(std::stoi(line)) );
  77     std::getline(myfile, line);
  78     hypParameters->SetRemeshAlgo( (GMSHPlugin_Hypothesis::RemeshAlgo)(std::stoi(line)) );
  79     std::getline(myfile, line);
  80     hypParameters->SetSmouthSteps( std::stod(line) );
  81     std::getline(myfile, line);
  82     hypParameters->SetSizeFactor( std::stod(line) );
  83     std::getline(myfile, line);
  84 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
  85     hypParameters->SetMeshCurvatureSize( std::stod(line) );
  86 #endif
  87     std::getline(myfile, line);
  88     hypParameters->SetMaxSize( std::stod(line) );
  89     std::getline(myfile, line);
  90     hypParameters->SetMinSize( std::stod(line) );
  91     std::getline(myfile, line);
  92     hypParameters->SetSecondOrder( std::stoi(line) );
  93     std::getline(myfile, line);
  94     hypParameters->SetUseIncomplElem( std::stoi(line) );
  95     std::getline(myfile, line); // Get Compound names comma separated
  96     if ( line != "0" ) // Mark no presence of groups
  97     {
  98       std::stringstream compoundLine (line);
  99       std::string compound;
 100       while ( std::getline (compoundLine, compound, ',') )
 101         hypParameters->SetCompoundOnEntry( compound );
 102     }
 103     _hypothesis = dynamic_cast< const GMSHPlugin_Hypothesis *> (hypParameters);
 104   }
 105   else
 106   {
 107     // in case the parameters are not defined the _hypothesis should remain undefined
 108     // so default parameters are used to mesh!
 109     // _hypothesis = NULL;
 110   }
 111 }
 112
 113 /**
 114  * @brief Write a binary file containing information on the elements/nodes
 115  *        created by the mesher
 116  *
 117  * @param nodeVec mapping between the mesh id and the netgen structure id
 118  * @param mesher gmhsplugin mesher
 119  * @param new_element_file Name of the output file
 120  * @return true if there are some error
 121  */
 122 void GMSHPlugin_GMSH_3D_SA::fillNewElementFile( std::vector< const SMDS_MeshNode* > &nodeVec, GMSHPlugin_Mesher &mesher, std::string new_element_file )
 123 {
 124   GModel* gModel = mesher.GetGModel();
 125   const int numberOfNodes = nodeVec.size() - 1;
 126   const int numOfVolumens = gModel->getNumMeshElements( 3 );
 127   int numberOfTotalNodes  = numberOfNodes;
 128   std::map< MVertex *,int> nodeMap;
 129   bool isOK = ( numOfVolumens > 0 );
 130   if ( isOK && !new_element_file.empty() )
 131   {
 132     MESSAGE("Writting new elements");
 133
 134     std::ofstream df(new_element_file, ios::out|ios::binary);
 135     double points[3];
 136     int    volumens[4];
 137
 138     df.write((char*) &numberOfNodes, sizeof(int));
 139
 140     // To get 3D elements we need to iterate in regions
 141     std::map<int,MVertex*> vertexIdToCoordinate;
 142     // Index nodes of new volumetric elements in order from the vol region
 143     for ( GModel::riter it = gModel->firstRegion(); it != gModel->lastRegion(); ++it)
 144     {
 145       GRegion *gRegion = *it;
 146       for( size_t i = 0; i < gRegion->mesh_vertices.size(); i++)
 147       {
 148         MVertex *v = gRegion->mesh_vertices[i];
 149         const SMDS_MeshNode * preMeshedNode =  mesher.PremeshedNode( v );
 150         if ( !preMeshedNode )
 151         {
 152           numberOfTotalNodes++;
 153           vertexIdToCoordinate[ numberOfTotalNodes ] = v;
 154           nodeMap.insert({ v, numberOfTotalNodes });
 155         }
 156       }
 157     }
 158
 159     df.write((char*) &numberOfTotalNodes, sizeof(int));
 160
 161     for (int nodeIndex = 1 ; nodeIndex <= numberOfNodes; ++nodeIndex )
 162     {
 163       //Id of the point
 164       int id = nodeVec.at(nodeIndex)->GetID();
 165       df.write((char*) &id, sizeof(int));
 166     }
 167
 168     // Writing info on new points
 169     for (int nodeIndex = numberOfNodes+1; nodeIndex <= numberOfTotalNodes; ++nodeIndex )
 170     {
 171       if ( vertexIdToCoordinate[nodeIndex] )
 172       {
 173         df.write((char *) &vertexIdToCoordinate[nodeIndex]->x(), sizeof(double) );
 174         df.write((char *) &vertexIdToCoordinate[nodeIndex]->y(), sizeof(double) );
 175         df.write((char *) &vertexIdToCoordinate[nodeIndex]->z(), sizeof(double) );
 176       }
 177     }
 178
 179     df.write((char*) &numOfVolumens, sizeof(int));
 180     for ( GModel::riter it = gModel->firstRegion(); it != gModel->lastRegion(); ++it)
 181     {
 182       GRegion *gRegion = *it;
 183       std::vector<MVertex*> verts;
 184
 185       for( size_t i = 0; i < gRegion->getNumMeshElements(); i++)
 186       {
 187         MElement *element = gRegion->getMeshElement(i);
 188         verts.clear();
 189         element->getVertices(verts);
 190         for( MVertex* v : verts )
 191         {
 192           const SMDS_MeshNode * node = mesher.PremeshedNode( v );
 193           auto it = find(nodeVec.begin(), nodeVec.end(), node );
 194           int nodeId = node ? (it - nodeVec.begin()) : nodeMap[ v ];
 195           df.write((char*) &nodeId, sizeof(int) );
 196         }
 197       }
 198     }
 199     df.close();
 200   }
 201 }
 202
 203 /**
 204  * @brief Fill the list of elements in order and associate the read orientation (if any) to then.
 205  *        Index the elements and orientation to the ordered map so we are sure the write orientation
 206  *        match the element from the load mesh. IF no orientation file was written, then all the faces
 207  *        of the read mesh match the face orientation of his original associated face
 208  * @param aMesh the loaded mesh
 209  * @param element_orientation_file name of the file where to read the orientation
 210  * @param listElements the ordered map of elements with his orientation associated
 211  * @return the listElements filled.
 212  */
 213 void GMSHPlugin_GMSH_3D_SA::fillListOfElementsOriented( SMESH_Mesh& aMesh, std::string element_orientation_file,
 214                                                         std::map<const SMDS_MeshElement*, bool, TIDCompare>& listElements )
 215 {
 216
 217   // fill the elements found in the surface to the listOfElements so it can be used by the mesher!
 218   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 219   std::map<vtkIdType, bool> elemOrientation;
 220   {
 221     // Setting all element orientation to false if there no element orientation file
 222     if(element_orientation_file.empty())
 223     {
 224       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 225       while ( iteratorElem->more() ) // loop on elements on a geom face
 226       {
 227         // check mesh face
 228         const SMDS_MeshElement* elem = iteratorElem->next();
 229         listElements[elem] = false;
 230       }
 231     }
 232     else
 233     {
 234       std::ifstream df(element_orientation_file, ios::binary|ios::in);
 235       int nbElement;
 236       bool orient;
 237
 238       vtkIdType id;
 239       df.read((char*)&nbElement, sizeof(int));
 240
 241       for(int ielem=0;ielem<nbElement;++ielem){
 242         df.read((char*) &id, sizeof(vtkIdType));
 243         df.read((char*) &orient, sizeof(bool));
 244         elemOrientation[id] = orient;
 245       }
 246       df.close();
 247       //
 248       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
 249       while ( iteratorElem->more() ) // loop on elements on a geom face
 250       {
 251         // check mesh face
 252         const SMDS_MeshElement* elem = iteratorElem->next();
 253         // only index registered elements
 254         bool isIn = elemOrientation.count(elem->GetID())==1;
 255         if(!isIn) continue;
 256         listElements[elem] = elemOrientation[elem->GetID()];
 257       }
 258     }
 259   }
 260 }
 261
 262
 263 /**
 264  * @brief Compute the mesh by first filling premeshed faces to gmsh and then calling the mesher for the upper dimension
 265  * @param aShape the loaded shape
 266  * @param aMesh the read Mesh (contain 2D elements covering the entire geometry)
 267  * @param new_element_file output file containing info the elements created by the mesher
 268  * @param element_orientation_file Binary file containing the orientation of surface elemnts
 269  * @param output_mesh whether or not write the created elements into the mesh
 270  * @return negation of mesh fail: true, false
 271  * */
 272 bool GMSHPlugin_GMSH_3D_SA::Compute( TopoDS_Shape &aShape, SMESH_Mesh& aMesh, std::string new_element_file,
 273                                       std::string element_orientation_file, bool output_mesh )
 274 {
 275   GMSHPlugin_Mesher mesher(&aMesh, aShape,/*2d=*/false, true);
 276   std::vector< const SMDS_MeshNode* > nodeVec; // to save premeshed elements
 277   mesher.SetParameters(dynamic_cast<const GMSHPlugin_Hypothesis*>(_hypothesis));
 278
 279   std::map<const SMDS_MeshElement*, bool, TIDCompare> listElements;
 280   fillListOfElementsOriented( aMesh, element_orientation_file, listElements );
 281   bool Compute = mesher.Compute3D( nodeVec, listElements, output_mesh );
 282
 283   fillNewElementFile( nodeVec, mesher, new_element_file );
 284   mesher.finalizeGModel();
 285   return Compute;
 286 }
 287
 288 /**
 289  * @brief Running the mesher on the given files
 290  *
 291  * @param input_mesh_file Mesh file (containing 2D elements)
 292  * @param shape_file Shape file (BREP or STEP format)
 293  * @param hypo_file Ascii file containing the gmsh parameters
 294  * @param element_orientation_file Binary file containing the orientation of surface elements
 295  * @param new_element_file output file containing info the elements created by the mesher
 296  * @param output_mesh_file output mesh file (if empty it will not be created)
 297  * @return int
 298  */
 299 int GMSHPlugin_GMSH_3D_SA::run(const std::string input_mesh_file,
 300                                 const std::string shape_file,
 301                                 const std::string hypo_file,
 302                                 const std::string element_orientation_file,
 303                                 const std::string new_element_file,
 304                                 const std::string output_mesh_file)
 305 {
 306   std::unique_ptr<SMESH_Mesh> myMesh(_gen->CreateMesh(false));
 307
 308   SMESH_DriverMesh::importMesh(input_mesh_file, *myMesh);
 309
 310   // Importing shape
 311   TopoDS_Shape myShape;
 312   SMESH_DriverShape::importShape(shape_file, myShape);
 313
 314   // Define _hypothesis to then be able to call SetParameters( hypothesis )
 315   importGMSHParameters(hypo_file);
 316
 317   MESSAGE("Meshing with gmsh3d");
 318   int ret = Compute(myShape, *myMesh, new_element_file, element_orientation_file, !output_mesh_file.empty());
 319
 320   if(ret){
 321     std::cerr << "Meshing failed" << std::endl;
 322     return ret;
 323   }
 324
 325   if(!output_mesh_file.empty()){
 326     std::string meshName = "MESH";
 327     SMESH_DriverMesh::exportMesh(output_mesh_file, *myMesh, meshName);
 328   }
 329
 330   return ret;
 331 }