X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FNETGENPlugin%2FNETGENPlugin_Mesher.cxx;h=433c7bf7327528ffb9ba9447d5513ee16b2a6ca3;hb=99cf7a6a3cf1daa95f13c1d91044a5fac89c83b1;hp=a15b44f2d1f460e945d97f2bd2f62040a8e289b4;hpb=f39fa7b8d8ca01e534c2e5f5bf3a002f23759a61;p=plugins%2Fnetgenplugin.git diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index a15b44f..433c7bf 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -83,11 +83,15 @@ namespace netgen { //extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); extern MeshingParameters mparam; extern volatile multithreadt multithread; + extern bool merge_solids; } #include #include +#ifdef WIN32 +#include +#endif using namespace nglib; using namespace std; @@ -99,10 +103,12 @@ using namespace std; #define NGPOINT_COORDS(p) p(0),p(1),p(2) +#ifdef _DEBUG_ // dump elements added to ng mesh //#define DUMP_SEGMENTS //#define DUMP_TRIANGLES //#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug AddIntVerticesInSolids() +#endif TopTools_IndexedMapOfShape ShapesWithLocalSize; std::map VertexId2LocalSize; @@ -181,22 +187,24 @@ void NETGENPlugin_Mesher::SetDefaultParameters() { netgen::MeshingParameters& mparams = netgen::mparam; // maximal mesh edge size - mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize(); - mparams.minh = 0; + mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize(); + mparams.minh = 0; // minimal number of segments per edge mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge(); // rate of growth of size between elements - mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate(); + mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate(); // safety factor for curvatures (elements per radius) mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius(); // create elements of second order - mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0; + mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder(); // quad-dominated surface meshing if (_isVolume) - mparams.quad = 0; + mparams.quad = 0; else - mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0; - _fineness = NETGENPlugin_Hypothesis::GetDefaultFineness(); + mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed(); + _fineness = NETGENPlugin_Hypothesis::GetDefaultFineness(); + mparams.uselocalh = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature(); + netgen::merge_solids = NETGENPlugin_Hypothesis::GetDefaultFuseEdges(); } //============================================================================= @@ -239,23 +247,25 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) netgen::MeshingParameters& mparams = netgen::mparam; // Initialize global NETGEN parameters: // maximal mesh segment size - mparams.maxh = hyp->GetMaxSize(); + mparams.maxh = hyp->GetMaxSize(); // maximal mesh element linear size - mparams.minh = hyp->GetMinSize(); + mparams.minh = hyp->GetMinSize(); // minimal number of segments per edge mparams.segmentsperedge = hyp->GetNbSegPerEdge(); // rate of growth of size between elements - mparams.grading = hyp->GetGrowthRate(); + mparams.grading = hyp->GetGrowthRate(); // safety factor for curvatures (elements per radius) mparams.curvaturesafety = hyp->GetNbSegPerRadius(); // create elements of second order - mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0; + mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0; // quad-dominated surface meshing // only triangles are allowed for volumic mesh (before realizing IMP 0021676) //if (!_isVolume) - mparams.quad = hyp->GetQuadAllowed() ? 1 : 0; - _optimize = hyp->GetOptimize(); - _fineness = hyp->GetFineness(); + mparams.quad = hyp->GetQuadAllowed() ? 1 : 0; + _optimize = hyp->GetOptimize(); + _fineness = hyp->GetFineness(); + mparams.uselocalh = hyp->GetSurfaceCurvature(); + netgen::merge_solids = hyp->GetFuseEdges(); _simpleHyp = NULL; SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); @@ -469,11 +479,8 @@ namespace // { // BRepTools::Clean (shape); try { -#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100 OCC_CATCH_SIGNALS; -#endif BRepMesh_IncrementalMesh e(shape, 0.01, true); - } catch (Standard_Failure) { @@ -484,6 +491,60 @@ namespace // } // } } + //================================================================================ + /*! + * \brief Returns a medium node either existing in SMESH of created by NETGEN + * \param [in] corner1 - corner node 1 + * \param [in] corner2 - corner node 2 + * \param [in] defaultMedium - the node created by NETGEN + * \param [in] helper - holder of medium nodes existing in SMESH + * \return const SMDS_MeshNode* - the result node + */ + //================================================================================ + + const SMDS_MeshNode* mediumNode( const SMDS_MeshNode* corner1, + const SMDS_MeshNode* corner2, + const SMDS_MeshNode* defaultMedium, + const SMESH_MesherHelper* helper) + { + if ( helper ) + { + TLinkNodeMap::const_iterator l2n = + helper->GetTLinkNodeMap().find( SMESH_TLink( corner1, corner2 )); + if ( l2n != helper->GetTLinkNodeMap().end() ) + defaultMedium = l2n->second; + } + return defaultMedium; + } + + //================================================================================ + /*! + * \brief Assure that mesh on given shapes is quadratic + */ + //================================================================================ + + void makeQuadratic( const TopTools_IndexedMapOfShape& shapes, + SMESH_Mesh* mesh ) + { + for ( int i = 1; i <= shapes.Extent(); ++i ) + { + SMESHDS_SubMesh* smDS = mesh->GetMeshDS()->MeshElements( shapes(i) ); + if ( !smDS ) continue; + SMDS_ElemIteratorPtr elemIt = smDS->GetElements(); + if ( !elemIt->more() ) continue; + const SMDS_MeshElement* e = elemIt->next(); + if ( !e || e->IsQuadratic() ) + continue; + + TIDSortedElemSet elems; + elems.insert( e ); + while ( elemIt->more() ) + elems.insert( elems.end(), elemIt->next() ); + + SMESH_MeshEditor( mesh ).ConvertToQuadratic( /*3d=*/false, elems, /*biQuad=*/false ); + } + } + } //================================================================================ @@ -518,8 +579,9 @@ void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo, // get root submeshes list< SMESH_subMesh* > rootSM; - if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) { - rootSM.push_back( sm ); + const int shapeID = mesh.GetMeshDS()->ShapeToIndex( shape ); + if ( shapeID > 0 ) { // SMESH_subMesh with ID 0 may exist, don't use it! + rootSM.push_back( mesh.GetSubMesh( shape )); } else { for ( TopoDS_Iterator it( shape ); it.More(); it.Next() ) @@ -630,9 +692,12 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom, */ //================================================================================ -void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& p, const double size) +void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, + const gp_XYZ& p, + const double size, + const bool overrideMinH) { - if ( netgen::mparam.minh > size ) + if ( overrideMinH && netgen::mparam.minh > size ) { ngMesh.SetMinimalH( size ); netgen::mparam.minh = size; @@ -651,6 +716,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, netgen::Mesh& ngMesh, vector& nodeVec, const list< SMESH_subMesh* > & meshedSM, + SMESH_MesherHelper* quadHelper, SMESH_ProxyMesh::Ptr proxyMesh) { TNode2IdMap nodeNgIdMap; @@ -824,6 +890,17 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, } } // loop on geomEdge ancestors + if ( quadHelper ) // remember medium nodes of sub-meshes + { + SMDS_ElemIteratorPtr edges = smDS->GetElements(); + while ( edges->more() ) + { + const SMDS_MeshElement* e = edges->next(); + if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e ))) + break; + } + } + break; } // case TopAbs_EDGE @@ -860,8 +937,8 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false )); // Problem with the second order in a quadrangular mesh remains. - // 1) All quadrangles geberated by NETGEN are moved to an inexistent face - // by FillSMesh() (find AddFaceDescriptor) + // 1) All quadrangles generated by NETGEN are moved to an inexistent face + // by FillSMesh() (find "AddFaceDescriptor") // 2) Temporary triangles generated by StdMeshers_QuadToTriaAdaptor // are on faces where quadrangles were. // Due to these 2 points, wrong geom faces are used while conversion to qudratic @@ -914,7 +991,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, if ( helper.IsSeamShape( shapeID )) if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->getshapeId() )) inFaceNode = f->GetNodeWrap( i-1 ); - else + else inFaceNode = f->GetNodeWrap( i+1 ); gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode ); @@ -934,10 +1011,22 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, swap( tri[1], tri[2] ); ngMesh.AddSurfaceElement (tri); #ifdef DUMP_TRIANGLES - cout << tri << endl; + cout << tri << endl; #endif } } + + if ( quadHelper ) // remember medium nodes of sub-meshes + { + SMDS_ElemIteratorPtr faces = smDS->GetElements(); + while ( faces->more() ) + { + const SMDS_MeshElement* f = faces->next(); + if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshFace*>( f ))) + break; + } + } + break; } // case TopAbs_FACE @@ -1146,7 +1235,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInFaces(const netgen::OCCGeometry& o nodeVec.push_back( nV ); // get node UV - bool uvOK = false; + bool uvOK = true; vData.uv = helper.GetNodeUV( face, nV, 0, &uvOK ); if ( !uvOK ) helper.CheckNodeUV( face, nV, vData.uv, BRep_Tool::Tolerance(V),/*force=*/1); @@ -1282,8 +1371,10 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry& #ifdef DUMP_TRIANGLES_SCRIPT // create a python script making a mesh containing triangles added for internal vertices ofstream py(DUMP_TRIANGLES_SCRIPT); - py << "from smesh import * "<< endl - << "m = Mesh(name='triangles')" << endl; + py << "import SMESH"<< endl + << "from salome.smesh import smeshBuilder"< & nodeVec) + vector< const SMDS_MeshNode* > & nodeVec, + const bool overrideMinH) { // ---------------------------- // Check wires and count nodes @@ -1662,11 +1754,14 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node ); // get an average size of adjacent segments to avoid sharp change of // element size (regression on issue 0020452, note 0010898) - int iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments ); - int iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments ); - double avgH = ( segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]) / 3; - - RestrictLocalSize( ngMesh, 0.5*(np1+np2), avgH ); + int iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments ); + int iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments ); + double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]; + int nbSeg = ( int( segLen[ iPrev ] > sumH / 100.) + + int( segLen[ i ] > sumH / 100.) + + int( segLen[ iNext ] > sumH / 100.)); + if ( nbSeg > 0 ) + RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH ); } if ( isInternalWire ) { @@ -1743,6 +1838,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh& ngMesh, * \param initState - bn of entities in netgen mesh before computing * \param sMesh - SMESH mesh to fill in * \param nodeVec - vector of nodes in which node index == netgen ID + * \param comment - returns problem description + * \param quadHelper - holder of medium nodes of sub-meshes * \retval int - error */ //================================================================================ @@ -1752,7 +1849,8 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, const NETGENPlugin_ngMeshInfo& initState, SMESH_Mesh& sMesh, std::vector& nodeVec, - SMESH_Comment& comment) + SMESH_Comment& comment, + SMESH_MesherHelper* quadHelper) { int nbNod = ngMesh.GetNP(); int nbSeg = ngMesh.GetNSeg(); @@ -1761,6 +1859,13 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, SMESHDS_Mesh* meshDS = sMesh.GetMeshDS(); + // quadHelper is used for either + // 1) making quadratic elements when a lower dimention mesh is loaded + // to SMESH before convertion to quadratic by NETGEN + // 2) sewing of quadratic elements with quadratic elements of sub-meshes + if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() ) + quadHelper = 0; + // ------------------------------------- // Create and insert nodes into nodeVec // ------------------------------------- @@ -1843,7 +1948,10 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, { if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]))) continue; - edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])); + if ( quadHelper ) // final mesh must be quadratic + edge = quadHelper->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])); + else + edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])); } else { @@ -1880,6 +1988,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, // from computation of 3D mesh ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0)); + vector nodes; for (i = nbInitFac+1; i <= nbFac; ++i ) { const netgen::Element2d& elem = ngMesh.SurfaceElement(i); @@ -1887,7 +1996,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, TopoDS_Face aFace; if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent()) aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd)); - vector nodes; + nodes.clear(); for (int j=1; j <= elem.GetNP(); ++j) { int pind = elem.PNum(j); @@ -1895,7 +2004,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, break; if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind)) { - nodes.push_back(node); + nodes.push_back( node ); if (!aFace.IsNull() && node->getshapeId() < 1) { const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j); @@ -1913,17 +2022,30 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, switch (elem.GetType()) { case netgen::TRIG: - face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]); + if ( quadHelper ) // final mesh must be quadratic + face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2]); + else + face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]); break; case netgen::QUAD: - face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); + if ( quadHelper ) // final mesh must be quadratic + face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); + else + face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); // exclude qudrangle elements from computation of 3D mesh const_cast< netgen::Element2d& >( elem ).SetIndex( quadFaceID ); break; case netgen::TRIG6: + nodes[5] = mediumNode( nodes[0],nodes[1],nodes[5], quadHelper ); + nodes[3] = mediumNode( nodes[1],nodes[2],nodes[3], quadHelper ); + nodes[4] = mediumNode( nodes[2],nodes[0],nodes[4], quadHelper ); face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]); break; case netgen::QUAD8: + nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper ); + nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper ); + nodes[5] = mediumNode( nodes[2],nodes[3],nodes[5], quadHelper ); + nodes[6] = mediumNode( nodes[3],nodes[0],nodes[6], quadHelper ); face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3], nodes[4],nodes[7],nodes[5],nodes[6]); // exclude qudrangle elements from computation of 3D mesh @@ -1955,7 +2077,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, TopoDS_Solid aSolid; if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent()) aSolid = TopoDS::Solid(occgeo.somap(aSolidInd)); - vector nodes; + nodes.clear(); for (int j=1; j <= elem.GetNP(); ++j) { int pind = elem.PNum(j); @@ -1981,6 +2103,12 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]); break; case netgen::TET10: + nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper ); + nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper ); + nodes[5] = mediumNode( nodes[2],nodes[0],nodes[5], quadHelper ); + nodes[6] = mediumNode( nodes[0],nodes[3],nodes[6], quadHelper ); + nodes[8] = mediumNode( nodes[1],nodes[3],nodes[8], quadHelper ); + nodes[9] = mediumNode( nodes[2],nodes[3],nodes[9], quadHelper ); vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3], nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]); break; @@ -2084,6 +2212,39 @@ namespace return str; } + //================================================================================ + /*! + * \brief Looks for triangles lying on a SOLID + */ + //================================================================================ + + bool hasBadElemOnSolid( const list& elems, + SMESH_subMesh* solidSM ) + { + TopTools_IndexedMapOfShape solidSubs; + TopExp::MapShapes( solidSM->GetSubShape(), solidSubs ); + SMESHDS_Mesh* mesh = solidSM->GetFather()->GetMeshDS(); + + list::const_iterator e = elems.begin(); + for ( ; e != elems.end(); ++e ) + { + const SMDS_MeshElement* elem = *e; + if ( elem->GetType() != SMDSAbs_Face ) + continue; + int nbNodesOnSolid = 0; + SMDS_NodeIteratorPtr nIt = elem->nodeIterator(); + while ( nIt->more() ) + { + const SMDS_MeshNode* n = nIt->next(); + const TopoDS_Shape& s = mesh->IndexToShape( n->getshapeId() ); + nbNodesOnSolid += ( !s.IsNull() && solidSubs.Contains( s )); + if ( nbNodesOnSolid > 2 ) + return true; + } + } + return false; + } + const double edgeMeshingTime = 0.001; const double faceMeshingTime = 0.019; const double edgeFaceMeshingTime = edgeMeshingTime + faceMeshingTime; @@ -2110,10 +2271,13 @@ bool NETGENPlugin_Mesher::Compute() " growth rate = " << mparams.grading << "\n" " elements per radius = " << mparams.curvaturesafety << "\n" " second order = " << mparams.secondorder << "\n" - " quad allowed = " << mparams.quad); - //cout << " quad allowed = " << mparams.quad<GetSubMesh( face.Current() )); meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); } - FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, proxyMesh); + FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh); } } // fill _ngMesh with faces of sub-meshes - err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ])); + err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper)); initState = NETGENPlugin_ngMeshInfo(_ngMesh); //toPython( _ngMesh, "/tmp/ngPython.py"); } @@ -2523,7 +2687,7 @@ bool NETGENPlugin_Mesher::Compute() { // store computed faces in SMESH in order not to create SMESH // faces for ng faces added here - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); + FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); // add ng faces to solids with internal vertices AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, internals ); // duplicate mesh faces on internal faces @@ -2595,8 +2759,26 @@ bool NETGENPlugin_Mesher::Compute() try { OCC_CATCH_SIGNALS; + if ( !meshedSM[ MeshDim_1D ].empty() ) + { + // remove segments not attached to geometry (IPAL0052479) + for (int i = 1; i <= _ngMesh->GetNSeg(); ++i) + { + const netgen::Segment & seg = _ngMesh->LineSegment (i); + if ( seg.epgeominfo[ 0 ].edgenr == 0 ) + _ngMesh->DeleteSegment( i ); + } + _ngMesh->Compress(); + } + // convert to quadratic netgen::OCCRefinementSurfaces ref (occgeo); ref.MakeSecondOrder (*_ngMesh); + + // care of elements already loaded to SMESH + // if ( initState._nbSegments > 0 ) + // makeQuadratic( occgeo.emap, _mesh ); + // if ( initState._nbFaces > 0 ) + // makeQuadratic( occgeo.fmap, _mesh ); } catch (Standard_Failure& ex) { @@ -2627,8 +2809,14 @@ bool NETGENPlugin_Mesher::Compute() // Feed back the SMESHDS with the generated Nodes and Elements if ( true /*isOK*/ ) // get whatever built - FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); //!< + { + 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 ); + } SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec); if ( readErr && !readErr->myBadElements.empty() ) error = readErr; @@ -2686,7 +2874,14 @@ bool NETGENPlugin_Mesher::Compute() { smError.reset( new SMESH_ComputeError( *error )); if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK ) + { smError->myName = COMPERR_WARNING; + } + else if ( !smError->myBadElements.empty() ) // bad surface mesh + { + if ( !hasBadElemOnSolid( smError->myBadElements, sm )) + smError.reset(); + } } pb3D = pb3D || ( smError && smError->IsKO() ); } @@ -2968,6 +3163,8 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder, { ((int&) _progressTic ) = *algoProgressTic + 1; + if ( !_occgeom ) return 0; + double progress = -1; if ( !_isVolume ) { @@ -3030,11 +3227,14 @@ double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder, void NETGENPlugin_Mesher::RemoveTmpFiles() { - if ( SMESH_File("test.out").remove() && netgen::testout) + bool rm = SMESH_File("test.out").remove() ; +#ifndef WIN32 + if (rm && netgen::testout) { delete netgen::testout; netgen::testout = 0; } +#endif SMESH_File("problemfaces").remove(); SMESH_File("occmesh.rep").remove(); } @@ -3118,7 +3318,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh, ofstream outfile(pyFile.c_str(), ios::out); if ( !outfile ) return; - outfile << "import smesh, SMESH" << endl + outfile << "import SMESH" << endl + << "from salome.smesh import smeshBuilder" << endl + << "smesh = smeshBuilder.New(salome.myStudy)" << endl << "mesh = smesh.Mesh()" << endl << endl; using namespace netgen; @@ -3589,13 +3791,22 @@ NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper() { Ng_Init(); - _isComputeOk = false; - _outputFileName = getOutputFileName(); - netgen::mycout = new ofstream ( _outputFileName.c_str() ); - netgen::myerr = netgen::mycout; + _isComputeOk = false; + _coutBuffer = NULL; + if ( !getenv( "KEEP_NETGEN_OUTPUT" )) + { + // redirect all netgen output (mycout,myerr,cout) to _outputFileName + _outputFileName = getOutputFileName(); + netgen::mycout = new ofstream ( _outputFileName.c_str() ); + netgen::myerr = netgen::mycout; + _coutBuffer = std::cout.rdbuf(); #ifdef _DEBUG_ - cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl; + cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl; +#else + std::cout.rdbuf( netgen::mycout->rdbuf() ); #endif + } + _ngMesh = Ng_NewMesh(); } @@ -3610,6 +3821,8 @@ NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper() Ng_DeleteMesh( _ngMesh ); Ng_Exit(); NETGENPlugin_Mesher::RemoveTmpFiles(); + if ( _coutBuffer ) + std::cout.rdbuf( _coutBuffer ); #ifdef _DEBUG_ if( _isComputeOk ) #endif @@ -3641,7 +3854,11 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName() TCollection_AsciiString aGenericName = (char*)aTmpDir.c_str(); aGenericName += "NETGEN_"; +#ifndef WIN32 aGenericName += getpid(); +#else + aGenericName += _getpid(); +#endif aGenericName += "_"; aGenericName += Abs((Standard_Integer)(long) aGenericName.ToCString()); aGenericName += ".out"; @@ -3657,20 +3874,20 @@ std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName() void NETGENPlugin_NetgenLibWrapper::removeOutputFile() { - string tmpDir = SALOMEDS_Tool::GetDirFromPath( _outputFileName ); - SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames; - aFiles->length(1); - std::string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out"; - aFiles[0] = aFileName.c_str(); - if ( netgen::mycout) + if ( !_outputFileName.empty() ) { - delete netgen::mycout; - netgen::mycout = 0; - netgen::myerr = 0; + if ( netgen::mycout ) + { + delete netgen::mycout; + netgen::mycout = 0; + netgen::myerr = 0; + } + 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::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true ); } - - SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true ); -#ifdef _DEBUG_ - cout << "NOTE: netgen output log was REMOVED " << _outputFileName << endl; -#endif }