extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
#endif
//extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
+#if defined(NETGEN_V5) && defined(WIN32)
+ DLL_HEADER
+#endif
extern MeshingParameters mparam;
+#if defined(NETGEN_V5) && defined(WIN32)
+ DLL_HEADER
+#endif
extern volatile multithreadt multithread;
+
+#if defined(NETGEN_V5) && defined(WIN32)
+ DLL_HEADER
+#endif
extern bool merge_solids;
// values used for occgeo.facemeshstatus
_progressTic(1),
_totalTime(1.0),
_simpleHyp(NULL),
+ _viscousLayersHyp(NULL),
_ptrToMe(NULL)
{
SetDefaultParameters();
SetDefaultParameters();
}
+//================================================================================
+/*!
+ * \brief Store a Viscous Layers hypothesis
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::SetParameters(const StdMeshers_ViscousLayers* hyp )
+{
+ _viscousLayersHyp = hyp;
+}
+
//=============================================================================
/*!
* Link - a pair of integer numbers
}
else
{
- minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine
+ minh = sqrt( minh ); // triangulation for visualization is rather fine
//cout << "TRIANGULATION minh = " <<minh << endl;
}
if ( minh > 0.5 * maxSize )
// get all nodes from connected <edges>
const bool isQuad = smDS->IsQuadratic();
- StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad );
+ StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad, &helper );
const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
if ( points.empty() )
return false; // invalid node params?
// Find solids the geomFace bounds
int solidID1 = 0, solidID2 = 0;
- StdMeshers_QuadToTriaAdaptor* quadAdaptor =
- dynamic_cast<StdMeshers_QuadToTriaAdaptor*>( proxyMesh.get() );
- if ( quadAdaptor )
- {
- solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() );
- }
- else
{
PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
while ( const TopoDS_Shape * solid = solidIt->next() )
else solidID1 = id;
}
}
+ if ( proxyMesh && proxyMesh->GetProxySubMesh( geomFace ))
+ {
+ // if a proxy sub-mesh contains temporary faces, then these faces
+ // should be used to mesh only one SOLID
+ bool hasTmp = false;
+ smDS = proxyMesh->GetSubMesh( geomFace );
+ SMDS_ElemIteratorPtr faces = smDS->GetElements();
+ while ( faces->more() )
+ {
+ const SMDS_MeshElement* f = faces->next();
+ if ( proxyMesh->IsTemporary( f ))
+ {
+ hasTmp = true;
+ std::vector<const SMDS_MeshNode*> fNodes( f->begin_nodes(), f->end_nodes() );
+ std::vector<const SMDS_MeshElement*> vols;
+ if ( _mesh->GetMeshDS()->GetElementsByNodes( fNodes, vols, SMDSAbs_Volume ) == 1 )
+ {
+ int geomID = vols[0]->getshapeId();
+ const TopoDS_Shape& solid = helper.GetMeshDS()->IndexToShape( geomID );
+ if ( !solid.IsNull() )
+ solidID1 = occgeom.somap.FindIndex ( solid );
+ solidID2 = 0;
+ break;
+ }
+ }
+ }
+ // exclude faces generated by NETGEN from computation of 3D mesh
+ const int fID = occgeom.fmap.FindIndex( geomFace );
+ if ( !hasTmp ) // shrunk mesh
+ {
+ // move netgen points according to moved nodes
+ SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
+ while ( smIt->more() )
+ {
+ SMESH_subMesh* sub = smIt->next();
+ if ( !sub->GetSubMeshDS() ) continue;
+ SMDS_NodeIteratorPtr nodeIt = sub->GetSubMeshDS()->GetNodes();
+ while ( nodeIt->more() )
+ {
+ const SMDS_MeshNode* n = nodeIt->next();
+ int ngID = ngNodeId( n, ngMesh, nodeNgIdMap );
+ netgen::MeshPoint& ngPoint = ngMesh.Point( ngID );
+ ngPoint(0) = n->X();
+ ngPoint(1) = n->Y();
+ ngPoint(2) = n->Z();
+ }
+ }
+ // remove faces near boundary to avoid their overlapping
+ // with shrunk faces
+ for ( int i = 1; i <= ngMesh.GetNSE(); ++i )
+ {
+ const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+ if ( elem.GetIndex() == fID )
+ {
+ for ( int iN = 0; iN < elem.GetNP(); ++iN )
+ if ( ngMesh[ elem[ iN ]].Type() != netgen::SURFACEPOINT )
+ {
+ ngMesh.DeleteSurfaceElement( i );
+ break;
+ }
+ }
+ }
+ }
+ //if ( hasTmp )
+ {
+ faceNgID++;
+ ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID,/*solid1=*/0,/*solid2=*/0,0 ));
+ for (int i = 1; i <= ngMesh.GetNSE(); ++i )
+ {
+ const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
+ if ( elem.GetIndex() == fID )
+ const_cast< netgen::Element2d& >( elem ).SetIndex( faceNgID );
+ }
+ }
+ }
// Add ng face descriptors of meshed faces
faceNgID++;
ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, solidID2, 0 ));
cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
<< " internal="<<isInternalFace << endl;
#endif
- if ( proxyMesh )
- smDS = proxyMesh->GetSubMesh( geomFace );
SMDS_ElemIteratorPtr faces = smDS->GetElements();
while ( faces->more() )
int ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
};
- inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
+ inline double dist2( const netgen::MeshPoint& p1, const netgen::MeshPoint& p2 )
{
return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
}
+
+ // inline double dist2(const netgen::MeshPoint& p, const SMDS_MeshNode* n )
+ // {
+ // return gp_Pnt( NGPOINT_COORDS(p)).SquareDistance( SMESH_NodeXYZ(n));
+ // }
}
//================================================================================
if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
quadHelper = 0;
+ int i, nbInitNod = initState._nbNodes;
+ if ( initState._elementsRemoved )
+ {
+ // PAL23427. Update nodeVec to track removal of netgen free points as a result
+ // of removal of faces in FillNgMesh() in the case of a shrunk sub-mesh
+ int ngID, nodeVecSize = nodeVec.size();
+ const double eps = std::numeric_limits<double>::min();
+ for ( ngID = i = 1; i < nodeVecSize; ++ngID, ++i )
+ {
+ gp_Pnt ngPnt( NGPOINT_COORDS( ngMesh.Point( ngID )));
+ gp_Pnt node ( SMESH_NodeXYZ ( nodeVec[ i ]));
+ if ( ngPnt.SquareDistance( node ) < eps )
+ {
+ nodeVec[ ngID ] = nodeVec[ i ];
+ }
+ else
+ {
+ --ngID;
+ }
+ }
+ nodeVec.resize( ngID );
+ nbInitNod = ngID - 1;
+ }
// -------------------------------------
// Create and insert nodes into nodeVec
// -------------------------------------
nodeVec.resize( nbNod + 1 );
- int i, nbInitNod = initState._nbNodes;
- for (i = nbInitNod+1; i <= nbNod; ++i )
+ for ( i = nbInitNod+1; i <= nbNod; ++i )
{
const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
SMDS_MeshNode* node = NULL;
SMESH_MesherHelper quadHelper( *_mesh );
quadHelper.SetIsQuadratic( mparams.secondorder );
- static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( _ngMesh, debugFile )
- while debugging netgen */
// -------------------------
// Prepare OCC geometry
// -------------------------
helper.SetSubShape( F );
TSideVector wires =
StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
- error, viscousMesh );
+ error, &helper, viscousMesh );
error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
if ( !error ) error = SMESH_ComputeError::New();
// generate volume mesh
// ---------------------
// Fill _ngMesh with nodes and faces of computed 2D submeshes
- if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
+ if ( !err && _isVolume &&
+ ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp ))
{
// load SMESH with computed segments and faces
FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
+ // compute prismatic boundary volumes
+ int nbQuad = _mesh->NbQuadrangles();
+ SMESH_ProxyMesh::Ptr viscousMesh;
+ if ( _viscousLayersHyp )
+ {
+ viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape );
+ if ( !viscousMesh )
+ return false;
+ }
// compute pyramids on quadrangles
- SMESH_ProxyMesh::Ptr proxyMesh;
- if ( _mesh->NbQuadrangles() > 0 )
+ vector<SMESH_ProxyMesh::Ptr> pyramidMeshes( occgeo.somap.Extent() );
+ if ( nbQuad > 0 )
for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
{
- StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
- proxyMesh.reset( Adaptor );
-
- int nbPyrams = _mesh->NbPyramids();
- Adaptor->Compute( *_mesh, occgeo.somap(iS) );
- if ( nbPyrams != _mesh->NbPyramids() )
+ StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor;
+ pyramidMeshes[ iS-1 ].reset( adaptor );
+ bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() );
+ if ( !ok )
+ return false;
+ }
+ // add proxy faces to NG mesh
+ list< SMESH_subMesh* > viscousSM;
+ for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
+ {
+ list< SMESH_subMesh* > quadFaceSM;
+ for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
+ if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() ))
{
- list< SMESH_subMesh* > quadFaceSM;
- for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
- if ( Adaptor->GetProxySubMesh( face.Current() ))
- {
- quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
- meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
- }
- FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh);
+ quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
+ meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
}
- }
+ else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() ))
+ {
+ viscousSM.push_back( _mesh->GetSubMesh( face.Current() ));
+ meshedSM[ MeshDim_2D ].remove( viscousSM.back() );
+ }
+ if ( !quadFaceSM.empty() )
+ FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]);
+ }
+ if ( !viscousSM.empty() )
+ FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh );
+
// fill _ngMesh with faces of sub-meshes
err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
- initState = NETGENPlugin_ngMeshInfo(_ngMesh);
- //toPython( _ngMesh, "/tmp/ngPython.py");
+ initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true);
+ // toPython( _ngMesh );
}
if (!err && _isVolume)
{
ofstream outfile( pyFile, ios::out );
if ( !outfile ) return;
- outfile << "import SMESH" << endl
+ outfile << "import salome, SMESH" << endl
<< "from salome.smesh import smeshBuilder" << endl
<< "smesh = smeshBuilder.New()" << endl
<< "mesh = smesh.Mesh()" << endl << endl;
*/
//================================================================================
-NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
- _copyOfLocalH(0)
+NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh,
+ bool checkRemovedElems):
+ _elementsRemoved( false ), _copyOfLocalH(0)
{
if ( ngMesh )
{
_nbSegments = ngMesh->GetNSeg();
_nbFaces = ngMesh->GetNSE();
_nbVolumes = ngMesh->GetNE();
+
+ if ( checkRemovedElems )
+ for ( int i = 1; i <= ngMesh->GetNSE() && !_elementsRemoved; ++i )
+ _elementsRemoved = ngMesh->SurfaceElement(i).IsDeleted();
}
else
{
