ForkMaster/netgenplugin
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

0020676: EDF 1212 GEOM: Partition operation creates vertices which causes mesh computation to fail with netgen

Browse Source 
* Solve problem with internal edges and faces
This commit is contained in:
eap 2010-03-18 09:51:47 +00:00
parent 6dbabf8a00
commit 1ebad126db
5 changed files with 810 additions and 388 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

792
src/NETGENPlugin/NETGENPlugin_Mesher.cxx
View File

File diff suppressed because it is too large Load Diff

85
src/NETGENPlugin/NETGENPlugin_Mesher.hxx
View File

@ -30,8 +30,15 @@
#include "NETGENPlugin_Defs.hxx"
#include "StdMeshers_FaceSide.hxx"
#include "SMDS_MeshElement.hxx"
namespace nglib {
#include <nglib.h>
}
#include <map>
#include <vector>
#include <set>
class SMESH_Mesh;
class SMESH_Comment;
@ -40,6 +47,7 @@ class TopoDS_Shape;
class TopTools_DataMapOfShapeShape;
class NETGENPlugin_Hypothesis;
class NETGENPlugin_SimpleHypothesis_2D;
class NETGENPlugin_Internals;
namespace netgen {
class OCCGeometry;
class Mesh;
@ -81,19 +89,20 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
const TopoDS_Shape& shape,
SMESH_Mesh& mesh,
std::list< SMESH_subMesh* > * meshedSM=0,
TopTools_DataMapOfShapeShape* internalE2F=0);
NETGENPlugin_Internals* internalShapes=0);
static int FillSMesh(const netgen::OCCGeometry& occgeom,
const netgen::Mesh& ngMesh,
const NETGENPlugin_ngMeshInfo& initState,
SMESH_Mesh& sMesh,
std::vector<SMDS_MeshNode*>& nodeVec,
std::vector<const SMDS_MeshNode*>& nodeVec,
SMESH_Comment& comment);
bool fillNgMesh(netgen::OCCGeometry& occgeom,
netgen::Mesh& ngMesh,
std::vector<SMDS_MeshNode*>& nodeVec,
const std::list< SMESH_subMesh* > & meshedSM);
std::vector<const SMDS_MeshNode*>& nodeVec,
const std::list< SMESH_subMesh* > & meshedSM,
NETGENPlugin_Internals* internalShapes=0);
void defaultParameters();
@ -111,4 +120,72 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher
std::map< int, std::pair<int,int> > _faceDescriptors;
};
//=============================================================================
/*!
* \brief Container of info needed to solve problems with internal shapes.
*
* Issue 0020676. It is made up as a class to be ready to extract from NETGEN
* and put in SMESH as soon as the same solution is needed somewhere else.
* The approach is to precompute internal edges in 2D and internal faces in 3D
* and put their mesh correctly (twice) into netgen mesh.
* In 2D, this class finds internal edges in faces and their vertices.
* In 3D, it additionally finds internal faces, their edges shared with other faces,
* and their vertices shared by several internal edges. Nodes built on the found
* shapes and mesh faces built on the found internal faces are to be doubled in
* netgen mesh to emulate a "crack"
*/
//=============================================================================
class NETGENPLUGIN_EXPORT NETGENPlugin_Internals
{
SMESH_Mesh& _mesh;
bool _is3D;
//2D
std::map<int,int> _ev2face; //!< edges and vertices in faces where they are TopAbs_INTERNAL
// 3D
std::set<int> _intShapes;
std::set<int> _borderFaces; //!< non-intrnal faces sharing the internal edge
public:
NETGENPlugin_Internals( SMESH_Mesh& mesh, const TopoDS_Shape& shape, bool is3D );
bool isShapeToPrecompute(const TopoDS_Shape& s);
// 2D
bool hasInternalEdges() const { return !_ev2face.empty(); }
bool isInternalEdge(int id ) const { return _ev2face.count( id ); }
bool isInternalVertex(int id ) const { return _ev2face.count( id ); }
const std::map<int,int>& getEdgesAndVerticesWithFaces() const { return _ev2face; }
void getInternalEdges( TopTools_IndexedMapOfShape& fmap,
TopTools_IndexedMapOfShape& emap,
TopTools_IndexedMapOfShape& vmap,
list< SMESH_subMesh* >& smToPrecompute);
// 3D
bool hasInternalFaces() const { return !_intShapes.empty(); }
bool isInternalShape(int id ) const { return _intShapes.count( id ); }
void findBorderElements( std::set< const SMDS_MeshElement*, TIDCompare > & borderElems );
bool isBorderFace(int faceID ) const { return _borderFaces.count( faceID ); }
void getInternalFaces( TopTools_IndexedMapOfShape& fmap,
TopTools_IndexedMapOfShape& emap,
list< SMESH_subMesh* >& facesSM,
list< SMESH_subMesh* >& boundarySM);
};
//================================================================================
/*!
* \brief It correctly initializes netgen library at constructor and
* correctly finishes using netgen library at destructor
*/
//================================================================================
struct NETGENPLUGIN_EXPORT NETGENPlugin_NetgenLibWrapper
{
nglib::Ng_Mesh * _ngMesh;
NETGENPlugin_NetgenLibWrapper();
~NETGENPlugin_NetgenLibWrapper();
void setMesh( nglib::Ng_Mesh* mesh );
};
#endif

