X-Git-Url: http://git.salome-platform.org/gitweb/?p=plugins%2Fnetgenplugin.git;a=blobdiff_plain;f=src%2FNETGENPlugin%2FNETGENPlugin_Mesher.cxx;h=d0979451989a9480e2e77a86e0c174b21bc2734f;hp=5d5c401002498e3878c19df461a38adc993caf63;hb=42a2e30144e2a98ea1fc906e80946b2421ad6906;hpb=e1a0d0a23f4ae76cef3888df02e448877ebb8d19 diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 5d5c401..d097945 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -52,21 +52,29 @@ #include +#include #include +#include +#include +#include #include #include +#include #include +#include #include #include #include #include #include +#include #include #include #include #include #include #include +#include // Netgen include files #ifndef OCCGEOMETRY @@ -82,8 +90,18 @@ namespace netgen { 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 @@ -141,12 +159,14 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh, _optimize(true), _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()), _isViscousLayers2D(false), + _chordalError(-1), // means disabled _ngMesh(NULL), _occgeom(NULL), _curShapeIndex(-1), _progressTic(1), _totalTime(1.0), _simpleHyp(NULL), + _viscousLayersHyp(NULL), _ptrToMe(NULL) { SetDefaultParameters(); @@ -161,7 +181,7 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh, //================================================================================ /*! - * Destuctor + * Destructor */ //================================================================================ @@ -282,17 +302,15 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) _fineness = hyp->GetFineness(); mparams.uselocalh = hyp->GetSurfaceCurvature(); netgen::merge_solids = hyp->GetFuseEdges(); + _chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.; _simpleHyp = NULL; // mesh size file mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str(); - SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); - CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager"); - SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject); - SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId()); - if ( !myStudy->_is_nil() ) + const NETGENPlugin_Hypothesis::TLocalSize& localSizes = hyp->GetLocalSizesAndEntries(); + if ( !localSizes.empty() ) { - const NETGENPlugin_Hypothesis::TLocalSize localSizes = hyp->GetLocalSizesAndEntries(); + SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin(); for ( ; it != localSizes.end() ; it++) { @@ -300,7 +318,7 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) double val = (*it).second; // -- GEOM::GEOM_Object_var aGeomObj; - SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() ); + SALOMEDS::SObject_var aSObj = SMESH_Gen_i::getStudyServant()->FindObjectID( entry.c_str() ); if ( !aSObj->_is_nil() ) { CORBA::Object_var obj = aSObj->GetObject(); aGeomObj = GEOM::GEOM_Object::_narrow(obj); @@ -326,6 +344,17 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* SetDefaultParameters(); } +//================================================================================ +/*! + * \brief Store a Viscous Layers hypothesis + */ +//================================================================================ + +void NETGENPlugin_Mesher::SetParameters(const StdMeshers_ViscousLayers* hyp ) +{ + _viscousLayersHyp = hyp; +} + //============================================================================= /*! * Link - a pair of integer numbers @@ -575,7 +604,8 @@ namespace void setLocalSize(const TopoDS_Edge& edge, double size, - netgen::Mesh& mesh) + netgen::Mesh& mesh, + const bool overrideMinH = true) { if ( size <= std::numeric_limits::min() ) return; @@ -586,7 +616,7 @@ namespace TopoDS_Iterator vIt( edge ); if ( !vIt.More() ) return; gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() )); - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH ); } else { @@ -596,15 +626,29 @@ namespace { Standard_Real u = u1 + delta*i; gp_Pnt p = curve->Value(u); - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH ); netgen::Point3d pi(p.X(), p.Y(), p.Z()); double resultSize = mesh.GetH(pi); if ( resultSize - size > 0.1*size ) // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136) - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201, overrideMinH ); } } } + + //================================================================================ + /*! + * \brief Return triangle size for a given chordalError and radius of curvature + */ + //================================================================================ + + double elemSizeForChordalError( double chordalError, double radius ) + { + if ( 2 * radius < chordalError ) + return 1.5 * radius; + return Sqrt( 3 ) * Sqrt( chordalError * ( 2 * radius - chordalError )); + } + } // namespace //================================================================================ @@ -614,16 +658,19 @@ namespace //================================================================================ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, - netgen::Mesh& ngMesh ) + netgen::Mesh& ngMesh) { - for(std::map::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++) + // edges + std::map::const_iterator it; + for( it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++) { int key = (*it).first; double hi = (*it).second; const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key); setLocalSize( TopoDS::Edge(shape), hi, ngMesh ); } - for(std::map::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++) + // vertices + for(it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++) { int key = (*it).first; double hi = (*it).second; @@ -631,7 +678,8 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex(shape) ); NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, p.XYZ(), hi ); } - for(map::const_iterator it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++) + // faces + for(it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++) { int key = (*it).first; double val = (*it).second; @@ -649,7 +697,8 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, ShapesWithControlPoints.insert( key ); } } - for(map::const_iterator it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++) + //solids + for(it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++) { int key = (*it).first; double val = (*it).second; @@ -666,6 +715,146 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, for ( size_t i = 0; i < ControlPoints.size(); ++i ) NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, ControlPoints[i].XYZ(), ControlPoints[i].Size() ); } + return; +} + +//================================================================================ +/*! + * \brief Restrict local size to achieve a required _chordalError + */ +//================================================================================ + +void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo, + netgen::Mesh& ngMesh) +{ + if ( _chordalError <= 0. ) + return; + + TopLoc_Location loc; + BRepLProp_SLProps surfProp( 2, 1e-6 ); + const double sizeCoef = 0.95; + + // find non-planar FACEs with non-constant curvature + std::vector fInd; + for ( int i = 1; i <= occgeo.fmap.Extent(); ++i ) + { + const TopoDS_Face& face = TopoDS::Face( occgeo.fmap( i )); + BRepAdaptor_Surface surfAd( face, false ); + switch ( surfAd.GetType() ) + { + case GeomAbs_Plane: + continue; + case GeomAbs_Cylinder: + case GeomAbs_Sphere: + case GeomAbs_Torus: // constant curvature + { + surfProp.SetSurface( surfAd ); + surfProp.SetParameters( 0, 0 ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + double size = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + occgeo.SetFaceMaxH( i, size * sizeCoef ); + // limit size one edges + TopTools_MapOfShape edgeMap; + for ( TopExp_Explorer eExp( face, TopAbs_EDGE ); eExp.More(); eExp.Next() ) + if ( edgeMap.Add( eExp.Current() )) + setLocalSize( TopoDS::Edge( eExp.Current() ), size, ngMesh, /*overrideMinH=*/false ); + break; + } + default: + Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc ); + if ( GeomLib_IsPlanarSurface( surf ).IsPlanar() ) + continue; + fInd.push_back( i ); + } + } + // set local size + if ( !fInd.empty() ) + { + BRep_Builder b; + TopoDS_Compound allFacesComp; + b.MakeCompound( allFacesComp ); + for ( size_t i = 0; i < fInd.size(); ++i ) + b.Add( allFacesComp, occgeo.fmap( fInd[i] )); + + // copy the shape to avoid spoiling its triangulation + TopoDS_Shape allFacesCompCopy = BRepBuilderAPI_Copy( allFacesComp ); + + // create triangulation with desired chordal error + BRepMesh_IncrementalMesh( allFacesCompCopy, + _chordalError, + /*isRelative = */Standard_False, + /*theAngDeflection = */ 0.5, + /*isInParallel = */Standard_True); + + // loop on FACEs + for ( TopExp_Explorer fExp( allFacesCompCopy, TopAbs_FACE ); fExp.More(); fExp.Next() ) + { + const TopoDS_Face& face = TopoDS::Face( fExp.Current() ); + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation ( face, loc ); + if ( triangulation.IsNull() ) continue; + + BRepAdaptor_Surface surf( face, false ); + surfProp.SetSurface( surf ); + + gp_XY uv[3]; + gp_XYZ p[3]; + double size[3]; + for ( int i = 1; i <= triangulation->NbTriangles(); ++i ) + { + Standard_Integer n1,n2,n3; + triangulation->Triangles()(i).Get( n1,n2,n3 ); + p [0] = triangulation->Nodes()(n1).Transformed(loc).XYZ(); + p [1] = triangulation->Nodes()(n2).Transformed(loc).XYZ(); + p [2] = triangulation->Nodes()(n3).Transformed(loc).XYZ(); + uv[0] = triangulation->UVNodes()(n1).XY(); + uv[1] = triangulation->UVNodes()(n2).XY(); + uv[2] = triangulation->UVNodes()(n3).XY(); + surfProp.SetParameters( uv[0].X(), uv[0].Y() ); + if ( !surfProp.IsCurvatureDefined() ) + break; + + for ( int n = 0; n < 3; ++n ) // get size at triangle nodes + { + surfProp.SetParameters( uv[n].X(), uv[n].Y() ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + size[n] = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + } + for ( int n1 = 0; n1 < 3; ++n1 ) // limit size along each triangle edge + { + int n2 = ( n1 + 1 ) % 3; + double minSize = size[n1], maxSize = size[n2]; + if ( size[n1] > size[n2] ) + minSize = size[n2], maxSize = size[n1]; + + if ( maxSize / minSize < 1.2 ) // netgen ignores size difference < 1.2 + { + ngMesh.RestrictLocalHLine ( netgen::Point3d( p[n1].X(), p[n1].Y(), p[n1].Z() ), + netgen::Point3d( p[n2].X(), p[n2].Y(), p[n2].Z() ), + sizeCoef * minSize ); + } + else + { + gp_XY uvVec( uv[n2] - uv[n1] ); + double len = ( p[n1] - p[n2] ).Modulus(); + int nb = int( len / minSize ) + 1; + for ( int j = 0; j <= nb; ++j ) + { + double r = double( j ) / nb; + gp_XY uvj = uv[n1] + r * uvVec; + + surfProp.SetParameters( uvj.X(), uvj.Y() ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + double h = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + + const gp_Pnt& pj = surfProp.Value(); + netgen::Point3d ngP( pj.X(), pj.Y(), pj.Z()); + ngMesh.RestrictLocalH( ngP, h * sizeCoef ); + } + } + } + } + } + } } //================================================================================ @@ -729,7 +918,7 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo, if ( shape.ShapeType() != TopAbs_VERTEX ) shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape if ( shape.Orientation() >= TopAbs_INTERNAL ) - shape.Orientation( TopAbs_FORWARD ); // isuue 0020676 + shape.Orientation( TopAbs_FORWARD ); // issue 0020676 switch ( shape.ShapeType() ) { case TopAbs_FACE : occgeo.fmap.Add( shape ); break; case TopAbs_EDGE : occgeo.emap.Add( shape ); break; @@ -798,7 +987,7 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom, } else { - minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine + minh = sqrt( minh ); // triangulation for visualization is rather fine //cout << "TRIANGULATION minh = " < 0.5 * maxSize ) @@ -916,7 +1105,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, // get all nodes from connected const bool isQuad = smDS->IsQuadratic(); - StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad ); + StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad, &helper ); const vector& points = fSide.GetUVPtStruct(); if ( points.empty() ) return false; // invalid node params? @@ -1044,13 +1233,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, // Find solids the geomFace bounds int solidID1 = 0, solidID2 = 0; - StdMeshers_QuadToTriaAdaptor* quadAdaptor = - dynamic_cast( 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() ) @@ -1060,6 +1242,81 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, 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 fNodes( f->begin_nodes(), f->end_nodes() ); + std::vector 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 )); @@ -1102,8 +1359,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace ) << " internal="<GetSubMesh( geomFace ); SMDS_ElemIteratorPtr faces = smDS->GetElements(); while ( faces->more() ) @@ -1114,9 +1369,11 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID); if ( const TopoDS_Shape * solid = solidIt->next() ) sm = _mesh->GetSubMesh( *solid ); - SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); - smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle sub-mesh")); - smError->myBadElements.push_back( f ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( helper.GetMeshDS(), COMPERR_BAD_INPUT_MESH, + "Not triangle sub-mesh"); + badElems->add( f ); + sm->GetComputeError().reset( badElems ); return false; } @@ -1402,7 +1659,7 @@ namespace struct TIntVData { gp_XY uv; //!< UV in face parametric space - int ngId; //!< ng id of corrsponding node + int ngId; //!< ng id of corresponding node gp_XY uvClose; //!< UV of closest boundary node int ngIdClose; //!< ng id of closest boundary node }; @@ -1420,10 +1677,15 @@ namespace 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)); + // } } //================================================================================ @@ -1631,7 +1893,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry& ofstream py(DUMP_TRIANGLES_SCRIPT); py << "import SMESH"<< endl << "from salome.smesh import smeshBuilder"<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::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; @@ -2493,8 +2777,6 @@ bool NETGENPlugin_Mesher::Compute() 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 // ------------------------- @@ -2611,6 +2893,7 @@ bool NETGENPlugin_Mesher::Compute() { // Local size on shapes SetLocalSize( occgeo, *_ngMesh ); + SetLocalSizeForChordalError( occgeo, *_ngMesh ); } // Precompute internal edges (issue 0020676) in order to @@ -2753,8 +3036,9 @@ bool NETGENPlugin_Mesher::Compute() } // Build viscous layers - if ( _isViscousLayers2D || - StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh )) + if (( _isViscousLayers2D ) || + ( !occgeo.fmap.IsEmpty() && + StdMeshers_ViscousLayers2D::HasProxyMesh( TopoDS::Face( occgeo.fmap(1) ), *_mesh ))) { if ( !internals.hasInternalVertexInFace() ) { FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); @@ -2781,7 +3065,7 @@ bool NETGENPlugin_Mesher::Compute() 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(); @@ -2825,37 +3109,58 @@ bool NETGENPlugin_Mesher::Compute() // 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 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) { @@ -2905,18 +3210,18 @@ bool NETGENPlugin_Mesher::Compute() if(netgen::multithread.terminate) return false; - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(err); } catch (Standard_Failure& ex) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(ex); err = 1; } catch (netgen::NgException exc) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(exc); err = 1; } @@ -2937,17 +3242,17 @@ bool NETGENPlugin_Mesher::Compute() if(netgen::multithread.terminate) return false; - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(err); } catch (Standard_Failure& ex) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(ex); } catch (netgen::NgException exc) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << text(exc); } } @@ -2980,12 +3285,12 @@ bool NETGENPlugin_Mesher::Compute() } catch (Standard_Failure& ex) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << "Exception in netgen at passing to 2nd order "; } catch (netgen::NgException exc) { - if ( comment.empty() ) // do not overwrite a previos error + if ( comment.empty() ) // do not overwrite a previous error comment << exc.What(); } } @@ -3005,12 +3310,22 @@ bool NETGENPlugin_Mesher::Compute() FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); if ( quadHelper.GetIsQuadratic() ) // remove free nodes + { for ( size_t i = 0; i < nodeVec.size(); ++i ) if ( nodeVec[i] && nodeVec[i]->NbInverseElements() == 0 ) + { _mesh->GetMeshDS()->RemoveFreeNode( nodeVec[i], 0, /*fromGroups=*/false ); + nodeVec[i]=0; + } + for ( size_t i = nodeVec.size()-1; i > 0; --i ) // remove trailing removed nodes + if ( !nodeVec[i] ) + nodeVec.resize( i ); + else + break; + } } SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec); - if ( readErr && !readErr->myBadElements.empty() ) + if ( readErr && readErr->HasBadElems() ) { error = readErr; if ( !comment.empty() && !readErr->myComment.empty() ) comment += "\n"; @@ -3072,9 +3387,10 @@ bool NETGENPlugin_Mesher::Compute() { smError->myName = COMPERR_WARNING; } - else if ( !smError->myBadElements.empty() ) // bad surface mesh + else if ( smError->HasBadElems() ) // bad surface mesh { - if ( !hasBadElemOnSolid( smError->myBadElements, sm )) + if ( !hasBadElemOnSolid + ( static_cast( smError.get() )->myBadElements, sm )) smError.reset(); } } @@ -3103,9 +3419,8 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) // Prepare OCC geometry // ------------------------- netgen::OCCGeometry occgeo; - list< SMESH_subMesh* > meshedSM[4]; // for 0-3 dimensions NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume ); - PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals ); + PrepareOCCgeometry( occgeo, _shape, *_mesh, 0, &internals ); bool tooManyElems = false; const int hugeNb = std::numeric_limits::max() / 100; @@ -3156,25 +3471,25 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED )); return false; } - if ( _simpleHyp ) - { - // Pass 1D simple parameters to NETGEN - // -------------------------------- - int nbSeg = _simpleHyp->GetNumberOfSegments(); - double segSize = _simpleHyp->GetLocalLength(); - for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) - { - const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); - if ( nbSeg ) - segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); - setLocalSize( e, segSize, *ngMesh ); - } - } - else // if ( ! _simpleHyp ) - { - // Local size on shapes - SetLocalSize( occgeo, *ngMesh ); - } + // if ( _simpleHyp ) + // { + // // Pass 1D simple parameters to NETGEN + // // -------------------------------- + // int nbSeg = _simpleHyp->GetNumberOfSegments(); + // double segSize = _simpleHyp->GetLocalLength(); + // for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) + // { + // const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); + // if ( nbSeg ) + // segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); + // setLocalSize( e, segSize, *ngMesh ); + // } + // } + // else // if ( ! _simpleHyp ) + // { + // // Local size on shapes + // SetLocalSize( occgeo, *ngMesh ); + // } // calculate total nb of segments and length of edges double fullLen = 0.0; int fullNbSeg = 0; @@ -3253,9 +3568,9 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) int nb1d = 0; if ( !tooManyElems ) { - TopTools_MapOfShape egdes; + TopTools_MapOfShape edges; for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next()) - if ( egdes.Add( exp1.Current() )) + if ( edges.Add( exp1.Current() )) nb1d += Edge2NbSeg.Find(exp1.Current()); } int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.))); @@ -3404,8 +3719,10 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder, SMESH_ComputeErrorPtr NETGENPlugin_Mesher::ReadErrors(const vector& nodeVec) { - SMESH_ComputeErrorPtr err = SMESH_ComputeError::New - (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh"); + if ( nodeVec.size() < 2 ) return SMESH_ComputeErrorPtr(); + SMESH_BadInputElements* err = + new SMESH_BadInputElements( nodeVec.back()->GetMesh(), COMPERR_BAD_INPUT_MESH, + "Some edges multiple times in surface mesh"); SMESH_File file("test.out"); vector two(2); vector three1(3), three2(3); @@ -3465,7 +3782,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector& nodeVec) if ( nbBadElems ) nbBadElems++; // avoid warning: variable set but not used #endif - return err; + return SMESH_ComputeErrorPtr( err ); } //================================================================================ @@ -3482,9 +3799,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) 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(salome.myStudy)" << endl + << "smesh = smeshBuilder.New()" << endl << "mesh = smesh.Mesh()" << endl << endl; using namespace netgen; @@ -3546,8 +3863,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) */ //================================================================================ -NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh): - _copyOfLocalH(0) +NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh, + bool checkRemovedElems): + _elementsRemoved( false ), _copyOfLocalH(0) { if ( ngMesh ) { @@ -3555,6 +3873,10 @@ NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* 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 { @@ -3653,7 +3975,7 @@ NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh& mesh, { _intShapes.insert( meshDS->ShapeToIndex( f.Current() )); - // egdes + // edges list< TopoDS_Shape > edges; for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next()) if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 ) @@ -4090,10 +4412,10 @@ void NETGENPlugin_NetgenLibWrapper::removeOutputFile() } string tmpDir = SALOMEDS_Tool::GetDirFromPath ( _outputFileName ); string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out"; - SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames; - aFiles->length(1); - aFiles[0] = aFileName.c_str(); + SALOMEDS_Tool::ListOfFiles aFiles; + aFiles.reserve(1); + aFiles.push_back(aFileName.c_str()); - SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true ); + SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles, true ); } }