Refactoring of runner

This commit is contained in:
Yoann Audouin 2022-09-19 15:13:31 +02:00
parent 9ca5a5ccb4
commit a52ccf7016
2 changed files with 248 additions and 371 deletions

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( bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
SMESH_Mesh& aMesh, SMESH_Mesh& aMesh,
@ -516,9 +580,8 @@ bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
// --------------------------------- // ---------------------------------
// Feed the Netgen with surface mesh // Feed the Netgen with surface mesh
// --------------------------------- // ---------------------------------
bool isRev=false;
TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); bool isInternalFace=false;
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 ( aParams._viscousLayersHyp ) if ( aParams._viscousLayersHyp )
@ -536,88 +599,56 @@ bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
proxyMesh.reset( Adaptor ); 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(); isRev = get<0>(info);
int faceID = meshDS->ShapeToIndex( aShapeFace ); isInternalFace = get<1>(info);
bool isInternalFace = internals.isInternalShape( faceID ); // Add nodes of triangles and triangles them-selves to netgen mesh
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 ); // add three nodes of triangle
if ( !aSubMeshDSFace ) continue; bool hasDegen = false;
for ( int iN = 0; iN < 3; ++iN )
SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
if ( aParams._quadraticMesh &&
dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
{ {
// add medium nodes of proxy triangles to helper (#16843) const SMDS_MeshNode* node = elem->GetNode( iN );
while ( iteratorElem->more() ) const int shapeID = node->getshapeId();
helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() )); if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
helper.IsDegenShape( shapeID ))
iteratorElem = aSubMeshDSFace->GetElements(); {
// ignore all nodes on degeneraged edge and use node on its vertex instead
TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
hasDegen = true;
}
int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
if ( ngID == invalid_ID )
{
ngID = ++Netgen_NbOfNodes;
Netgen_point [ 0 ] = node->X();
Netgen_point [ 1 ] = node->Y();
Netgen_point [ 2 ] = node->Z();
Ng_AddPoint(Netgen_mesh, Netgen_point);
}
Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
} }
while ( iteratorElem->more() ) // loop on elements on a geom face // add triangle
if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
Netgen_triangle[0] == Netgen_triangle[2] ||
Netgen_triangle[2] == Netgen_triangle[1] ))
continue;
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
{ {
// check mesh face swap( Netgen_triangle[1], Netgen_triangle[2] );
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;
}
// Add nodes of triangles and triangles them-selves to netgen mesh
// add three nodes of triangle
bool hasDegen = false;
for ( int iN = 0; iN < 3; ++iN )
{
const SMDS_MeshNode* node = elem->GetNode( iN );
const int shapeID = node->getshapeId();
if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
helper.IsDegenShape( shapeID ))
{
// ignore all nodes on degeneraged edge and use node on its vertex instead
TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
hasDegen = true;
}
int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
if ( ngID == invalid_ID )
{
ngID = ++Netgen_NbOfNodes;
Netgen_point [ 0 ] = node->X();
Netgen_point [ 1 ] = node->Y();
Netgen_point [ 2 ] = node->Z();
Ng_AddPoint(Netgen_mesh, Netgen_point);
}
Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
}
// add triangle
if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
Netgen_triangle[0] == Netgen_triangle[2] ||
Netgen_triangle[2] == Netgen_triangle[1] ))
continue;
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
}
if ( isInternalFace && !proxyMesh->IsTemporary( elem )) } // loop on elements on a face
{
swap( Netgen_triangle[1], Netgen_triangle[2] );
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 // insert old nodes into nodeVec
nodeVec.resize( nodeToNetgenID.size() + 1, 0 ); nodeVec.resize( nodeToNetgenID.size() + 1, 0 );

View File

@ -27,7 +27,7 @@
#include "NETGENPlugin_Runner.hxx" #include "NETGENPlugin_Runner.hxx"
#include "NETGENPlugin_NETGEN_3D.hxx"
#include "NETGENPlugin_DriverParam.hxx" #include "NETGENPlugin_DriverParam.hxx"
#include <fstream> #include <fstream>
@ -228,6 +228,84 @@ int netgen3d(const std::string input_mesh_file,
return ret; return ret;
} }
bool getSurfaceElements(
SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape,
SMESH_ProxyMesh::Ptr proxyMesh,
NETGENPlugin_Internals &internals,
SMESH_MesherHelper &helper,
netgen_params &aParams,
std::string element_orientation_file,
std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
)
{
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()){
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 {
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 isInternalFace = false;
bool isIn;
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;
}
// 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( bool mycomputeFillNgMesh(
SMESH_Mesh& aMesh, SMESH_Mesh& aMesh,
const TopoDS_Shape& aShape, const TopoDS_Shape& aShape,
@ -238,14 +316,16 @@ bool mycomputeFillNgMesh(
std::string element_orientation_file, std::string element_orientation_file,
int &Netgen_NbOfNodes) int &Netgen_NbOfNodes)
{ {
netgen::multithread.terminate = 0; netgen::multithread.terminate = 0;
netgen::multithread.task = "Volume meshing"; netgen::multithread.task = "Volume meshing";
aParams._progressByTic = -1.;
SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape); aParams._quadraticMesh = helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true ); helper.SetElementsOnShape( true );
Netgen_NbOfNodes = 0;
double Netgen_point[3]; double Netgen_point[3];
int Netgen_triangle[3]; int Netgen_triangle[3];
@ -268,9 +348,8 @@ bool mycomputeFillNgMesh(
// --------------------------------- // ---------------------------------
// Feed the Netgen with surface mesh // Feed the Netgen with surface mesh
// --------------------------------- // ---------------------------------
bool isRev=false;
TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); bool isInternalFace=false;
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 ( aParams._viscousLayersHyp ) if ( aParams._viscousLayersHyp )
@ -288,106 +367,56 @@ bool mycomputeFillNgMesh(
proxyMesh.reset( Adaptor ); proxyMesh.reset( Adaptor );
} }
// Get list of elements + their orientation from element_orientation file std::map<const SMDS_MeshElement*, tuple<bool, bool>> listElements;
std::map<vtkIdType, bool> elemOrientation; 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
{ {
// Setting all element orientation to false if there no element orientation file isRev = get<0>(info);
if(element_orientation_file.empty()){ isInternalFace = get<1>(info);
SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face); // Add nodes of triangles and triangles them-selves to netgen mesh
while ( iteratorElem->more() ) // loop on elements on a geom face
{
// check mesh face
const SMDS_MeshElement* elem = iteratorElem->next();
elemOrientation[elem->GetID()] = false;
}
} else {
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 // add three nodes of triangle
// Sizeof was the same but how he othered the type was different bool hasDegen = false;
// Maybe using another type (uint64_t) instead would be better for ( int iN = 0; iN < 3; ++iN )
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 isInternalFace = false;
bool isIn;
while ( iteratorElem->more() ) // loop on elements on a geom face
{ {
// check mesh face const SMDS_MeshNode* node = elem->GetNode( iN );
const SMDS_MeshElement* elem = iteratorElem->next(); const int shapeID = node->getshapeId();
if ( !elem ) if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); helper.IsDegenShape( shapeID ))
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()];
// Add nodes of triangles and triangles them-selves to netgen mesh
// add three nodes of triangle
bool hasDegen = false;
for ( int iN = 0; iN < 3; ++iN )
{ {
const SMDS_MeshNode* node = elem->GetNode( iN ); // ignore all nodes on degeneraged edge and use node on its vertex instead
const int shapeID = node->getshapeId(); TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
helper.IsDegenShape( shapeID )) hasDegen = true;
{
// ignore all nodes on degeneraged edge and use node on its vertex instead
TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
hasDegen = true;
}
int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
if ( ngID == invalid_ID )
{
ngID = ++Netgen_NbOfNodes;
Netgen_point [ 0 ] = node->X();
Netgen_point [ 1 ] = node->Y();
Netgen_point [ 2 ] = node->Z();
Ng_AddPoint(Netgen_mesh, Netgen_point);
}
Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
} }
// add triangle int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] || if ( ngID == invalid_ID )
Netgen_triangle[0] == Netgen_triangle[2] ||
Netgen_triangle[2] == Netgen_triangle[1] ))
continue;
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] ); ngID = ++Netgen_NbOfNodes;
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); Netgen_point [ 0 ] = node->X();
Netgen_point [ 1 ] = node->Y();
Netgen_point [ 2 ] = node->Z();
Ng_AddPoint(Netgen_mesh, Netgen_point);
} }
Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
} }
// add triangle
if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
Netgen_triangle[0] == Netgen_triangle[2] ||
Netgen_triangle[2] == Netgen_triangle[1] ))
continue;
Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
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 // insert old nodes into nodeVec
nodeVec.resize( nodeToNetgenID.size() + 1, 0 ); nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
@ -396,143 +425,19 @@ bool mycomputeFillNgMesh(
nodeVec[ n_id->second ] = n_id->first; nodeVec[ n_id->second ] = n_id->first;
nodeToNetgenID.clear(); nodeToNetgenID.clear();
// TODO: Handle internal vertex if ( internals.hasInternalVertexInSolid() )
//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::OCCGeometry occgeo;
{ NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
netgen::Point3d pmin, pmax; (netgen::Mesh&) *Netgen_mesh,
ngMesh->GetBox( pmin, pmax, 0 ); nodeVec,
ngMesh->SetLocalH( pmin, pmax, aParams.grading ); internals);
}
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_NbOfNodes = Ng_GetNP( Netgen_mesh );
{
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; 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( bool mycomputeFillNewElementFile(
std::vector< const SMDS_MeshNode* > &nodeVec, std::vector< const SMDS_MeshNode* > &nodeVec,
NETGENPlugin_NetgenLibWrapper &ngLib, NETGENPlugin_NetgenLibWrapper &ngLib,
@ -545,7 +450,7 @@ bool mycomputeFillNewElementFile(
int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
int Netgen_NbOfTetra = Ng_GetNE(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() ) if ( isOK && !new_element_file.empty() )
{ {
std::ofstream df(new_element_file, ios::out|ios::binary); std::ofstream df(new_element_file, ios::out|ios::binary);
@ -583,58 +488,6 @@ bool mycomputeFillNewElementFile(
return false; 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 * @brief Compute aShape within aMesh using netgen3d
* *
@ -651,7 +504,7 @@ int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aPa
std::string new_element_file, std::string element_orientation_file, std::string new_element_file, std::string element_orientation_file,
bool output_mesh) bool output_mesh)
{ {
// vector of nodes in which node index == netgen ID // vector of nodes in which node index == netgen ID
vector< const SMDS_MeshNode* > nodeVec; vector< const SMDS_MeshNode* > nodeVec;
NETGENPlugin_NetgenLibWrapper ngLib; NETGENPlugin_NetgenLibWrapper ngLib;
SMESH_MesherHelper helper(aMesh); SMESH_MesherHelper helper(aMesh);
@ -660,17 +513,14 @@ int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aPa
int Netgen_NbOfNodes=0; int Netgen_NbOfNodes=0;
bool ret; bool ret;
std::cout << "mycomputeFillNgMesh" << std::endl;
ret = mycomputeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, element_orientation_file, Netgen_NbOfNodes); ret = mycomputeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, aParams, element_orientation_file, Netgen_NbOfNodes);
if(ret) if(ret)
return error( aParams._error, aParams._comment); return error( aParams._error, aParams._comment);
std::cout << "mycomputePrepareParam" << std::endl;
netgen::OCCGeometry occgeo; netgen::OCCGeometry occgeo;
mycomputePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith); NETGENPlugin_NETGEN_3D::computePrepareParam(aMesh, ngLib, occgeo, helper, aParams, endWith);
std::cout << "mycomputeRunMesher" << std::endl;
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(ret){
if(aParams._error) if(aParams._error)
return error(aParams._error, aParams._comment); return error(aParams._error, aParams._comment);
@ -678,15 +528,11 @@ int netgen3dInternal(TopoDS_Shape &aShape, SMESH_Mesh& aMesh, netgen_params& aPa
error(aParams._comment); error(aParams._comment);
return true; return true;
} }
std::cout << "mycomputeFillNewElementFile" << std::endl;
mycomputeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes); mycomputeFillNewElementFile(nodeVec, ngLib, new_element_file, Netgen_NbOfNodes);
std::cout << "mycomputeFillMesh" << std::endl;
if(output_mesh) if(output_mesh)
mycomputeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes); NETGENPlugin_NETGEN_3D::computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
std::cout << "Done" << std::endl;
return false; return false;
} }