1 // Copyright (C) 2005 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License.
9 // This library is distributed in the hope that it will be useful
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File: SMESH_MesherHelper.cxx
21 // Created: 15.02.06 15:22:41
22 // Author: Sergey KUUL
23 // Copyright: Open CASCADE 2006
26 #include "SMESH_MesherHelper.hxx"
28 #include "SMDS_FacePosition.hxx"
29 #include "SMDS_EdgePosition.hxx"
30 #include "SMESH_MeshEditor.hxx"
32 #include <BRepAdaptor_Surface.hxx>
33 #include <BRepTools.hxx>
34 #include <BRep_Tool.hxx>
35 #include <BRepTools_WireExplorer.hxx>
36 #include <Geom2d_Curve.hxx>
37 #include <Geom_Curve.hxx>
38 #include <Geom_Surface.hxx>
39 #include <ShapeAnalysis.hxx>
41 #include <TopExp_Explorer.hxx>
42 #include <TopTools_ListIteratorOfListOfShape.hxx>
43 #include <TopTools_MapOfShape.hxx>
45 #include <gp_Pnt2d.hxx>
47 #include <Standard_Failure.hxx>
48 #include <Standard_ErrorHandler.hxx>
50 #include <utilities.h>
52 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
54 //================================================================================
58 //================================================================================
60 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
61 : myMesh(&theMesh), myShapeID(-1), myCreateQuadratic(false)
63 mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
66 //=======================================================================
67 //function : CheckShape
69 //=======================================================================
71 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
73 SMESHDS_Mesh* meshDS = GetMeshDS();
74 // we can create quadratic elements only if all elements
75 // created on subshapes of given shape are quadratic
76 // also we have to fill myNLinkNodeMap
77 myCreateQuadratic = true;
78 mySeamShapeIds.clear();
79 myDegenShapeIds.clear();
80 TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
81 SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
83 TopExp_Explorer exp( aSh, subType );
84 for (; exp.More() && myCreateQuadratic; exp.Next()) {
85 if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
86 if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
88 const SMDS_MeshElement* e = it->next();
89 if ( e->GetType() != elemType || !e->IsQuadratic() ) {
90 myCreateQuadratic = false;
95 switch ( e->NbNodes() ) {
97 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
99 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
100 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
101 AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
103 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
104 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
105 AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
106 AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
109 myCreateQuadratic = false;
118 if(!myCreateQuadratic) {
119 myNLinkNodeMap.clear();
123 return myCreateQuadratic;
126 //================================================================================
128 * \brief Set geomerty to make elements on
129 * \param aSh - geomertic shape
131 //================================================================================
133 void SMESH_MesherHelper::SetSubShape(const int aShID)
135 if ( aShID == myShapeID )
138 SetSubShape( GetMeshDS()->IndexToShape( aShID ));
140 SetSubShape( TopoDS_Shape() );
143 //================================================================================
145 * \brief Set geomerty to make elements on
146 * \param aSh - geomertic shape
148 //================================================================================
150 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
152 if ( myShape.IsSame( aSh ))
156 mySeamShapeIds.clear();
157 myDegenShapeIds.clear();
159 if ( myShape.IsNull() ) {
163 SMESHDS_Mesh* meshDS = GetMeshDS();
164 myShapeID = meshDS->ShapeToIndex(aSh);
166 // treatment of periodic faces
167 for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
169 const TopoDS_Face& face = TopoDS::Face( eF.Current() );
170 BRepAdaptor_Surface surface( face );
171 if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
173 for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
175 // look for a seam edge
176 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
177 bool isClosed = BRep_Tool::IsClosed( edge, face );
178 // BEGIN: jfa for bug 0019943
180 MESSAGE("$$$ CLOSED 1 $$$")
182 for (TopExp_Explorer expw (face, TopAbs_WIRE); expw.More() && !isClosed; expw.Next()) {
183 const TopoDS_Wire& wire = TopoDS::Wire(expw.Current());
185 for (BRepTools_WireExplorer we (wire, face); we.More() && !isClosed; we.Next()) {
186 if (we.Current().IsSame(edge)) {
188 if (nbe == 2) isClosed = true;
193 // END: jfa for bug 0019943
195 MESSAGE("$$$ CLOSED 2 $$$")
196 // initialize myPar1, myPar2 and myParIndex
197 if ( mySeamShapeIds.empty() ) {
199 BRep_Tool::UVPoints( edge, face, uv1, uv2 );
200 if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
202 myParIndex = 1; // U periodic
203 myPar1 = surface.FirstUParameter();
204 myPar2 = surface.LastUParameter();
207 myParIndex = 2; // V periodic
208 myPar1 = surface.FirstVParameter();
209 myPar2 = surface.LastVParameter();
212 // store shapes indices
213 mySeamShapeIds.insert( meshDS->ShapeToIndex( edge ));
214 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
215 mySeamShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
218 // look for a degenerated edge
219 if ( BRep_Tool::Degenerated( edge )) {
220 myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
221 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
222 myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
229 //================================================================================
231 * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
232 * \param F - the face
233 * \retval bool - return true if the face is periodic
235 //================================================================================
237 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
239 if ( F.IsNull() ) return !mySeamShapeIds.empty();
241 if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
242 return !mySeamShapeIds.empty();
244 Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F );
245 if ( !aSurface.IsNull() )
246 return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
251 //=======================================================================
252 //function : IsMedium
254 //=======================================================================
256 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode* node,
257 const SMDSAbs_ElementType typeToCheck)
259 return SMESH_MeshEditor::IsMedium( node, typeToCheck );
262 //=======================================================================
263 //function : AddNLinkNode
265 //=======================================================================
267 * Auxilary function for filling myNLinkNodeMap
269 void SMESH_MesherHelper::AddNLinkNode(const SMDS_MeshNode* n1,
270 const SMDS_MeshNode* n2,
271 const SMDS_MeshNode* n12)
273 NLink link( n1, n2 );
274 if ( n1 > n2 ) link = NLink( n2, n1 );
275 // add new record to map
276 myNLinkNodeMap.insert( make_pair(link,n12));
279 //=======================================================================
281 * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
282 * \param uv1 - UV on the seam
283 * \param uv2 - UV within a face
284 * \retval gp_Pnt2d - selected UV
286 //=======================================================================
288 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
290 double p1 = uv1.Coord( myParIndex );
291 double p2 = uv2.Coord( myParIndex );
292 double p3 = ( Abs( p1 - myPar1 ) < Abs( p1 - myPar2 )) ? myPar2 : myPar1;
293 if ( Abs( p2 - p1 ) > Abs( p2 - p3 ))
295 gp_Pnt2d result = uv1;
296 result.SetCoord( myParIndex, p1 );
300 //=======================================================================
302 * \brief Return node UV on face
303 * \param F - the face
304 * \param n - the node
305 * \param n2 - a node of element being created located inside a face
306 * \retval gp_XY - resulting UV
308 * Auxilary function called form GetMediumNode()
310 //=======================================================================
312 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F,
313 const SMDS_MeshNode* n,
314 const SMDS_MeshNode* n2) const
316 gp_Pnt2d uv( 1e100, 1e100 );
317 const SMDS_PositionPtr Pos = n->GetPosition();
318 if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
320 // node has position on face
321 const SMDS_FacePosition* fpos =
322 static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
323 uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
325 else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
327 // node has position on edge => it is needed to find
328 // corresponding edge from face, get pcurve for this
329 // edge and recieve value from this pcurve
330 const SMDS_EdgePosition* epos =
331 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
332 SMESHDS_Mesh* meshDS = GetMeshDS();
333 int edgeID = Pos->GetShapeId();
334 TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
337 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
338 uv = C2d->Value( epos->GetUParameter() );
339 // for a node on a seam edge select one of UVs on 2 pcurves
340 if ( n2 && mySeamShapeIds.find( edgeID ) != mySeamShapeIds.end() )
341 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
343 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
345 if ( int vertexID = n->GetPosition()->GetShapeId() ) {
347 const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
349 uv = BRep_Tool::Parameters( V, F );
351 catch (Standard_Failure& exc) {
355 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() )
356 ok = ( V == vert.Current() );
359 MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
360 << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
362 // get UV of a vertex closest to the node
364 gp_Pnt pn ( n->X(),n->Y(),n->Z() );
365 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() ) {
366 TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
367 gp_Pnt p = BRep_Tool::Pnt( curV );
368 double curDist = p.SquareDistance( pn );
369 if ( curDist < dist ) {
371 uv = BRep_Tool::Parameters( curV, F );
372 if ( dist < DBL_MIN ) break;
377 TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
378 for ( ; it.More(); it.Next() ) {
379 if ( it.Value().ShapeType() == TopAbs_EDGE ) {
380 const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
382 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
383 if ( !C2d.IsNull() ) {
384 double u = ( V == TopExp::FirstVertex( edge ) ) ? f : l;
385 uv = C2d->Value( u );
392 if ( n2 && mySeamShapeIds.find( vertexID ) != mySeamShapeIds.end() )
393 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
399 //=======================================================================
401 * \brief Return node U on edge
402 * \param E - the Edge
403 * \param n - the node
404 * \retval double - resulting U
406 * Auxilary function called form GetMediumNode()
408 //=======================================================================
410 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E,
411 const SMDS_MeshNode* n)
414 const SMDS_PositionPtr Pos = n->GetPosition();
415 if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
416 const SMDS_EdgePosition* epos =
417 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
418 param = epos->GetUParameter();
420 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
421 SMESHDS_Mesh * meshDS = GetMeshDS();
422 int vertexID = n->GetPosition()->GetShapeId();
423 const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
424 param = BRep_Tool::Parameter( V, E );
429 //=======================================================================
430 //function : GetMediumNode
432 //=======================================================================
434 * Special function for search or creation medium node
436 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
437 const SMDS_MeshNode* n2,
440 TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
442 NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
443 ItNLinkNode itLN = myNLinkNodeMap.find( link );
444 if ( itLN != myNLinkNodeMap.end() ) {
445 return (*itLN).second;
448 // create medium node
450 SMESHDS_Mesh* meshDS = GetMeshDS();
451 int faceID = -1, edgeID = -1;
452 const SMDS_PositionPtr Pos1 = n1->GetPosition();
453 const SMDS_PositionPtr Pos2 = n2->GetPosition();
455 if( myShape.IsNull() )
457 if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
458 faceID = Pos1->GetShapeId();
460 else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
461 faceID = Pos2->GetShapeId();
464 if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
465 edgeID = Pos1->GetShapeId();
467 if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
468 edgeID = Pos2->GetShapeId();
473 // we try to create medium node using UV parameters of
474 // nodes, else - medium between corresponding 3d points
475 if(faceID>-1 || shapeType == TopAbs_FACE) {
476 // obtaining a face and 2d points for nodes
478 if( myShape.IsNull() )
479 F = TopoDS::Face(meshDS->IndexToShape(faceID));
481 F = TopoDS::Face(myShape);
485 gp_XY p1 = GetNodeUV(F,n1,n2);
486 gp_XY p2 = GetNodeUV(F,n2,n1);
488 //checking if surface is periodic
489 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
490 Standard_Real UF,UL,VF,VL;
491 S->Bounds(UF,UL,VF,VL);
494 Standard_Boolean isUPeriodic = S->IsUPeriodic();
496 Standard_Real UPeriod = S->UPeriod();
497 Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
498 Standard_Real pmid = (p1.X()+p2x)/2.;
499 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
502 u= (p1.X()+p2.X())/2.;
504 Standard_Boolean isVPeriodic = S->IsVPeriodic();
506 Standard_Real VPeriod = S->VPeriod();
507 Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
508 Standard_Real pmid = (p1.Y()+p2y)/2.;
509 v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
512 v = (p1.Y()+p2.Y())/2.;
514 gp_Pnt P = S->Value(u, v);
515 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
516 meshDS->SetNodeOnFace(n12, faceID, u, v);
517 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
520 if (edgeID>-1 || shapeType == TopAbs_EDGE) {
523 if( myShape.IsNull() )
524 E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
526 E = TopoDS::Edge(myShape);
530 double p1 = GetNodeU(E,n1);
531 double p2 = GetNodeU(E,n2);
534 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
537 Standard_Boolean isPeriodic = C->IsPeriodic();
540 Standard_Real Period = C->Period();
541 Standard_Real p = p2+ShapeAnalysis::AdjustByPeriod(p2,p1,Period);
542 Standard_Real pmid = (p1+p)/2.;
543 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
548 gp_Pnt P = C->Value( u );
549 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
550 meshDS->SetNodeOnEdge(n12, edgeID, u);
551 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
557 double x = ( n1->X() + n2->X() )/2.;
558 double y = ( n1->Y() + n2->Y() )/2.;
559 double z = ( n1->Z() + n2->Z() )/2.;
560 n12 = meshDS->AddNode(x,y,z);
562 meshDS->SetNodeOnEdge(n12, edgeID);
564 meshDS->SetNodeOnFace(n12, faceID);
566 meshDS->SetNodeInVolume(n12, myShapeID);
567 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
572 //=======================================================================
576 //=======================================================================
578 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
580 SMESHDS_Mesh * meshDS = GetMeshDS();
581 SMDS_MeshNode* node = 0;
583 node = meshDS->AddNodeWithID( x, y, z, ID );
585 node = meshDS->AddNode( x, y, z );
586 if ( mySetElemOnShape && myShapeID > 0 ) {
587 switch ( myShape.ShapeType() ) {
588 case TopAbs_SOLID: meshDS->SetNodeInVolume( node, myShapeID); break;
589 case TopAbs_SHELL: meshDS->SetNodeInVolume( node, myShapeID); break;
590 case TopAbs_FACE: meshDS->SetNodeOnFace( node, myShapeID); break;
591 case TopAbs_EDGE: meshDS->SetNodeOnEdge( node, myShapeID); break;
592 case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
599 //=======================================================================
601 * Creates quadratic or linear edge
603 //=======================================================================
605 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
606 const SMDS_MeshNode* n2,
610 SMESHDS_Mesh * meshDS = GetMeshDS();
612 SMDS_MeshEdge* edge = 0;
613 if (myCreateQuadratic) {
614 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
616 edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
618 edge = meshDS->AddEdge(n1, n2, n12);
622 edge = meshDS->AddEdgeWithID(n1, n2, id);
624 edge = meshDS->AddEdge(n1, n2);
627 if ( mySetElemOnShape && myShapeID > 0 )
628 meshDS->SetMeshElementOnShape( edge, myShapeID );
633 //=======================================================================
635 * Creates quadratic or linear triangle
637 //=======================================================================
639 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
640 const SMDS_MeshNode* n2,
641 const SMDS_MeshNode* n3,
645 SMESHDS_Mesh * meshDS = GetMeshDS();
646 SMDS_MeshFace* elem = 0;
647 if(!myCreateQuadratic) {
649 elem = meshDS->AddFaceWithID(n1, n2, n3, id);
651 elem = meshDS->AddFace(n1, n2, n3);
654 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
655 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
656 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
659 elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
661 elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
663 if ( mySetElemOnShape && myShapeID > 0 )
664 meshDS->SetMeshElementOnShape( elem, myShapeID );
669 //=======================================================================
671 * Creates quadratic or linear quadrangle
673 //=======================================================================
675 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
676 const SMDS_MeshNode* n2,
677 const SMDS_MeshNode* n3,
678 const SMDS_MeshNode* n4,
682 SMESHDS_Mesh * meshDS = GetMeshDS();
683 SMDS_MeshFace* elem = 0;
684 if(!myCreateQuadratic) {
686 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
688 elem = meshDS->AddFace(n1, n2, n3, n4);
691 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
692 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
693 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
694 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
697 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
699 elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
701 if ( mySetElemOnShape && myShapeID > 0 )
702 meshDS->SetMeshElementOnShape( elem, myShapeID );
707 //=======================================================================
709 * Creates quadratic or linear volume
711 //=======================================================================
713 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
714 const SMDS_MeshNode* n2,
715 const SMDS_MeshNode* n3,
716 const SMDS_MeshNode* n4,
717 const SMDS_MeshNode* n5,
718 const SMDS_MeshNode* n6,
722 SMESHDS_Mesh * meshDS = GetMeshDS();
723 SMDS_MeshVolume* elem = 0;
724 if(!myCreateQuadratic) {
726 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
728 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
731 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
732 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
733 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
735 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
736 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
737 const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
739 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
740 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
741 const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
744 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
745 n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
747 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
748 n12, n23, n31, n45, n56, n64, n14, n25, n36);
750 if ( mySetElemOnShape && myShapeID > 0 )
751 meshDS->SetMeshElementOnShape( elem, myShapeID );
756 //=======================================================================
758 * Creates quadratic or linear volume
760 //=======================================================================
762 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
763 const SMDS_MeshNode* n2,
764 const SMDS_MeshNode* n3,
765 const SMDS_MeshNode* n4,
769 SMESHDS_Mesh * meshDS = GetMeshDS();
770 SMDS_MeshVolume* elem = 0;
771 if(!myCreateQuadratic) {
773 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
775 elem = meshDS->AddVolume(n1, n2, n3, n4);
778 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
779 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
780 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
782 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
783 const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
784 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
787 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
789 elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
791 if ( mySetElemOnShape && myShapeID > 0 )
792 meshDS->SetMeshElementOnShape( elem, myShapeID );
797 //=======================================================================
799 * Creates quadratic or linear pyramid
801 //=======================================================================
803 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
804 const SMDS_MeshNode* n2,
805 const SMDS_MeshNode* n3,
806 const SMDS_MeshNode* n4,
807 const SMDS_MeshNode* n5,
811 SMDS_MeshVolume* elem = 0;
812 if(!myCreateQuadratic) {
814 elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
816 elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
819 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
820 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
821 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
822 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
824 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
825 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
826 const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
827 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
830 elem = GetMeshDS()->AddVolumeWithID ( n1, n2, n3, n4, n5,
835 elem = GetMeshDS()->AddVolume( n1, n2, n3, n4, n5,
839 if ( mySetElemOnShape && myShapeID > 0 )
840 GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
845 //=======================================================================
847 * Creates quadratic or linear hexahedron
849 //=======================================================================
851 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
852 const SMDS_MeshNode* n2,
853 const SMDS_MeshNode* n3,
854 const SMDS_MeshNode* n4,
855 const SMDS_MeshNode* n5,
856 const SMDS_MeshNode* n6,
857 const SMDS_MeshNode* n7,
858 const SMDS_MeshNode* n8,
862 SMESHDS_Mesh * meshDS = GetMeshDS();
863 SMDS_MeshVolume* elem = 0;
864 if(!myCreateQuadratic) {
866 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
868 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
871 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
872 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
873 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
874 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
876 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
877 const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
878 const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
879 const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
881 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
882 const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
883 const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
884 const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
887 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
888 n12, n23, n34, n41, n56, n67,
889 n78, n85, n15, n26, n37, n48, id);
891 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
892 n12, n23, n34, n41, n56, n67,
893 n78, n85, n15, n26, n37, n48);
895 if ( mySetElemOnShape && myShapeID > 0 )
896 meshDS->SetMeshElementOnShape( elem, myShapeID );
901 //=======================================================================
903 * \brief Load nodes bound to face into a map of node columns
904 * \param theParam2ColumnMap - map of node columns to fill
905 * \param theFace - the face on which nodes are searched for
906 * \param theBaseEdge - the edge nodes of which are columns' bases
907 * \param theMesh - the mesh containing nodes
908 * \retval bool - false if something is wrong
910 * The key of the map is a normalized parameter of each
911 * base node on theBaseEdge.
912 * This method works in supposition that nodes on the face
913 * forms a rectangular grid and elements can be quardrangles or triangles
915 //=======================================================================
917 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
918 const TopoDS_Face& theFace,
919 const TopoDS_Edge& theBaseEdge,
920 SMESHDS_Mesh* theMesh)
922 // get vertices of theBaseEdge
923 TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
924 TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
925 TopExp::Vertices( eFrw, vfb, vlb );
927 // find the other edges of theFace and orientation of e1
928 TopoDS_Edge e1, e2, eTop;
929 bool rev1, CumOri = false;
930 TopExp_Explorer exp( theFace, TopAbs_EDGE );
932 for ( ; exp.More(); exp.Next() ) {
933 if ( ++nbEdges > 4 ) {
934 return false; // more than 4 edges in theFace
936 TopoDS_Edge e = TopoDS::Edge( exp.Current() );
937 if ( theBaseEdge.IsSame( e ))
939 TopoDS_Vertex vCommon;
940 if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
942 else if ( vCommon.IsSame( vfb )) {
944 vft = TopExp::LastVertex( e1, CumOri );
945 rev1 = vfb.IsSame( vft );
947 vft = TopExp::FirstVertex( e1, CumOri );
953 return false; // less than 4 edges in theFace
955 if ( e2.IsNull() && vfb.IsSame( vlb ))
958 // submeshes corresponding to shapes
959 SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
960 SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
961 SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
962 SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
963 SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
964 SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
965 SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
966 SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
967 if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
968 RETURN_BAD_RESULT( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
969 sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
971 if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
972 RETURN_BAD_RESULT(" Diff nb of nodes on opposite edges" );
974 if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
975 RETURN_BAD_RESULT("Empty submesh of vertex");
977 // define whether mesh is quadratic
978 bool isQuadraticMesh = false;
979 SMDS_ElemIteratorPtr eIt = smFace->GetElements();
980 if ( !eIt->more() ) {
981 RETURN_BAD_RESULT("No elements on the face");
983 const SMDS_MeshElement* e = eIt->next();
984 isQuadraticMesh = e->IsQuadratic();
986 if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
987 // check quadratic case
988 if ( isQuadraticMesh ) {
989 // what if there are quadrangles and triangles mixed?
990 // int n1 = sm1->NbNodes()/2;
991 // int n2 = smb->NbNodes()/2;
992 // int n3 = sm1->NbNodes() - n1;
993 // int n4 = smb->NbNodes() - n2;
994 // int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
995 // if( nf != smFace->NbNodes() ) {
996 // MESSAGE( "Wrong nb face nodes: " <<
997 // sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
1002 RETURN_BAD_RESULT( "Wrong nb face nodes: " <<
1003 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
1007 int vsize = sm1->NbNodes() + 2;
1008 int hsize = smb->NbNodes() + 2;
1009 if(isQuadraticMesh) {
1010 vsize = vsize - sm1->NbNodes()/2 -1;
1011 hsize = hsize - smb->NbNodes()/2 -1;
1014 // load nodes from theBaseEdge
1016 set<const SMDS_MeshNode*> loadedNodes;
1017 const SMDS_MeshNode* nullNode = 0;
1019 vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
1020 nVecf.resize( vsize, nullNode );
1021 loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
1023 vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
1024 nVecl.resize( vsize, nullNode );
1025 loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
1028 BRep_Tool::Range( eFrw, f, l );
1029 double range = l - f;
1030 SMDS_NodeIteratorPtr nIt = smb->GetNodes();
1031 const SMDS_MeshNode* node;
1032 while ( nIt->more() ) {
1034 if(IsMedium(node, SMDSAbs_Edge))
1036 const SMDS_EdgePosition* pos =
1037 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1041 double u = ( pos->GetUParameter() - f ) / range;
1042 vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
1043 nVec.resize( vsize, nullNode );
1044 loadedNodes.insert( nVec[ 0 ] = node );
1046 if ( theParam2ColumnMap.size() != hsize ) {
1047 RETURN_BAD_RESULT( "Wrong node positions on theBaseEdge" );
1050 // load nodes from e1
1052 map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
1053 nIt = sm1->GetNodes();
1054 while ( nIt->more() ) {
1058 const SMDS_EdgePosition* pos =
1059 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1063 sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
1065 loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
1066 map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
1067 int row = rev1 ? vsize - 1 : 0;
1068 int dRow = rev1 ? -1 : +1;
1069 for ( ; u_n != sortedNodes.end(); u_n++ ) {
1071 loadedNodes.insert( nVecf[ row ] = u_n->second );
1074 // try to load the rest nodes
1076 // get all faces from theFace
1077 TIDSortedElemSet allFaces, foundFaces;
1078 eIt = smFace->GetElements();
1079 while ( eIt->more() ) {
1080 const SMDS_MeshElement* e = eIt->next();
1081 if ( e->GetType() == SMDSAbs_Face )
1082 allFaces.insert( e );
1084 // Starting from 2 neighbour nodes on theBaseEdge, look for a face
1085 // the nodes belong to, and between the nodes of the found face,
1086 // look for a not loaded node considering this node to be the next
1087 // in a column of the starting second node. Repeat, starting
1088 // from nodes next to the previous starting nodes in their columns,
1089 // and so on while a face can be found. Then go the the next pair
1090 // of nodes on theBaseEdge.
1091 TParam2ColumnMap::iterator par_nVec_1 = theParam2ColumnMap.begin();
1092 TParam2ColumnMap::iterator par_nVec_2 = par_nVec_1;
1095 for ( par_nVec_2++; par_nVec_2 != theParam2ColumnMap.end(); par_nVec_1++, par_nVec_2++ ) {
1098 const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
1099 const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
1100 const SMDS_MeshElement* face = 0;
1101 bool lastColOnClosedFace = ( nVecf[ row ] == n2 );
1103 // look for a face by 2 nodes
1104 face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
1106 int nbFaceNodes = face->NbNodes();
1107 if ( face->IsQuadratic() )
1109 if ( nbFaceNodes>4 ) {
1110 RETURN_BAD_RESULT(" Too many nodes in a face: " << nbFaceNodes );
1112 // look for a not loaded node of the <face>
1114 const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
1115 for ( int i = 0; i < nbFaceNodes && !found; ++i ) {
1116 node = face->GetNode( i );
1117 found = loadedNodes.insert( node ).second;
1118 if ( !found && node != n1 && node != n2 )
1121 if ( lastColOnClosedFace && row + 1 < vsize ) {
1122 node = nVecf[ row + 1 ];
1123 found = ( face->GetNodeIndex( node ) >= 0 );
1126 if ( ++row > vsize - 1 ) {
1127 RETURN_BAD_RESULT( "Too many nodes in column "<< col <<": "<< row+1);
1129 par_nVec_2->second[ row ] = node;
1130 foundFaces.insert( face );
1132 if ( nbFaceNodes==4 ) {
1133 n1 = par_nVec_1->second[ row ];
1136 else if ( nbFaceNodes==3 && n3 == par_nVec_1->second[ row + 1 ] ) {
1140 RETURN_BAD_RESULT( "Not quad mesh, column "<< col );
1144 while ( face && n1 && n2 );
1146 if ( row < vsize - 1 ) {
1147 MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
1148 MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
1149 MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
1150 if ( n1 ) { MESSAGE( "Current node 1: "<< n1); }
1151 else { MESSAGE( "Current node 1: NULL"); }
1152 if ( n2 ) { MESSAGE( "Current node 2: "<< n2); }
1153 else { MESSAGE( "Current node 2: NULL"); }
1154 MESSAGE( "first base node: "<< theParam2ColumnMap.begin()->second[0]);
1155 MESSAGE( "last base node: "<< theParam2ColumnMap.rbegin()->second[0]);
1158 } // loop on columns
1164 * Check mesh without geometry for: if all elements on this shape are quadratic,
1165 * quadratic elements will be created.
1166 * Used then generated 3D mesh without geometry.
1168 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
1170 int NbAllEdgsAndFaces=0;
1171 int NbQuadFacesAndEdgs=0;
1172 int NbFacesAndEdges=0;
1173 //All faces and edges
1174 NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
1176 //Quadratic faces and edges
1177 NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
1179 //Linear faces and edges
1180 NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
1182 if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
1184 return SMESH_MesherHelper::QUADRATIC;
1186 else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
1188 return SMESH_MesherHelper::LINEAR;
1191 //Mesh with both type of elements
1192 return SMESH_MesherHelper::COMP;