netgenplugin/src/NETGENPlugin/NETGENPlugin_NETGEN_3D_SA.cxx

// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//

//=============================================================================
// File      : NETGENPlugin_NETGEN_3D_SA.cxx
// Created   : lundi 19 Septembre 2022
// Author    : Yoann AUDOUIN (CEA)
// Project   : SALOME
//=============================================================================
//
//
#include "NETGENPlugin_NETGEN_3D_SA.hxx"

#include "NETGENPlugin_DriverParam.hxx"
#include "NETGENPlugin_Hypothesis.hxx"
#include "StdMeshers_MaxElementVolume.hxx"

#include <SMESH_Gen.hxx>
#include <SMESH_Mesh.hxx>
#include <SMESH_MesherHelper.hxx>
#include <SMESH_DriverShape.hxx>
#include <SMESH_DriverMesh.hxx>
#include <SMESHDS_Mesh.hxx>


#include <boost/filesystem.hpp>
namespace fs = boost::filesystem;

/*
  Netgen include files
*/

#ifndef OCCGEOMETRY
#define OCCGEOMETRY
#endif
#include <occgeom.hpp>

#ifdef NETGEN_V5
#include <ngexception.hpp>
#endif
#ifdef NETGEN_V6
#include <core/exception.hpp>
#endif

namespace nglib {
#include <nglib.h>
}
namespace netgen {

  NETGENPLUGIN_DLL_HEADER
  extern MeshingParameters mparam;

  NETGENPLUGIN_DLL_HEADER
  extern volatile multithreadt multithread;
}
using namespace nglib;

//=============================================================================
/*!
 * Constructor
 */
//=============================================================================

NETGENPlugin_NETGEN_3D_SA::NETGENPlugin_NETGEN_3D_SA()
  :  NETGENPlugin_NETGEN_3D(0, new SMESH_Gen())
{
  _name = "NETGEN_3D_SA";
}

//=============================================================================
/*!
 * Destructor
 */
//=============================================================================

NETGENPlugin_NETGEN_3D_SA::~NETGENPlugin_NETGEN_3D_SA()
{
  if(_gen)
    delete _gen;
}

/**
 * @brief fill plugin hypothesis from the netgen_params structure
 *
 * @param aParams the structure
 * @param gen SMESH_Gen associate with the SA
 */
void NETGENPlugin_NETGEN_3D_SA::fillHyp(netgen_params aParams)
{
  if(_gen)
    std::cout << "_gen is set"  << std::endl;
  if(aParams.has_netgen_param){
    NETGENPlugin_Hypothesis * hypParameters = new NETGENPlugin_Hypothesis(0, GetGen());

    hypParameters->SetMaxSize(aParams.maxh);
    hypParameters->SetMinSize(aParams.minh);
    hypParameters->SetNbSegPerEdge(aParams.segmentsperedge);
    hypParameters->SetGrowthRate(aParams.grading);
    hypParameters->SetNbSegPerRadius(aParams.curvaturesafety);
    hypParameters->SetSecondOrder(aParams.secondorder);
    hypParameters->SetQuadAllowed(aParams.quad);
    hypParameters->SetOptimize(aParams.optimize);
    hypParameters->SetFineness((NETGENPlugin_Hypothesis::Fineness)aParams.fineness);
    hypParameters->SetSurfaceCurvature(aParams.uselocalh);
    hypParameters->SetFuseEdges(aParams.merge_solids);
    hypParameters->SetChordalErrorEnabled(aParams.chordalError);
    if(aParams.optimize){
      hypParameters->SetNbSurfOptSteps(aParams.optsteps2d);
      hypParameters->SetNbVolOptSteps(aParams.optsteps3d);
    }
    hypParameters->SetElemSizeWeight(aParams.elsizeweight);
    hypParameters->SetWorstElemMeasure(aParams.opterrpow);
    hypParameters->SetUseDelauney(aParams.delaunay);
    hypParameters->SetCheckOverlapping(aParams.checkoverlap);
    hypParameters->SetCheckChartBoundary(aParams.checkchartboundary);
    hypParameters->SetMeshSizeFile(aParams.meshsizefilename);

    _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis *> (hypParameters);
  }
  if(aParams.has_maxelementvolume_hyp){
    _hypMaxElementVolume = new StdMeshers_MaxElementVolume(1, GetGen());
    _maxElementVolume = aParams.maxElementVolume;
  }
  // TODO: Handle viscous layer
}

