ForkMaster/netgenplugin
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Refactoring of normal Compute

Browse Source
This commit is contained in:
Yoann Audouin 2022-09-19 11:18:08 +02:00
parent f48857dd37
commit 7b58079a8f
2 changed files with 283 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

348
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
View File

@ -315,11 +315,12 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat"); fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat");
fs::path new_element_file=tmp_folder / fs::path("new_elements.dat"); fs::path new_element_file=tmp_folder / fs::path("new_elements.dat");
fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med"); fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med");
// TODO: Remove that file we do not use it // Not used kept for debug
fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med"); //fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med");
fs::path shape_file=tmp_folder / fs::path("shape.brep"); fs::path shape_file=tmp_folder / fs::path("shape.brep");
fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt"); fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt");
fs::path log_file=tmp_folder / fs::path("run.log"); fs::path log_file=tmp_folder / fs::path("run.log");
fs::path cmd_file=tmp_folder / fs::path("cmd.log");
//TODO: Handle variable mesh_name //TODO: Handle variable mesh_name
std::string mesh_name = "Maillage_1"; std::string mesh_name = "Maillage_1";
@ -364,7 +365,7 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
+ "NONE"; + "NONE";
// Writing command in log // Writing command in log
{ {
std::ofstream flog(log_file.string()); std::ofstream flog(cmd_file.string());
flog << cmd << endl; flog << cmd << endl;
flog << endl; flog << endl;
} }
@ -395,7 +396,6 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4); elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time5-time4);
std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl; std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl;
// TODO: better error handling (display log ?)
if(ret != 0){ if(ret != 0){
// Run crahed // Run crahed
std::cout << "Issue with command: " << std::endl; std::cout << "Issue with command: " << std::endl;
@ -416,8 +416,8 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
int nodeID; int nodeID;
SMESH_MesherHelper helper(aMesh); SMESH_MesherHelper helper(aMesh);
// This function // This function is mandatory for setElementsOnShape to work
int _quadraticMesh = helper.IsQuadraticSubMesh(aShape); helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true ); helper.SetElementsOnShape( true );
// Number of nodes in intial mesh // Number of nodes in intial mesh
@ -474,32 +474,31 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
*/ */
//============================================================================= //=============================================================================
bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape)
{
if(aMesh.IsParallel())
return RemoteCompute(aMesh, aShape);
auto time0 = std::chrono::high_resolution_clock::now();
bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_MesherHelper &helper,
netgen_params &aParams,
int &Netgen_NbOfNodes)
{
netgen::multithread.terminate = 0; netgen::multithread.terminate = 0;
netgen::multithread.task = "Volume meshing"; netgen::multithread.task = "Volume meshing";
_progressByTic = -1.; aParams._progressByTic = -1.;
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
SMESH_MesherHelper helper(aMesh); aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
_quadraticMesh = helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true ); helper.SetElementsOnShape( true );
int Netgen_NbOfNodes = 0; Netgen_NbOfNodes = 0;
double Netgen_point[3]; double Netgen_point[3];
int Netgen_triangle[3]; int Netgen_triangle[3];
NETGENPlugin_NetgenLibWrapper ngLib;
Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh; Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
// vector of nodes in which node index == netgen ID
vector< const SMDS_MeshNode* > nodeVec;
{ {
const int invalid_ID = -1; const int invalid_ID = -1;
@ -522,10 +521,10 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
if ( _viscousLayersHyp ) if ( aParams._viscousLayersHyp )
{ {
netgen::multithread.percent = 3; netgen::multithread.percent = 3;
proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape ); proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
if ( !proxyMesh ) if ( !proxyMesh )
return false; return false;
} }
@ -553,7 +552,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
if ( !aSubMeshDSFace ) continue; if ( !aSubMeshDSFace ) continue;
SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
if ( _quadraticMesh && if ( aParams._quadraticMesh &&
dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace )) dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
{ {
// add medium nodes of proxy triangles to helper (#16843) // add medium nodes of proxy triangles to helper (#16843)
@ -566,10 +565,16 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
{ {
// check mesh face // check mesh face
const SMDS_MeshElement* elem = iteratorElem->next(); const SMDS_MeshElement* elem = iteratorElem->next();
if ( !elem ) if ( !elem ){
return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); aParams._error = COMPERR_BAD_INPUT_MESH;
if ( elem->NbCornerNodes() != 3 ) aParams._comment = "Null element encounters";
return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); return true;
}
if ( elem->NbCornerNodes() != 3 ){
aParams._error = COMPERR_BAD_INPUT_MESH;
aParams._comment = "Not triangle element encounters";
return true;
}
// Add nodes of triangles and triangles them-selves to netgen mesh // Add nodes of triangles and triangles them-selves to netgen mesh
@ -630,88 +635,231 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
internals); internals);
} }
} }
Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
// ------------------------- return false;
// Generate the volume mesh
// -------------------------
auto time1 = std::chrono::high_resolution_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(time1-time0);
std::cout << "Time for seq:fill_in_ngmesh: " << elapsed.count() * 1e-9 << std::endl;
return (ngLib._isComputeOk = compute( aMesh, helper, nodeVec, ngLib ));
} }
// namespace bool NETGENPlugin_NETGEN_3D::computePrepareParam(
// { SMESH_Mesh& aMesh,
// void limitVolumeSize( netgen::Mesh* ngMesh, NETGENPlugin_NetgenLibWrapper &ngLib,
// double maxh ) netgen::OCCGeometry &occgeo,
// { SMESH_MesherHelper &helper,
// // get average h of faces netgen_params &aParams,
// double faceh = 0; int &endWith)
// int nbh = 0;
// for (int i = 1; i <= ngMesh->GetNSE(); i++)
// {
// const netgen::Element2d& face = ngMesh->SurfaceElement(i);
// for (int j=1; j <= face.GetNP(); ++j)
// {
// const netgen::PointIndex & i1 = face.PNumMod(j);
// const netgen::PointIndex & i2 = face.PNumMod(j+1);
// if ( i1 < i2 )
// {
// const netgen::Point3d & p1 = ngMesh->Point( i1 );
// const netgen::Point3d & p2 = ngMesh->Point( i2 );
// faceh += netgen::Dist2( p1, p2 );
// nbh++;
// }
// }
// }
// faceh = Sqrt( faceh / nbh );
// double compareh; {
// if ( faceh < 0.5 * maxh ) compareh = -1; netgen::multithread.terminate = 0;
// else if ( faceh > 1.5 * maxh ) compareh = 1;
// else compareh = 0;
// // cerr << "faceh " << faceh << endl;
// // cerr << "init maxh " << maxh << endl;
// // cerr << "compareh " << compareh << endl;
// if ( compareh > 0 ) netgen::Mesh* ngMesh = ngLib._ngMesh;
// maxh *= 1.2;
// else
// maxh *= 0.8;
// // cerr << "maxh " << maxh << endl;
// // get bnd box NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
// netgen::Point3d pmin, pmax;
// ngMesh->GetBox( pmin, pmax, 0 );
// const double dx = pmax.X() - pmin.X();
// const double dy = pmax.Y() - pmin.Y();
// const double dz = pmax.Z() - pmin.Z();
// if ( ! & ngMesh->LocalHFunction() )
// ngMesh->SetLocalH( pmin, pmax, compareh <= 0 ? 0.1 : 0.5 );
// // adjusted by SALOME_TESTS/Grids/smesh/bugs_08/I8 if ( aParams._hypParameters )
// const int nbX = Max( 2, int( dx / maxh * 2 )); {
// const int nbY = Max( 2, int( dy / maxh * 2 )); aMesher.SetParameters( aParams._hypParameters );
// const int nbZ = Max( 2, int( dz / maxh * 2 ));
// netgen::Point3d p; if ( !aParams._hypParameters->GetLocalSizesAndEntries().empty() ||
// for ( int i = 0; i <= nbX; ++i ) !aParams._hypParameters->GetMeshSizeFile().empty() )
// { {
// p.X() = pmin.X() + i * dx / nbX; if ( ! &ngMesh->LocalHFunction() )
// for ( int j = 0; j <= nbY; ++j ) {
// { netgen::Point3d pmin, pmax;
// p.Y() = pmin.Y() + j * dy / nbY; ngMesh->GetBox( pmin, pmax, 0 );
// for ( int k = 0; k <= nbZ; ++k ) ngMesh->SetLocalH( pmin, pmax, aParams._hypParameters->GetGrowthRate() );
// { }
// p.Z() = pmin.Z() + k * dz / nbZ; aMesher.SetLocalSize( occgeo, *ngMesh );
// ngMesh->RestrictLocalH( p, maxh );
// } try {
// } ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
// } } catch (netgen::NgException & ex) {
// } aParams._error = COMPERR_BAD_PARMETERS;
// } aParams._comment = ex.What();
return false;
}
}
if ( !aParams._hypParameters->GetOptimize() )
endWith = netgen::MESHCONST_MESHVOLUME;
}
else if ( aParams._hypMaxElementVolume )
{
netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
// limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
}
else if ( aMesh.HasShapeToMesh() )
{
aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
}
else
{
netgen::Point3d pmin, pmax;
ngMesh->GetBox (pmin, pmax);
netgen::mparam.maxh = Dist(pmin, pmax)/2;
}
if ( !aParams._hypParameters && aMesh.HasShapeToMesh() )
{
netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
}
return false;
}
bool NETGENPlugin_NETGEN_3D::computeRunMesher(
netgen::OCCGeometry &occgeo,
vector< const SMDS_MeshNode* > &nodeVec,
netgen::Mesh* ngMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
netgen_params &aParams,
int &startWith, int &endWith)
{
int err = 1;
try
{
OCC_CATCH_SIGNALS;
ngLib.CalcLocalH(ngMesh);
err = ngLib.GenerateMesh(occgeo, startWith, endWith);
if(netgen::multithread.terminate)
return false;
if ( err ){
aParams._comment = SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task;
return true;
}
}
catch (Standard_Failure& ex)
{
SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
str << " at " << netgen::multithread.task
<< ": " << ex.DynamicType()->Name();
if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
str << ": " << ex.GetMessageString();
aParams._comment = str;
return true;
}
catch (netgen::NgException& exc)
{
SMESH_Comment str("NgException");
if ( strlen( netgen::multithread.task ) > 0 )
str << " at " << netgen::multithread.task;
str << ": " << exc.What();
aParams._comment = str;
return true;
}
catch (...)
{
SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
if ( strlen( netgen::multithread.task ) > 0 )
str << " at " << netgen::multithread.task;
aParams._comment = str;
return true;
}
if ( err )
{
SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
if ( ce && ce->HasBadElems() ){
aParams._error = ce->myName;
aParams._comment = ce->myComment;
return true;
}
}
return false;
}
bool NETGENPlugin_NETGEN_3D::computeFillMesh(
vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_MesherHelper &helper,
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 = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
if ( isOK )
{
double Netgen_point[3];
int Netgen_tetrahedron[4];
// create and insert new nodes into nodeVec
nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
int nodeIndex = Netgen_NbOfNodes + 1;
for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
{
Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
}
// create tetrahedrons
for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
{
Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
try
{
helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
nodeVec.at( Netgen_tetrahedron[1] ),
nodeVec.at( Netgen_tetrahedron[2] ),
nodeVec.at( Netgen_tetrahedron[3] ));
}
catch (...)
{
}
}
}
return false;
}
bool NETGENPlugin_NETGEN_3D::Compute(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape)
{
if(aMesh.IsParallel())
return RemoteCompute(aMesh, aShape);
// 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;
netgen_params aParams;
aParams._hypParameters = const_cast<NETGENPlugin_Hypothesis*>(_hypParameters);
aParams._hypMaxElementVolume = const_cast<StdMeshers_MaxElementVolume*>(_hypMaxElementVolume);
aParams.maxElementVolume = _maxElementVolume;
aParams._progressByTic = _progressByTic;
aParams._quadraticMesh = _quadraticMesh;
aParams._viscousLayersHyp = const_cast<StdMeshers_ViscousLayers*>(_viscousLayersHyp);
bool ret;
ret = computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, Netgen_NbOfNodes);
if(ret)
return error( aParams._error, aParams._comment);
netgen::OCCGeometry occgeo;
computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
ret = computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith);
if(ret){
if(aParams._error)
return error(aParams._error, aParams._comment);
error(aParams._comment);
return true;
}
computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
return false;
}
//================================================================================ //================================================================================
/*! /*!

35
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.hxx
View File

@ -37,11 +37,14 @@
#include "SMESH_Algo.hxx" #include "SMESH_Algo.hxx"
#include "Utils_SALOME_Exception.hxx" #include "Utils_SALOME_Exception.hxx"
#include <vector>
class StdMeshers_ViscousLayers; class StdMeshers_ViscousLayers;
class StdMeshers_MaxElementVolume; class StdMeshers_MaxElementVolume;
class NETGENPlugin_Hypothesis; class NETGENPlugin_Hypothesis;
class NETGENPlugin_NetgenLibWrapper; class NETGENPlugin_NetgenLibWrapper;
class netgen_params; class netgen_params;
class SMDS_MeshNode;
class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo
{ {
@ -67,6 +70,37 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo
const TopoDS_Shape& aShape, const TopoDS_Shape& aShape,
MapShapeNbElems& aResMap); MapShapeNbElems& aResMap);
static bool computeFillNgMesh(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_MesherHelper &helper,
netgen_params &aParams,
int &Netgen_NbOfNodes);
static bool computePrepareParam(
SMESH_Mesh& aMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
netgen::OCCGeometry &occgeo,
SMESH_MesherHelper &helper,
netgen_params &aParams,
int &endWith);
static bool computeRunMesher(
netgen::OCCGeometry &occgeo,
std::vector< const SMDS_MeshNode* > &nodeVec,
netgen::Mesh* ngMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
netgen_params &aParams,
int &startWith, int &endWith);
static bool computeFillMesh(
std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_MesherHelper &helper,
int &Netgen_NbOfNodes);
protected: protected:
void exportElementOrientation(SMESH_Mesh& aMesh, void exportElementOrientation(SMESH_Mesh& aMesh,
@ -81,6 +115,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo
const TopoDS_Shape& aShape); const TopoDS_Shape& aShape);
bool compute(SMESH_Mesh& mesh, bool compute(SMESH_Mesh& mesh,
SMESH_MesherHelper& helper, SMESH_MesherHelper& helper,
std::vector< const SMDS_MeshNode* >& nodeVec, std::vector< const SMDS_MeshNode* >& nodeVec,