90
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx
View File

@ -174,6 +174,7 @@ static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
// Check wires and count nodes
// ----------------------------
int nbNodes = 0;
double totalLength = 0;
for ( int iW = 0; iW < wires.size(); ++iW )
{
StdMeshers_FaceSidePtr wire = wires[ iW ];
@ -187,7 +188,8 @@ static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
(new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,
SMESH_Comment("Unexpected nb of points on wire ") << iW
<< ": " << uvPtVec.size()<<" != "<<wire->NbPoints()));
nbNodes += wire->NbPoints();
nbNodes += wire->NbPoints();
totalLength += wire->Length();
}
nodeVec.reserve( nbNodes );
@ -200,7 +202,8 @@ static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
// ngMesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); // set grading
const int faceID = 1, solidID = 0;
ngMesh.AddFaceDescriptor (FaceDescriptor(faceID, solidID, solidID, 0));
if ( ngMesh.GetNFD() < 1 )
ngMesh.AddFaceDescriptor (FaceDescriptor(faceID, solidID, solidID, 0));
for ( int iW = 0; iW < wires.size(); ++iW )
{
@ -231,8 +234,8 @@ static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
Segment seg;
#ifdef NETGEN_NEW
seg.pnums[0] = ngMesh.GetNP(); // ng node id
seg.pnums[1] = seg.pnums[0] + 1; // ng node id
seg.pnums[0] = ngMesh.GetNP(); // ng node id
seg.pnums[1] = seg.pnums[0] + 1; // ng node id
#else
seg.p1 = ngMesh.GetNP(); // ng node id
seg.p2 = seg.p1 + 1; // ng node id
@ -313,7 +316,71 @@ static TError AddSegmentsToMesh(netgen::Mesh& ngMesh,
} // loop on wires of a face
ngMesh.CalcSurfacesOfNode();
// add a segment instead of internal vertex
// const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
// for ( TopoDS_Iterator sh (face); sh.More(); sh.Next())
// {
// if ( sh.Value().ShapeType() != TopAbs_VERTEX ) continue;
// const TopoDS_Vertex V = TopoDS::Vertex( sh.Value() );
// SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
// sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
// const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, helper.GetMeshDS() );
// if ( !nV ) continue;
// double segLen = totalLength / ngMesh.GetNSeg() / 2;
// bool uvOK = false;
// gp_XY uvV = helper.GetNodeUV( face, nV, 0, &uvOK );
// if ( !uvOK ) helper.CheckNodeUV( face, nV, uvV, BRep_Tool::Tolerance( V ),/*force=*/true);
// gp_XY uvP( uvV.X() + segLen, uvV.Y() );
// TopLoc_Location loc;
// Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
// gp_Pnt P = surf->Value( uvP.X(), uvP.Y() ).Transformed( loc );
// MeshPoint mpV( Point<3> (nV->X(), nV->Y(), nV->Z()) );
// MeshPoint mpP( Point<3> (P.X(), P.Y(), P.Z()));
// ngMesh.AddPoint ( mpV, 1, EDGEPOINT );
// ngMesh.AddPoint ( mpP, 1, EDGEPOINT );
// nodeVec.push_back( nV );
// // Add the segment
// Segment seg;
// #ifdef NETGEN_NEW
// seg.pnums[0] = ngMesh.GetNP()-1; // ng node id
// seg.pnums[1] = ngMesh.GetNP(); // ng node id
// #else
// seg.p1 = ngMesh.GetNP()-1; // ng node id
// seg.p2 = ngMesh.GetNP(); // ng node id
// #endif
// seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
// seg.si = faceID; // = geom.fmap.FindIndex (face);
// seg.epgeominfo[ 0 ].dist = 0; // param on curve
// seg.epgeominfo[ 0 ].u = uvV.X();
// seg.epgeominfo[ 0 ].v = uvV.Y();
// seg.epgeominfo[ 1 ].dist = segLen; // param on curve
// seg.epgeominfo[ 1 ].u = uvP.X();
// seg.epgeominfo[ 1 ].v = uvP.Y();
// seg.epgeominfo[ 0 ].edgenr = 10; // = geom.emap.FindIndex(edge);
// seg.epgeominfo[ 1 ].edgenr = 10; // = geom.emap.FindIndex(edge);
// ngMesh.AddSegment (seg);
// // add reverse segment
// #ifdef NETGEN_NEW
// swap (seg.pnums[0], seg.pnums[1]);
// #else
// swap (seg.p1, seg.p2);
// #endif
// swap( seg.epgeominfo[0], seg.epgeominfo[1] );
// seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
// ngMesh.AddSegment (seg);
// }
ngMesh.CalcSurfacesOfNode();
return TError();
}
@ -356,8 +423,8 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh,
// Make input netgen mesh
// -------------------------
Ng_Init();
netgen::Mesh * ngMesh = new netgen::Mesh ();
NETGENPlugin_NetgenLibWrapper ngLib;
netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh;
netgen::OCCGeometry occgeo;
NETGENPlugin_Mesher::PrepareOCCgeometry( occgeo, F, aMesh );
@ -366,10 +433,8 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh,
vector< const SMDS_MeshNode* > nodeVec;
problem = AddSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
if ( problem && !problem->IsOK() ) {
delete ngMesh; Ng_Exit();
if ( problem && !problem->IsOK() )
return error( problem );
}
// --------------------
// compute edge length
@ -480,11 +545,6 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh,
face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
}
Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
Ng_Exit();
NETGENPlugin_Mesher::RemoveTmpFiles();
return !err;
}

