X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FNETGENPlugin%2FNETGENPlugin_NETGEN_2D_ONLY.cxx;h=9c122c3f236351ad432cb26abcc3e745d1419b2b;hb=5bf2f4cc8d6f7bfe308be651a6f8b3c4cd7e5852;hp=08b09509c371f86997f576abb946e9883054a7bf;hpb=fbceed67d3744002f1e438129b3b3183dd6ac861;p=plugins%2Fnetgenplugin.git diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index 08b0950..9c122c3 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2022 CEA/DEN, EDF R&D, OPEN CASCADE // // 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 @@ -63,34 +63,31 @@ namespace nglib { #include //#include namespace netgen { + NETGENPLUGIN_DLL_HEADER + extern MeshingParameters mparam; #ifdef NETGEN_V5 - extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int); -#else - extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*); + extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh); #endif - extern MeshingParameters mparam; } using namespace std; using namespace netgen; using namespace nglib; -//#define DUMP_SEGMENTS - //============================================================================= /*! - * + * */ //============================================================================= -NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, int studyId, +NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, SMESH_Gen* gen) - : SMESH_2D_Algo(hypId, studyId, gen) + : SMESH_2D_Algo(hypId, gen) { - MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY"); _name = "NETGEN_2D_ONLY"; - + _shapeType = (1 << TopAbs_FACE);// 1 bit /shape type + _onlyUnaryInput = false; // treat all FACEs at once _compatibleHypothesis.push_back("MaxElementArea"); _compatibleHypothesis.push_back("LengthFromEdges"); @@ -98,26 +95,26 @@ NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, int studyId, _compatibleHypothesis.push_back("NETGEN_Parameters_2D"); _compatibleHypothesis.push_back("ViscousLayers2D"); - _hypMaxElementArea = 0; - _hypLengthFromEdges = 0; + _hypMaxElementArea = 0; + _hypLengthFromEdges = 0; _hypQuadranglePreference = 0; - _hypParameters = 0; + _hypParameters = 0; } //============================================================================= /*! - * + * */ //============================================================================= NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY() { - MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY"); + //MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY"); } //============================================================================= /*! - * + * */ //============================================================================= @@ -128,6 +125,8 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, _hypMaxElementArea = 0; _hypLengthFromEdges = 0; _hypQuadranglePreference = 0; + _hypParameters = 0; + _progressByTic = -1; const list& hyps = GetUsedHypothesis(aMesh, aShape, false); @@ -139,6 +138,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, aStatus = HYP_MISSING; + bool hasVL = false; list::const_iterator ith; for (ith = hyps.begin(); ith != hyps.end(); ++ith ) { @@ -155,7 +155,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, else if ( hypName == "NETGEN_Parameters_2D" ) _hypParameters = static_cast(hyp); else if ( hypName == StdMeshers_ViscousLayers2D::GetHypType() ) - continue; + hasVL = true; else { aStatus = HYP_INCOMPATIBLE; return false; @@ -164,13 +164,59 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, int nbHyps = bool(_hypMaxElementArea) + bool(_hypLengthFromEdges) + bool(_hypParameters ); if ( nbHyps > 1 ) - aStatus = HYP_CONCURENT; + aStatus = HYP_CONCURRENT; + else if ( hasVL ) + error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus )); else aStatus = HYP_OK; + if ( aStatus == HYP_OK && _hypParameters && _hypQuadranglePreference ) + { + aStatus = HYP_INCOMPAT_HYPS; + return error(SMESH_Comment("\"") << _hypQuadranglePreference->GetName() + << "\" and \"" << _hypParameters->GetName() + << "\" are incompatible hypotheses"); + } + return ( aStatus == HYP_OK ); } +// namespace +// { +// void limitSize( netgen::Mesh* ngMesh, +// const double maxh ) +// { +// // get bnd box +// netgen::Point3d pmin, pmax; +// ngMesh->GetBox( pmin, pmax, 0 ); +// const double dx = pmax.X() - pmin.X(); +// const double dy = pmax.Y() - pmin.Y(); +// const double dz = pmax.Z() - pmin.Z(); + +// const int nbX = Max( 2, int( dx / maxh * 3 )); +// const int nbY = Max( 2, int( dy / maxh * 3 )); +// const int nbZ = Max( 2, int( dz / maxh * 3 )); + +// if ( ! & ngMesh->LocalHFunction() ) +// ngMesh->SetLocalH( pmin, pmax, 0.1 ); + +// netgen::Point3d p; +// for ( int i = 0; i <= nbX; ++i ) +// { +// p.X() = pmin.X() + i * dx / nbX; +// for ( int j = 0; j <= nbY; ++j ) +// { +// p.Y() = pmin.Y() + j * dy / nbY; +// for ( int k = 0; k <= nbZ; ++k ) +// { +// p.Z() = pmin.Z() + k * dz / nbZ; +// ngMesh->RestrictLocalH( p, maxh ); +// } +// } +// } +// } +// } + //============================================================================= /*! *Here we are going to use the NETGEN mesher @@ -180,203 +226,416 @@ bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh& aMesh, bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { -#ifdef WITH_SMESH_CANCEL_COMPUTE netgen::multithread.terminate = 0; -#endif - MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::Compute()"); + //netgen::multithread.task = "Surface meshing"; SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - int faceID = meshDS->ShapeToIndex( aShape ); - SMESH_MesherHelper helper(aMesh); - _quadraticMesh = helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); - const bool ignoreMediumNodes = _quadraticMesh; - - // build viscous layers if required - const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD )); - SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); - if ( !proxyMesh ) - return false; - // ------------------------ - // get all edges of a face - // ------------------------ - TError problem; - TSideVector wires = - StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem, proxyMesh ); - if ( problem && !problem->IsOK() ) - return error( problem ); - int nbWires = wires.size(); - if ( nbWires == 0 ) - return error( "Problem in StdMeshers_FaceSide::GetFaceWires()"); - if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments - return error(COMPERR_BAD_INPUT_MESH, - SMESH_Comment("Too few segments: ")<NbSegments()); - - // -------------------- - // compute edge length - // -------------------- + NETGENPlugin_NetgenLibWrapper ngLib; + ngLib._isComputeOk = false; + + netgen::Mesh ngMeshNoLocSize; + netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh, & ngMeshNoLocSize }; + netgen::OCCGeometry occgeoComm; + + // min / max sizes are set as follows: + // if ( _hypParameters ) + // min and max are defined by the user + // else if ( _hypLengthFromEdges ) + // min = aMesher.GetDefaultMinSize() + // max = average segment len of a FACE + // else if ( _hypMaxElementArea ) + // min = aMesher.GetDefaultMinSize() + // max = f( _hypMaxElementArea ) + // else + // min = aMesher.GetDefaultMinSize() + // max = max segment len of a FACE NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false); - netgen::OCCGeometry occgeo; - aMesher.PrepareOCCgeometry( occgeo, F, aMesh ); - occgeo.fmap.Clear(); // face can be reversed, which is wrong in this case (issue 19978) - occgeo.fmap.Add( F ); - - if ( _hypParameters ) + aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam + const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true; + if ( _hypMaxElementArea ) { - aMesher.SetParameters(_hypParameters); + netgen::mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) ); } - else + if ( _hypQuadranglePreference ) + netgen::mparam.quad = true; + + // local size is common for all FACEs in aShape? + const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && netgen::mparam.uselocalh ); + const bool isDefaultHyp = ( !_hypLengthFromEdges && !_hypMaxElementArea && !_hypParameters ); + + if ( isCommonLocalSize ) // compute common local size in ngMeshes[0] { - double edgeLength = 0; - if (_hypLengthFromEdges || (!_hypLengthFromEdges && !_hypMaxElementArea)) + //list< SMESH_subMesh* > meshedSM[4]; --> all sub-shapes are added to occgeoComm + aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM ); + + // local size set at MESHCONST_ANALYSE step depends on + // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user + if ( !_hypParameters || netgen::mparam.minh < DBL_MIN ) { - int nbSegments = 0; - for ( int iW = 0; iW < nbWires; ++iW ) + if ( !_hypParameters ) + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.; + netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh ); + } + // set local size depending on curvature and NOT closeness of EDGEs +#ifdef NETGEN_V6 + const double factor = 2; //netgen::occparam.resthcloseedgefac; +#else + const double factor = netgen::occparam.resthcloseedgefac; + netgen::occparam.resthcloseedgeenable = false; + netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading; +#endif + occgeoComm.face_maxh = netgen::mparam.maxh; +#ifdef NETGEN_V6 + netgen::OCCParameters occparam; + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0], netgen::mparam, occparam ); +#else + netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes[0] ); +#endif + occgeoComm.emap.Clear(); + occgeoComm.vmap.Clear(); + + // set local size according to size of existing segments + TopTools_IndexedMapOfShape edgeMap; + TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap ); + for ( int iE = 1; iE <= edgeMap.Extent(); ++iE ) + { + const TopoDS_Shape& edge = edgeMap( iE ); + if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge ))) + continue; + SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge ); + if ( !smDS ) continue; + SMDS_ElemIteratorPtr segIt = smDS->GetElements(); + while ( segIt->more() ) { - edgeLength += wires[ iW ]->Length(); - nbSegments += wires[ iW ]->NbSegments(); + const SMDS_MeshElement* seg = segIt->next(); + SMESH_TNodeXYZ n1 = seg->GetNode(0); + SMESH_TNodeXYZ n2 = seg->GetNode(1); + gp_XYZ p = 0.5 * ( n1 + n2 ); + netgen::Point3d pi(p.X(), p.Y(), p.Z()); + ngMeshes[0]->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() ); } - if ( nbSegments ) - edgeLength /= nbSegments; } - if ( _hypMaxElementArea ) + + // set local size defined on shapes + aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); + try { + ngMeshes[0]->LoadLocalMeshSize( mparam.meshsizefilename ); + } catch (NgException & ex) { + return error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } + netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization + + // ================== + // Loop on all FACEs + // ================== + + vector< const SMDS_MeshNode* > nodeVec; + + TopExp_Explorer fExp( aShape, TopAbs_FACE ); + for ( int iF = 0; fExp.More(); fExp.Next(), ++iF ) + { + TopoDS_Face F = TopoDS::Face( fExp.Current() /*.Oriented( TopAbs_FORWARD )*/); + int faceID = meshDS->ShapeToIndex( F ); + SMESH_ComputeErrorPtr& faceErr = aMesh.GetSubMesh( F )->GetComputeError(); + + _quadraticMesh = helper.IsQuadraticSubMesh( F ); + const bool ignoreMediumNodes = _quadraticMesh; + + // build viscous layers if required + if ( F.Orientation() != TopAbs_FORWARD && + F.Orientation() != TopAbs_REVERSED ) + F.Orientation( TopAbs_FORWARD ); // avoid pb with TopAbs_INTERNAL + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !proxyMesh ) + continue; + + // ------------------------ + // get all EDGEs of a FACE + // ------------------------ + TSideVector wires = + StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, faceErr, &helper, proxyMesh ); + if ( faceErr && !faceErr->IsOK() ) + continue; + size_t nbWires = wires.size(); + if ( nbWires == 0 ) { - double maxArea = _hypMaxElementArea->GetMaxArea(); - edgeLength = sqrt(2. * maxArea/sqrt(3.0)); + faceErr.reset + ( new SMESH_ComputeError + ( COMPERR_ALGO_FAILED, "Problem in StdMeshers_FaceSide::GetFaceWires()" )); + continue; + } + if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments + { + faceErr.reset + ( new SMESH_ComputeError + ( COMPERR_BAD_INPUT_MESH, SMESH_Comment("Too few segments: ")<NbSegments()) ); + continue; } - if ( edgeLength < DBL_MIN ) - edgeLength = occgeo.GetBoundingBox().Diam(); - netgen::mparam.maxh = edgeLength; - netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh ); - netgen::mparam.quad = _hypQuadranglePreference ? 1 : 0; - netgen::mparam.grading = 0.7; // very coarse mesh by default - } - occgeo.face_maxh = netgen::mparam.maxh; + // ---------------------- + // compute maxh of a FACE + // ---------------------- - // ------------------------- - // Make input netgen mesh - // ------------------------- + if ( !_hypParameters ) + { + double edgeLength = 0; + if (_hypLengthFromEdges ) + { + // compute edgeLength as an average segment length + smIdType nbSegments = 0; + for ( size_t iW = 0; iW < nbWires; ++iW ) + { + edgeLength += wires[ iW ]->Length(); + nbSegments += wires[ iW ]->NbSegments(); + } + if ( nbSegments ) + edgeLength /= double( nbSegments ); + netgen::mparam.maxh = edgeLength; + } + else if ( isDefaultHyp ) + { + // set edgeLength by a longest segment + double maxSeg2 = 0; + for ( size_t iW = 0; iW < nbWires; ++iW ) + { + const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct(); + if ( points.empty() ) + return error( COMPERR_BAD_INPUT_MESH ); + gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node ); + for ( size_t i = 1; i < points.size(); ++i ) + { + gp_Pnt p = SMESH_TNodeXYZ( points[i].node ); + maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev )); + pPrev = p; + } + } + edgeLength = sqrt( maxSeg2 ) * 1.05; + netgen::mparam.maxh = edgeLength; + } + if ( netgen::mparam.maxh < DBL_MIN ) + netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam(); - NETGENPlugin_NetgenLibWrapper ngLib; - netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh; + if ( !isCommonLocalSize ) + { + netgen::mparam.minh = aMesher.GetDefaultMinSize( F, netgen::mparam.maxh ); + } + } - Box<3> bb = occgeo.GetBoundingBox(); - bb.Increase (bb.Diam()/10); - ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading); - ngMesh->SetGlobalH (netgen::mparam.maxh); + // prepare occgeom + netgen::OCCGeometry occgeom; + occgeom.shape = F; + occgeom.fmap.Add( F ); + occgeom.CalcBoundingBox(); + occgeom.facemeshstatus.SetSize(1); + occgeom.facemeshstatus = 0; + occgeom.face_maxh_modified.SetSize(1); + occgeom.face_maxh_modified = 0; + occgeom.face_maxh.SetSize(1); + occgeom.face_maxh = netgen::mparam.maxh; + + // ------------------------- + // Fill netgen mesh + // ------------------------- + + // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt + // w/o MESHCONST_ANALYSE at the second loop + int err = 0; + enum { LOC_SIZE, NO_LOC_SIZE }; + int iLoop = isCommonLocalSize ? 0 : 1; + for ( ; iLoop < 2; iLoop++ ) + { + //bool isMESHCONST_ANALYSE = false; + InitComputeError(); - vector< const SMDS_MeshNode* > nodeVec; - problem = aMesher.AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec ); - if ( problem && !problem->IsOK() ) - return error( problem ); + netgen::Mesh * ngMesh = ngMeshes[ iLoop ]; + ngMesh->DeleteMesh(); - // ------------------------- - // Generate surface mesh - // ------------------------- + if ( iLoop == NO_LOC_SIZE ) + { + ngMesh->SetGlobalH ( mparam.maxh ); + ngMesh->SetMinimalH( mparam.minh ); + Box<3> bb = occgeom.GetBoundingBox(); + bb.Increase (bb.Diam()/10); + ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); + aMesher.SetLocalSize( occgeom, *ngMesh ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh ); + try { + ngMesh->LoadLocalMeshSize( mparam.meshsizefilename ); + } catch (NgException & ex) { + return error( COMPERR_BAD_PARMETERS, ex.What() ); + } + } -#ifndef NETGEN_V5 - char *optstr = 0; -#endif - int startWith = MESHCONST_MESHSURFACE; - int endWith = MESHCONST_OPTSURFACE; - int err = 1; + nodeVec.clear(); + faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec, + /*overrideMinH=*/!_hypParameters); + if ( faceErr && !faceErr->IsOK() ) + break; - try { -#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100 - OCC_CATCH_SIGNALS; -#endif -#ifdef NETGEN_V5 - err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith); -#else - err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr); -#endif -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(netgen::multithread.terminate) - return false; -#endif - if ( err ) - 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(); - error(str); - } - catch (...) { - SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); - str << " at " << netgen::multithread.task; - error(str); - } + //if ( !isCommonLocalSize ) + //limitSize( ngMesh, mparam.maxh * 0.8); - // ---------------------------------------------------- - // Fill the SMESHDS with the generated nodes and faces - // ---------------------------------------------------- + // ------------------------- + // Generate surface mesh + // ------------------------- - int nbNodes = ngMesh->GetNP(); - int nbFaces = ngMesh->GetNSE(); + const int startWith = MESHCONST_MESHSURFACE; + const int endWith = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE; - int nbInputNodes = nodeVec.size()-1; - nodeVec.resize( nbNodes+1, 0 ); + SMESH_Comment str; + try { + OCC_CATCH_SIGNALS; - // add nodes - for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) - { - const MeshPoint& ngPoint = ngMesh->Point( ngID ); - SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); - nodeVec[ ngID ] = node; - } + err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh); - // create faces - bool reverse = ( aShape.Orientation() == TopAbs_REVERSED ); - int i,j; - for ( i = 1; i <= nbFaces ; ++i ) - { - const Element2d& elem = ngMesh->SurfaceElement(i); - vector nodes( elem.GetNP() ); - for (j=1; j <= elem.GetNP(); ++j) - { - int pind = elem.PNum(j); - if ( pind < 1 ) - break; - const SMDS_MeshNode* node = nodeVec[ pind ]; - if ( reverse ) - nodes[ nodes.size()-j ] = node; - else - nodes[ j-1 ] = node; - if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) + if ( netgen::multithread.terminate ) + return false; + if ( err ) + str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task; + } + catch (Standard_Failure& ex) { - const PointGeomInfo& pgi = elem.GeomInfoPi(j); - meshDS->SetNodeOnFace((SMDS_MeshNode*)node, faceID, pgi.u, pgi.v); + err = 1; + str << "Exception in netgen::OCCGenerateMesh()" + << " at " << netgen::multithread.task + << ": " << ex.DynamicType()->Name(); + if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) + str << ": " << ex.GetMessageString(); + } + catch (...) { + err = 1; + str << "Exception in netgen::OCCGenerateMesh()" + << " at " << netgen::multithread.task; + } + if ( err ) + { + if ( aMesher.FixFaceMesh( occgeom, *ngMesh, 1 )) + break; + if ( iLoop == LOC_SIZE ) + { + netgen::mparam.minh = netgen::mparam.maxh; + netgen::mparam.maxh = 0; + for ( size_t iW = 0; iW < wires.size(); ++iW ) + { + StdMeshers_FaceSidePtr wire = wires[ iW ]; + const vector& uvPtVec = wire->GetUVPtStruct(); + for ( size_t iP = 1; iP < uvPtVec.size(); ++iP ) + { + SMESH_TNodeXYZ p( uvPtVec[ iP ].node ); + netgen::Point3d np( p.X(),p.Y(),p.Z()); + double segLen = p.Distance( uvPtVec[ iP-1 ].node ); + double size = ngMesh->GetH( np ); + netgen::mparam.minh = Min( netgen::mparam.minh, size ); + netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen ); + } + } + //cerr << "min " << netgen::mparam.minh << " max " << netgen::mparam.maxh << endl; + netgen::mparam.minh *= 0.9; + netgen::mparam.maxh *= 1.1; + continue; + } + else + { + faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str )); + } + } + + + // ---------------------------------------------------- + // Fill the SMESHDS with the generated nodes and faces + // ---------------------------------------------------- + + int nbNodes = ngMesh->GetNP(); + int nbFaces = ngMesh->GetNSE(); + + int nbInputNodes = (int) nodeVec.size()-1; + nodeVec.resize( nbNodes+1, 0 ); + + // add nodes + for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID ) + { + const MeshPoint& ngPoint = ngMesh->Point( ngID ); + SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2)); + nodeVec[ ngID ] = node; + } + + // create faces + int i,j; + vector nodes; + for ( i = 1; i <= nbFaces ; ++i ) + { + const Element2d& elem = ngMesh->SurfaceElement(i); + nodes.resize( elem.GetNP() ); + for (j=1; j <= elem.GetNP(); ++j) + { + int pind = elem.PNum(j); + if ( pind < 1 ) + break; + nodes[ j-1 ] = nodeVec[ pind ]; + if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) + { + const PointGeomInfo& pgi = elem.GeomInfoPi(j); + meshDS->SetNodeOnFace( nodes[ j-1 ], faceID, pgi.u, pgi.v); + } + } + if ( j > elem.GetNP() ) + { + if ( elem.GetType() == TRIG ) + helper.AddFace(nodes[0],nodes[1],nodes[2]); + else + helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); + } } - } - if ( j > elem.GetNP() ) - { - SMDS_MeshFace* face = 0; - if ( elem.GetType() == TRIG ) - face = helper.AddFace(nodes[0],nodes[1],nodes[2]); - else - face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]); - } - } - return !err; + break; + } // two attempts + } // loop on FACEs + + return true; } -#ifdef WITH_SMESH_CANCEL_COMPUTE void NETGENPlugin_NETGEN_2D_ONLY::CancelCompute() { SMESH_Algo::CancelCompute(); netgen::multithread.terminate = 1; } -#endif + +//================================================================================ +/*! + * \brief Return progress of Compute() [0.,1] + */ +//================================================================================ + +double NETGENPlugin_NETGEN_2D_ONLY::GetProgress() const +{ + return -1; + // const char* task1 = "Surface meshing"; + // //const char* task2 = "Optimizing surface"; + // double& progress = const_cast( this )->_progress; + // if ( _progressByTic < 0. && + // strncmp( netgen::multithread.task, task1, 3 ) == 0 ) + // { + // progress = Min( 0.25, SMESH_Algo::GetProgressByTic() ); // [0, 0.25] + // } + // else //if ( strncmp( netgen::multithread.task, task2, 3 ) == 0) + // { + // if ( _progressByTic < 0 ) + // { + // NETGENPlugin_NETGEN_2D_ONLY* me = (NETGENPlugin_NETGEN_2D_ONLY*) this; + // me->_progressByTic = 0.25 / (_progressTic+1); + // } + // const_cast( this )->_progressTic++; + // progress = Max( progress, _progressByTic * _progressTic ); + // } + // //cout << netgen::multithread.task << " " << _progressTic << endl; + // return Min( progress, 0.99 ); +} //============================================================================= /*! @@ -393,7 +652,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, return false; // collect info from edges - int nb0d = 0, nb1d = 0; + smIdType nb0d = 0, nb1d = 0; bool IsQuadratic = false; bool IsFirst = true; double fullLen = 0.0; @@ -411,9 +670,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this)); return false; } - std::vector aVec = (*anIt).second; + std::vector aVec = (*anIt).second; nb0d += aVec[SMDSEntity_Node]; - nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); + nb1d += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); double aLen = SMESH_Algo::EdgeLength(E); fullLen += aLen; if(IsFirst) { @@ -425,9 +684,9 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, // compute edge length double ELen = 0; - if (_hypLengthFromEdges || !_hypLengthFromEdges && !_hypMaxElementArea) { + if (( _hypLengthFromEdges ) || ( !_hypLengthFromEdges && !_hypMaxElementArea )) { if ( nb1d > 0 ) - ELen = fullLen / nb1d; + ELen = fullLen / double( nb1d ); } if ( _hypMaxElementArea ) { double maxArea = _hypMaxElementArea->GetMaxArea(); @@ -445,10 +704,10 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh, smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this)); return false; } - int nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) ); - int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 ); - std::vector aVec(SMDSEntity_Last); - for(int i=SMDSEntity_Node; i aVec(SMDSEntity_Last); + for(smIdType i=SMDSEntity_Node; i