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 0de7ea3ca0
commit 3c4eb540f8
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 new_element_file=tmp_folder / fs::path("new_elements.dat");
fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med");
// TODO: Remove that file we do not use it
fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med");
// Not used kept for debug
//fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med");
fs::path shape_file=tmp_folder / fs::path("shape.brep");
fs::path param_file=tmp_folder / fs::path("netgen3d_param.txt");
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
std::string mesh_name = "Maillage_1";
@ -364,7 +365,7 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
+ "NONE";
// Writing command in log
{
std::ofstream flog(log_file.string());
std::ofstream flog(cmd_file.string());
flog << cmd << 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);
std::cout << "Time for exec of run_mesher: " << elapsed.count() * 1e-9 << std::endl;
// TODO: better error handling (display log ?)
if(ret != 0){
// Run crahed
std::cout << "Issue with command: " << std::endl;
@ -416,8 +416,8 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
int nodeID;
SMESH_MesherHelper helper(aMesh);
// This function
int _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
// This function is mandatory for setElementsOnShape to work
helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true );
// 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.task = "Volume meshing";
_progressByTic = -1.;
aParams._progressByTic = -1.;
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
SMESH_MesherHelper helper(aMesh);
_quadraticMesh = helper.IsQuadraticSubMesh(aShape);
aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true );
int Netgen_NbOfNodes = 0;
Netgen_NbOfNodes = 0;
double Netgen_point[3];
int Netgen_triangle[3];
NETGENPlugin_NetgenLibWrapper ngLib;
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;
@ -522,10 +521,10 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
if ( _viscousLayersHyp )
if ( aParams._viscousLayersHyp )
{
netgen::multithread.percent = 3;
proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
if ( !proxyMesh )
return false;
}
@ -553,7 +552,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
if ( !aSubMeshDSFace ) continue;
SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
if ( _quadraticMesh &&
if ( aParams._quadraticMesh &&
dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
{
// add medium nodes of proxy triangles to helper (#16843)
@ -566,10 +565,16 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
{
// 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");
if ( !elem ){
aParams._error = COMPERR_BAD_INPUT_MESH;
aParams._comment = "Null 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
@ -630,88 +635,231 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
internals);
}
}
// -------------------------
// 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 ));
Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
return false;
}
// namespace
// {
// void limitVolumeSize( netgen::Mesh* ngMesh,
// double maxh )
// {
// // get average h of faces
// double faceh = 0;
// 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 );
bool NETGENPlugin_NETGEN_3D::computePrepareParam(
SMESH_Mesh& aMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
netgen::OCCGeometry &occgeo,
SMESH_MesherHelper &helper,
netgen_params &aParams,
int &endWith)
// double compareh;
// if ( faceh < 0.5 * maxh ) compareh = -1;
// else if ( faceh > 1.5 * maxh ) compareh = 1;
// else compareh = 0;
// // cerr << "faceh " << faceh << endl;
// // cerr << "init maxh " << maxh << endl;
// // cerr << "compareh " << compareh << endl;
{
netgen::multithread.terminate = 0;
// if ( compareh > 0 )
// maxh *= 1.2;
// else
// maxh *= 0.8;
// // cerr << "maxh " << maxh << endl;
netgen::Mesh* ngMesh = ngLib._ngMesh;
// // get bnd box
// 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();
NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
// if ( ! & ngMesh->LocalHFunction() )
// ngMesh->SetLocalH( pmin, pmax, compareh <= 0 ? 0.1 : 0.5 );
// // adjusted by SALOME_TESTS/Grids/smesh/bugs_08/I8
// const int nbX = Max( 2, int( dx / maxh * 2 ));
// const int nbY = Max( 2, int( dy / maxh * 2 ));
// const int nbZ = Max( 2, int( dz / maxh * 2 ));
if ( aParams._hypParameters )
{
aMesher.SetParameters( aParams._hypParameters );
// netgen::Point3d p;
// for ( int i = 0; i <= nbX; ++i )
// {
// p.X() = pmin.X() + i * dx / nbX;
// for ( int j = 0; j <= nbY; ++j )
// {
// p.Y() = pmin.Y() + j * dy / nbY;
// for ( int k = 0; k <= nbZ; ++k )
// {
// p.Z() = pmin.Z() + k * dz / nbZ;
// ngMesh->RestrictLocalH( p, maxh );
// }
// }
// }
// }
// }
if ( !aParams._hypParameters->GetLocalSizesAndEntries().empty() ||
!aParams._hypParameters->GetMeshSizeFile().empty() )
{
if ( ! &ngMesh->LocalHFunction() )
{
netgen::Point3d pmin, pmax;
ngMesh->GetBox( pmin, pmax, 0 );
ngMesh->SetLocalH( pmin, pmax, aParams._hypParameters->GetGrowthRate() );
}
aMesher.SetLocalSize( occgeo, *ngMesh );
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 "Utils_SALOME_Exception.hxx"
#include <vector>
class StdMeshers_ViscousLayers;
class StdMeshers_MaxElementVolume;
class NETGENPlugin_Hypothesis;
class NETGENPlugin_NetgenLibWrapper;
class netgen_params;
class SMDS_MeshNode;
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,
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:
void exportElementOrientation(SMESH_Mesh& aMesh,
@ -81,6 +115,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_3D: public SMESH_3D_Algo
const TopoDS_Shape& aShape);
bool compute(SMESH_Mesh& mesh,
SMESH_MesherHelper& helper,
std::vector< const SMDS_MeshNode* >& nodeVec,