X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FNETGENPlugin%2FNETGENPlugin_NETGEN_3D.cxx;h=c3fa281b91dcbe64a3a0c8388abcad7202a48633;hb=229f703ec60fcc5c0b5e97956e540193273d80e9;hp=bb34c62e9570824a4c3b876ab7bfb8c81b447c2c;hpb=14f4b7bcf8ce41078b583005dbeae402dc6a8584;p=plugins%2Fnetgenplugin.git diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx index bb34c62..c3fa281 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx @@ -1,30 +1,63 @@ +// Copyright (C) 2007-2013 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 +// +// 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. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + //============================================================================= // File : NETGENPlugin_NETGEN_3D.cxx // Moved here from SMESH_NETGEN_3D.cxx // Created : lundi 27 Janvier 2003 // Author : Nadir BOUHAMOU (CEA) // Project : SALOME -// Copyright : CEA 2003 -// $Header$ //============================================================================= -using namespace std; - +// #include "NETGENPlugin_NETGEN_3D.hxx" -#include "SMESH_Gen.hxx" -#include "SMESH_Mesh.hxx" -#include "SMESH_ControlsDef.hxx" -#include "SMESHDS_Mesh.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "StdMeshers_Helper.hxx" - +#include "NETGENPlugin_Hypothesis.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include #include +#include #include #include +#include #include -#include "utilities.h" +#include +#include + +#include #include #include @@ -34,7 +67,24 @@ using namespace std; Netgen include files */ -#include "nglib.h" +#ifndef OCCGEOMETRY +#define OCCGEOMETRY +#endif +#include +namespace nglib { +#include +} +namespace netgen { +#ifdef NETGEN_V5 + extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int); +#else + extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*); +#endif + extern MeshingParameters mparam; + extern volatile multithreadt multithread; +} +using namespace nglib; +using namespace std; //============================================================================= /*! @@ -43,17 +93,23 @@ using namespace std; //============================================================================= NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId, - SMESH_Gen* gen) + SMESH_Gen* gen) : SMESH_3D_Algo(hypId, studyId, gen) { MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D"); _name = "NETGEN_3D"; _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type _compatibleHypothesis.push_back("MaxElementVolume"); + _compatibleHypothesis.push_back("NETGEN_Parameters"); + _compatibleHypothesis.push_back("ViscousLayers"); _maxElementVolume = 0.; _hypMaxElementVolume = NULL; + _hypParameters = NULL; + _viscousLayersHyp = NULL; + + _requireShape = false; // can work without shape } //============================================================================= @@ -73,47 +129,56 @@ NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D() */ //============================================================================= -bool NETGENPlugin_NETGEN_3D::CheckHypothesis - (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - SMESH_Hypothesis::Hypothesis_Status& aStatus) +bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus) { MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis"); _hypMaxElementVolume = NULL; + _hypParameters = NULL; + _viscousLayersHyp = NULL; _maxElementVolume = DBL_MAX; + // for correct work of GetProgress(): + netgen::multithread.percent = 0.; + netgen::multithread.task = "Volume meshing"; + _optimizationStarted = false; + list::const_iterator itl; const SMESHDS_Hypothesis* theHyp; - const list& hyps = GetUsedHypothesis(aMesh, aShape); - int nbHyp = hyps.size(); - if (!nbHyp) + const list& hyps = + GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false); + list ::const_iterator h = hyps.begin(); + if ( h == hyps.end()) { aStatus = SMESH_Hypothesis::HYP_OK; - //aStatus = SMESH_Hypothesis::HYP_MISSING; return true; // can work with no hypothesis } - itl = hyps.begin(); - theHyp = (*itl); // use only the first hypothesis - - string hypName = theHyp->GetName(); - - bool isOk = false; - - if (hypName == "MaxElementVolume") + aStatus = HYP_OK; + for ( ; h != hyps.end(); ++h ) { - _hypMaxElementVolume = static_cast (theHyp); - ASSERT(_hypMaxElementVolume); - _maxElementVolume = _hypMaxElementVolume->GetMaxVolume(); - isOk =true; - aStatus = SMESH_Hypothesis::HYP_OK; + if ( !_hypMaxElementVolume ) + _hypMaxElementVolume = dynamic_cast< const StdMeshers_MaxElementVolume*> ( *h ); + if ( !_viscousLayersHyp ) + _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h ); + if ( ! _hypParameters ) + _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis*> ( *h ); + + if ( *h != _hypMaxElementVolume && + *h != _viscousLayersHyp && + *h != _hypParameters) + aStatus = HYP_INCOMPATIBLE; } - else - aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE; + if ( _hypMaxElementVolume && _hypParameters ) + aStatus = HYP_INCOMPATIBLE; + + if ( _hypMaxElementVolume ) + _maxElementVolume = _hypMaxElementVolume->GetMaxVolume(); - return isOk; + return aStatus == HYP_OK; } //============================================================================= @@ -125,183 +190,248 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) { - MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume); - + netgen::multithread.terminate = 0; + SMESHDS_Mesh* meshDS = aMesh.GetMeshDS(); - const int invalid_ID = -1; + SMESH_MesherHelper helper(aMesh); + bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + helper.SetElementsOnShape( true ); - SMESH::Controls::Area areaControl; - SMESH::Controls::TSequenceOfXYZ nodesCoords; + int Netgen_NbOfNodes = 0; + double Netgen_point[3]; + int Netgen_triangle[3]; - // ------------------------------------------------------------------- - // get triangles on aShell and make a map of nodes to Netgen node IDs - // ------------------------------------------------------------------- + NETGENPlugin_NetgenLibWrapper ngLib; + Ng_Mesh * Netgen_mesh = ngLib._ngMesh; - StdMeshers_Helper* myTool = new StdMeshers_Helper(aMesh); - bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); + // vector of nodes in which node index == netgen ID + vector< const SMDS_MeshNode* > nodeVec; + { + const int invalid_ID = -1; - typedef map< const SMDS_MeshNode*, int> TNodeToIDMap; - TNodeToIDMap nodeToNetgenID; - list< const SMDS_MeshElement* > triangles; - list< bool > isReversed; // orientation of triangles + 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; + // maps nodes to ng ID + typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap; + typedef TNodeToIDMap::value_type TN2ID; + TNodeToIDMap nodeToNetgenID; - 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 ) + // find internal shapes + NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true ); + + // --------------------------------- + // Feed the Netgen with surface mesh + // --------------------------------- + + TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType(); + bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID ); + + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); + if ( _viscousLayersHyp ) { - bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); + proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape ); + if ( !proxyMesh ) + return false; + } + if ( aMesh.NbQuadrangles() > 0 ) + { + StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; + Adaptor->Compute(aMesh,aShape,proxyMesh.get()); + proxyMesh.reset( Adaptor ); + } + for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next()) + { + const TopoDS_Shape& aShapeFace = exFa.Current(); + int faceID = meshDS->ShapeToIndex( aShapeFace ); + bool isInternalFace = internals.isInternalShape( faceID ); + bool isRev = false; + 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 = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace )); + + const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( 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 || - !( elem->NbNodes()==3 || ( _quadraticMesh && elem->NbNodes()==6) ) ) { - INFOS( "NETGENPlugin_NETGEN_3D::Compute(), bad mesh"); - delete myTool; myTool = 0; - return false; - } - // 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 )); - } -#ifdef _DEBUG_ - // check if a trainge is degenerated - areaControl.GetPoints( elem, nodesCoords ); - double area = areaControl.GetValue( nodesCoords ); - if ( area <= DBL_MIN ) { - 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 - // --------------------------------- + if ( !elem ) + return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + if ( elem->NbCornerNodes() != 3 ) + return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); - int Netgen_NbOfNodes = 0; - int Netgen_param2ndOrder = 0; - double Netgen_paramFine = 1.; - double Netgen_paramSize = _maxElementVolume; + // Add nodes of triangles and triangles them-selves to netgen mesh - double Netgen_point[3]; - int Netgen_triangle[3]; - int Netgen_tetrahedron[4]; + // add three nodes of triangle + bool hasDegen = false; + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = elem->GetNode( iN ); + const int shapeID = node->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; + } + int& ngID = nodeToNetgenID.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[ isRev ? 2-iN : iN ] = ngID; + } + // add triangle + if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] || + Netgen_triangle[0] == Netgen_triangle[2] || + Netgen_triangle[2] == Netgen_triangle[1] )) + continue; - Ng_Init(); + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); - Ng_Mesh * Netgen_mesh = Ng_NewMesh(); + if ( isInternalFace && !proxyMesh->IsTemporary( elem )) + { + swap( Netgen_triangle[1], Netgen_triangle[2] ); + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + } + } // loop on elements on a face + } // loop on faces of a SOLID or SHELL - // 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; - } + // insert old nodes into nodeVec + nodeVec.resize( nodeToNetgenID.size() + 1, 0 ); + TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); + for ( ; n_id != nodeToNetgenID.end(); ++n_id ) + nodeVec[ n_id->second ] = n_id->first; + nodeToNetgenID.clear(); - // 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 ( internals.hasInternalVertexInSolid() ) + { + netgen::OCCGeometry occgeo; + NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo, + (netgen::Mesh&) *Netgen_mesh, + nodeVec, + internals); } - 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); } // ------------------------- // Generate the volume mesh // ------------------------- - Ng_Meshing_Parameters Netgen_param; + return ( ngLib._isComputeOk = compute( aMesh, helper, nodeVec, Netgen_mesh)); +} - Netgen_param.secondorder = Netgen_param2ndOrder; - Netgen_param.fineness = Netgen_paramFine; - Netgen_param.maxh = Netgen_paramSize; +//================================================================================ +/*! + * \brief set parameters and generate the volume mesh + */ +//================================================================================ + +bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh& aMesh, + SMESH_MesherHelper& helper, + vector< const SMDS_MeshNode* >& nodeVec, + Ng_Mesh * Netgen_mesh) +{ + netgen::multithread.terminate = 0; - Ng_Result status; + netgen::Mesh* ngMesh = (netgen::Mesh*)Netgen_mesh; + int Netgen_NbOfNodes = Ng_GetNP(Netgen_mesh); - try { - status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param); +#ifndef NETGEN_V5 + char *optstr = 0; +#endif + int startWith = netgen::MESHCONST_MESHVOLUME; + int endWith = netgen::MESHCONST_OPTVOLUME; + int err = 1; + + NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true ); + netgen::OCCGeometry occgeo; + + if ( _hypParameters ) + { + aMesher.SetParameters( _hypParameters ); + if ( !_hypParameters->GetOptimize() ) + endWith = netgen::MESHCONST_MESHVOLUME; } - catch (...) { - MESSAGE("An exception has been caught during the Volume Mesh Generation ..."); - status = NG_VOLUME_FAILURE; + else if ( _hypMaxElementVolume ) + { + netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); + } + else if ( aMesh.HasShapeToMesh() ) + { + aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh ); + netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2; + } + else + { + netgen::Point3d pmin, pmax; + ngMesh->GetBox (pmin, pmax); + netgen::mparam.maxh = Dist(pmin, pmax)/2; } - int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); + if ( !_hypParameters && aMesh.HasShapeToMesh() ) + { + netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh ); + } - int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); + try + { +#if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100 + OCC_CATCH_SIGNALS; +#endif +#ifdef NETGEN_V5 + ngMesh->CalcLocalH(netgen::mparam.grading); + err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith); +#else + ngMesh->CalcLocalH(); + err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr); +#endif + if(netgen::multithread.terminate) + return false; + 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 (netgen::NgException exc) + { + SMESH_Comment str("NgException"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + str << ": " << exc.What(); + error(str); + } + catch (...) + { + SMESH_Comment str("Exception in netgen::OCCGenerateMesh()"); + if ( strlen( netgen::multithread.task ) > 0 ) + str << " at " << netgen::multithread.task; + error(str); + } - MESSAGE("End of Volume Mesh Generation. status=" << status << + int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh); + int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh); + + MESSAGE("End of Volume Mesh Generation. err=" << err << ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes << ", nb tetra: " << Netgen_NbOfTetra); @@ -309,88 +439,258 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, // Feed back the SMESHDS with the generated Nodes and Volume Elements // ------------------------------------------------------------------- - bool isOK = ( status == NG_OK && Netgen_NbOfTetra > 0 ); + if ( err ) + { + SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec); + if ( ce && !ce->myBadElements.empty() ) + error( ce ); + } + + 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; - } + double Netgen_point[3]; + int Netgen_tetrahedron[4]; + // create and insert new nodes into nodeVec + nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 ); int nodeIndex = Netgen_NbOfNodes + 1; - int shapeID = meshDS->ShapeToIndex( aShape ); for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex ) { Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point ); - SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0], - Netgen_point[1], - Netgen_point[2]); - meshDS->SetNodeInVolume(node, shapeID); - nodeVec.at(nodeIndex) = node; + nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]); } // create tetrahedrons 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] )); - meshDS->SetMeshElementOnShape(elt, shapeID ); + try + { + helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ), + nodeVec.at( Netgen_tetrahedron[1] ), + nodeVec.at( Netgen_tetrahedron[2] ), + nodeVec.at( Netgen_tetrahedron[3] )); + } + catch (...) + { + } } } - Ng_DeleteMesh(Netgen_mesh); - Ng_Exit(); - - delete myTool; myTool = 0; - - return isOK; + return !err; } -//============================================================================= +//================================================================================ /*! - * + * \brief Compute tetrahedral mesh from 2D mesh without geometry */ -//============================================================================= +//================================================================================ -ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save) +bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh, + SMESH_MesherHelper* aHelper) { - return save; -} + const int invalid_ID = -1; -//============================================================================= -/*! - * - */ -//============================================================================= + netgen::multithread.terminate = 0; -istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load) + SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh(); + if ( MeshType == SMESH_MesherHelper::COMP ) + return error( COMPERR_BAD_INPUT_MESH, + SMESH_Comment("Mesh with linear and quadratic elements given")); + + aHelper->SetIsQuadratic( MeshType == SMESH_MesherHelper::QUADRATIC ); + + // --------------------------------- + // Feed the Netgen with surface mesh + // --------------------------------- + + 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]; + + NETGENPlugin_NetgenLibWrapper ngLib; + Ng_Mesh * Netgen_mesh = ngLib._ngMesh; + + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh )); + if ( aMesh.NbQuadrangles() > 0 ) + { + StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; + Adaptor->Compute(aMesh); + proxyMesh.reset( Adaptor ); + } + + // maps nodes to ng ID + typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap; + typedef TNodeToIDMap::value_type TN2ID; + TNodeToIDMap nodeToNetgenID; + + SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces(); + while( fIt->more()) + { + // check element + const SMDS_MeshElement* elem = fIt->next(); + if ( !elem ) + return error( COMPERR_BAD_INPUT_MESH, "Null element encounters"); + if ( elem->NbCornerNodes() != 3 ) + return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters"); + + // add three nodes of triangle + for ( int iN = 0; iN < 3; ++iN ) + { + const SMDS_MeshNode* node = elem->GetNode( iN ); + int& ngID = nodeToNetgenID.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; + } + Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle); + } + proxyMesh.reset(); // delete tmp faces + + // vector of nodes in which node index == netgen ID + vector< const SMDS_MeshNode* > nodeVec ( nodeToNetgenID.size() + 1 ); + // insert old nodes into nodeVec + TNodeToIDMap::iterator n_id = nodeToNetgenID.begin(); + for ( ; n_id != nodeToNetgenID.end(); ++n_id ) + nodeVec.at( n_id->second ) = n_id->first; + nodeToNetgenID.clear(); + + // ------------------------- + // Generate the volume mesh + // ------------------------- + + return ( ngLib._isComputeOk = compute( aMesh, *aHelper, nodeVec, Netgen_mesh)); +} + +void NETGENPlugin_NETGEN_3D::CancelCompute() { - return load; + SMESH_Algo::CancelCompute(); + netgen::multithread.terminate = 1; } -//============================================================================= +//================================================================================ /*! - * + * \brief Return Compute progress */ -//============================================================================= +//================================================================================ -ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp) +double NETGENPlugin_NETGEN_3D::GetProgress() const { - return hyp.SaveTo( save ); + double res; + const char* volMeshing = "Volume meshing"; + const char* dlnMeshing = "Delaunay meshing"; + if ( !_optimizationStarted && + ( strncmp( netgen::multithread.task, dlnMeshing, 3 ) == 0 || + strncmp( netgen::multithread.task, volMeshing, 3 ) == 0 )) + { + res = 0.001 + 0.5 * netgen::multithread.percent / 100.; // [0., 0.5] + } + else // different otimizations + { + ((NETGENPlugin_NETGEN_3D*)this)->_optimizationStarted = true; + res = 0.5 + 0.5 * SMESH_Algo::GetProgressByTic(); // [0.5, 1.] + } + return res; } //============================================================================= /*! - * + * */ //============================================================================= -istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp) +bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + MapShapeNbElems& aResMap) { - return hyp.LoadFrom( load ); + int nbtri = 0, nbqua = 0; + double fullArea = 0.0; + 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() ) { + SMESH_ComputeErrorPtr& smError = sm->GetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this)); + return false; + } + std::vector aVec = (*anIt).second; + nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]); + nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]); + GProp_GProps G; + BRepGProp::SurfaceProperties(F,G); + double anArea = G.Mass(); + fullArea += anArea; + } + + // collect info from edges + int nb0d_e = 0, nb1d_e = 0; + bool IsQuadratic = false; + bool IsFirst = true; + TopTools_MapOfShape tmpMap; + 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(expF.Current()); + MapShapeNbElemsItr anIt = aResMap.find(aSubMesh); + if( anIt==aResMap.end() ) { + SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED, + "Submesh can not be evaluated",this)); + return false; + } + std::vector aVec = (*anIt).second; + nb0d_e += aVec[SMDSEntity_Node]; + nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); + if(IsFirst) { + IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]); + IsFirst = false; + } + } + tmpMap.Clear(); + + double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) ); + double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. ); + double ELen = Min(ELen_vol,ELen_face*2); + + GProp_GProps G; + BRepGProp::VolumeProperties(aShape,G); + double aVolume = G.Mass(); + double tetrVol = 0.1179*ELen*ELen*ELen; + double CoeffQuality = 0.9; + int nbVols = int( aVolume/tetrVol/CoeffQuality ); + int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2; + int nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5; + std::vector aVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; i