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 <Geom2d_Curve.hxx>
36 #include <Geom_Curve.hxx>
37 #include <Geom_Surface.hxx>
38 #include <TopExp_Explorer.hxx>
39 #include <TopTools_MapOfShape.hxx>
40 #include <gp_Pnt2d.hxx>
41 #include <ShapeAnalysis.hxx>
43 #include <utilities.h>
45 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
47 //=======================================================================
48 //function : CheckShape
50 //=======================================================================
52 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
54 SMESHDS_Mesh* meshDS = GetMeshDS();
55 // we can create quadratic elements only if all elements
56 // created on subshapes of given shape are quadratic
57 // also we have to fill myNLinkNodeMap
58 myCreateQuadratic = true;
59 mySeamShapeIds.clear();
60 TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
61 SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
63 TopExp_Explorer exp( aSh, subType );
64 for (; exp.More() && myCreateQuadratic; exp.Next()) {
65 if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
66 if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
68 const SMDS_MeshElement* e = it->next();
69 if ( e->GetType() != elemType || !e->IsQuadratic() ) {
70 myCreateQuadratic = false;
75 switch ( e->NbNodes() ) {
77 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
79 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
80 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
81 AddNLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
83 AddNLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
84 AddNLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
85 AddNLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
86 AddNLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
89 myCreateQuadratic = false;
98 if(!myCreateQuadratic) {
99 myNLinkNodeMap.clear();
103 return myCreateQuadratic;
106 //================================================================================
108 * \brief Set geomerty to make elements on
109 * \param aSh - geomertic shape
111 //================================================================================
113 void SMESH_MesherHelper::SetSubShape(const int aShID)
115 if ( aShID == myShapeID )
118 SetSubShape( GetMesh()->GetMeshDS()->IndexToShape( aShID ));
120 SetSubShape( TopoDS_Shape() );
123 //================================================================================
125 * \brief Set geomerty to make elements on
126 * \param aSh - geomertic shape
128 //================================================================================
130 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
132 if ( myShape.IsSame( aSh ))
136 mySeamShapeIds.clear();
138 if ( myShape.IsNull() ) {
142 SMESHDS_Mesh* meshDS = GetMeshDS();
143 myShapeID = meshDS->ShapeToIndex(aSh);
145 // treatment of periodic faces
146 for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
148 const TopoDS_Face& face = TopoDS::Face( eF.Current() );
149 BRepAdaptor_Surface surface( face );
150 if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
152 // look for a seam edge
153 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) {
154 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
155 if ( BRep_Tool::IsClosed( edge, face )) {
156 // initialize myPar1, myPar2 and myParIndex
157 if ( mySeamShapeIds.empty() ) {
159 BRep_Tool::UVPoints( edge, face, uv1, uv2 );
160 if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
162 myParIndex = 1; // U periodic
163 myPar1 = surface.FirstUParameter();
164 myPar2 = surface.LastUParameter();
167 myParIndex = 2; // V periodic
168 myPar1 = surface.FirstVParameter();
169 myPar2 = surface.LastVParameter();
172 // store shapes indices
173 mySeamShapeIds.insert( meshDS->ShapeToIndex( exp.Current() ));
174 for ( TopExp_Explorer v( exp.Current(), TopAbs_VERTEX ); v.More(); v.Next() )
175 mySeamShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
182 //================================================================================
184 * \brief Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
185 * \param F - the face
186 * \retval bool - return true if the face is periodic
188 //================================================================================
190 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
192 if ( F.IsNull() ) return !mySeamShapeIds.empty();
194 if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
195 return !mySeamShapeIds.empty();
197 Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F );
198 if ( !aSurface.IsNull() )
199 return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
204 //=======================================================================
205 //function : IsMedium
207 //=======================================================================
209 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode* node,
210 const SMDSAbs_ElementType typeToCheck)
212 return SMESH_MeshEditor::IsMedium( node, typeToCheck );
215 //=======================================================================
216 //function : AddNLinkNode
218 //=======================================================================
220 * Auxilary function for filling myNLinkNodeMap
222 void SMESH_MesherHelper::AddNLinkNode(const SMDS_MeshNode* n1,
223 const SMDS_MeshNode* n2,
224 const SMDS_MeshNode* n12)
226 NLink link( n1, n2 );
227 if ( n1 > n2 ) link = NLink( n2, n1 );
228 // add new record to map
229 myNLinkNodeMap.insert( make_pair(link,n12));
232 //=======================================================================
234 * \brief Select UV on either of 2 pcurves of a seam edge, closest to the given UV
235 * \param uv1 - UV on the seam
236 * \param uv2 - UV within a face
237 * \retval gp_Pnt2d - selected UV
239 //=======================================================================
241 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
243 double p1 = uv1.Coord( myParIndex );
244 double p2 = uv2.Coord( myParIndex );
245 double p3 = ( Abs( p1 - myPar1 ) < Abs( p1 - myPar2 )) ? myPar2 : myPar1;
246 if ( Abs( p2 - p1 ) > Abs( p2 - p3 ))
248 gp_Pnt2d result = uv1;
249 result.SetCoord( myParIndex, p1 );
253 //=======================================================================
255 * \brief Return node UV on face
256 * \param F - the face
257 * \param n - the node
258 * \param n2 - a node of element being created located inside a face
259 * \retval gp_XY - resulting UV
261 * Auxilary function called form GetMediumNode()
263 //=======================================================================
265 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F,
266 const SMDS_MeshNode* n,
267 const SMDS_MeshNode* n2) const
269 gp_Pnt2d uv( 1e100, 1e100 );
270 const SMDS_PositionPtr Pos = n->GetPosition();
271 if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
273 // node has position on face
274 const SMDS_FacePosition* fpos =
275 static_cast<const SMDS_FacePosition*>(n->GetPosition().get());
276 uv = gp_Pnt2d(fpos->GetUParameter(),fpos->GetVParameter());
278 else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
280 // node has position on edge => it is needed to find
281 // corresponding edge from face, get pcurve for this
282 // edge and recieve value from this pcurve
283 const SMDS_EdgePosition* epos =
284 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
285 SMESHDS_Mesh* meshDS = GetMesh()->GetMeshDS();
286 int edgeID = Pos->GetShapeId();
287 TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
290 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
291 uv = C2d->Value( epos->GetUParameter() );
292 // for a node on a seam edge select one of UVs on 2 pcurves
293 if ( n2 && mySeamShapeIds.find( edgeID ) != mySeamShapeIds.end() )
294 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
296 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
298 int vertexID = n->GetPosition()->GetShapeId();
299 const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
300 uv = BRep_Tool::Parameters( V, F );
301 if ( n2 && mySeamShapeIds.find( vertexID ) != mySeamShapeIds.end() )
302 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
307 //=======================================================================
309 * \brief Return node U on edge
310 * \param E - the Edge
311 * \param n - the node
312 * \retval double - resulting U
314 * Auxilary function called form GetMediumNode()
316 //=======================================================================
318 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E,
319 const SMDS_MeshNode* n)
322 const SMDS_PositionPtr Pos = n->GetPosition();
323 if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE) {
324 const SMDS_EdgePosition* epos =
325 static_cast<const SMDS_EdgePosition*>(n->GetPosition().get());
326 param = epos->GetUParameter();
328 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX) {
329 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
330 int vertexID = n->GetPosition()->GetShapeId();
331 const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
332 param = BRep_Tool::Parameter( V, E );
337 //=======================================================================
338 //function : GetMediumNode
340 //=======================================================================
342 * Special function for search or creation medium node
344 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
345 const SMDS_MeshNode* n2,
348 TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
350 NLink link(( n1 < n2 ? n1 : n2 ), ( n1 < n2 ? n2 : n1 ));
351 ItNLinkNode itLN = myNLinkNodeMap.find( link );
352 if ( itLN != myNLinkNodeMap.end() ) {
353 return (*itLN).second;
356 // create medium node
358 SMESHDS_Mesh* meshDS = GetMeshDS();
359 int faceID = -1, edgeID = -1;
360 const SMDS_PositionPtr Pos1 = n1->GetPosition();
361 const SMDS_PositionPtr Pos2 = n2->GetPosition();
363 if( myShape.IsNull() )
365 if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
366 faceID = Pos1->GetShapeId();
368 else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
369 faceID = Pos2->GetShapeId();
372 if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
373 edgeID = Pos1->GetShapeId();
375 if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
376 edgeID = Pos2->GetShapeId();
381 // we try to create medium node using UV parameters of
382 // nodes, else - medium between corresponding 3d points
383 if(faceID>-1 || shapeType == TopAbs_FACE) {
384 // obtaining a face and 2d points for nodes
386 if( myShape.IsNull() )
387 F = TopoDS::Face(meshDS->IndexToShape(faceID));
389 F = TopoDS::Face(myShape);
393 gp_XY p1 = GetNodeUV(F,n1,n2);
394 gp_XY p2 = GetNodeUV(F,n2,n1);
396 //checking if surface is periodic
397 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
398 Standard_Real UF,UL,VF,VL;
399 S->Bounds(UF,UL,VF,VL);
402 Standard_Boolean isUPeriodic = S->IsUPeriodic();
404 Standard_Real UPeriod = S->UPeriod();
405 Standard_Real p2x = p2.X()+ShapeAnalysis::AdjustByPeriod(p2.X(),p1.X(),UPeriod);
406 Standard_Real pmid = (p1.X()+p2x)/2.;
407 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,UF,UL);
410 u= (p1.X()+p2.X())/2.;
412 Standard_Boolean isVPeriodic = S->IsVPeriodic();
414 Standard_Real VPeriod = S->VPeriod();
415 Standard_Real p2y = p2.Y()+ShapeAnalysis::AdjustByPeriod(p2.Y(),p1.Y(),VPeriod);
416 Standard_Real pmid = (p1.Y()+p2y)/2.;
417 v = pmid+ShapeAnalysis::AdjustToPeriod(pmid,VF,VL);
420 v = (p1.Y()+p2.Y())/2.;
422 gp_Pnt P = S->Value(u, v);
423 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
424 meshDS->SetNodeOnFace(n12, faceID, u, v);
425 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
428 if (edgeID>-1 || shapeType == TopAbs_EDGE) {
431 if( myShape.IsNull() )
432 E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
434 E = TopoDS::Edge(myShape);
438 double p1 = GetNodeU(E,n1);
439 double p2 = GetNodeU(E,n2);
442 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
445 Standard_Boolean isPeriodic = C->IsPeriodic();
448 Standard_Real Period = C->Period();
449 Standard_Real p = p2+ShapeAnalysis::AdjustByPeriod(p2,p1,Period);
450 Standard_Real pmid = (p1+p)/2.;
451 u = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
456 gp_Pnt P = C->Value( u );
457 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
458 meshDS->SetNodeOnEdge(n12, edgeID, u);
459 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
465 double x = ( n1->X() + n2->X() )/2.;
466 double y = ( n1->Y() + n2->Y() )/2.;
467 double z = ( n1->Z() + n2->Z() )/2.;
468 n12 = meshDS->AddNode(x,y,z);
470 meshDS->SetNodeOnEdge(n12, edgeID);
472 meshDS->SetNodeOnFace(n12, faceID);
474 meshDS->SetNodeInVolume(n12, myShapeID);
475 myNLinkNodeMap.insert(NLinkNodeMap::value_type(link,n12));
480 //=======================================================================
481 //function : AddQuadraticEdge
483 //=======================================================================
485 * Special function for creation quadratic edge
487 SMDS_QuadraticEdge* SMESH_MesherHelper::AddQuadraticEdge(const SMDS_MeshNode* n1,
488 const SMDS_MeshNode* n2,
492 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
494 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
496 myCreateQuadratic = true;
499 return (SMDS_QuadraticEdge*)(meshDS->AddEdgeWithID(n1, n2, n12, id));
501 return (SMDS_QuadraticEdge*)(meshDS->AddEdge(n1, n2, n12));
504 //=======================================================================
507 //=======================================================================
509 * Special function for creation quadratic triangle
511 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
512 const SMDS_MeshNode* n2,
513 const SMDS_MeshNode* n3,
517 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
518 if(!myCreateQuadratic) {
520 return meshDS->AddFaceWithID(n1, n2, n3, id);
522 return meshDS->AddFace(n1, n2, n3);
525 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
526 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
527 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
530 return meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
532 return meshDS->AddFace(n1, n2, n3, n12, n23, n31);
536 //=======================================================================
539 //=======================================================================
541 * Special function for creation quadratic quadrangle
543 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
544 const SMDS_MeshNode* n2,
545 const SMDS_MeshNode* n3,
546 const SMDS_MeshNode* n4,
550 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
551 if(!myCreateQuadratic) {
553 return meshDS->AddFaceWithID(n1, n2, n3, n4, id);
555 return meshDS->AddFace(n1, n2, n3, n4);
558 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
559 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
560 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
561 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
564 return meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
566 return meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
570 //=======================================================================
571 //function : AddVolume
573 //=======================================================================
575 * Special function for creation quadratic volume
577 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
578 const SMDS_MeshNode* n2,
579 const SMDS_MeshNode* n3,
580 const SMDS_MeshNode* n4,
581 const SMDS_MeshNode* n5,
582 const SMDS_MeshNode* n6,
586 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
587 if(!myCreateQuadratic) {
589 return meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
591 return meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
594 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
595 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
596 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
598 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
599 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
600 const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
602 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
603 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
604 const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
607 return meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
608 n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
610 return meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
611 n12, n23, n31, n45, n56, n64, n14, n25, n36);
614 //=======================================================================
616 * Special function for creation quadratic volume
618 //=======================================================================
620 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
621 const SMDS_MeshNode* n2,
622 const SMDS_MeshNode* n3,
623 const SMDS_MeshNode* n4,
627 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
628 if(!myCreateQuadratic) {
630 return meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
632 return meshDS->AddVolume(n1, n2, n3, n4);
635 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
636 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
637 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
639 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
640 const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
641 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
644 return meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
646 return meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
649 //=======================================================================
651 * Special function for creation quadratic pyramid
653 //=======================================================================
655 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
656 const SMDS_MeshNode* n2,
657 const SMDS_MeshNode* n3,
658 const SMDS_MeshNode* n4,
659 const SMDS_MeshNode* n5,
663 if(!myCreateQuadratic) {
665 return GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
667 return GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
670 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
671 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
672 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
673 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
675 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
676 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
677 const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
678 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
681 return GetMeshDS()->AddVolumeWithID ( n1, n2, n3, n4, n5,
686 return GetMeshDS()->AddVolume( n1, n2, n3, n4, n5,
691 //=======================================================================
693 * Special function for creation of quadratic hexahedron
695 //=======================================================================
697 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
698 const SMDS_MeshNode* n2,
699 const SMDS_MeshNode* n3,
700 const SMDS_MeshNode* n4,
701 const SMDS_MeshNode* n5,
702 const SMDS_MeshNode* n6,
703 const SMDS_MeshNode* n7,
704 const SMDS_MeshNode* n8,
708 SMESHDS_Mesh * meshDS = GetMesh()->GetMeshDS();
709 if(!myCreateQuadratic) {
711 return meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
713 return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
716 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
717 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
718 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
719 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
721 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
722 const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
723 const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
724 const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
726 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
727 const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
728 const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
729 const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
732 return meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
733 n12, n23, n34, n41, n56, n67,
734 n78, n85, n15, n26, n37, n48, id);
736 return meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
737 n12, n23, n34, n41, n56, n67,
738 n78, n85, n15, n26, n37, n48);
741 //=======================================================================
743 * \brief Load nodes bound to face into a map of node columns
744 * \param theParam2ColumnMap - map of node columns to fill
745 * \param theFace - the face on which nodes are searched for
746 * \param theBaseEdge - the edge nodes of which are columns' bases
747 * \param theMesh - the mesh containing nodes
748 * \retval bool - false if something is wrong
750 * The key of the map is a normalized parameter of each
751 * base node on theBaseEdge.
752 * This method works in supposition that nodes on the face
753 * forms a rectangular grid and elements can be quardrangles or triangles
755 //=======================================================================
757 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
758 const TopoDS_Face& theFace,
759 const TopoDS_Edge& theBaseEdge,
760 SMESHDS_Mesh* theMesh)
762 // get vertices of theBaseEdge
763 TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
764 TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
765 TopExp::Vertices( eFrw, vfb, vlb );
767 // find the other edges of theFace and orientation of e1
768 TopoDS_Edge e1, e2, eTop;
769 bool rev1, CumOri = false;
770 TopExp_Explorer exp( theFace, TopAbs_EDGE );
772 for ( ; exp.More(); exp.Next() ) {
773 if ( ++nbEdges > 4 ) {
774 return false; // more than 4 edges in theFace
776 TopoDS_Edge e = TopoDS::Edge( exp.Current() );
777 if ( theBaseEdge.IsSame( e ))
779 TopoDS_Vertex vCommon;
780 if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
782 else if ( vCommon.IsSame( vfb )) {
784 vft = TopExp::LastVertex( e1, CumOri );
785 rev1 = vfb.IsSame( vft );
787 vft = TopExp::FirstVertex( e1, CumOri );
793 return false; // less than 4 edges in theFace
795 if ( e2.IsNull() && vfb.IsSame( vlb ))
798 // submeshes corresponding to shapes
799 SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
800 SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
801 SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
802 SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
803 SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
804 SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
805 SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
806 SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
807 if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
808 RETURN_BAD_RESULT( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
809 sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
811 if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
812 RETURN_BAD_RESULT(" Diff nb of nodes on opposite edges" );
814 if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
815 RETURN_BAD_RESULT("Empty submesh of vertex");
817 // define whether mesh is quadratic
818 bool isQuadraticMesh = false;
819 SMDS_ElemIteratorPtr eIt = smFace->GetElements();
820 if ( !eIt->more() ) {
821 RETURN_BAD_RESULT("No elements on the face");
823 const SMDS_MeshElement* e = eIt->next();
824 isQuadraticMesh = e->IsQuadratic();
826 if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
827 // check quadratic case
828 if ( isQuadraticMesh ) {
829 // what if there are quadrangles and triangles mixed?
830 // int n1 = sm1->NbNodes()/2;
831 // int n2 = smb->NbNodes()/2;
832 // int n3 = sm1->NbNodes() - n1;
833 // int n4 = smb->NbNodes() - n2;
834 // int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
835 // if( nf != smFace->NbNodes() ) {
836 // MESSAGE( "Wrong nb face nodes: " <<
837 // sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
842 RETURN_BAD_RESULT( "Wrong nb face nodes: " <<
843 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
847 int vsize = sm1->NbNodes() + 2;
848 int hsize = smb->NbNodes() + 2;
849 if(isQuadraticMesh) {
850 vsize = vsize - sm1->NbNodes()/2 -1;
851 hsize = hsize - smb->NbNodes()/2 -1;
854 // load nodes from theBaseEdge
856 set<const SMDS_MeshNode*> loadedNodes;
857 const SMDS_MeshNode* nullNode = 0;
859 vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
860 nVecf.resize( vsize, nullNode );
861 loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
863 vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
864 nVecl.resize( vsize, nullNode );
865 loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
868 BRep_Tool::Range( eFrw, f, l );
869 double range = l - f;
870 SMDS_NodeIteratorPtr nIt = smb->GetNodes();
871 const SMDS_MeshNode* node;
872 while ( nIt->more() ) {
876 const SMDS_EdgePosition* pos =
877 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
881 double u = ( pos->GetUParameter() - f ) / range;
882 vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
883 nVec.resize( vsize, nullNode );
884 loadedNodes.insert( nVec[ 0 ] = node );
886 if ( theParam2ColumnMap.size() != hsize ) {
887 RETURN_BAD_RESULT( "Wrong node positions on theBaseEdge" );
890 // load nodes from e1
892 map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
893 nIt = sm1->GetNodes();
894 while ( nIt->more() ) {
898 const SMDS_EdgePosition* pos =
899 dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
903 sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
905 loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
906 map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
907 int row = rev1 ? vsize - 1 : 0;
908 int dRow = rev1 ? -1 : +1;
909 for ( ; u_n != sortedNodes.end(); u_n++ ) {
911 loadedNodes.insert( nVecf[ row ] = u_n->second );
914 // try to load the rest nodes
916 // get all faces from theFace
917 TIDSortedElemSet allFaces, foundFaces;
918 eIt = smFace->GetElements();
919 while ( eIt->more() ) {
920 const SMDS_MeshElement* e = eIt->next();
921 if ( e->GetType() == SMDSAbs_Face )
922 allFaces.insert( e );
924 // Starting from 2 neighbour nodes on theBaseEdge, look for a face
925 // the nodes belong to, and between the nodes of the found face,
926 // look for a not loaded node considering this node to be the next
927 // in a column of the starting second node. Repeat, starting
928 // from nodes next to the previous starting nodes in their columns,
929 // and so on while a face can be found. Then go the the next pair
930 // of nodes on theBaseEdge.
931 TParam2ColumnMap::iterator par_nVec_1 = theParam2ColumnMap.begin();
932 TParam2ColumnMap::iterator par_nVec_2 = par_nVec_1;
935 for ( par_nVec_2++; par_nVec_2 != theParam2ColumnMap.end(); par_nVec_1++, par_nVec_2++ ) {
938 const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
939 const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
940 const SMDS_MeshElement* face = 0;
941 bool lastColOnClosedFace = ( nVecf[ row ] == n2 );
943 // look for a face by 2 nodes
944 face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
946 int nbFaceNodes = face->NbNodes();
947 if ( face->IsQuadratic() )
949 if ( nbFaceNodes>4 ) {
950 RETURN_BAD_RESULT(" Too many nodes in a face: " << nbFaceNodes );
952 // look for a not loaded node of the <face>
954 const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
955 for ( int i = 0; i < nbFaceNodes && !found; ++i ) {
956 node = face->GetNode( i );
957 found = loadedNodes.insert( node ).second;
958 if ( !found && node != n1 && node != n2 )
961 if ( lastColOnClosedFace && row + 1 < vsize ) {
962 node = nVecf[ row + 1 ];
963 found = ( face->GetNodeIndex( node ) >= 0 );
966 if ( ++row > vsize - 1 ) {
967 RETURN_BAD_RESULT( "Too many nodes in column "<< col <<": "<< row+1);
969 par_nVec_2->second[ row ] = node;
970 foundFaces.insert( face );
972 if ( nbFaceNodes==4 ) {
973 n1 = par_nVec_1->second[ row ];
976 else if ( nbFaceNodes==3 && n3 == par_nVec_1->second[ row + 1 ] ) {
980 RETURN_BAD_RESULT( "Not quad mesh, column "<< col );
984 while ( face && n1 && n2 );
986 if ( row < vsize - 1 ) {
987 MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
988 MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
989 MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
990 if ( n1 ) { MESSAGE( "Current node 1: "<< n1); }
991 else { MESSAGE( "Current node 1: NULL"); }
992 if ( n2 ) { MESSAGE( "Current node 2: "<< n2); }
993 else { MESSAGE( "Current node 2: NULL"); }
994 MESSAGE( "first base node: "<< theParam2ColumnMap.begin()->second[0]);
995 MESSAGE( "last base node: "<< theParam2ColumnMap.rbegin()->second[0]);