16
src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.hxx
View File

@ -29,25 +29,9 @@
#include "SMESH_2D_Algo.hxx"
#include "SMESH_Mesh.hxx"
/*#define OCCGEOMETRY
#include <occgeom.hpp>
#include <meshing.hpp>//amv*/
class StdMeshers_MaxElementArea;
class StdMeshers_LengthFromEdges;
class StdMeshers_QuadranglePreference;
//class NETGENPlugin_Hypothesis;
/*namespace netgen {
class OCCGeometry;
}*/
/*namespace netgen {
class OCCGeometry;
extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
extern MeshingParameters mparam;
}*/
//using namespace netgen;
/*!
* \brief Mesher generating 2D elements on a geometrical face taking

215
src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx
View File

@ -151,158 +151,6 @@ 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<int>& internalShapeIds,
SMESH_MesherHelper& helper,
TIDSortedElemSet & borderElems,
set<int> & borderFaceIds )
{
SMESH_Mesh* mesh = helper.GetMesh();
SMESHDS_Mesh* meshDS = helper.GetMeshDS();
// loop on internal geom edges
set<int>::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 <s>
int intFace = 0;
set<int>::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 <s>
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 );
}
}
}
}
// <suspectFaces> having link with same orientation as mesh faces on
// the internal geom face are <borderElems>.
// Some of them have only one node on edge s, we collect them to decide
// later by links of found <borderElems>
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
@ -328,7 +176,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
int Netgen_triangle[3];
int Netgen_tetrahedron[4];
TNetgenLibWrapper ngLib;
NETGENPlugin_NetgenLibWrapper ngLib;
Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
// maps of 1) ordinary nodes and 2) doubled nodes on internal shapes
@ -342,60 +190,17 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
SMESH::Controls::Area areaControl;
SMESH::Controls::TSequenceOfXYZ nodesCoords;
// --------------------------------------------------------------------
// 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".
// --------------------------------------------------------------------
set<int> internalShapeIds;
set<int> borderFaceIds; //!< non-internal geom faces sharing internal edge
NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
// mesh faces on <borderFaceIds>, having nodes on internal edge that are
// to be replaced by doubled nodes
// mesh faces on non-internal geom faces sharing internal edge, whose some nodes
// are on internal edge and 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<int>::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<int> 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 );
}
internals.findBorderElements( borderElems );
// ---------------------------------
// Feed the Netgen with surface mesh
@ -408,12 +213,12 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
if ( aMesh.NbQuadrangles() > 0 )
Adaptor.Compute(aMesh,aShape);
for ( exFa.ReInit(); exFa.More(); exFa.Next())
for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); 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 isInternalFace = internals.isInternalShape( faceID );
bool isBorderFace = internals.isBorderFace( faceID );
bool isRev = false;
if ( checkReverse && !isInternalFace &&
helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
@ -471,7 +276,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
hasDegen = true;
}
bool isDblN = isDblF && internalShapeIds.count( shapeID );
bool isDblN = isDblF && internals.isInternalShape( shapeID );
int& ngID = nodeToNetgenID[isDblN].insert(TN2ID( node, invalid_ID )).first->second;
if ( ngID == invalid_ID )
{
@ -681,7 +486,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
int Netgen_triangle[3];
int Netgen_tetrahedron[4];
TNetgenLibWrapper ngLib;
NETGENPlugin_NetgenLibWrapper ngLib;
Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
// set nodes and remember thier netgen IDs