/**
 * @brief Write a binary file containing information on the elements/nodes
 *        created by the mesher
 *
 * @param nodeVec mapping between the mesh id and the netgen structure id
 * @param ngLib Wrapper on netgen library
 * @param new_element_file Name of the output file
 * @param Netgen_NbOfNodes Number of nodes in the netgen structure
 * @return true if there are some error
 */
bool NETGENPlugin_NETGEN_3D_SA::computeFillNewElementFile(
    std::vector< const SMDS_MeshNode* > &nodeVec,
    NETGENPlugin_NetgenLibWrapper &ngLib,
    std::string new_element_file,
    int &Netgen_NbOfNodes
)
{
  Ng_Mesh* Netgen_mesh = ngLib.ngMesh();

  int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
  int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);

  bool isOK = ( Netgen_NbOfTetra > 0 );
  if ( isOK && !new_element_file.empty() )
  {
    std::ofstream df(new_element_file, ios::out|ios::binary);

    double Netgen_point[3];
    int    Netgen_tetrahedron[4];

    // Writing nodevec (correspondence netgen numbering mesh numbering)
    // Number of nodes
    df.write((char*) &Netgen_NbOfNodes, sizeof(int));
    df.write((char*) &Netgen_NbOfNodesNew, sizeof(int));
    for (int nodeIndex = 1 ; nodeIndex <= Netgen_NbOfNodes; ++nodeIndex )
    {
      //Id of the point
      int id = nodeVec.at(nodeIndex)->GetID();
      df.write((char*) &id, sizeof(int));
    }

    // Writing info on new points
    for (int nodeIndex = Netgen_NbOfNodes +1 ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
    {
      Ng_GetPoint(Netgen_mesh, nodeIndex, Netgen_point );
      // Coordinates of the point
      df.write((char *) &Netgen_point, sizeof(double)*3);
    }

    // create tetrahedrons
    df.write((char*) &Netgen_NbOfTetra, sizeof(int));
    for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
    {
      Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
      df.write((char*) &Netgen_tetrahedron, sizeof(int)*4);
    }
  }
  return false;
}

/**
 * @brief Compute mesh associated to shape
 *
 * @param aShape the shape
 * @param aMesh the mesh
 * @param aParams netgen_params structure
 * @param new_element_file Name of the file containing new element
 * @param output_mesh Name of the output mesh (if empty it will not be written)
 * @return true if there are some error
 */
bool NETGENPlugin_NETGEN_3D_SA::Compute(
        TopoDS_Shape &aShape,
        SMESH_Mesh& aMesh,
        netgen_params& aParams,
        std::string new_element_file,
        bool output_mesh)
{
  // vector of nodes in which node index == netgen ID
  vector< const SMDS_MeshNode* > nodeVec;
  NETGENPlugin_NetgenLibWrapper ngLib;
  SMESH_MesherHelper helper(aMesh);
  int startWith = netgen::MESHCONST_MESHVOLUME;
  int endWith   = netgen::MESHCONST_OPTVOLUME;
  int Netgen_NbOfNodes=0;

  // Changing netgen log_file putting it next to new_element_file
  fs::path netgen_log_file = fs::path(new_element_file).remove_filename() / fs::path("NETGEN.out");

  ngLib.setOutputFile(netgen_log_file.string());

  NETGENPlugin_NETGEN_3D::computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, Netgen_NbOfNodes);

  netgen::OCCGeometry occgeo;
  NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, endWith);

  NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);

  computeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);

  if(output_mesh)
    NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);

  return false;
}


/**
 * @brief Running the mesher on the given files
 *
 * @param input_mesh_file Mesh file (containing 2D elements)
 * @param shape_file Shape file (BREP or STEP format)
 * @param hypo_file Ascii file containing the netgen parameters
 * @param element_orientation_file Binary file containing the orientation of surface elemnts
 * @param new_element_file output file containing info the elements created by the mesher
 * @param output_mesh_file output mesh file (if empty it will not be created)
 * @return int
 */
