-// 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
// 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
#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<int,double> VertexId2LocalSize;
// }
// }
}
+ //================================================================================
+ /*!
+ * \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 );
+ }
+ }
+
}
//================================================================================
*/
//================================================================================
-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;
netgen::Mesh& ngMesh,
vector<const SMDS_MeshNode*>& nodeVec,
const list< SMESH_subMesh* > & meshedSM,
+ SMESH_MesherHelper* quadHelper,
SMESH_ProxyMesh::Ptr proxyMesh)
{
TNode2IdMap nodeNgIdMap;
}
} // 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
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
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 );
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
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);
netgen::OCCGeometry& geom,
const TSideVector& wires,
SMESH_MesherHelper& helper,
- vector< const SMDS_MeshNode* > & nodeVec)
+ vector< const SMDS_MeshNode* > & nodeVec,
+ const bool overrideMinH)
{
// ----------------------------
// Check wires and count nodes
// 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 sunH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
- int nbSeg = ( int( segLen[ iPrev ] > sunH / 100.) +
- int( segLen[ i ] > sunH / 100.) +
- int( segLen[ iNext ] > sunH / 100.));
+ 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), sunH / nbSeg );
+ RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg, overrideMinH );
}
if ( isInternalWire )
{
* \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
*/
//================================================================================
const NETGENPlugin_ngMeshInfo& initState,
SMESH_Mesh& sMesh,
std::vector<const SMDS_MeshNode*>& nodeVec,
- SMESH_Comment& comment)
+ SMESH_Comment& comment,
+ SMESH_MesherHelper* quadHelper)
{
int nbNod = ngMesh.GetNP();
int nbSeg = ngMesh.GetNSeg();
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
// -------------------------------------
{
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
{
// from computation of 3D mesh
ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
+ vector<const SMDS_MeshNode*> nodes;
for (i = nbInitFac+1; i <= nbFac; ++i )
{
const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
TopoDS_Face aFace;
if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
- vector<SMDS_MeshNode*> nodes;
+ nodes.clear();
for (int j=1; j <= elem.GetNP(); ++j)
{
int pind = elem.PNum(j);
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);
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
TopoDS_Solid aSolid;
if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
- vector<SMDS_MeshNode*> nodes;
+ nodes.clear();
for (int j=1; j <= elem.GetNP(); ++j)
{
int pind = elem.PNum(j);
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;
return str;
}
+ //================================================================================
+ /*!
+ * \brief Looks for triangles lying on a SOLID
+ */
+ //================================================================================
+
+ bool hasBadElemOnSolid( const list<const SMDS_MeshElement*>& elems,
+ SMESH_subMesh* solidSM )
+ {
+ TopTools_IndexedMapOfShape solidSubs;
+ TopExp::MapShapes( solidSM->GetSubShape(), solidSubs );
+ SMESHDS_Mesh* mesh = solidSM->GetFather()->GetMeshDS();
+
+ list<const SMDS_MeshElement*>::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;
" fuse edges = " << netgen::merge_solids);
SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
+ SMESH_MesherHelper quadHelper( *_mesh );
+ quadHelper.SetIsQuadratic( mparams.secondorder );
static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( ngMesh, debugFile )
while debugging netgen */
if ( !err )
{
err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
- FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
+ FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper));
}
initState = NETGENPlugin_ngMeshInfo(_ngMesh);
if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
{
// load SMESH with computed segments and faces
- FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
+ FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
// compute pyramids on quadrangles
SMESH_ProxyMesh::Ptr proxyMesh;
quadFaceSM.push_back( _mesh->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");
}
{
// 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
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)
{
// 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;
{
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() );
}
{
Ng_Init();
- // redirect all netgen output (mycout,myerr,cout) to _outputFileName
- _isComputeOk = false;
- _outputFileName = getOutputFileName();
- netgen::mycout = new ofstream ( _outputFileName.c_str() );
- netgen::myerr = netgen::mycout;
- _coutBuffer = std::cout.rdbuf();
+ _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() );
+ std::cout.rdbuf( netgen::mycout->rdbuf() );
#endif
+ }
_ngMesh = Ng_NewMesh();
}
Ng_DeleteMesh( _ngMesh );
Ng_Exit();
NETGENPlugin_Mesher::RemoveTmpFiles();
- std::cout.rdbuf( _coutBuffer );
+ if ( _coutBuffer )
+ std::cout.rdbuf( _coutBuffer );
#ifdef _DEBUG_
if( _isComputeOk )
#endif
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
}