From: eap Date: Fri, 12 Mar 2010 08:36:16 +0000 (+0000) Subject: 0020682: EDF 1222 SMESH: 3D mesh from a skin mesh and with volumic cells X-Git-Tag: V5_1_4a1~8 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=1d3a5a5454bdc92417e1ce928f645503270064ff;p=plugins%2Fnetgenplugin.git 0020682: EDF 1222 SMESH: 3D mesh from a skin mesh and with volumic cells * Redesign to solve pb with internal face (StudyFiss_bugNetgen3D.hdf) --- diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index 08fc384..5345129 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -42,11 +42,12 @@ #include "SMESH_MeshEditor.hxx" #include "StdMeshers_QuadToTriaAdaptor.hxx" +#include #include #include -#include #include #include +#include #include #include @@ -150,6 +151,158 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis return isOk; } +namespace +{ + //================================================================================ + /*! + * \brief It correctly initializes netgen library at constructor and + * correctly finishes using netgen library at destructor + */ + //================================================================================ + + struct TNetgenLibWrapper + { + Ng_Mesh* _ngMesh; + TNetgenLibWrapper() + { + Ng_Init(); + _ngMesh = Ng_NewMesh(); + } + ~TNetgenLibWrapper() + { + Ng_DeleteMesh( _ngMesh ); + Ng_Exit(); + NETGENPlugin_Mesher::RemoveTmpFiles(); + } + }; + + //================================================================================ + /*! + * \brief Find mesh faces on non-internal geom faces sharing internal edge + * some nodes of which are to be doubled to make the second border of the "crack" + */ + //================================================================================ + + void findBorders( const set& internalShapeIds, + SMESH_MesherHelper& helper, + TIDSortedElemSet & borderElems, + set & borderFaceIds ) + { + SMESH_Mesh* mesh = helper.GetMesh(); + SMESHDS_Mesh* meshDS = helper.GetMeshDS(); + + // loop on internal geom edges + set::const_iterator intShapeId = internalShapeIds.begin(); + for ( ; intShapeId != internalShapeIds.end(); ++intShapeId ) + { + const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId ); + if ( s.ShapeType() != TopAbs_EDGE ) continue; + + // get internal and non-internal geom faces sharing the internal edge + int intFace = 0; + set::iterator bordFace = borderFaceIds.end(); + TopTools_ListIteratorOfListOfShape ancIt( mesh->GetAncestors( s )); + for ( ; ancIt.More(); ancIt.Next() ) + { + if ( ancIt.Value().ShapeType() != TopAbs_FACE ) continue; + int faceID = meshDS->ShapeToIndex( ancIt.Value() ); + if ( internalShapeIds.count( faceID )) + intFace = faceID; + else + bordFace = borderFaceIds.insert( faceID ).first; + } + if ( bordFace == borderFaceIds.end() || !intFace ) continue; + + // get all links of mesh faces on internal geom face sharing nodes on edge + set< SMESH_OrientedLink > links; //!< links of faces on internal geom face + TIDSortedElemSet suspectFaces; //!< mesh faces on border geom faces + SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace ); + if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue; + SMESH_subMeshIteratorPtr smIt = mesh->GetSubMesh( s )->getDependsOnIterator(true,true); + while ( smIt->more() ) + { + SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS(); + if ( !sm ) continue; + SMDS_NodeIteratorPtr nIt = sm->GetNodes(); + while ( nIt->more() ) + { + const SMDS_MeshNode* nOnEdge = nIt->next(); + SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( intFaceSM->Contains( f )) + { + int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 ); + for ( int i = 0; i < nbNodes; ++i ) + links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes))); + } + else + { + suspectFaces.insert( f ); + } + } + } + } + // having link with same orientation as mesh faces on + // the internal geom face are . + // Some of them have only one node on edge s, we collect them to decide + // later by links of found + TIDSortedElemSet posponedFaces; + set< SMESH_OrientedLink > borderLinks; + TIDSortedElemSet::iterator fIt = suspectFaces.begin(); + for ( ; fIt != suspectFaces.end(); ++fIt ) + { + const SMDS_MeshElement* f = *fIt; + bool linkFound = false, isBorder = false; + list< SMESH_OrientedLink > faceLinks; + int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 ); + for ( int i = 0; i < nbNodes; ++i ) + { + SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes)); + faceLinks.push_back( link ); + if ( !linkFound ) + { + set< SMESH_OrientedLink >::iterator foundLink = links.find( link ); + if ( foundLink != links.end() ) + { + linkFound= true; + isBorder = ( foundLink->_reversed == link._reversed ); + if ( !isBorder ) break; + } + } + } + if ( isBorder ) + { + borderElems.insert( f ); + borderLinks.insert( faceLinks.begin(), faceLinks.end() ); + } + else if ( !linkFound ) + { + posponedFaces.insert( f ); + } + } + // decide on posponedFaces + for ( fIt = posponedFaces.begin(); fIt != posponedFaces.end(); ++fIt ) + { + const SMDS_MeshElement* f = *fIt; + int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 ); + for ( int i = 0; i < nbNodes; ++i ) + { + SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes)); + set< SMESH_OrientedLink >::iterator foundLink = borderLinks.find( link ); + if ( foundLink != borderLinks.end() ) + { + if ( foundLink->_reversed != link._reversed ) + borderElems.insert( f ); + break; + } + } + } + } + } +} + //============================================================================= /*! *Here we are going to use the NETGEN mesher @@ -163,90 +316,182 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - const int invalid_ID = -1; + SMESH_MesherHelper helper(aMesh); + bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape); - SMESH::Controls::Area areaControl; - SMESH::Controls::TSequenceOfXYZ nodesCoords; + int Netgen_NbOfNodes = 0; + int Netgen_param2ndOrder = 0; + double Netgen_paramFine = 1.; + double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); - // ------------------------------------------------------------------- - // get triangles on aShell and make a map of nodes to Netgen node IDs - // ------------------------------------------------------------------- + double Netgen_point[3]; + int Netgen_triangle[3]; + int Netgen_tetrahedron[4]; - SMESH_MesherHelper helper(aMesh); - SMESH_MesherHelper* myTool = &helper; - bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + TNetgenLibWrapper ngLib; + Ng_Mesh * Netgen_mesh = ngLib._ngMesh; + // maps of 1) ordinary nodes and 2) doubled nodes on internal shapes typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap; - TNodeToIDMap nodeToNetgenID; - list< const SMDS_MeshElement* > triangles; - list< bool > isReversed; // orientation of triangles + typedef TNodeToIDMap::value_type TN2ID; + TNodeToIDMap nodeToNetgenID[2]; + + { + const int invalid_ID = -1; - TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); - bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); + SMESH::Controls::Area areaControl; + SMESH::Controls::TSequenceOfXYZ nodesCoords; - // for the degeneraged edge: ignore all but one node on it; - // map storing ids of degen edges and vertices and their netgen id: - map< int, int* > degenShapeIdToPtrNgId; - map< int, int* >::iterator shId_ngId; - list< int > degenNgIds; + // -------------------------------------------------------------------- + // Issue 0020676 (StudyFiss_bugNetgen3D.hdf). Pb with internal face. + // Find internal geom faces, edges and vertices. + // Nodes and faces built on the found internal shapes + // will be doubled in Netgen input to make two borders of the "crack". + // -------------------------------------------------------------------- - StdMeshers_QuadToTriaAdaptor Adaptor; - Adaptor.Compute(aMesh,aShape); + set internalShapeIds; + set borderFaceIds; //!< non-internal geom faces sharing internal edge - for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) - { - const TopoDS_Shape& aShapeFace = exp.Current(); - const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); - if ( aSubMeshDSFace ) + // mesh faces on , having nodes on internal edge that are + // to be replaced by doubled nodes + TIDSortedElemSet borderElems; + + // find "internal" faces and edges + TopExp_Explorer exFa(aShape,TopAbs_FACE), exEd, exVe; + for ( ; exFa.More(); exFa.Next()) + { + if ( exFa.Current().Orientation() == TopAbs_INTERNAL ) + { + internalShapeIds.insert( meshDS->ShapeToIndex( exFa.Current() )); + for ( exEd.Init( exFa.Current(), TopAbs_EDGE ); exEd.More(); exEd.Next()) + if ( helper.NbAncestors( exEd.Current(), aMesh, TopAbs_FACE ) > 1 ) + internalShapeIds.insert( meshDS->ShapeToIndex( exEd.Current() )); + } + } + if ( !internalShapeIds.empty() ) + { + // find internal vertices, + // we consider vertex internal if it is shared by more than one internal edge + TopTools_ListIteratorOfListOfShape ancIt; + set::iterator intShapeId = internalShapeIds.begin(); + for ( ; intShapeId != internalShapeIds.end(); ++intShapeId ) + { + const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId ); + if ( s.ShapeType() != TopAbs_EDGE ) continue; + for ( exVe.Init( s, TopAbs_VERTEX ); exVe.More(); exVe.Next()) + { + set internalEdges; + for ( ancIt.Initialize( aMesh.GetAncestors( exVe.Current() )); + ancIt.More(); ancIt.Next() ) + { + if ( ancIt.Value().ShapeType() != TopAbs_EDGE ) continue; + int edgeID = meshDS->ShapeToIndex( ancIt.Value() ); + if ( internalShapeIds.count( edgeID )) + internalEdges.insert( edgeID ); + } + if ( internalEdges.size() > 1 ) + internalShapeIds.insert( meshDS->ShapeToIndex( exVe.Current() )); + } + } + // find border shapes and mesh elements + findBorders( internalShapeIds, helper, borderElems, borderFaceIds ); + } + + // --------------------------------- + // Feed the Netgen with surface mesh + // --------------------------------- + + TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); + bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); + + StdMeshers_QuadToTriaAdaptor Adaptor; + if ( aMesh.NbQuadrangles() > 0 ) + Adaptor.Compute(aMesh,aShape); + + for ( exFa.ReInit(); exFa.More(); exFa.Next()) { + const TopoDS_Shape& aShapeFace = exFa.Current(); + int faceID = meshDS->ShapeToIndex( aShapeFace ); + bool isInternalFace = internalShapeIds.count( faceID ); + bool isBorderFace = borderFaceIds.count( faceID ); bool isRev = false; - if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 ) + if ( checkReverse && !isInternalFace && + helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 ) // IsReversedSubMesh() can work wrong on strongly curved faces, // so we use it as less as possible isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); + const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + if ( !aSubMeshDSFace ) continue; SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); - while ( iteratorElem->more() ) // loop on elements on a face + while ( iteratorElem->more() ) // loop on elements on a geom face { - // check element + // check mesh face const SMDS_MeshElement* elem = iteratorElem->next(); if ( !elem ) return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + vector< const SMDS_MeshElement* > trias; bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 )); - if ( !isTraingle ) { - // using adaptor + if ( !isTraingle ) + { + // use adaptor to convert quadrangle face into triangles const list* faces = Adaptor.GetTriangles(elem); - if(faces==0) { + if(faces==0) return error( COMPERR_BAD_INPUT_MESH, SMESH_Comment("No triangles in adaptor for element ")<GetID()); - } - list::const_iterator itf = faces->begin(); - for(; itf!=faces->end(); itf++ ) { - triangles.push_back( (*itf) ); - isReversed.push_back( isRev ); - // put triange's nodes to nodeToNetgenID map - SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator(); - while ( triangleNodesIt->more() ) { - const SMDS_MeshNode * node = - static_cast(triangleNodesIt->next()); - if(myTool->IsMedium(node)) - continue; - nodeToNetgenID.insert( make_pair( node, invalid_ID )); - } - } + trias.assign( faces->begin(), faces->end() ); + } + else + { + trias.push_back( elem ); } - else { - // keep a triangle - triangles.push_back( elem ); - isReversed.push_back( isRev ); - // put elem nodes to nodeToNetgenID map - SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator(); - while ( triangleNodesIt->more() ) { - const SMDS_MeshNode * node = - static_cast(triangleNodesIt->next()); - if(myTool->IsMedium(node)) - continue; - nodeToNetgenID.insert( make_pair( node, invalid_ID )); + // Add nodes of triangles and triangles them-selves to netgen mesh + + // a triangle on internal face is added twice, + // on border face, once but with doubled nodes + bool isBorder = ( isBorderFace && borderElems.count( elem )); + int nbDblLoops = ( isInternalFace && isTraingle || isBorder ) ? 2 : 1; + + for ( int i = 0; i < trias.size(); ++i ) + { + bool reverse = isRev; + for ( int isDblF = isBorder; isDblF < nbDblLoops; ++isDblF, reverse = !reverse ) + { + // add three nodes of triangle + bool hasDegen = false; + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = trias[i]->GetNode( iN ); + int shapeID = node->GetPosition()->GetShapeId(); + if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && + helper.IsDegenShape( shapeID )) + { + // ignore all nodes on degeneraged edge and use node on its vertex instead + TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value(); + node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS ); + hasDegen = true; + } + bool isDblN = isDblF && internalShapeIds.count( shapeID ); + int& ngID = nodeToNetgenID[isDblN].insert(TN2ID( node, invalid_ID )).first->second; + if ( ngID == invalid_ID ) + { + ngID = ++Netgen_NbOfNodes; + Netgen_point [ 0 ] = node->X(); + Netgen_point [ 1 ] = node->Y(); + Netgen_point [ 2 ] = node->Z(); + Ng_AddPoint(Netgen_mesh, Netgen_point); + } + Netgen_triangle[ iN ] = ngID; + } + // add triangle + if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] || + Netgen_triangle[0] == Netgen_triangle[2] || + Netgen_triangle[2] == Netgen_triangle[1] )) + break; + if ( reverse ) + swap( Netgen_triangle[1], Netgen_triangle[2] ); + + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); } } #ifdef _DEBUG_ @@ -257,92 +502,9 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, MESSAGE( "Warning: Degenerated " << elem ); } #endif - } - // look for degeneraged edges and vetices - for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next()) - { - TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() ); - if ( BRep_Tool::Degenerated( aShapeEdge )) - { - degenNgIds.push_back( invalid_ID ); - int* ptrIdOnEdge = & degenNgIds.back(); - // remember edge id - int edgeID = meshDS->ShapeToIndex( aShapeEdge ); - degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge )); - // remember vertex id - int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge )); - degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge )); - } - } - } - } - // --------------------------------- - // Feed the Netgen with surface mesh - // --------------------------------- + } // loop on elements on a face + } // loop on faces of a SOLID or SHELL - int Netgen_NbOfNodes = 0; - int Netgen_param2ndOrder = 0; - double Netgen_paramFine = 1.; - double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); - - double Netgen_point[3]; - int Netgen_triangle[3]; - int Netgen_tetrahedron[4]; - - Ng_Init(); - - Ng_Mesh * Netgen_mesh = Ng_NewMesh(); - - // set nodes and remember thier netgen IDs - bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty(); - TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); - for ( ; n_id != nodeToNetgenID.end(); ++n_id ) - { - const SMDS_MeshNode* node = n_id->first; - - // ignore nodes on degenerated edge - if ( hasDegen ) { - int shapeId = node->GetPosition()->GetShapeId(); - shId_ngId = degenShapeIdToPtrNgId.find( shapeId ); - isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() ); - if ( isDegen && *(shId_ngId->second) != invalid_ID ) { - n_id->second = *(shId_ngId->second); - continue; - } - } - Netgen_point [ 0 ] = node->X(); - Netgen_point [ 1 ] = node->Y(); - Netgen_point [ 2 ] = node->Z(); - Ng_AddPoint(Netgen_mesh, Netgen_point); - n_id->second = ++Netgen_NbOfNodes; // set netgen ID - - if ( isDegen ) // all nodes on a degen edge get one netgen ID - *(shId_ngId->second) = n_id->second; - } - - // set triangles - list< const SMDS_MeshElement* >::iterator tria = triangles.begin(); - list< bool >::iterator reverse = isReversed.begin(); - for ( ; tria != triangles.end(); ++tria, ++reverse ) - { - int i = 0; - SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator(); - while ( triangleNodesIt->more() ) { - const SMDS_MeshNode * node = - static_cast(triangleNodesIt->next()); - if(myTool->IsMedium(node)) - continue; - Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ]; - ++i; - } - if ( !hasDegen || - // ignore degenerated triangles, they have 2 or 3 same ids - (Netgen_triangle[0] != Netgen_triangle[1] && - Netgen_triangle[0] != Netgen_triangle[2] && - Netgen_triangle[2] != Netgen_triangle[1] )) - { - Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); - } } // ------------------------- @@ -352,8 +514,8 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, Ng_Meshing_Parameters Netgen_param; Netgen_param.secondorder = Netgen_param2ndOrder; - Netgen_param.fineness = Netgen_paramFine; - Netgen_param.maxh = Netgen_paramSize; + Netgen_param.fineness = Netgen_paramFine; + Netgen_param.maxh = Netgen_paramSize; Ng_Result status; @@ -393,15 +555,26 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, // Feed back the SMESHDS with the generated Nodes and Volume Elements // ------------------------------------------------------------------- + // vector of nodes in which node index == netgen ID + vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 ); + // insert old nodes into nodeVec + for ( int isDbl = 0; isDbl < 2; ++isDbl ) + { + TNodeToIDMap::iterator n_id = nodeToNetgenID[isDbl].begin(); + for ( ; n_id != nodeToNetgenID[isDbl].end(); ++n_id ) + nodeVec[ n_id->second ] = n_id->first; + nodeToNetgenID[isDbl].clear(); + } + if ( status == NG_VOLUME_FAILURE ) + { + SMESH_ComputeErrorPtr err = NETGENPlugin_Mesher::readErrors(nodeVec); + if ( err && !err->myBadElements.empty() ) + error( err ); + } + bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built if ( isOK ) { - // vector of nodes in which node index == netgen ID - vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 ); - // insert old nodes into nodeVec - for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) { - nodeVec.at( n_id->second ) = n_id->first; - } // create and insert new nodes into nodeVec int nodeIndex = Netgen_NbOfNodes + 1; int shapeID = meshDS->ShapeToIndex( aShape ); @@ -419,23 +592,24 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex ) { Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron); - SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), - nodeVec.at( Netgen_tetrahedron[1] ), - nodeVec.at( Netgen_tetrahedron[2] ), - nodeVec.at( Netgen_tetrahedron[3] )); + SMDS_MeshVolume * elt = helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); meshDS->SetMeshElementOnShape(elt, shapeID ); } } - Ng_DeleteMesh(Netgen_mesh); - Ng_Exit(); - - NETGENPlugin_Mesher::RemoveTmpFiles(); - return (status == NG_OK); } -bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, +//================================================================================ +/*! + * \brief Compute tetrahedral mesh from 2D mesh without geometry + */ +//================================================================================ + +bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, SMESH_MesherHelper* aHelper) { MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume); @@ -462,42 +636,33 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, while( fIt->more()) sortedFaces.insert( fIt->next() ); TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end(); - for ( ; itFace != fEnd; ++itFace ) { + for ( ; itFace != fEnd; ++itFace ) + { // check element const SMDS_MeshElement* elem = *itFace; if ( !elem ) return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + + vector< const SMDS_MeshElement* > trias; bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 )); if ( !isTraingle ) { // using adaptor const list* faces = Adaptor.GetTriangles(elem); - if(faces==0) { - continue; // Issue 0020682. There already can be 3d mesh - } - list::const_iterator itf = faces->begin(); - for(; itf!=faces->end(); itf++ ) { - triangles.push_back( (*itf) ); - // put triange's nodes to nodeToNetgenID map - SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator(); - while ( triangleNodesIt->more() ) { - const SMDS_MeshNode * node = - static_cast(triangleNodesIt->next()); - if(aHelper->IsMedium(node)) - continue; - nodeToNetgenID.insert( make_pair( node, invalid_ID )); - } - } + if(faces==0) + return error( COMPERR_BAD_INPUT_MESH, + SMESH_Comment("No triangles in adaptor for element ")<GetID()); + trias.assign( faces->begin(), faces->end() ); } else { - // keep a triangle - triangles.push_back( elem ); - // put elem nodes to nodeToNetgenID map - SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator(); - while ( triangleNodesIt->more() ) { - const SMDS_MeshNode * node = - static_cast(triangleNodesIt->next()); - if(aHelper->IsMedium(node)) - continue; + trias.push_back( elem ); + } + for ( int i = 0; i < trias.size(); ++i ) + { + triangles.push_back( trias[i] ); + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = trias[i]->GetNode( iN ); + // put elem nodes to nodeToNetgenID map nodeToNetgenID.insert( make_pair( node, invalid_ID )); } } @@ -516,11 +681,10 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, int Netgen_triangle[3]; int Netgen_tetrahedron[4]; - Ng_Init(); - - Ng_Mesh * Netgen_mesh = Ng_NewMesh(); + TNetgenLibWrapper ngLib; + Ng_Mesh * Netgen_mesh = ngLib._ngMesh; - // set nodes and remember thier netgen IDs + // set nodes and remember thier netgen IDs TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); for ( ; n_id != nodeToNetgenID.end(); ++n_id ) @@ -625,11 +789,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, } } - Ng_DeleteMesh(Netgen_mesh); - Ng_Exit(); - - NETGENPlugin_Mesher::RemoveTmpFiles(); - return (status == NG_OK); } @@ -646,8 +805,8 @@ bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh, { int nbtri = 0, nbqua = 0; double fullArea = 0.0; - for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) { - TopoDS_Face F = TopoDS::Face( exp.Current() ); + for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) { + TopoDS_Face F = TopoDS::Face( expF.Current() ); SMESH_subMesh *sm = aMesh.GetSubMesh(F); MapShapeNbElemsItr anIt = aResMap.find(sm); if( anIt==aResMap.end() ) { @@ -669,12 +828,12 @@ bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh, bool IsQuadratic = false; bool IsFirst = true; TopTools_MapOfShape tmpMap; - for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) { - TopoDS_Edge E = TopoDS::Edge(exp.Current()); + for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) { + TopoDS_Edge E = TopoDS::Edge(expF.Current()); if( tmpMap.Contains(E) ) continue; tmpMap.Add(E); - SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current()); + SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current()); MapShapeNbElemsItr anIt = aResMap.find(aSubMesh); if( anIt==aResMap.end() ) { SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();