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

Refactoring of runner

Browse Source
This commit is contained in:
Yoann Audouin 2022-09-19 15:13:31 +02:00
parent 1b4859d135
commit d9a37018c5
2 changed files with 248 additions and 371 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

117
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
View File

@ -474,6 +474,70 @@ int NETGENPlugin_NETGEN_3D::RemoteCompute(SMESH_Mesh& aMesh,
*/
//=============================================================================
/**
* @brief Get an iterator on the Surface element with their orientation
*
*/
bool getSurfaceElements(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
SMESH_ProxyMesh::Ptr proxyMesh,
NETGENPlugin_Internals &internals,
SMESH_MesherHelper &helper,
netgen_params &aParams,
std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
)
{
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
{
const TopoDS_Shape& aShapeFace = exFa.Current();
int faceID = meshDS->ShapeToIndex( aShapeFace );
bool isInternalFace = internals.isInternalShape( faceID );
bool isRev = false;
if ( checkReverse && !isInternalFace &&
helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
// IsReversedSubMesh() can work wrong on strongly curved faces,
// so we use it as less as possible
isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
if ( !aSubMeshDSFace ) continue;
SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
if ( aParams._quadraticMesh &&
dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
{
// add medium nodes of proxy triangles to helper (#16843)
while ( iteratorElem->more() )
helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
iteratorElem = aSubMeshDSFace->GetElements();
}
while(iteratorElem->more()){
const SMDS_MeshElement* elem = iteratorElem->next();
// check mesh face
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;
}
listElements[elem] = tuple(isRev, isInternalFace);
}
}
return false;
}
bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
SMESH_Mesh& aMesh,
@ -516,9 +580,8 @@ bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
// ---------------------------------
// Feed the Netgen with surface mesh
// ---------------------------------
TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
bool isRev=false;
bool isInternalFace=false;
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
if ( aParams._viscousLayersHyp )
@ -536,46 +599,15 @@ bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
proxyMesh.reset( Adaptor );
}
for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
std::map<const SMDS_MeshElement*, tuple<bool, bool>> listElements;
bool ret = getSurfaceElements(aMesh, aShape, proxyMesh, internals, helper, aParams, listElements);
if(ret)
return ret;
for ( auto const& [elem, info] : listElements ) // loop on elements on a geom face
{
const TopoDS_Shape& aShapeFace = exFa.Current();
int faceID = meshDS->ShapeToIndex( aShapeFace );
bool isInternalFace = internals.isInternalShape( faceID );
bool isRev = false;
if ( checkReverse && !isInternalFace &&
helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
// IsReversedSubMesh() can work wrong on strongly curved faces,
// so we use it as less as possible
isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
if ( !aSubMeshDSFace ) continue;
SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
if ( aParams._quadraticMesh &&
dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
{
// add medium nodes of proxy triangles to helper (#16843)
while ( iteratorElem->more() )
helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
iteratorElem = aSubMeshDSFace->GetElements();
}
while ( iteratorElem->more() ) // loop on elements on a geom face
{
// check mesh face
const SMDS_MeshElement* elem = iteratorElem->next();
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;
}
isRev = get<0>(info);
isInternalFace = get<1>(info);
// Add nodes of triangles and triangles them-selves to netgen mesh
// add three nodes of triangle
@ -617,7 +649,6 @@ bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
}
} // loop on elements on a face
} // loop on faces of a SOLID or SHELL
// insert old nodes into nodeVec
nodeVec.resize( nodeToNetgenID.size() + 1, 0 );

362
src/NETGENPlugin/NETGENPlugin_Runner.cxx
View File

@ -27,7 +27,7 @@
#include "NETGENPlugin_Runner.hxx"
#include "NETGENPlugin_NETGEN_3D.hxx"
#include "NETGENPlugin_DriverParam.hxx"
#include <fstream>
@ -228,66 +228,19 @@ int netgen3d(const std::string input_mesh_file,
return ret;
}
bool mycomputeFillNgMesh(
bool getSurfaceElements(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_ProxyMesh::Ptr proxyMesh,
NETGENPlugin_Internals &internals,
SMESH_MesherHelper &helper,
netgen_params &aParams,
std::string element_orientation_file,
int &Netgen_NbOfNodes)
std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
)
{
netgen::multithread.terminate = 0;
netgen::multithread.task = "Volume meshing";
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true );
double Netgen_point[3];
int Netgen_triangle[3];
Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
{
const int invalid_ID = -1;
SMESH::Controls::Area areaControl;
SMESH::Controls::TSequenceOfXYZ nodesCoords;
// maps nodes to ng ID
typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
typedef TNodeToIDMap::value_type TN2ID;
TNodeToIDMap nodeToNetgenID;
// find internal shapes
NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
// ---------------------------------
// Feed the Netgen with surface mesh
// ---------------------------------
TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
if ( aParams._viscousLayersHyp )
{
netgen::multithread.percent = 3;
proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
if ( !proxyMesh )
return false;
}
if ( aMesh.NbQuadrangles() > 0 )
{
netgen::multithread.percent = 6;
StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
Adaptor->Compute(aMesh,aShape,proxyMesh.get());
proxyMesh.reset( Adaptor );
}
// Get list of elements + their orientation from element_orientation file
std::map<vtkIdType, bool> elemOrientation;
{
@ -330,21 +283,99 @@ bool mycomputeFillNgMesh(
{
// 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;
}
// 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;
}
bool mycomputeFillNgMesh(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
SMESH_MesherHelper &helper,
netgen_params &aParams,
std::string element_orientation_file,
int &Netgen_NbOfNodes)
{
netgen::multithread.terminate = 0;
netgen::multithread.task = "Volume meshing";
aParams._progressByTic = -1.;
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true );
Netgen_NbOfNodes = 0;
double Netgen_point[3];
int Netgen_triangle[3];
Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
{
const int invalid_ID = -1;
SMESH::Controls::Area areaControl;
SMESH::Controls::TSequenceOfXYZ nodesCoords;
// maps nodes to ng ID
typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
typedef TNodeToIDMap::value_type TN2ID;
TNodeToIDMap nodeToNetgenID;
// find internal shapes
NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
// ---------------------------------
// Feed the Netgen with surface mesh
// ---------------------------------
bool isRev=false;
bool isInternalFace=false;
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
if ( aParams._viscousLayersHyp )
{
netgen::multithread.percent = 3;
proxyMesh = aParams._viscousLayersHyp->Compute( aMesh, aShape );
if ( !proxyMesh )
return false;
}
if ( aMesh.NbQuadrangles() > 0 )
{
netgen::multithread.percent = 6;
StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
Adaptor->Compute(aMesh,aShape,proxyMesh.get());
proxyMesh.reset( Adaptor );
}
std::map<const SMDS_MeshElement*, tuple<bool, bool>> listElements;
bool ret = getSurfaceElements(aMesh, aShape, proxyMesh, internals, helper, aParams, element_orientation_file, listElements);
if(ret)
return ret;
for ( auto const& [elem, info] : listElements ) // loop on elements on a geom face
{
isRev = get<0>(info);
isInternalFace = get<1>(info);
// Add nodes of triangles and triangles them-selves to netgen mesh
// add three nodes of triangle
@ -370,7 +401,6 @@ bool mycomputeFillNgMesh(
Netgen_point [ 2 ] = node->Z();
Ng_AddPoint(Netgen_mesh, Netgen_point);
}
Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
}
// add triangle
@ -381,13 +411,12 @@ bool mycomputeFillNgMesh(
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
// TODO: Handle that case (quadrangle 2D) (isInternal is set to false)
if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
{
swap( Netgen_triangle[1], Netgen_triangle[2] );
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
}
}
} // loop on elements on a face
// insert old nodes into nodeVec
nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
@ -396,143 +425,19 @@ bool mycomputeFillNgMesh(
nodeVec[ n_id->second ] = n_id->first;
nodeToNetgenID.clear();
// TODO: Handle internal vertex
//if ( internals.hasInternalVertexInSolid() )
//{
// netgen::OCCGeometry occgeo;
// NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
// (netgen::Mesh&) *Netgen_mesh,
// nodeVec,
// internals);
//}
if ( internals.hasInternalVertexInSolid() )
{
netgen::OCCGeometry occgeo;
NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
(netgen::Mesh&) *Netgen_mesh,
nodeVec,
internals);
}
}
Netgen_mesh = ngLib.ngMesh();
Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
return false;
}
bool mycomputePrepareParam(
SMESH_Mesh& aMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
netgen::OCCGeometry &occgeo,
SMESH_MesherHelper &helper,
netgen_params &aParams,
int &endWith)
{
// -------------------------
// Generate the volume mesh
// -------------------------
netgen::multithread.terminate = 0;
netgen::Mesh* ngMesh = ngLib._ngMesh;
NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
set_netgen_parameters(aParams);
if ( aParams.has_netgen_param )
{
if ( aParams.has_local_size)
{
if ( ! &ngMesh->LocalHFunction() )
{
netgen::Point3d pmin, pmax;
ngMesh->GetBox( pmin, pmax, 0 );
ngMesh->SetLocalH( pmin, pmax, aParams.grading );
}
aMesher.SetLocalSize( occgeo, *ngMesh );
try {
ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
} catch (netgen::NgException & ex) {
return error( COMPERR_BAD_PARMETERS, ex.What() );
}
}
if ( !aParams.optimize )
endWith = netgen::MESHCONST_MESHVOLUME;
}
else if ( aParams.has_maxelementvolume_hyp )
{
netgen::mparam.maxh = pow( 72, 1/6. ) * pow( aParams.maxElementVolume, 1/3. );
// limitVolumeSize( ngMesh, netgen::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.has_netgen_param && aMesh.HasShapeToMesh() )
{
netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
}
return false;
}
bool mycomputeRunMesher(
netgen::OCCGeometry &occgeo,
std::vector< const SMDS_MeshNode* > &nodeVec,
netgen::Mesh* ngMesh,
NETGENPlugin_NetgenLibWrapper &ngLib,
int &startWith, int &endWith)
{
int err = 1;
try
{
OCC_CATCH_SIGNALS;
ngLib.CalcLocalH(ngMesh);
err = ngLib.GenerateMesh(occgeo, startWith, endWith, ngMesh);
if(netgen::multithread.terminate)
return false;
if ( err )
return error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
}
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();
return error(str);
}
catch (netgen::NgException& exc)
{
SMESH_Comment str("NgException");
if ( strlen( netgen::multithread.task ) > 0 )
str << " at " << netgen::multithread.task;
str << ": " << exc.What();
return error(str);
}
catch (...)
{
SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
if ( strlen( netgen::multithread.task ) > 0 )
str << " at " << netgen::multithread.task;
return error(str);
}
if ( err )
{
SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
if ( ce && ce->HasBadElems() )
return error( ce );
}
return false;
}
bool mycomputeFillNewElementFile(
std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib,
@ -545,7 +450,7 @@ bool mycomputeFillNewElementFile(
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
bool isOK = ( Netgen_NbOfTetra > 0 );
if ( isOK && !new_element_file.empty() )
{
std::ofstream df(new_element_file, ios::out|ios::binary);
@ -583,58 +488,6 @@ bool mycomputeFillNewElementFile(
return false;
}
bool mycomputeFillMesh(
std::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);
// -------------------------------------------------------------------
// Feed back the SMESHDS with the generated Nodes and Volume Elements
// -------------------------------------------------------------------
// Adding new elements in aMesh as well
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;
}
/**
* @brief Compute aShape within aMesh using netgen3d
*
@ -660,17 +513,14 @@ int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aPa
int Netgen_NbOfNodes=0;
bool ret;
std::cout << "mycomputeFillNgMesh" << std::endl;
ret = mycomputeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, element_orientation_file, Netgen_NbOfNodes);
if(ret)
return error( aParams._error, aParams._comment);
std::cout << "mycomputePrepareParam" << std::endl;
netgen::OCCGeometry occgeo;
mycomputePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
std::cout << "mycomputeRunMesher" << std::endl;
NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
ret = mycomputeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);
ret = NETGENPlugin_NETGEN_3D::computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, aParams, startWith, endWith);
if(ret){
if(aParams._error)
return error(aParams._error, aParams._comment);
@ -678,15 +528,11 @@ int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aPa
error(aParams._comment);
return true;
}
std::cout << "mycomputeFillNewElementFile" << std::endl;
mycomputeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
std::cout << "mycomputeFillMesh" << std::endl;
if(output_mesh)
mycomputeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
std::cout << "Done" << std::endl;
NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
return false;
}