#include <BRepLProp_SLProps.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRep_Builder.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
#include <BRep_Tool.hxx>
#include <Bnd_B3d.hxx>
#include <GeomLib_IsPlanarSurface.hxx>
else if (meshedSM)
{
const int dim = SMESH_Gen::GetShapeDim( shape );
- meshedSM[ dim ].push_back( sm );
+ meshedSM[ dim ].push_back( sm );
}
}
}
occgeo.BuildFMapFromPrefilled();
}
+//================================================================================
+/*!
+ * \brief Initialize netgen::OCCGeometry with OCCT shape
+ */
+//================================================================================
+
+void NETGENPlugin_Mesher::PrepareOCCgeometry2(netgen::OCCGeometry& occgeo,
+ const TopoDS_Shape& shape,
+ SMESH_Mesh& mesh,
+ list< SMESH_subMesh* > * meshedSM,
+ NETGENPlugin_Internals* intern)
+{
+ updateTriangulation( shape );
+
+ Bnd_Box bb;
+ // BRepBndLib::Add (shape, bb);
+ GEOMUtils::PreciseBoundingBox(shape,bb);
+ double x1,y1,z1,x2,y2,z2;
+ bb.Get (x1,y1,z1,x2,y2,z2);
+ netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
+ netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
+ occgeo.boundingbox = netgen::Box<3> (p1,p2);
+
+ occgeo.shape = shape;
+ occgeo.changed = 1;
+
+ // fill maps of shapes of occgeo with not yet meshed subshapes
+
+ // get root submeshes
+ list< SMESH_subMesh* > rootSM;
+ const int shapeID = mesh.GetMeshDS()->ShapeToIndex( shape );
+ if ( shapeID > 0 ) { // SMESH_subMesh with ID 0 may exist, don't use it!
+ rootSM.push_back( mesh.GetSubMesh( shape ));
+ }
+ else {
+ for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
+ rootSM.push_back( mesh.GetSubMesh( it.Value() ));
+ }
+
+ int totNbFaces = 0;
+
+ // add subshapes of empty submeshes
+ list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
+ for ( ; rootIt != rootEnd; ++rootIt ) {
+ SMESH_subMesh * root = *rootIt;
+ SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
+ /*complexShapeFirst=*/true);
+ // to find a right orientation of subshapes (PAL20462)
+ TopTools_IndexedMapOfShape subShapes;
+ TopExp::MapShapes(root->GetSubShape(), subShapes);
+ while ( smIt->more() )
+ {
+ SMESH_subMesh* sm = smIt->next();
+ TopoDS_Shape shape = sm->GetSubShape();
+ totNbFaces += ( shape.ShapeType() == TopAbs_FACE );
+ if ( intern && intern->isShapeToPrecompute( shape ))
+ continue;
+
+ if ( shape.ShapeType() != TopAbs_VERTEX )
+ shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape
+ if ( shape.Orientation() >= TopAbs_INTERNAL )
+ shape.Orientation( TopAbs_FORWARD ); // issue 0020676
+
+ switch ( shape.ShapeType() ) {
+ case TopAbs_FACE : occgeo.fmap.Add( shape ); break;
+ case TopAbs_EDGE : occgeo.emap.Add( shape ); break;
+ case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
+ case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
+ default:;
+ }
+
+ // collect submeshes of meshed shapes
+ if ( meshedSM && !sm->IsEmpty())
+ {
+ const int dim = SMESH_Gen::GetShapeDim( shape );
+ meshedSM[ dim ].push_back( sm );
+
+ TopoDS_Shape shape = sm->GetSubShape();
+ switch ( shape.ShapeType() ) {
+ /* Fill map with vertex of premeshed elements*/
+ case TopAbs_EDGE : {
+ /* Rules to create geometric elements see TopoDS_Builder.cxx */
+ /* Vertex can be added to edges!*/
+ SMDS_NodeIteratorPtr nodeIt = sm->GetSubMeshDS()->GetNodes();
+
+ while ( nodeIt->more() )
+ {
+ printf("Including geom vertex associated to mesh nodes!\n\n");
+ const SMDS_MeshNode* node = nodeIt->next();
+ gp_Pnt point(node->X(),node->Y(),node->Z());
+ TopoDS_Vertex v = BRepBuilderAPI_MakeVertex( point );
+ /* Optionaly we could associate the edge to the edge itself too!*/
+ occgeo.vmap.Add( v );
+ }
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ occgeo.facemeshstatus.SetSize (totNbFaces);
+ occgeo.facemeshstatus = 0;
+ occgeo.face_maxh_modified.SetSize(totNbFaces);
+ occgeo.face_maxh_modified = 0;
+ occgeo.face_maxh.SetSize(totNbFaces);
+ occgeo.face_maxh = netgen::mparam.maxh;
+
+ //
+ occgeo.BuildFMapFromPrefilled();
+}
+
//================================================================================
/*!
* \brief Return a default min size value suitable for the given geometry.
switch ( sm->GetSubShape().ShapeType() )
{
case TopAbs_EDGE: { // EDGE
-
+ // ----------------------
TopoDS_Edge geomEdge = TopoDS::Edge( sm->GetSubShape() );
if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
// get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
helper.SetSubShape( face );
list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
- visitedEdgeSM2Faces );
+ visitedEdgeSM2Faces );
if ( edges.empty() )
continue; // wrong ancestor?
if ( points.empty() )
return false; // invalid node params?
smIdType i, nbSeg = fSide.NbSegments();
- for ( i = 0; i < nbSeg; ++i )
- {
- const SMDS_MeshNode * n0 = points[ i ].node;
- const SMDS_MeshNode * n1 = points[ i+1 ].node;
- double size = SMESH_TNodeXYZ( n0 ).Distance( n1 );
- netgen::Point3d p0(n0->X(), n0->Y(), n0->Z());
- netgen::Point3d p1(n1->X(), n1->Y(), n1->Z());
- ngMesh.RestrictLocalHLine( p0, p1, size, true );
- }
-
- }
-///////////////////////////////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////////////////////////////
- {
- /*
- // ----------------------
- TopoDS_Edge geomEdge = TopoDS::Edge( sm->GetSubShape() );
- if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
- geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
-
- // Add ng segments for each not meshed FACE the EDGE bounds
- PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE );
- while ( const TopoDS_Shape * anc = fIt->next() )
- {
- faceNgID = occgeom.fmap.FindIndex( *anc );
- if ( faceNgID < 1 )
- continue; // meshed face
-
- int faceSMDSId = meshDS->ShapeToIndex( *anc );
- if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId ))
- continue; // already treated EDGE
-
- TopoDS_Face face = TopoDS::Face( occgeom.fmap( faceNgID ));
- if ( face.Orientation() >= TopAbs_INTERNAL )
- face.Orientation( TopAbs_FORWARD ); // issue 0020676
-
- // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
- helper.SetSubShape( face );
- list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
- visitedEdgeSM2Faces );
- if ( edges.empty() )
- continue; // wrong ancestor?
-
- // find out orientation of <edges> within <face>
- TopoDS_Edge eNotSeam = edges.front();
- if ( helper.HasSeam() )
- {
- list< TopoDS_Edge >::iterator eIt = edges.begin();
- while ( helper.IsRealSeam( *eIt )) ++eIt;
- if ( eIt != edges.end() )
- eNotSeam = *eIt;
- }
- TopAbs_Orientation fOri = helper.GetSubShapeOri( face, eNotSeam );
- bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL );
- // get all nodes from connected <edges>
- const bool skipMedium = netgen::mparam.secondorder;//smDS->IsQuadratic();
- StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, skipMedium, &helper );
- const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
- if ( points.empty() )
- return false; // invalid node params?
- smIdType i, nbSeg = fSide.NbSegments();
+ // remember EDGEs of fSide to treat only once
+ for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
+ visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
- // remember EDGEs of fSide to treat only once
- for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
- visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
+ double otherSeamParam = 0;
+ bool isSeam = false;
- double otherSeamParam = 0;
- bool isSeam = false;
+ // add segments
- // add segments
+ int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
- int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
+ for ( i = 0; i < nbSeg; ++i )
+ {
+ const UVPtStruct& p1 = points[ i ];
+ const UVPtStruct& p2 = points[ i+1 ];
- for ( i = 0; i < nbSeg; ++i )
+ if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
+ {
+ isSeam = false;
+ if ( helper.IsRealSeam( p1.node->GetShapeID() ))
{
- const UVPtStruct& p1 = points[ i ];
- const UVPtStruct& p2 = points[ i+1 ];
-
- if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
- {
- isSeam = false;
- if ( helper.IsRealSeam( p1.node->GetShapeID() ))
- {
- TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
- isSeam = helper.IsRealSeam( e );
- if ( isSeam )
- {
- otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
- }
- }
- }
- netgen::Segment seg;
- // ng node ids
- seg[0] = prevNgId;
- seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
- seg.edgenr = countEdgeId+1;
- seg.si = countEdgeId+1;
- seg.epgeominfo[0].edgenr = countEdgeId;
- seg.epgeominfo[1].edgenr = countEdgeId;
- // node param on curve
- seg.epgeominfo[ 0 ].dist = p1.param;
- seg.epgeominfo[ 1 ].dist = p2.param;
- // uv on face
- seg.epgeominfo[ 0 ].u = p1.u;
- seg.epgeominfo[ 0 ].v = p1.v;
- seg.epgeominfo[ 1 ].u = p2.u;
- seg.epgeominfo[ 1 ].v = p2.v;
-
- //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
- //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge );
-
- //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
- // seg.si = faceNgID; // = geom.fmap.FindIndex (face);
- // seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
- ngMesh.AddSegment (seg);
-
- SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
- RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
-
- #ifdef DUMP_SEGMENTS
- cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl
- << "\tface index: " << seg.si << endl
- << "\tp1: " << seg[0] << endl
- << "\tp2: " << seg[1] << endl
- << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
- << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
- //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
- << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
- << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
- //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
- #endif
+ TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
+ isSeam = helper.IsRealSeam( e );
if ( isSeam )
{
- if ( helper.GetPeriodicIndex() && 1 ) {
- seg.epgeominfo[ 0 ].u = otherSeamParam;
- seg.epgeominfo[ 1 ].u = otherSeamParam;
- swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
- } else {
- seg.epgeominfo[ 0 ].v = otherSeamParam;
- seg.epgeominfo[ 1 ].v = otherSeamParam;
- swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
- }
- swap( seg[0], seg[1] );
- swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist );
- seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
- ngMesh.AddSegment( seg );
- #ifdef DUMP_SEGMENTS
- cout << "Segment: " << seg.edgenr << endl
- << "\t is SEAM (reverse) of the previous. "
- << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
- << " = " << otherSeamParam << endl;
- #endif
- }
- else if ( fOri == TopAbs_INTERNAL )
- {
- swap( seg[0], seg[1] );
- swap( seg.epgeominfo[0], seg.epgeominfo[1] );
- seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
- ngMesh.AddSegment( seg );
- #ifdef DUMP_SEGMENTS
- cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
- #endif
+ otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
}
}
- } // loop on geomEdge ancestors
+ }
+ netgen::Segment seg;
+ // ng node ids
+ seg[0] = prevNgId;
+ seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
+ // node param on curve
+ seg.epgeominfo[ 0 ].dist = p1.param;
+ seg.epgeominfo[ 1 ].dist = p2.param;
+ // uv on face
+ seg.epgeominfo[ 0 ].u = p1.u;
+ seg.epgeominfo[ 0 ].v = p1.v;
+ seg.epgeominfo[ 1 ].u = p2.u;
+ seg.epgeominfo[ 1 ].v = p2.v;
+
+ geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
+ int edgeID = occgeom.GetEdge(geomEdge).nr;
+
+ seg.epgeominfo[ 0 ].edgenr = edgeID;
+ seg.epgeominfo[ 1 ].edgenr = edgeID;
+
+ //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
+ // seg.si = faceNgID; // = geom.fmap.FindIndex (face);
+ // seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+
+ seg.si = edgeID+1; // = geom.fmap.FindIndex (face);
+ seg.edgenr = edgeID+1; // segment id
- if ( quadHelper ) // remember medium nodes of sub-meshes
+ ngMesh.AddSegment (seg);
+
+ SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
+ RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
+
+#ifdef DUMP_SEGMENTS
+ cout << "Segment: " << seg.edgenr << " on SMESH face " << meshDS->ShapeToIndex( face ) << endl
+ << "\tface index: " << seg.si << endl
+ << "\tp1: " << seg[0] << endl
+ << "\tp2: " << seg[1] << endl
+ << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
+ << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
+ //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
+ << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
+ << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
+ //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
+#endif
+ if ( isSeam )
{
- SMDS_ElemIteratorPtr edges = smDS->GetElements();
- while ( edges->more() )
- {
- const SMDS_MeshElement* e = edges->next();
- if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e )))
- break;
+ if ( helper.GetPeriodicIndex() && 1 ) {
+ seg.epgeominfo[ 0 ].u = otherSeamParam;
+ seg.epgeominfo[ 1 ].u = otherSeamParam;
+ swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
+ } else {
+ seg.epgeominfo[ 0 ].v = otherSeamParam;
+ seg.epgeominfo[ 1 ].v = otherSeamParam;
+ swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
}
+ swap( seg[0], seg[1] );
+ swap( seg.epgeominfo[0].dist, seg.epgeominfo[1].dist );
+ seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
+ ngMesh.AddSegment( seg );
+#ifdef DUMP_SEGMENTS
+ cout << "Segment: " << seg.edgenr << endl
+ << "\t is SEAM (reverse) of the previous. "
+ << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
+ << " = " << otherSeamParam << endl;
+#endif
}
-
- countEdgeId++; // Assume edge is iterate incrementally as done by the emap!
- break;
-
- */
+ else if ( fOri == TopAbs_INTERNAL )
+ {
+ swap( seg[0], seg[1] );
+ swap( seg.epgeominfo[0], seg.epgeominfo[1] );
+ seg.edgenr = edgeID+1;//ngMesh.GetNSeg() + 1; // segment id
+ ngMesh.AddSegment( seg );
+#ifdef DUMP_SEGMENTS
+ cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
+#endif
+ }
+ }
+ } // loop on geomEdge ancestors
- } /* scope of original code*/
- } // case TopAbs_EDGE
+ 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
case TopAbs_FACE: { // FACE
// ----------------------
const TopoDS_Face& geomFace = TopoDS::Face( sm->GetSubShape() );
}
+// Basicatly try to reproduce what is done in MeshFace method in basegeom.cpp class of netgen
+// where Meshing2 class is instantiated with the face, paramters, bounding box
+// then the mesh points are added to netgen and then the connectivity of the clossed chaing of nodes is feed into Meshing2 class
+SMESH_ComputeErrorPtr
+NETGENPlugin_Mesher::AddSegmentsToMesh3(netgen::Mesh& ngMesh,
+ netgen::OCCGeometry& geom,
+ const TSideVector& wires,
+ SMESH_MesherHelper& helper,
+ vector< const SMDS_MeshNode* > & nodeVec,
+ netgen::MeshingParameters & mparams,
+ const bool overrideMinH )
+{
+ // ----------------------------
+ // Check wires and count nodes
+ // ----------------------------
+ smIdType nbNodes = 0;
+ for ( size_t iW = 0; iW < wires.size(); ++iW )
+ {
+ StdMeshers_FaceSidePtr wire = wires[ iW ];
+ const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+ if ((int) uvPtVec.size() != wire->NbPoints() )
+ return SMESH_ComputeError::New(COMPERR_BAD_INPUT_MESH,
+ SMESH_Comment("Unexpected nb of points on wire ") << iW
+ << ": " << uvPtVec.size()<<" != "<<wire->NbPoints());
+ nbNodes += wire->NbPoints();
+ }
+ nodeVec.reserve( nodeVec.size() + nbNodes + 1 );
+ if ( nodeVec.empty() )
+ nodeVec.push_back( 0 );
+
+ const int predefinedNodes = ngMesh.GetNP();
+ const bool wasNgMeshEmpty = ( ngMesh.GetNP() < 1 );
+
+ map<const SMDS_MeshNode*, int > node2ngID;
+ // was empty part is ignored for now!
+ netgen::Box<3> bb = geom.GetBoundingBox();
+ bb.Increase(bb.Diam()/10);
+ netgen::Meshing2 meshing(geom, mparams, bb);
+ for ( size_t iW = 0; iW < wires.size(); ++iW )
+ {
+ StdMeshers_FaceSidePtr wire = wires[ iW ];
+ const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
+ const smIdType nbSegments = wire->NbPoints() - 1;
+
+ // assure the 1st node to be in node2ngID, which is needed to correctly
+ // "close chain of segments" (see below) in case if the 1st node is not
+ // onVertex because it is on a Viscous layer
+ node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 - predefinedNodes ));
+ // MANDATORY TO ADD FACE DESCRIPTO. OTHERWISE THE Meshing2 fails
+ ngMesh.AddFaceDescriptor(netgen::FaceDescriptor(1, 1, 1, 1));
+ // ngMesh.SetBCName(0, "Myname");
+ printf("wire id, nbSegments: %d, %d\n", iW, nbSegments );
+ for ( int i = 0; i < nbSegments; ++i ) // loop on segments
+ {
+ // Add the first point of a segment
+
+ const SMDS_MeshNode * n = uvPtVec[ i ].node;
+ const int posShapeID = n->GetShapeID();
+ bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
+ bool onEdge = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE );
+
+ // skip nodes on degenerated edges
+ if ( helper.IsDegenShape( posShapeID ) &&
+ helper.IsDegenShape( uvPtVec[ i+1 ].node->GetShapeID() ))
+ continue;
+
+ int ngID1 = ngMesh.GetNP() + 1 - predefinedNodes, ngID2 = ngID1+1;
+
+ if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID ) )
+ ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
+
+ ngID1 = node2ngID.insert( make_pair( uvPtVec[ i ].node, ngID1 )).first->second;
+
+ if ( ngID1 > ngMesh.GetNP() - predefinedNodes )
+ {
+ netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) );
+ auto pIndex = ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
+ // add point to mesh2 class
+ int pointId = meshing.AddPoint( ngMesh.Point(pIndex), ngID1 /*id of the added point*/ );
+ // printf("Point id of added point:%d %d\n", pIndex, pointId);
+ nodeVec.push_back( n );
+ }
+ }
+
+ for ( int i = 0; i < nbSegments; ++i ) // loop on segments
+ {
+ const UVPtStruct& pnt0 = uvPtVec[ i ];
+ const UVPtStruct& pnt1 = uvPtVec[ i + 1 ];
+ netgen::PointGeomInfo gi0, gi1;
+ gi0.trignum = gi1.trignum = 1; /*face index + 1*/
+ gi0.u = pnt0.u;
+ gi0.v = pnt0.v;
+ gi1.u = pnt1.u;
+ gi1.v = pnt1.v;
+ meshing.AddBoundaryElement(node2ngID[uvPtVec[ i ].node],node2ngID[uvPtVec[i+1].node], gi0, gi1);
+ }
+
+ }
+
+ auto noldsurfels = ngMesh.GetNSE();
+ // printf( "noldsurfels %d\n", noldsurfels );
+ // printf("Before Generate the mesh: %d,%d,%d \n", ngMesh.GetNP(), ngMesh.GetNSeg(), ngMesh.GetNSE() );
+ // printf("Maxh %0.1lf\n", mparams.maxh );
+ netgen::MESHING2_RESULT result = meshing.GenerateMesh(ngMesh,mparams,mparams.maxh,1);
+ printf("Mesh Generated: %d,%d,%d \n", ngMesh.GetNP(), ngMesh.GetNSeg(), ngMesh.GetNSE() );
+ ngMesh.CalcSurfacesOfNode();
+
+ return TError();
+}
+
+
// Iterate in all edges of the geometry, and use RestrictLocalHLine function to set the value of allow edges in 1D mesh
SMESH_ComputeErrorPtr
NETGENPlugin_Mesher::AddSegmentsToMesh2(netgen::Mesh& ngMesh,
{
const SMDS_MeshNode * n0 = uvPtVec[ i ].node;
const SMDS_MeshNode * n1 = uvPtVec[ i + 1 ].node;
+
+ bool onEdge = n0->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE && n1->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE;
+
const int posShapeID = n0->GetShapeID();
// skip nodes on degenerated edges
if ( helper.IsDegenShape( posShapeID ) && helper.IsDegenShape( n1->GetShapeID() ))
netgen::Point3d p0(n0->X(), n0->Y(), n0->Z());
netgen::Point3d p1(n1->X(), n1->Y(), n1->Z());
- // ngMesh.RestrictLocalHLine( p0, p1, size, overrideMinH );
- ngMesh.RestrictLocalH( p0, size, overrideMinH );
-
- // if ( node2ngID.count(n0) == 0 )
- // {
- // netgen::MeshPoint mp0( p0 );
- // ngMesh.AddPoint ( mp0, 1, netgen::FIXEDPOINT );
- // node2ngID.insert( make_pair( n0, ngMesh.GetNP() + 1 ) );
- // nodeVec.push_back( n0 );
- // }
- // if ( node2ngID.count(n1) == 0 )
- // {
- // netgen::MeshPoint mp1( p1 );
- // ngMesh.AddPoint ( mp1, 1, netgen::FIXEDPOINT );
- // node2ngID.insert( make_pair( n1, ngMesh.GetNP() + 1 ) );
- // nodeVec.push_back( n1 );
- // }
-
+ if ( onEdge && size > 1e-6 )
+ {
+ ngMesh.RestrictLocalHLine( p0, p1, size, overrideMinH );
+ // ngMesh.RestrictLocalH( p0, size, overrideMinH ); /*similar behavior is observed for this case*/
+
+ }
}
}
if ( onVertex || ( !wasNgMeshEmpty && onEdge ) || helper.IsRealSeam( posShapeID ) )
{
- ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
-
- // printf("OnVertex: %d \n", ngID1 );
- // if ( node2ngID.count( n ) == 0 )
- // node2ngID.insert( make_pair( n, ngID1 ));
- // ngID1 = node2ngID[ n ];
-
- // printf("OnVertex: %d, %0.1lf, %0.1lf, %0.1lf\n", ngID1, n->X(), n->Y(), n->Z() );
+ ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
}
if ( ngID1 > ngMesh.GetNP() )
{
if (aEdge.IsNull())
{
- int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
+ // int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
+ int aGeomEdgeInd = seg.epgeominfo[j].edgenr + 1; /*+1 is needed for netgen6!!!!!*/
+
if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
}
// add a faces descriptor to exclude qudrangle elements generated by NETGEN
// from computation of 3D mesh
ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
+
+ printf("nbInitFac and nbFac %d,%d:\n", nbInitFac, nbFac );
vector<const SMDS_MeshNode*> nodes;
for ( int i = nbInitFac+1; i <= nbFac; ++i )
{
const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
const int aGeomFaceInd = elem.GetIndex();
+
+ printf("Face index %d:\n", aGeomFaceInd );
+
TopoDS_Face aFace;
if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
}
bool NETGENPlugin_Mesher::Fill2DViscousLayer( netgen::OCCGeometry& occgeo, vector< const SMDS_MeshNode* >& nodeVec,
- NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState )
+ NETGENPlugin_Internals* internals, NETGENPlugin_ngMeshInfo& initState,
+ netgen::MeshingParameters &mparams, bool & wasAlreadyComputed )
{
SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
SMESH_Comment comment;
FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
initState = NETGENPlugin_ngMeshInfo(_ngMesh);
}
+
+ // Move init state here BEFORE calling AddSegmentsToMesh3
+ initState = NETGENPlugin_ngMeshInfo(_ngMesh);
+
SMESH_ProxyMesh::Ptr viscousMesh;
SMESH_MesherHelper helper( *_mesh );
for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID )
// add new segments to _ngMesh instead of excluded ones
helper.SetSubShape( F );
TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true, error, &helper, viscousMesh );
- error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
+ error = AddSegmentsToMesh3( *_ngMesh, occgeo, wires, helper, nodeVec, mparams );
+ // error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
if ( !error ) error = SMESH_ComputeError::New();
+
+ wasAlreadyComputed = toCompute;
}
- initState = NETGENPlugin_ngMeshInfo(_ngMesh);
+
+ // initState = NETGENPlugin_ngMeshInfo(_ngMesh);
}
return true;
}
return err;
}
+
void NETGENPlugin_Mesher::InitialSetup( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo,
list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals,
SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState,
}
+
+void NETGENPlugin_Mesher::InitialSetup2( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo,
+ list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals,
+ SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState,
+ netgen::MeshingParameters &mparams )
+{
+ // Init occ geometry maps for non meshed object and fill meshedSM with premeshed objects
+ PrepareOCCgeometry2( occgeo, _shape, *_mesh, meshedSM, internals );
+ _occgeom = &occgeo;
+ _ngMesh = NULL;
+
+ // Fill basic mesh param and
+ // define _ngMesh
+ SetBasicMeshParameters( ngLib, mparams, occgeo );
+
+ // Mesh internal faces and edges with netgen and then fill smesh with those entries!
+ // if ( internals->hasInternalEdges() )
+ // FillInternalElements( ngLib, *internals, occgeo, initState, quadHelper, meshedSM );
+
+}
+
void NETGENPlugin_Mesher::InitialSetupSA( NETGENPlugin_NetgenLibWrapper& ngLib, netgen::OCCGeometry& occgeo,
list< SMESH_subMesh* >* meshedSM, NETGENPlugin_Internals* internals,
SMESH_MesherHelper &quadHelper, NETGENPlugin_ngMeshInfo& initState,
vector< const SMDS_MeshNode* >& nodeVec, SMESH_Comment& comment, DIM dim )
{
SetBasicMeshParametersFor2D(occgeo, nodeVec, mparams, internals, initState );
- if ( !Fill2DViscousLayer(occgeo, nodeVec, internals, initState ) )
+ bool wasAlreadyComputed = false;
+ if ( !Fill2DViscousLayer(occgeo, nodeVec, internals, initState, mparams, wasAlreadyComputed ) )
return false;
mparams.uselocalh = true; // restore as it is used at surface optimization
- int err = CallNetgenMeshFaces( ngLib, occgeo, comment );
+
+ int err = 0;
+ if ( !wasAlreadyComputed )
+ err = CallNetgenMeshFaces( ngLib, occgeo, comment );
if ( !err && dim == DIM::D2 /* if mesh is 3D then second order is defined after volumens are computed*/ )
MakeSecondOrder( mparams, occgeo, meshedSM, initState, comment );
// Generate the mesh
// -------------------------
+ /* Original version
InitialSetup( ngLib, occgeo, meshedSM, &internals, quadHelper, initState, mparams );
err = Fill0D1DElements( occgeo, nodeVec, meshedSM, quadHelper );
+ */
+
+ /* Alternative version where premeshed 1D elements are passed as geometrical vertex of the geometry */
+ InitialSetup2( ngLib, occgeo, meshedSM, &internals, quadHelper, initState, mparams );
+ // err = Fill0D1DElements( occgeo, nodeVec, meshedSM, quadHelper );
+
initState = NETGENPlugin_ngMeshInfo(_ngMesh);
err = CallNetgenMeshEdges( ngLib, occgeo );
SetBasicMeshParametersFor2D(occgeo, nodeVec, mparams, &internals, initState );
- if ( !Fill2DViscousLayer(occgeo, nodeVec, &internals, initState ) )
+ bool wasAlreadyComputed = false;
+ if ( !Fill2DViscousLayer(occgeo, nodeVec, &internals, initState, mparams, wasAlreadyComputed ) )
return false;
SMESH_Comment comment;
- {
-
- if ( _isVolume )
- _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic;
+ if ( _isVolume )
+ _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic;
- mparams.uselocalh = true; // restore as it is used at surface optimization
+ mparams.uselocalh = true; // restore as it is used at surface optimization
+ if ( !wasAlreadyComputed )
err = CallNetgenMeshFaces( ngLib, occgeo, comment );
-
- if ( _isVolume )
- {
- doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 );
- _ticTime = doneTime / _totalTime / _progressTic;
- }
-
- if ( !err )
- if ( !Fill3DViscousLayerAndQuadAdaptor( occgeo, nodeVec, mparams, initState, meshedSM, quadHelper, err ) )
- return false;
+
+ if ( _isVolume )
+ {
+ doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 );
+ _ticTime = doneTime / _totalTime / _progressTic;
+ }
- if (!err && _isVolume)
- {
- SetBasicMeshParametersFor3D( ngLib, occgeo, nodeVec, mparams, &internals, initState, quadHelper, comment );
- err = CallNetgenMeshVolumens( ngLib, occgeo, comment );
- }
+ if ( !err )
+ if ( !Fill3DViscousLayerAndQuadAdaptor( occgeo, nodeVec, mparams, initState, meshedSM, quadHelper, err ) )
+ return false;
- if (!err )
- MakeSecondOrder( mparams, occgeo, meshedSM, initState, comment );
+ if (!err && _isVolume)
+ {
+ SetBasicMeshParametersFor3D( ngLib, occgeo, nodeVec, mparams, &internals, initState, quadHelper, comment );
+ err = CallNetgenMeshVolumens( ngLib, occgeo, comment );
}
+ if (!err )
+ MakeSecondOrder( mparams, occgeo, meshedSM, initState, comment );
+
_ticTime = 0.98 / _progressTic;
//int nbNod = _ngMesh->GetNP();
// netgen::mparam.Print(std::cerr);
#ifdef NETGEN_V6
-
ngMesh->SetGeometry( shared_ptr<netgen::NetgenGeometry>( &occgeo, &NOOP_Deleter ));
netgen::mparam.perfstepsstart = startWith;
netgen::mparam.perfstepsend = endWith;
std::shared_ptr<netgen::Mesh> meshPtr( ngMesh, &NOOP_Deleter );
err = occgeo.GenerateMesh( meshPtr, netgen::mparam );
-
#else
- #ifdef NETGEN_V5
-
- err = netgen::OCCGenerateMesh(occgeo, ngMesh, netgen::mparam, startWith, endWith);
-
- #else
-
char *optstr = 0;
err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
-
- #endif
#endif
return err;
}
