1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File: SMESH_MesherHelper.cxx
23 // Created: 15.02.06 15:22:41
24 // Author: Sergey KUUL
26 #include "SMESH_MesherHelper.hxx"
28 #include "SMDS_FacePosition.hxx"
29 #include "SMDS_EdgePosition.hxx"
31 #include <BRepAdaptor_Surface.hxx>
32 #include <BRepTools.hxx>
33 #include <BRep_Tool.hxx>
34 #include <BRepTools_WireExplorer.hxx>
35 #include <Geom2d_Curve.hxx>
36 #include <Geom_Curve.hxx>
37 #include <Geom_Surface.hxx>
38 #include <ShapeAnalysis.hxx>
40 #include <TopExp_Explorer.hxx>
41 #include <TopTools_ListIteratorOfListOfShape.hxx>
42 #include <TopTools_MapOfShape.hxx>
44 #include <gp_Pnt2d.hxx>
46 #include <Standard_Failure.hxx>
47 #include <Standard_ErrorHandler.hxx>
49 #include <utilities.h>
51 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
53 //================================================================================
57 //================================================================================
59 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
60 : myMesh(&theMesh), myShapeID(-1), myCreateQuadratic(false)
62 mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
65 //=======================================================================
66 //function : CheckShape
68 //=======================================================================
70 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
72 SMESHDS_Mesh* meshDS = GetMeshDS();
73 // we can create quadratic elements only if all elements
74 // created on subshapes of given shape are quadratic
75 // also we have to fill myNLinkNodeMap
76 myCreateQuadratic = true;
77 mySeamShapeIds.clear();
78 myDegenShapeIds.clear();
79 TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
80 SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
82 int nbOldLinks = myNLinkNodeMap.size();
84 TopExp_Explorer exp( aSh, subType );
85 for (; exp.More() && myCreateQuadratic; exp.Next()) {
86 if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
87 if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
89 const SMDS_MeshElement* e = it->next();
90 if ( e->GetType() != elemType || !e->IsQuadratic() ) {
91 myCreateQuadratic = false;
96 switch ( e->NbNodes() ) {
98 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
100 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
101 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
102 AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
104 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
105 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
106 AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
107 AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
110 myCreateQuadratic = false;
119 if ( nbOldLinks == myNLinkNodeMap.size() )
120 myCreateQuadratic = false;
122 if(!myCreateQuadratic) {
123 myNLinkNodeMap.clear();
127 return myCreateQuadratic;
130 //================================================================================
132 * \brief Set geomerty to make elements on
133 * \param aSh - geomertic shape
135 //================================================================================
137 void SMESH_MesherHelper::SetSubShape(const int aShID)
139 if ( aShID == myShapeID )
142 SetSubShape( GetMeshDS()->IndexToShape( aShID ));
144 SetSubShape( TopoDS_Shape() );
147 //================================================================================
149 * \brief Set geomerty to make elements on
150 * \param aSh - geomertic shape
152 //================================================================================
154 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
156 if ( myShape.IsSame( aSh ))
160 mySeamShapeIds.clear();
161 myDegenShapeIds.clear();
163 if ( myShape.IsNull() ) {
167 SMESHDS_Mesh* meshDS = GetMeshDS();
168 myShapeID = meshDS->ShapeToIndex(aSh);
170 // treatment of periodic faces
171 for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
173 const TopoDS_Face& face = TopoDS::Face( eF.Current() );
174 BRepAdaptor_Surface surface( face );
175 if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
177 for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
179 // look for a seam edge
180 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
181 if ( BRep_Tool::IsClosed( edge, face )) {
182 // initialize myPar1, myPar2 and myParIndex
183 if ( mySeamShapeIds.empty() ) {
185 BRep_Tool::UVPoints( edge, face, uv1, uv2 );
186 if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
188 myParIndex = 1; // U periodic
189 myPar1 = surface.FirstUParameter();
190 myPar2 = surface.LastUParameter();
193 myParIndex = 2; // V periodic
194 myPar1 = surface.FirstVParameter();
195 myPar2 = surface.LastVParameter();
198 // store seam shape indices, negative if shape encounters twice
199 int edgeID = meshDS->ShapeToIndex( edge );
200 mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
201 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
202 int vertexID = meshDS->ShapeToIndex( v.Current() );
203 mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
207 // look for a degenerated edge
208 if ( BRep_Tool::Degenerated( edge )) {
209 myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
210 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
211 myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
218 //================================================================================
220 * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
221 * \param F - the face
222 * \retval bool - return true if the face is periodic
224 //================================================================================
226 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
228 if ( F.IsNull() ) return !mySeamShapeIds.empty();
230 if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
231 return !mySeamShapeIds.empty();
233 Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F );
234 if ( !aSurface.IsNull() )
235 return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
240 //=======================================================================
241 //function : IsMedium
243 //=======================================================================
245 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode* node,
246 const SMDSAbs_ElementType typeToCheck)
248 return SMESH_MeshEditor::IsMedium( node, typeToCheck );
251 //=======================================================================
252 //function : AddNLinkNode
254 //=======================================================================
256 * Auxilary function for filling myNLinkNodeMap
258 void SMESH_MesherHelper::AddNLinkNode(const SMDS_MeshNode* n1,
259 const SMDS_MeshNode* n2,
260 const SMDS_MeshNode* n12)
262 NLink link( n1, n2 );
263 if ( n1 > n2 ) link = NLink( n2, n1 );
264 // add new record to map
265 myNLinkNodeMap.insert( make_pair(link,n12));
268 //=======================================================================
270 * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
271 * \param uv1 - UV on the seam
272 * \param uv2 - UV within a face
273 * \retval gp_Pnt2d - selected UV
275 //=======================================================================
277 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
279 double p1 = uv1.Coord( myParIndex );
280 double p2 = uv2.Coord( myParIndex );
281 double p3 = ( Abs( p1 - myPar1 ) < Abs( p1 - myPar2 )) ? myPar2 : myPar1;
282 if ( Abs( p2 - p1 ) > Abs( p2 - p3 ))
284 gp_Pnt2d result = uv1;
285 result.SetCoord( myParIndex, p1 );
289 //=======================================================================
291 * \brief Return node UV on face
292 * \param F - the face
293 * \param n - the node
294 * \param n2 - a node of element being created located inside a face
295 * \retval gp_XY - resulting UV
297 * Auxilary function called form GetMediumNode()
299 //=======================================================================
301 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F,
302 const SMDS_MeshNode* n,
303 const SMDS_MeshNode* n2) const
305 gp_Pnt2d uv( 1e100, 1e100 );
306 const SMDS_PositionPtr Pos = n->GetPosition();
307 if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
309 // node has position on face
310 const SMDS_FacePosition* fpos =
311 static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
312 uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
314 else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
316 // node has position on edge => it is needed to find
317 // corresponding edge from face, get pcurve for this
318 // edge and recieve value from this pcurve
319 const SMDS_EdgePosition* epos =
320 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
321 SMESHDS_Mesh* meshDS = GetMeshDS();
322 int edgeID = Pos->GetShapeId();
323 TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
325 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
326 uv = C2d->Value( epos->GetUParameter() );
327 // for a node on a seam edge select one of UVs on 2 pcurves
328 if ( n2 && IsSeamShape( edgeID ) )
329 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
331 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
333 if ( int vertexID = n->GetPosition()->GetShapeId() ) {
335 const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
337 uv = BRep_Tool::Parameters( V, F );
339 catch (Standard_Failure& exc) {
343 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() )
344 ok = ( V == vert.Current() );
347 MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
348 << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
350 // get UV of a vertex closest to the node
352 gp_Pnt pn ( n->X(),n->Y(),n->Z() );
353 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() ) {
354 TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
355 gp_Pnt p = BRep_Tool::Pnt( curV );
356 double curDist = p.SquareDistance( pn );
357 if ( curDist < dist ) {
359 uv = BRep_Tool::Parameters( curV, F );
360 if ( dist < DBL_MIN ) break;
365 TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
366 for ( ; it.More(); it.Next() ) {
367 if ( it.Value().ShapeType() == TopAbs_EDGE ) {
368 const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
370 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
371 if ( !C2d.IsNull() ) {
372 double u = ( V == TopExp::FirstVertex( edge ) ) ? f : l;
373 uv = C2d->Value( u );
380 if ( n2 && IsSeamShape( vertexID ) )
381 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
387 //=======================================================================
389 * \brief Return node U on edge
390 * \param E - the Edge
391 * \param n - the node
392 * \retval double - resulting U
394 * Auxilary function called form GetMediumNode()
396 //=======================================================================
398 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E,
399 const SMDS_MeshNode* n)
402 const SMDS_PositionPtr Pos = n->GetPosition();
403 if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
404 const SMDS_EdgePosition* epos =
405 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
406 param = epos->GetUParameter();
408 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
409 SMESHDS_Mesh * meshDS = GetMeshDS();
410 int vertexID = n->GetPosition()->GetShapeId();
411 const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
412 param = BRep_Tool::Parameter( V, E );
417 //=======================================================================
418 //function : GetMediumNode
420 //=======================================================================
422 * Special function for search or creation medium node
424 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
425 const SMDS_MeshNode* n2,
428 TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
430 NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
431 ItNLinkNode itLN = myNLinkNodeMap.find( link );
432 if ( itLN != myNLinkNodeMap.end() ) {
433 return (*itLN).second;
436 // create medium node
438 SMESHDS_Mesh* meshDS = GetMeshDS();
439 int faceID = -1, edgeID = -1;
440 const SMDS_PositionPtr Pos1 = n1->GetPosition();
441 const SMDS_PositionPtr Pos2 = n2->GetPosition();
443 if( myShape.IsNull() )
445 if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
446 faceID = Pos1->GetShapeId();
448 else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
449 faceID = Pos2->GetShapeId();
452 if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
453 edgeID = Pos1->GetShapeId();
455 if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
456 edgeID = Pos2->GetShapeId();
461 // we try to create medium node using UV parameters of
462 // nodes, else - medium between corresponding 3d points
463 if(faceID>-1 || shapeType == TopAbs_FACE) {
464 // obtaining a face and 2d points for nodes
466 if( myShape.IsNull() )
467 F = TopoDS::Face(meshDS->IndexToShape(faceID));
469 F = TopoDS::Face(myShape);
473 gp_XY p1 = GetNodeUV(F,n1,n2);
474 gp_XY p2 = GetNodeUV(F,n2,n1);
476 if ( IsDegenShape( Pos1->GetShapeId() ))
477 p1.SetCoord( myParIndex, p2.Coord( myParIndex ));
478 else if ( IsDegenShape( Pos2->GetShapeId() ))
479 p2.SetCoord( myParIndex, p1.Coord( myParIndex ));
481 //checking if surface is periodic
482 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
483 Standard_Real UF,UL,VF,VL;
484 S->Bounds(UF,UL,VF,VL);
487 Standard_Boolean isUPeriodic = S->IsUPeriodic();
489 Standard_Real UPeriod = S->UPeriod();
490 Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
491 Standard_Real pmid = (p1.X()+p2x)/2.;
492 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
495 u= (p1.X()+p2.X())/2.;
497 Standard_Boolean isVPeriodic = S->IsVPeriodic();
499 Standard_Real VPeriod = S->VPeriod();
500 Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
501 Standard_Real pmid = (p1.Y()+p2y)/2.;
502 v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
505 v = (p1.Y()+p2.Y())/2.;
507 gp_Pnt P = S->Value(u, v);
508 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
509 meshDS->SetNodeOnFace(n12, faceID, u, v);
510 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
513 if (edgeID>-1 || shapeType == TopAbs_EDGE) {
516 if( myShape.IsNull() )
517 E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
519 E = TopoDS::Edge(myShape);
523 double p1 = GetNodeU(E,n1);
524 double p2 = GetNodeU(E,n2);
527 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
530 Standard_Boolean isPeriodic = C->IsPeriodic();
533 Standard_Real Period = C->Period();
534 Standard_Real p = p2+ShapeAnalysis::AdjustByPeriod(p2,p1,Period);
535 Standard_Real pmid = (p1+p)/2.;
536 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
541 gp_Pnt P = C->Value( u );
542 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
543 meshDS->SetNodeOnEdge(n12, edgeID, u);
544 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
550 double x = ( n1->X() + n2->X() )/2.;
551 double y = ( n1->Y() + n2->Y() )/2.;
552 double z = ( n1->Z() + n2->Z() )/2.;
553 n12 = meshDS->AddNode(x,y,z);
555 meshDS->SetNodeOnEdge(n12, edgeID);
557 meshDS->SetNodeOnFace(n12, faceID);
559 meshDS->SetNodeInVolume(n12, myShapeID);
560 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
565 //=======================================================================
569 //=======================================================================
571 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
573 SMESHDS_Mesh * meshDS = GetMeshDS();
574 SMDS_MeshNode* node = 0;
576 node = meshDS->AddNodeWithID( x, y, z, ID );
578 node = meshDS->AddNode( x, y, z );
579 if ( mySetElemOnShape && myShapeID > 0 ) {
580 switch ( myShape.ShapeType() ) {
581 case TopAbs_SOLID: meshDS->SetNodeInVolume( node, myShapeID); break;
582 case TopAbs_SHELL: meshDS->SetNodeInVolume( node, myShapeID); break;
583 case TopAbs_FACE: meshDS->SetNodeOnFace( node, myShapeID); break;
584 case TopAbs_EDGE: meshDS->SetNodeOnEdge( node, myShapeID); break;
585 case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
592 //=======================================================================
594 * Creates quadratic or linear edge
596 //=======================================================================
598 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
599 const SMDS_MeshNode* n2,
603 SMESHDS_Mesh * meshDS = GetMeshDS();
605 SMDS_MeshEdge* edge = 0;
606 if (myCreateQuadratic) {
607 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
609 edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
611 edge = meshDS->AddEdge(n1, n2, n12);
615 edge = meshDS->AddEdgeWithID(n1, n2, id);
617 edge = meshDS->AddEdge(n1, n2);
620 if ( mySetElemOnShape && myShapeID > 0 )
621 meshDS->SetMeshElementOnShape( edge, myShapeID );
626 //=======================================================================
628 * Creates quadratic or linear triangle
630 //=======================================================================
632 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
633 const SMDS_MeshNode* n2,
634 const SMDS_MeshNode* n3,
638 SMESHDS_Mesh * meshDS = GetMeshDS();
639 SMDS_MeshFace* elem = 0;
640 if(!myCreateQuadratic) {
642 elem = meshDS->AddFaceWithID(n1, n2, n3, id);
644 elem = meshDS->AddFace(n1, n2, n3);
647 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
648 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
649 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
652 elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
654 elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
656 if ( mySetElemOnShape && myShapeID > 0 )
657 meshDS->SetMeshElementOnShape( elem, myShapeID );
662 //=======================================================================
664 * Creates quadratic or linear quadrangle
666 //=======================================================================
668 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
669 const SMDS_MeshNode* n2,
670 const SMDS_MeshNode* n3,
671 const SMDS_MeshNode* n4,
675 SMESHDS_Mesh * meshDS = GetMeshDS();
676 SMDS_MeshFace* elem = 0;
677 if(!myCreateQuadratic) {
679 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
681 elem = meshDS->AddFace(n1, n2, n3, n4);
684 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
685 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
686 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
687 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
690 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
692 elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
694 if ( mySetElemOnShape && myShapeID > 0 )
695 meshDS->SetMeshElementOnShape( elem, myShapeID );
700 //=======================================================================
702 * Creates quadratic or linear volume
704 //=======================================================================
706 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
707 const SMDS_MeshNode* n2,
708 const SMDS_MeshNode* n3,
709 const SMDS_MeshNode* n4,
710 const SMDS_MeshNode* n5,
711 const SMDS_MeshNode* n6,
715 SMESHDS_Mesh * meshDS = GetMeshDS();
716 SMDS_MeshVolume* elem = 0;
717 if(!myCreateQuadratic) {
719 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
721 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
724 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
725 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
726 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
728 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
729 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
730 const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
732 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
733 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
734 const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
737 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
738 n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
740 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
741 n12, n23, n31, n45, n56, n64, n14, n25, n36);
743 if ( mySetElemOnShape && myShapeID > 0 )
744 meshDS->SetMeshElementOnShape( elem, myShapeID );
749 //=======================================================================
751 * Creates quadratic or linear volume
753 //=======================================================================
755 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
756 const SMDS_MeshNode* n2,
757 const SMDS_MeshNode* n3,
758 const SMDS_MeshNode* n4,
762 SMESHDS_Mesh * meshDS = GetMeshDS();
763 SMDS_MeshVolume* elem = 0;
764 if(!myCreateQuadratic) {
766 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
768 elem = meshDS->AddVolume(n1, n2, n3, n4);
771 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
772 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
773 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
775 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
776 const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
777 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
780 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
782 elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
784 if ( mySetElemOnShape && myShapeID > 0 )
785 meshDS->SetMeshElementOnShape( elem, myShapeID );
790 //=======================================================================
792 * Creates quadratic or linear pyramid
794 //=======================================================================
796 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
797 const SMDS_MeshNode* n2,
798 const SMDS_MeshNode* n3,
799 const SMDS_MeshNode* n4,
800 const SMDS_MeshNode* n5,
804 SMDS_MeshVolume* elem = 0;
805 if(!myCreateQuadratic) {
807 elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
809 elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
812 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
813 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
814 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
815 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
817 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
818 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
819 const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
820 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
823 elem = GetMeshDS()->AddVolumeWithID ( n1, n2, n3, n4, n5,
828 elem = GetMeshDS()->AddVolume( n1, n2, n3, n4, n5,
832 if ( mySetElemOnShape && myShapeID > 0 )
833 GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
838 //=======================================================================
840 * Creates quadratic or linear hexahedron
842 //=======================================================================
844 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
845 const SMDS_MeshNode* n2,
846 const SMDS_MeshNode* n3,
847 const SMDS_MeshNode* n4,
848 const SMDS_MeshNode* n5,
849 const SMDS_MeshNode* n6,
850 const SMDS_MeshNode* n7,
851 const SMDS_MeshNode* n8,
855 SMESHDS_Mesh * meshDS = GetMeshDS();
856 SMDS_MeshVolume* elem = 0;
857 if(!myCreateQuadratic) {
859 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
861 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
864 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
865 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
866 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
867 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
869 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
870 const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
871 const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
872 const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
874 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
875 const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
876 const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
877 const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
880 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
881 n12, n23, n34, n41, n56, n67,
882 n78, n85, n15, n26, n37, n48, id);
884 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
885 n12, n23, n34, n41, n56, n67,
886 n78, n85, n15, n26, n37, n48);
888 if ( mySetElemOnShape && myShapeID > 0 )
889 meshDS->SetMeshElementOnShape( elem, myShapeID );
894 //=======================================================================
896 * \brief Load nodes bound to face into a map of node columns
897 * \param theParam2ColumnMap - map of node columns to fill
898 * \param theFace - the face on which nodes are searched for
899 * \param theBaseEdge - the edge nodes of which are columns' bases
900 * \param theMesh - the mesh containing nodes
901 * \retval bool - false if something is wrong
903 * The key of the map is a normalized parameter of each
904 * base node on theBaseEdge.
905 * This method works in supposition that nodes on the face
906 * forms a rectangular grid and elements can be quardrangles or triangles
908 //=======================================================================
910 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
911 const TopoDS_Face& theFace,
912 const TopoDS_Edge& theBaseEdge,
913 SMESHDS_Mesh* theMesh)
915 // get vertices of theBaseEdge
916 TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
917 TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
918 TopExp::Vertices( eFrw, vfb, vlb );
920 // find the other edges of theFace and orientation of e1
921 TopoDS_Edge e1, e2, eTop;
922 bool rev1, CumOri = false;
923 TopExp_Explorer exp( theFace, TopAbs_EDGE );
925 for ( ; exp.More(); exp.Next() ) {
926 if ( ++nbEdges > 4 ) {
927 return false; // more than 4 edges in theFace
929 TopoDS_Edge e = TopoDS::Edge( exp.Current() );
930 if ( theBaseEdge.IsSame( e ))
932 TopoDS_Vertex vCommon;
933 if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
935 else if ( vCommon.IsSame( vfb )) {
937 vft = TopExp::LastVertex( e1, CumOri );
938 rev1 = vfb.IsSame( vft );
940 vft = TopExp::FirstVertex( e1, CumOri );
946 return false; // less than 4 edges in theFace
948 if ( e2.IsNull() && vfb.IsSame( vlb ))
951 // submeshes corresponding to shapes
952 SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
953 SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
954 SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
955 SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
956 SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
957 SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
958 SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
959 SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
960 if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
961 RETURN_BAD_RESULT( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
962 sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
964 if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
965 RETURN_BAD_RESULT(" Diff nb of nodes on opposite edges" );
967 if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
968 RETURN_BAD_RESULT("Empty submesh of vertex");
970 // define whether mesh is quadratic
971 bool isQuadraticMesh = false;
972 SMDS_ElemIteratorPtr eIt = smFace->GetElements();
973 if ( !eIt->more() ) {
974 RETURN_BAD_RESULT("No elements on the face");
976 const SMDS_MeshElement* e = eIt->next();
977 isQuadraticMesh = e->IsQuadratic();
979 if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
980 // check quadratic case
981 if ( isQuadraticMesh ) {
982 // what if there are quadrangles and triangles mixed?
983 // int n1 = sm1->NbNodes()/2;
984 // int n2 = smb->NbNodes()/2;
985 // int n3 = sm1->NbNodes() - n1;
986 // int n4 = smb->NbNodes() - n2;
987 // int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
988 // if( nf != smFace->NbNodes() ) {
989 // MESSAGE( "Wrong nb face nodes: " <<
990 // sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
995 RETURN_BAD_RESULT( "Wrong nb face nodes: " <<
996 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
1000 int vsize = sm1->NbNodes() + 2;
1001 int hsize = smb->NbNodes() + 2;
1002 if(isQuadraticMesh) {
1003 vsize = vsize - sm1->NbNodes()/2 -1;
1004 hsize = hsize - smb->NbNodes()/2 -1;
1007 // load nodes from theBaseEdge
1009 std::set<const SMDS_MeshNode*> loadedNodes;
1010 const SMDS_MeshNode* nullNode = 0;
1012 std::vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
1013 nVecf.resize( vsize, nullNode );
1014 loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
1016 std::vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
1017 nVecl.resize( vsize, nullNode );
1018 loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
1021 BRep_Tool::Range( eFrw, f, l );
1022 double range = l - f;
1023 SMDS_NodeIteratorPtr nIt = smb->GetNodes();
1024 const SMDS_MeshNode* node;
1025 while ( nIt->more() ) {
1027 if(IsMedium(node, SMDSAbs_Edge))
1029 const SMDS_EdgePosition* pos =
1030 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1034 double u = ( pos->GetUParameter() - f ) / range;
1035 std::vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
1036 nVec.resize( vsize, nullNode );
1037 loadedNodes.insert( nVec[ 0 ] = node );
1039 if ( theParam2ColumnMap.size() != hsize ) {
1040 RETURN_BAD_RESULT( "Wrong node positions on theBaseEdge" );
1043 // load nodes from e1
1045 std::map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
1046 nIt = sm1->GetNodes();
1047 while ( nIt->more() ) {
1051 const SMDS_EdgePosition* pos =
1052 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1056 sortedNodes.insert( std::make_pair( pos->GetUParameter(), node ));
1058 loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
1059 std::map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
1060 int row = rev1 ? vsize - 1 : 0;
1061 int dRow = rev1 ? -1 : +1;
1062 for ( ; u_n != sortedNodes.end(); u_n++ ) {
1064 loadedNodes.insert( nVecf[ row ] = u_n->second );
1067 // try to load the rest nodes
1069 // get all faces from theFace
1070 TIDSortedElemSet allFaces, foundFaces;
1071 eIt = smFace->GetElements();
1072 while ( eIt->more() ) {
1073 const SMDS_MeshElement* e = eIt->next();
1074 if ( e->GetType() == SMDSAbs_Face )
1075 allFaces.insert( e );
1077 // Starting from 2 neighbour nodes on theBaseEdge, look for a face
1078 // the nodes belong to, and between the nodes of the found face,
1079 // look for a not loaded node considering this node to be the next
1080 // in a column of the starting second node. Repeat, starting
1081 // from nodes next to the previous starting nodes in their columns,
1082 // and so on while a face can be found. Then go the the next pair
1083 // of nodes on theBaseEdge.
1084 TParam2ColumnMap::iterator par_nVec_1 = theParam2ColumnMap.begin();
1085 TParam2ColumnMap::iterator par_nVec_2 = par_nVec_1;
1088 for ( par_nVec_2++; par_nVec_2 != theParam2ColumnMap.end(); par_nVec_1++, par_nVec_2++ ) {
1091 const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
1092 const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
1093 const SMDS_MeshElement* face = 0;
1094 bool lastColOnClosedFace = ( nVecf[ row ] == n2 );
1096 // look for a face by 2 nodes
1097 face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
1099 int nbFaceNodes = face->NbNodes();
1100 if ( face->IsQuadratic() )
1102 if ( nbFaceNodes>4 ) {
1103 RETURN_BAD_RESULT(" Too many nodes in a face: " << nbFaceNodes );
1105 // look for a not loaded node of the <face>
1107 const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
1108 for ( int i = 0; i < nbFaceNodes && !found; ++i ) {
1109 node = face->GetNode( i );
1110 found = loadedNodes.insert( node ).second;
1111 if ( !found && node != n1 && node != n2 )
1114 if ( lastColOnClosedFace && row + 1 < vsize ) {
1115 node = nVecf[ row + 1 ];
1116 found = ( face->GetNodeIndex( node ) >= 0 );
1119 if ( ++row > vsize - 1 ) {
1120 RETURN_BAD_RESULT( "Too many nodes in column "<< col <<": "<< row+1);
1122 par_nVec_2->second[ row ] = node;
1123 foundFaces.insert( face );
1125 if ( nbFaceNodes==4 ) {
1126 n1 = par_nVec_1->second[ row ];
1129 else if ( nbFaceNodes==3 && n3 == par_nVec_1->second[ row + 1 ] ) {
1133 RETURN_BAD_RESULT( "Not quad mesh, column "<< col );
1137 while ( face && n1 && n2 );
1139 if ( row < vsize - 1 ) {
1140 MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
1141 MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
1142 MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
1143 if ( n1 ) { MESSAGE( "Current node 1: "<< n1); }
1144 else { MESSAGE( "Current node 1: NULL"); }
1145 if ( n2 ) { MESSAGE( "Current node 2: "<< n2); }
1146 else { MESSAGE( "Current node 2: NULL"); }
1147 MESSAGE( "first base node: "<< theParam2ColumnMap.begin()->second[0]);
1148 MESSAGE( "last base node: "<< theParam2ColumnMap.rbegin()->second[0]);
1151 } // loop on columns
1156 //=======================================================================
1158 * \brief Return number of unique ancestors of the shape
1160 //=======================================================================
1162 int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
1163 const SMESH_Mesh& mesh,
1164 TopAbs_ShapeEnum ancestorType/*=TopAbs_SHAPE*/)
1166 TopTools_MapOfShape ancestors;
1167 TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
1168 for ( ; ansIt.More(); ansIt.Next() ) {
1169 if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
1170 ancestors.Add( ansIt.Value() );
1172 return ancestors.Extent();
1175 //=======================================================================
1177 * Check mesh without geometry for: if all elements on this shape are quadratic,
1178 * quadratic elements will be created.
1179 * Used then generated 3D mesh without geometry.
1181 //=======================================================================
1183 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
1185 int NbAllEdgsAndFaces=0;
1186 int NbQuadFacesAndEdgs=0;
1187 int NbFacesAndEdges=0;
1188 //All faces and edges
1189 NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
1191 //Quadratic faces and edges
1192 NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
1194 //Linear faces and edges
1195 NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
1197 if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
1199 return SMESH_MesherHelper::QUADRATIC;
1201 else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
1203 return SMESH_MesherHelper::LINEAR;
1206 //Mesh with both type of elements
1207 return SMESH_MesherHelper::COMP;
1210 //=======================================================================
1212 * \brief Return an alternative parameter for a node on seam
1214 //=======================================================================
1216 double SMESH_MesherHelper::GetOtherParam(const double param) const
1218 return fabs(param-myPar1) < fabs(param-myPar2) ? myPar2 : myPar1;