int NETGENPlugin_NETGEN_3D_SA::run(const std::string input_mesh_file,
          const std::string shape_file,
          const std::string hypo_file,
          const std::string element_orientation_file,
          const std::string new_element_file,
          const std::string output_mesh_file)
{

  _element_orientation_file = element_orientation_file;

  std::unique_ptr<SMESH_Mesh> myMesh(_gen->CreateMesh(false));

  importMesh(input_mesh_file, *myMesh);

  // Importing shape
  TopoDS_Shape myShape;
  importShape(shape_file, myShape);

  // Importing hypothesis
  netgen_params myParams;

  importNetgenParams(hypo_file, myParams);
  fillHyp(myParams);

  MESSAGE("Meshing with netgen3d");
  int ret = Compute(myShape, *myMesh, myParams,
                      new_element_file,
                      !output_mesh_file.empty());


  if(ret){
    std::cerr << "Meshing failed" << std::endl;
    return ret;
  }

  if(!output_mesh_file.empty()){
    std::string meshName = "MESH";
    exportMesh(output_mesh_file, *myMesh, meshName);
  }

  return ret;
}

/**
 * @brief Compute the list of already meshed Surface elements and info
 *        on their orientation and if they are internal
 *
 * @param aMesh Global Mesh
 * @param aShape Shape associated to the mesh
 * @param proxyMesh pointer to mesh used fo find the elements
 * @param internals information on internal sub shapes
 * @param helper helper associated to the mesh
 * @param listElements map of surface element associated with
 *                     their orientation and internal status
 * @return true if their was some error
 */
bool NETGENPlugin_NETGEN_3D_SA::getSurfaceElements(
    SMESH_Mesh&         aMesh,
    const TopoDS_Shape& aShape,
    SMESH_ProxyMesh::Ptr proxyMesh,
    NETGENPlugin_Internals &internals,
    SMESH_MesherHelper &helper,
    std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
    )
{
  // To remove compilation warnings
  (void) aShape;
  (void) proxyMesh;
  (void) internals;
  (void) helper;
  SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();

  // Get list of elements + their orientation from element_orientation file
  std::map<vtkIdType, bool> elemOrientation;
  {
    // Setting all element orientation to false if there no element orientation file
    if(_element_orientation_file.empty()){
      MESSAGE("No element orientation file");

      SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
      while ( iteratorElem->more() ) // loop on elements on a geom face
        {
          // check mesh face
          const SMDS_MeshElement* elem = iteratorElem->next();
          elemOrientation[elem->GetID()] = false;
        }
    } else {
      MESSAGE("Reading from elements from file: " << _element_orientation_file);
      std::ifstream df(_element_orientation_file, ios::binary|ios::in);
      int nbElement;
      bool orient;

      // Warning of the use of vtkIdType (I had issue when run_mesher was compiled with internal vtk) and salome not
      // Sizeof was the same but how he othered the type was different
      // Maybe using another type (uint64_t) instead would be better
      vtkIdType id;
      df.read((char*)&nbElement, sizeof(int));

      for(int ielem=0;ielem<nbElement;++ielem){
        df.read((char*) &id, sizeof(vtkIdType));
        df.read((char*) &orient, sizeof(bool));
        elemOrientation[id] = orient;
      }
    }
  }

  // Adding elements from Mesh
  SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
  bool isRev;

  bool isIn;

  while ( iteratorElem->more() ) // loop on elements on a geom face
  {
    // check mesh face
    const SMDS_MeshElement* elem = iteratorElem->next();
    if ( !elem ){
      return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
    }
    if ( elem->NbCornerNodes() != 3 ){
      return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
    }
    // Keeping only element that are in the element orientation file
    isIn = elemOrientation.count(elem->GetID())==1;
    if(!isIn)
      continue;
    // Get orientation
    // Netgen requires that all the triangle point outside
    isRev = elemOrientation[elem->GetID()];
    listElements[elem] = tuple(isRev, false);
  }

  return false;
}