00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "SMESH_MesherHelper.hxx"
00028
00029 #include "SMDS_FacePosition.hxx"
00030 #include "SMDS_EdgePosition.hxx"
00031 #include "SMDS_VolumeTool.hxx"
00032 #include "SMESH_subMesh.hxx"
00033 #include "SMESH_ProxyMesh.hxx"
00034
00035 #include <BRepAdaptor_Surface.hxx>
00036 #include <BRepTools.hxx>
00037 #include <BRepTools_WireExplorer.hxx>
00038 #include <BRep_Tool.hxx>
00039 #include <Geom2d_Curve.hxx>
00040 #include <GeomAPI_ProjectPointOnCurve.hxx>
00041 #include <GeomAPI_ProjectPointOnSurf.hxx>
00042 #include <Geom_Curve.hxx>
00043 #include <Geom_RectangularTrimmedSurface.hxx>
00044 #include <Geom_Surface.hxx>
00045 #include <ShapeAnalysis.hxx>
00046 #include <TopExp.hxx>
00047 #include <TopExp_Explorer.hxx>
00048 #include <TopTools_ListIteratorOfListOfShape.hxx>
00049 #include <TopTools_MapIteratorOfMapOfShape.hxx>
00050 #include <TopTools_MapOfShape.hxx>
00051 #include <TopoDS.hxx>
00052 #include <gp_Ax3.hxx>
00053 #include <gp_Pnt2d.hxx>
00054 #include <gp_Trsf.hxx>
00055
00056 #include <Standard_Failure.hxx>
00057 #include <Standard_ErrorHandler.hxx>
00058
00059 #include <utilities.h>
00060
00061 #include <limits>
00062
00063 using namespace std;
00064
00065 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
00066
00067 namespace {
00068
00069 gp_XYZ XYZ(const SMDS_MeshNode* n) { return gp_XYZ(n->X(), n->Y(), n->Z()); }
00070
00071 enum { U_periodic = 1, V_periodic = 2 };
00072 }
00073
00074
00078
00079
00080 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
00081 : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false)
00082 {
00083 myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0;
00084 mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
00085 }
00086
00087
00088
00089
00090
00091
00092 SMESH_MesherHelper::~SMESH_MesherHelper()
00093 {
00094 {
00095 TID2ProjectorOnSurf::iterator i_proj = myFace2Projector.begin();
00096 for ( ; i_proj != myFace2Projector.end(); ++i_proj )
00097 delete i_proj->second;
00098 }
00099 {
00100 TID2ProjectorOnCurve::iterator i_proj = myEdge2Projector.begin();
00101 for ( ; i_proj != myEdge2Projector.end(); ++i_proj )
00102 delete i_proj->second;
00103 }
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
00114 {
00115 SMESHDS_Mesh* meshDS = GetMeshDS();
00116
00117
00118
00119 myCreateQuadratic = true;
00120 mySeamShapeIds.clear();
00121 myDegenShapeIds.clear();
00122 TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
00123 SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
00124
00125 int nbOldLinks = myTLinkNodeMap.size();
00126
00127 if ( !myMesh->HasShapeToMesh() )
00128 {
00129 if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC )))
00130 {
00131 SMDS_FaceIteratorPtr fIt = meshDS->facesIterator();
00132 while ( fIt->more() )
00133 AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
00134 }
00135 }
00136 else
00137 {
00138 TopExp_Explorer exp( aSh, subType );
00139 TopTools_MapOfShape checkedSubShapes;
00140 for (; exp.More() && myCreateQuadratic; exp.Next()) {
00141 if ( !checkedSubShapes.Add( exp.Current() ))
00142 continue;
00143 if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
00144 if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
00145 while(it->more()) {
00146 const SMDS_MeshElement* e = it->next();
00147 if ( e->GetType() != elemType || !e->IsQuadratic() ) {
00148 myCreateQuadratic = false;
00149 break;
00150 }
00151 else {
00152
00153 switch ( e->NbNodes() ) {
00154 case 3:
00155 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
00156 case 6:
00157 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
00158 AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
00159 AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
00160 case 8:
00161 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
00162 AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
00163 AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
00164 AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
00165 break;
00166 default:
00167 myCreateQuadratic = false;
00168 break;
00169 }
00170 }
00171 }
00172 }
00173 }
00174 }
00175 }
00176
00177 if ( nbOldLinks == myTLinkNodeMap.size() )
00178 myCreateQuadratic = false;
00179
00180 if(!myCreateQuadratic) {
00181 myTLinkNodeMap.clear();
00182 }
00183 SetSubShape( aSh );
00184
00185 return myCreateQuadratic;
00186 }
00187
00188
00189
00190
00191
00192
00193 void SMESH_MesherHelper::SetSubShape(const int aShID)
00194 {
00195 if ( aShID == myShapeID )
00196 return;
00197 if ( aShID > 0 )
00198 SetSubShape( GetMeshDS()->IndexToShape( aShID ));
00199 else
00200 SetSubShape( TopoDS_Shape() );
00201 }
00202
00203
00204
00205
00206
00207
00208 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
00209 {
00210 if ( myShape.IsSame( aSh ))
00211 return;
00212
00213 myShape = aSh;
00214 mySeamShapeIds.clear();
00215 myDegenShapeIds.clear();
00216
00217 if ( myShape.IsNull() ) {
00218 myShapeID = 0;
00219 return;
00220 }
00221 SMESHDS_Mesh* meshDS = GetMeshDS();
00222 myShapeID = meshDS->ShapeToIndex(aSh);
00223 myParIndex = 0;
00224
00225
00226 for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
00227 {
00228 const TopoDS_Face& face = TopoDS::Face( eF.Current() );
00229 TopLoc_Location loc;
00230 Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
00231
00232 if ( surface->IsUPeriodic() || surface->IsVPeriodic() )
00233 {
00234
00235
00236 GeomAdaptor_Surface surf( surface );
00237
00238 for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
00239 {
00240
00241 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
00242 if ( BRep_Tool::IsClosed( edge, face )) {
00243
00244 gp_Pnt2d uv1, uv2;
00245 BRep_Tool::UVPoints( edge, face, uv1, uv2 );
00246 if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
00247 {
00248 myParIndex |= U_periodic;
00249 myPar1[0] = surf.FirstUParameter();
00250 myPar2[0] = surf.LastUParameter();
00251 }
00252 else {
00253 myParIndex |= V_periodic;
00254 myPar1[1] = surf.FirstVParameter();
00255 myPar2[1] = surf.LastVParameter();
00256 }
00257
00258 int edgeID = meshDS->ShapeToIndex( edge );
00259 mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
00260 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
00261 int vertexID = meshDS->ShapeToIndex( v.Current() );
00262 mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
00263 }
00264 }
00265
00266
00267 if ( BRep_Tool::Degenerated( edge )) {
00268 myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
00269 for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
00270 myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
00271 }
00272 }
00273 }
00274 }
00275 }
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
00286 {
00287 if ( F.IsNull() ) return !mySeamShapeIds.empty();
00288
00289 if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
00290 return !mySeamShapeIds.empty();
00291
00292 TopLoc_Location loc;
00293 Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F,loc );
00294 if ( !aSurface.IsNull() )
00295 return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
00296
00297 return false;
00298 }
00299
00300
00301
00302
00303
00304
00305 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode* node,
00306 const SMDSAbs_ElementType typeToCheck)
00307 {
00308 return SMESH_MeshEditor::IsMedium( node, typeToCheck );
00309 }
00310
00311
00312
00313
00314
00315
00316 TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node,
00317 const SMESHDS_Mesh* meshDS)
00318 {
00319 int shapeID = node->getshapeId();
00320 if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() )
00321 return meshDS->IndexToShape( shapeID );
00322 else
00323 return TopoDS_Shape();
00324 }
00325
00326
00327
00328
00329
00330
00331
00332 void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
00333 const SMDS_MeshNode* n2,
00334 const SMDS_MeshNode* n12)
00335 {
00336
00337 SMESH_TLink link( n1, n2 );
00338 myTLinkNodeMap.insert( make_pair(link,n12));
00339 }
00340
00341
00345
00346
00347 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
00348 {
00349 if ( edge->IsQuadratic() )
00350 AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
00351 }
00352
00353
00357
00358
00359 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
00360 {
00361 if ( !f->IsPoly() )
00362 switch ( f->NbNodes() ) {
00363 case 6:
00364 AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(3));
00365 AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(4));
00366 AddTLinkNode(f->GetNode(2),f->GetNode(0),f->GetNode(5)); break;
00367 case 8:
00368 AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(4));
00369 AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(5));
00370 AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
00371 AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7));
00372 default:;
00373 }
00374 }
00375
00376
00380
00381
00382 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
00383 {
00384 if ( volume->IsQuadratic() )
00385 {
00386 SMDS_VolumeTool vTool( volume );
00387 const SMDS_MeshNode** nodes = vTool.GetNodes();
00388 set<int> addedLinks;
00389 for ( int iF = 1; iF < vTool.NbFaces(); ++iF )
00390 {
00391 const int nbN = vTool.NbFaceNodes( iF );
00392 const int* iNodes = vTool.GetFaceNodesIndices( iF );
00393 for ( int i = 0; i < nbN; )
00394 {
00395 int iN1 = iNodes[i++];
00396 int iN12 = iNodes[i++];
00397 int iN2 = iNodes[i++];
00398 if ( iN1 > iN2 ) std::swap( iN1, iN2 );
00399 int linkID = iN1 * vTool.NbNodes() + iN2;
00400 pair< set<int>::iterator, bool > it_isNew = addedLinks.insert( linkID );
00401 if ( it_isNew.second )
00402 AddTLinkNode( nodes[iN1], nodes[iN2], nodes[iN12] );
00403 else
00404 addedLinks.erase( it_isNew.first );
00405 }
00406 }
00407 }
00408 }
00409
00410
00415
00416
00417 bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const
00418 {
00419 map< int,bool >::const_iterator id_ok = myNodePosShapesValidity.find( shapeID );
00420 return ( id_ok == myNodePosShapesValidity.end() || !id_ok->second );
00421 }
00422
00423
00428
00429
00430 void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
00431 {
00432 ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok));
00433 }
00434
00435
00436
00437
00438
00439
00440 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
00441 {
00442 gp_Pnt2d result = uv1;
00443 for ( int i = U_periodic; i <= V_periodic ; ++i )
00444 {
00445 if ( myParIndex & i )
00446 {
00447 double p1 = uv1.Coord( i );
00448 double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]);
00449 if ( myParIndex == i ||
00450 dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ||
00451 dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. )
00452 {
00453 double p2 = uv2.Coord( i );
00454 double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1];
00455 if ( Abs( p2 - p1 ) > Abs( p2 - p1Alt ))
00456 result.SetCoord( i, p1Alt );
00457 }
00458 }
00459 }
00460 return result;
00461 }
00462
00463
00464
00465
00466
00467
00468 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F,
00469 const SMDS_MeshNode* n,
00470 const SMDS_MeshNode* n2,
00471 bool* check) const
00472 {
00473 gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() );
00474 const SMDS_PositionPtr Pos = n->GetPosition();
00475 bool uvOK = false;
00476 if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
00477 {
00478
00479 const SMDS_FacePosition* fpos =
00480 static_cast<const SMDS_FacePosition*>(n->GetPosition());
00481 uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
00482 if ( check )
00483 uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
00484 }
00485 else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
00486 {
00487
00488
00489
00490 const SMDS_EdgePosition* epos =
00491 static_cast<const SMDS_EdgePosition*>(n->GetPosition());
00492 int edgeID = n->getshapeId();
00493 TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
00494 double f, l, u = epos->GetUParameter();
00495 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
00496 bool validU = ( f < u && u < l );
00497 if ( validU )
00498 uv = C2d->Value( u );
00499 else
00500 uv.SetCoord( Precision::Infinite(),0.);
00501 if ( check || !validU )
00502 uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ), !validU );
00503
00504
00505 if ( n2 && IsSeamShape( edgeID ) )
00506 {
00507 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0, check ));
00508 }
00509 else
00510 {
00511 TopLoc_Location loc;
00512 Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
00513 Standard_Boolean isUPeriodic = S->IsUPeriodic();
00514 Standard_Boolean isVPeriodic = S->IsVPeriodic();
00515 if ( isUPeriodic || isVPeriodic ) {
00516 Standard_Real UF,UL,VF,VL;
00517 S->Bounds(UF,UL,VF,VL);
00518 if(isUPeriodic)
00519 uv.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod(uv.X(),UF,UL));
00520 if(isVPeriodic)
00521 uv.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod(uv.Y(),VF,VL));
00522 }
00523 }
00524 }
00525 else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
00526 {
00527 if ( int vertexID = n->getshapeId() ) {
00528 const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
00529 try {
00530 uv = BRep_Tool::Parameters( V, F );
00531 uvOK = true;
00532 }
00533 catch (Standard_Failure& exc) {
00534 }
00535 if ( !uvOK ) {
00536 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() )
00537 uvOK = ( V == vert.Current() );
00538 if ( !uvOK ) {
00539 #ifdef _DEBUG_
00540 MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
00541 << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
00542 #endif
00543
00544 double dist = 1e100;
00545 gp_Pnt pn = XYZ( n );
00546 for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) {
00547 TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
00548 gp_Pnt p = BRep_Tool::Pnt( curV );
00549 double curDist = p.SquareDistance( pn );
00550 if ( curDist < dist ) {
00551 dist = curDist;
00552 uv = BRep_Tool::Parameters( curV, F );
00553 uvOK = ( dist < DBL_MIN );
00554 }
00555 }
00556 }
00557 else {
00558 uvOK = false;
00559 TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
00560 for ( ; it.More(); it.Next() ) {
00561 if ( it.Value().ShapeType() == TopAbs_EDGE ) {
00562 const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
00563 double f,l;
00564 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
00565 if ( !C2d.IsNull() ) {
00566 double u = ( V == TopExp::FirstVertex( edge ) ) ? f : l;
00567 uv = C2d->Value( u );
00568 uvOK = true;
00569 break;
00570 }
00571 }
00572 }
00573 }
00574 }
00575 if ( n2 && IsSeamShape( vertexID ) )
00576 uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
00577 }
00578 }
00579 else
00580 {
00581 uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
00582 }
00583
00584 if ( check )
00585 *check = uvOK;
00586
00587 return uv.XY();
00588 }
00589
00590
00591
00592
00593
00594
00595 bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F,
00596 const SMDS_MeshNode* n,
00597 gp_XY& uv,
00598 const double tol,
00599 const bool force,
00600 double distXYZ[4]) const
00601 {
00602 int shapeID = n->getshapeId();
00603 bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
00604 if ( force || toCheckPosOnShape( shapeID ) || infinit )
00605 {
00606
00607 TopLoc_Location loc;
00608 Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
00609 gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
00610 double dist = 0;
00611 if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
00612 if ( infinit ||
00613 (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol )
00614 {
00615 setPosOnShapeValidity( shapeID, false );
00616 if ( !infinit && distXYZ ) {
00617 surfPnt.Transform( loc );
00618 distXYZ[0] = dist;
00619 distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
00620 }
00621
00622 GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
00623 projector.Perform( nodePnt );
00624 if ( !projector.IsDone() || projector.NbPoints() < 1 )
00625 {
00626 MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
00627 return false;
00628 }
00629 Quantity_Parameter U,V;
00630 projector.LowerDistanceParameters(U,V);
00631 uv.SetCoord( U,V );
00632 surfPnt = surface->Value( U, V );
00633 dist = nodePnt.Distance( surfPnt );
00634 if ( distXYZ ) {
00635 surfPnt.Transform( loc );
00636 distXYZ[0] = dist;
00637 distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
00638 }
00639 if ( dist > tol )
00640 {
00641 MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
00642 return false;
00643 }
00644
00645 if ( myShape.IsSame(F) && shapeID == myShapeID )
00646 const_cast<SMDS_MeshNode*>(n)->SetPosition
00647 ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
00648 }
00649 else if ( uv.Modulus() > numeric_limits<double>::min() )
00650 {
00651 setPosOnShapeValidity( shapeID, true );
00652 }
00653 }
00654 return true;
00655 }
00656
00657
00658
00659
00660
00661
00662 GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
00663 TopLoc_Location& loc,
00664 double tol ) const
00665 {
00666 Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
00667 int faceID = GetMeshDS()->ShapeToIndex( F );
00668 TID2ProjectorOnSurf& i2proj = const_cast< TID2ProjectorOnSurf&>( myFace2Projector );
00669 TID2ProjectorOnSurf::iterator i_proj = i2proj.find( faceID );
00670 if ( i_proj == i2proj.end() )
00671 {
00672 if ( tol == 0 ) tol = BRep_Tool::Tolerance( F );
00673 double U1, U2, V1, V2;
00674 surface->Bounds(U1, U2, V1, V2);
00675 GeomAPI_ProjectPointOnSurf* proj = new GeomAPI_ProjectPointOnSurf();
00676 proj->Init( surface, U1, U2, V1, V2, tol );
00677 i_proj = i2proj.insert( make_pair( faceID, proj )).first;
00678 }
00679 return *( i_proj->second );
00680 }
00681
00682 namespace
00683 {
00684 gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; }
00685 gp_XY_FunPtr(Added);
00686 gp_XY_FunPtr(Subtracted);
00687 }
00688
00689
00690
00691
00692
00693
00694
00695
00696 gp_XY SMESH_MesherHelper::applyIn2D(const Handle(Geom_Surface)& surface,
00697 const gp_XY& uv1,
00698 const gp_XY& uv2,
00699 xyFunPtr fun,
00700 const bool resultInPeriod)
00701 {
00702 Standard_Boolean isUPeriodic = surface.IsNull() ? false : surface->IsUPeriodic();
00703 Standard_Boolean isVPeriodic = surface.IsNull() ? false : surface->IsVPeriodic();
00704 if ( !isUPeriodic && !isVPeriodic )
00705 return fun(uv1,uv2);
00706
00707
00708 double u2 =
00709 uv2.X()+(isUPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.X(),uv1.X(),surface->UPeriod()) :0);
00710 double v2 =
00711 uv2.Y()+(isVPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.Y(),uv1.Y(),surface->VPeriod()) :0);
00712
00713
00714 gp_XY res = fun( uv1, gp_XY(u2,v2) );
00715
00716
00717 if ( resultInPeriod )
00718 {
00719 Standard_Real UF,UL,VF,VL;
00720 surface->Bounds(UF,UL,VF,VL);
00721 if ( isUPeriodic )
00722 res.SetX( res.X() + ShapeAnalysis::AdjustToPeriod(res.X(),UF,UL));
00723 if ( isVPeriodic )
00724 res.SetY( res.Y() + ShapeAnalysis::AdjustToPeriod(res.Y(),VF,VL));
00725 }
00726
00727 return res;
00728 }
00729
00730
00731
00732
00733
00734 gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
00735 const gp_XY& p1,
00736 const gp_XY& p2)
00737 {
00738
00739
00740
00741 Handle(Geom_Surface) surf = surface;
00742 while ( !surf.IsNull() && surf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
00743 surf = Handle(Geom_RectangularTrimmedSurface)::DownCast( surf )->BasisSurface();
00744
00745 return applyIn2D( surf, p1, p2, & AverageUV );
00746 }
00747
00748
00749
00750
00751
00752
00753 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E,
00754 const SMDS_MeshNode* n,
00755 const SMDS_MeshNode* inEdgeNode,
00756 bool* check)
00757 {
00758 double param = 0;
00759 const SMDS_PositionPtr pos = n->GetPosition();
00760 if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
00761 {
00762 const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
00763 param = epos->GetUParameter();
00764 }
00765 else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
00766 {
00767 if ( inEdgeNode && TopExp::FirstVertex( E ).IsSame( TopExp::LastVertex( E )))
00768 {
00769 Standard_Real f,l;
00770 BRep_Tool::Range( E, f,l );
00771 double uInEdge = GetNodeU( E, inEdgeNode );
00772 param = ( fabs( uInEdge - f ) < fabs( l - uInEdge )) ? f : l;
00773 }
00774 else
00775 {
00776 SMESHDS_Mesh * meshDS = GetMeshDS();
00777 int vertexID = n->getshapeId();
00778 const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
00779 param = BRep_Tool::Parameter( V, E );
00780 }
00781 }
00782 if ( check )
00783 {
00784 double tol = BRep_Tool::Tolerance( E );
00785 double f,l; BRep_Tool::Range( E, f,l );
00786 bool force = ( param < f-tol || param > l+tol );
00787 if ( !force && pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
00788 force = ( GetMeshDS()->ShapeToIndex( E ) != n->getshapeId() );
00789
00790 *check = CheckNodeU( E, n, param, 2*tol, force );
00791 }
00792 return param;
00793 }
00794
00795
00796
00797
00798
00799
00800
00801 bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E,
00802 const SMDS_MeshNode* n,
00803 double& u,
00804 const double tol,
00805 const bool force,
00806 double distXYZ[4]) const
00807 {
00808 int shapeID = n->getshapeId();
00809 if ( force || toCheckPosOnShape( shapeID ))
00810 {
00811 TopLoc_Location loc; double f,l;
00812 Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
00813 if ( curve.IsNull() )
00814 {
00815 if ( u+tol < f || u-tol > l )
00816 {
00817 double r = Max( 0.5, 1 - tol*n->GetID());
00818 u = f*r + l*(1-r);
00819 }
00820 }
00821 else
00822 {
00823 gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
00824 if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
00825 gp_Pnt curvPnt = curve->Value( u );
00826 double dist = nodePnt.Distance( curvPnt );
00827 if ( distXYZ ) {
00828 curvPnt.Transform( loc );
00829 distXYZ[0] = dist;
00830 distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
00831 }
00832 if ( dist > tol )
00833 {
00834 setPosOnShapeValidity( shapeID, false );
00835
00836 int edgeID = GetMeshDS()->ShapeToIndex( E );
00837 TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector );
00838 TID2ProjectorOnCurve::iterator i_proj =
00839 i2proj.insert( make_pair( edgeID, (GeomAPI_ProjectPointOnCurve*) 0 )).first;
00840 if ( !i_proj->second )
00841 {
00842 i_proj->second = new GeomAPI_ProjectPointOnCurve();
00843 i_proj->second->Init( curve, f, l );
00844 }
00845 GeomAPI_ProjectPointOnCurve* projector = i_proj->second;
00846 projector->Perform( nodePnt );
00847 if ( projector->NbPoints() < 1 )
00848 {
00849 MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
00850 return false;
00851 }
00852 Quantity_Parameter U = projector->LowerDistanceParameter();
00853 u = double( U );
00854 curvPnt = curve->Value( u );
00855 dist = nodePnt.Distance( curvPnt );
00856 if ( distXYZ ) {
00857 curvPnt.Transform( loc );
00858 distXYZ[0] = dist;
00859 distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
00860 }
00861 if ( dist > tol )
00862 {
00863 MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
00864 MESSAGE("distance " << dist << " " << tol );
00865 return false;
00866 }
00867
00868 if ( myShape.IsSame(E) && shapeID == myShapeID )
00869 const_cast<SMDS_MeshNode*>(n)->SetPosition
00870 ( SMDS_PositionPtr( new SMDS_EdgePosition( U )));
00871 }
00872 else if ( fabs( u ) > numeric_limits<double>::min() )
00873 {
00874 setPosOnShapeValidity( shapeID, true );
00875 }
00876 if (( u < f-tol || u > l+tol ) && force )
00877 {
00878
00879 try
00880 {
00881
00882 double period = curve->Period();
00883 u = ( u < f ) ? u + period : u - period;
00884 }
00885 catch (Standard_Failure& exc)
00886 {
00887 return false;
00888 }
00889 }
00890 }
00891 }
00892 return true;
00893 }
00894
00895
00896
00897
00898
00899
00900 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
00901 const SMDS_MeshNode* n2,
00902 bool force3d)
00903 {
00904
00905
00906 SMESH_TLink link(n1,n2);
00907 ItTLinkNode itLN = myTLinkNodeMap.find( link );
00908 if ( itLN != myTLinkNodeMap.end() ) {
00909 return (*itLN).second;
00910 }
00911
00912
00913
00914 SMDS_MeshNode* n12;
00915 SMESHDS_Mesh* meshDS = GetMeshDS();
00916
00917 if ( IsSeamShape( n1->getshapeId() ))
00918
00919 std::swap( n1, n2 );
00920
00921
00922 int faceID = -1, edgeID = -1;
00923 const SMDS_PositionPtr Pos1 = n1->GetPosition();
00924 const SMDS_PositionPtr Pos2 = n2->GetPosition();
00925
00926 TopoDS_Edge E; double u [2];
00927 TopoDS_Face F; gp_XY uv[2];
00928 bool uvOK[2] = { false, false };
00929
00930 if( myShape.IsNull() )
00931 {
00932 if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
00933 faceID = n1->getshapeId();
00934 }
00935 else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
00936 faceID = n2->getshapeId();
00937 }
00938
00939 if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
00940 edgeID = n1->getshapeId();
00941 }
00942 if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
00943 edgeID = n2->getshapeId();
00944 }
00945 }
00946
00947 TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
00948 if ( faceID>0 || shapeType == TopAbs_FACE)
00949 {
00950 if( myShape.IsNull() )
00951 F = TopoDS::Face(meshDS->IndexToShape(faceID));
00952 else {
00953 F = TopoDS::Face(myShape);
00954 faceID = myShapeID;
00955 }
00956 uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
00957 uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
00958 }
00959 else if (edgeID>0 || shapeType == TopAbs_EDGE)
00960 {
00961 if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
00962 Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
00963 n1->getshapeId() != n2->getshapeId() )
00964 return getMediumNodeOnComposedWire(n1,n2,force3d);
00965
00966 if( myShape.IsNull() )
00967 E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
00968 else {
00969 E = TopoDS::Edge(myShape);
00970 edgeID = myShapeID;
00971 }
00972 u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
00973 u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
00974 }
00975 if(!force3d)
00976 {
00977
00978
00979 if( ! F.IsNull() )
00980 {
00981 if ( uvOK[0] && uvOK[1] )
00982 {
00983 if ( IsDegenShape( n1->getshapeId() )) {
00984 if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
00985 else uv[0].SetCoord( 2, uv[1].Coord( 2 ));
00986 }
00987 else if ( IsDegenShape( n2->getshapeId() )) {
00988 if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
00989 else uv[1].SetCoord( 2, uv[0].Coord( 2 ));
00990 }
00991
00992 TopLoc_Location loc;
00993 Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
00994 gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
00995 gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
00996 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
00997 meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
00998 myTLinkNodeMap.insert(make_pair(link,n12));
00999 return n12;
01000 }
01001 }
01002 else if ( !E.IsNull() )
01003 {
01004 double f,l;
01005 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
01006 if(!C.IsNull())
01007 {
01008 Standard_Boolean isPeriodic = C->IsPeriodic();
01009 double U;
01010 if(isPeriodic) {
01011 Standard_Real Period = C->Period();
01012 Standard_Real p = u[1]+ShapeAnalysis::AdjustByPeriod(u[1],u[0],Period);
01013 Standard_Real pmid = (u[0]+p)/2.;
01014 U = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
01015 }
01016 else
01017 U = (u[0]+u[1])/2.;
01018
01019 gp_Pnt P = C->Value( U );
01020 n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
01021 meshDS->SetNodeOnEdge(n12, edgeID, U);
01022 myTLinkNodeMap.insert(make_pair(link,n12));
01023 return n12;
01024 }
01025 }
01026 }
01027
01028
01029 double x = ( n1->X() + n2->X() )/2.;
01030 double y = ( n1->Y() + n2->Y() )/2.;
01031 double z = ( n1->Z() + n2->Z() )/2.;
01032 n12 = meshDS->AddNode(x,y,z);
01033
01034 if ( !F.IsNull() )
01035 {
01036 gp_XY UV = ( uv[0] + uv[1] ) / 2.;
01037 CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), true);
01038 meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() );
01039 }
01040 else if ( !E.IsNull() )
01041 {
01042 double U = ( u[0] + u[1] ) / 2.;
01043 CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), true);
01044 meshDS->SetNodeOnEdge(n12, edgeID, U);
01045 }
01046 else if ( myShapeID > 0 )
01047 {
01048 meshDS->SetNodeInVolume(n12, myShapeID);
01049 }
01050
01051 myTLinkNodeMap.insert( make_pair( link, n12 ));
01052 return n12;
01053 }
01054
01055
01059
01060
01061 const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_MeshNode* n1,
01062 const SMDS_MeshNode* n2,
01063 bool force3d)
01064 {
01065 gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2);
01066 SMDS_MeshNode* n12 = AddNode( middle.X(), middle.Y(), middle.Z() );
01067
01068
01069
01070
01071 double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4];
01072 int iOkEdge = 0;
01073 TopoDS_Edge edges[2];
01074 for ( int is2nd = 0; is2nd < 2; ++is2nd )
01075 {
01076
01077 const SMDS_MeshNode* n = is2nd ? n2 : n1;
01078 TopoDS_Shape shape = GetSubShapeByNode( n, GetMeshDS() );
01079 if ( shape.IsNull() || shape.ShapeType() != TopAbs_EDGE )
01080 continue;
01081
01082
01083 TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
01084 double node2MiddleDist = middle.Distance( XYZ(n) );
01085 double foundU = GetNodeU( edge, n );
01086 CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), true, distXYZ );
01087 if ( distXYZ[0] < node2MiddleDist )
01088 {
01089 distMiddleProj = distXYZ[0];
01090 u = foundU;
01091 iOkEdge = is2nd;
01092 }
01093 }
01094 if ( Precision::IsInfinite( distMiddleProj ))
01095 {
01096
01097 TopoDS_Vertex vCommon;
01098 if ( TopExp::CommonVertex( edges[0], edges[1], vCommon ))
01099 u = BRep_Tool::Parameter( vCommon, edges[0] );
01100 else
01101 {
01102 double f,l, u0 = GetNodeU( edges[0], n1 );
01103 BRep_Tool::Range( edges[0],f,l );
01104 u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l;
01105 }
01106 iOkEdge = 0;
01107 distMiddleProj = 0;
01108 }
01109
01110
01111 double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
01112 if ( !force3d && distMiddleProj > 2*tol )
01113 {
01114 TopLoc_Location loc; double f,l;
01115 Handle(Geom_Curve) curve = BRep_Tool::Curve( edges[iOkEdge],loc,f,l );
01116 gp_Pnt p = curve->Value( u );
01117 GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
01118 }
01119
01120 GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u);
01121
01122 myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 ));
01123
01124 return n12;
01125 }
01126
01127
01128
01129
01130
01131
01132 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
01133 {
01134 SMESHDS_Mesh * meshDS = GetMeshDS();
01135 SMDS_MeshNode* node = 0;
01136 if ( ID )
01137 node = meshDS->AddNodeWithID( x, y, z, ID );
01138 else
01139 node = meshDS->AddNode( x, y, z );
01140 if ( mySetElemOnShape && myShapeID > 0 ) {
01141 switch ( myShape.ShapeType() ) {
01142 case TopAbs_SOLID: meshDS->SetNodeInVolume( node, myShapeID); break;
01143 case TopAbs_SHELL: meshDS->SetNodeInVolume( node, myShapeID); break;
01144 case TopAbs_FACE: meshDS->SetNodeOnFace( node, myShapeID); break;
01145 case TopAbs_EDGE: meshDS->SetNodeOnEdge( node, myShapeID); break;
01146 case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
01147 default: ;
01148 }
01149 }
01150 return node;
01151 }
01152
01153
01154
01155
01156
01157
01158 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
01159 const SMDS_MeshNode* n2,
01160 const int id,
01161 const bool force3d)
01162 {
01163 SMESHDS_Mesh * meshDS = GetMeshDS();
01164
01165 SMDS_MeshEdge* edge = 0;
01166 if (myCreateQuadratic) {
01167 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01168 if(id)
01169 edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
01170 else
01171 edge = meshDS->AddEdge(n1, n2, n12);
01172 }
01173 else {
01174 if(id)
01175 edge = meshDS->AddEdgeWithID(n1, n2, id);
01176 else
01177 edge = meshDS->AddEdge(n1, n2);
01178 }
01179
01180 if ( mySetElemOnShape && myShapeID > 0 )
01181 meshDS->SetMeshElementOnShape( edge, myShapeID );
01182
01183 return edge;
01184 }
01185
01186
01187
01188
01189
01190
01191 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
01192 const SMDS_MeshNode* n2,
01193 const SMDS_MeshNode* n3,
01194 const int id,
01195 const bool force3d)
01196 {
01197 SMESHDS_Mesh * meshDS = GetMeshDS();
01198 SMDS_MeshFace* elem = 0;
01199
01200 if( n1==n2 || n2==n3 || n3==n1 )
01201 return elem;
01202
01203 if(!myCreateQuadratic) {
01204 if(id)
01205 elem = meshDS->AddFaceWithID(n1, n2, n3, id);
01206 else
01207 elem = meshDS->AddFace(n1, n2, n3);
01208 }
01209 else {
01210 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01211 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01212 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01213
01214 if(id)
01215 elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
01216 else
01217 elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
01218 }
01219 if ( mySetElemOnShape && myShapeID > 0 )
01220 meshDS->SetMeshElementOnShape( elem, myShapeID );
01221
01222 return elem;
01223 }
01224
01225
01226
01227
01228
01229
01230 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
01231 const SMDS_MeshNode* n2,
01232 const SMDS_MeshNode* n3,
01233 const SMDS_MeshNode* n4,
01234 const int id,
01235 const bool force3d)
01236 {
01237 SMESHDS_Mesh * meshDS = GetMeshDS();
01238 SMDS_MeshFace* elem = 0;
01239
01240 if( n1==n2 ) {
01241 return AddFace(n1,n3,n4,id,force3d);
01242 }
01243 if( n1==n3 ) {
01244 return AddFace(n1,n2,n4,id,force3d);
01245 }
01246 if( n1==n4 ) {
01247 return AddFace(n1,n2,n3,id,force3d);
01248 }
01249 if( n2==n3 ) {
01250 return AddFace(n1,n2,n4,id,force3d);
01251 }
01252 if( n2==n4 ) {
01253 return AddFace(n1,n2,n3,id,force3d);
01254 }
01255 if( n3==n4 ) {
01256 return AddFace(n1,n2,n3,id,force3d);
01257 }
01258
01259 if(!myCreateQuadratic) {
01260 if(id)
01261 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
01262 else
01263 elem = meshDS->AddFace(n1, n2, n3, n4);
01264 }
01265 else {
01266 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01267 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01268 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01269 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01270
01271 if(id)
01272 elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
01273 else
01274 elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
01275 }
01276 if ( mySetElemOnShape && myShapeID > 0 )
01277 meshDS->SetMeshElementOnShape( elem, myShapeID );
01278
01279 return elem;
01280 }
01281
01282
01283
01284
01285
01286
01287 SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_MeshNode*>& nodes,
01288 const int id,
01289 const bool force3d)
01290 {
01291 SMESHDS_Mesh * meshDS = GetMeshDS();
01292 SMDS_MeshFace* elem = 0;
01293
01294 if(!myCreateQuadratic) {
01295 if(id)
01296 elem = meshDS->AddPolygonalFaceWithID(nodes, id);
01297 else
01298 elem = meshDS->AddPolygonalFace(nodes);
01299 }
01300 else {
01301 vector<const SMDS_MeshNode*> newNodes;
01302 for ( int i = 0; i < nodes.size(); ++i )
01303 {
01304 const SMDS_MeshNode* n1 = nodes[i];
01305 const SMDS_MeshNode* n2 = nodes[(i+1)/nodes.size()];
01306 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01307 newNodes.push_back( n1 );
01308 newNodes.push_back( n12 );
01309 }
01310 if(id)
01311 elem = meshDS->AddPolygonalFaceWithID(newNodes, id);
01312 else
01313 elem = meshDS->AddPolygonalFace(newNodes);
01314 }
01315 if ( mySetElemOnShape && myShapeID > 0 )
01316 meshDS->SetMeshElementOnShape( elem, myShapeID );
01317
01318 return elem;
01319 }
01320
01321
01322
01323
01324
01325
01326 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01327 const SMDS_MeshNode* n2,
01328 const SMDS_MeshNode* n3,
01329 const SMDS_MeshNode* n4,
01330 const SMDS_MeshNode* n5,
01331 const SMDS_MeshNode* n6,
01332 const int id,
01333 const bool force3d)
01334 {
01335 SMESHDS_Mesh * meshDS = GetMeshDS();
01336 SMDS_MeshVolume* elem = 0;
01337 if(!myCreateQuadratic) {
01338 if(id)
01339 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
01340 else
01341 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
01342 }
01343 else {
01344 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01345 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01346 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01347
01348 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
01349 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
01350 const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
01351
01352 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
01353 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
01354 const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
01355
01356 if(id)
01357 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
01358 n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
01359 else
01360 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
01361 n12, n23, n31, n45, n56, n64, n14, n25, n36);
01362 }
01363 if ( mySetElemOnShape && myShapeID > 0 )
01364 meshDS->SetMeshElementOnShape( elem, myShapeID );
01365
01366 return elem;
01367 }
01368
01369
01370
01371
01372
01373
01374 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01375 const SMDS_MeshNode* n2,
01376 const SMDS_MeshNode* n3,
01377 const SMDS_MeshNode* n4,
01378 const int id,
01379 const bool force3d)
01380 {
01381 SMESHDS_Mesh * meshDS = GetMeshDS();
01382 SMDS_MeshVolume* elem = 0;
01383 if(!myCreateQuadratic) {
01384 if(id)
01385 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
01386 else
01387 elem = meshDS->AddVolume(n1, n2, n3, n4);
01388 }
01389 else {
01390 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01391 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01392 const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01393
01394 const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
01395 const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
01396 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01397
01398 if(id)
01399 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
01400 else
01401 elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
01402 }
01403 if ( mySetElemOnShape && myShapeID > 0 )
01404 meshDS->SetMeshElementOnShape( elem, myShapeID );
01405
01406 return elem;
01407 }
01408
01409
01410
01411
01412
01413
01414 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01415 const SMDS_MeshNode* n2,
01416 const SMDS_MeshNode* n3,
01417 const SMDS_MeshNode* n4,
01418 const SMDS_MeshNode* n5,
01419 const int id,
01420 const bool force3d)
01421 {
01422 SMDS_MeshVolume* elem = 0;
01423 if(!myCreateQuadratic) {
01424 if(id)
01425 elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
01426 else
01427 elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
01428 }
01429 else {
01430 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01431 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01432 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01433 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01434
01435 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
01436 const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
01437 const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
01438 const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
01439
01440 if(id)
01441 elem = GetMeshDS()->AddVolumeWithID ( n1, n2, n3, n4, n5,
01442 n12, n23, n34, n41,
01443 n15, n25, n35, n45,
01444 id);
01445 else
01446 elem = GetMeshDS()->AddVolume( n1, n2, n3, n4, n5,
01447 n12, n23, n34, n41,
01448 n15, n25, n35, n45);
01449 }
01450 if ( mySetElemOnShape && myShapeID > 0 )
01451 GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
01452
01453 return elem;
01454 }
01455
01456
01457
01458
01459
01460
01461 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01462 const SMDS_MeshNode* n2,
01463 const SMDS_MeshNode* n3,
01464 const SMDS_MeshNode* n4,
01465 const SMDS_MeshNode* n5,
01466 const SMDS_MeshNode* n6,
01467 const SMDS_MeshNode* n7,
01468 const SMDS_MeshNode* n8,
01469 const int id,
01470 const bool force3d)
01471 {
01472 SMESHDS_Mesh * meshDS = GetMeshDS();
01473 SMDS_MeshVolume* elem = 0;
01474 if(!myCreateQuadratic) {
01475 if(id)
01476 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
01477 else
01478 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
01479 }
01480 else {
01481 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01482 const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01483 const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01484 const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01485
01486 const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
01487 const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
01488 const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
01489 const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
01490
01491 const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
01492 const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
01493 const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
01494 const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
01495
01496 if(id)
01497 elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
01498 n12, n23, n34, n41, n56, n67,
01499 n78, n85, n15, n26, n37, n48, id);
01500 else
01501 elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
01502 n12, n23, n34, n41, n56, n67,
01503 n78, n85, n15, n26, n37, n48);
01504 }
01505 if ( mySetElemOnShape && myShapeID > 0 )
01506 meshDS->SetMeshElementOnShape( elem, myShapeID );
01507
01508 return elem;
01509 }
01510
01511
01512
01513
01514
01515
01516 SMDS_MeshVolume*
01517 SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
01518 const std::vector<int>& quantities,
01519 const int id,
01520 const bool force3d)
01521 {
01522 SMESHDS_Mesh * meshDS = GetMeshDS();
01523 SMDS_MeshVolume* elem = 0;
01524 if(!myCreateQuadratic)
01525 {
01526 if(id)
01527 elem = meshDS->AddPolyhedralVolumeWithID(nodes, quantities, id);
01528 else
01529 elem = meshDS->AddPolyhedralVolume(nodes, quantities);
01530 }
01531 else
01532 {
01533 vector<const SMDS_MeshNode*> newNodes;
01534 vector<int> newQuantities;
01535 for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace)
01536 {
01537 int nbNodesInFace = quantities[iFace];
01538 newQuantities.push_back(0);
01539 for ( int i = 0; i < nbNodesInFace; ++i )
01540 {
01541 const SMDS_MeshNode* n1 = nodes[ iN + i ];
01542 newNodes.push_back( n1 );
01543 newQuantities.back()++;
01544
01545 const SMDS_MeshNode* n2 = nodes[ iN + ( i+1==nbNodesInFace ? 0 : i+1 )];
01546
01547
01548 {
01549 const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01550 newNodes.push_back( n12 );
01551 newQuantities.back()++;
01552 }
01553 }
01554 iN += nbNodesInFace;
01555 }
01556 if(id)
01557 elem = meshDS->AddPolyhedralVolumeWithID( newNodes, newQuantities, id );
01558 else
01559 elem = meshDS->AddPolyhedralVolume( newNodes, newQuantities );
01560 }
01561 if ( mySetElemOnShape && myShapeID > 0 )
01562 meshDS->SetMeshElementOnShape( elem, myShapeID );
01563
01564 return elem;
01565 }
01566
01567
01568
01569
01570
01571
01572 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
01573 const TopoDS_Face& theFace,
01574 const TopoDS_Edge& theBaseEdge,
01575 SMESHDS_Mesh* theMesh,
01576 SMESH_ProxyMesh* theProxyMesh)
01577 {
01578 const SMESHDS_SubMesh* faceSubMesh = 0;
01579 if ( theProxyMesh )
01580 {
01581 faceSubMesh = theProxyMesh->GetSubMesh( theFace );
01582 if ( !faceSubMesh ||
01583 faceSubMesh->NbElements() == 0 ||
01584 theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() ))
01585 {
01586
01587 faceSubMesh = 0;
01588 theProxyMesh = 0;
01589 }
01590 }
01591 if ( !faceSubMesh )
01592 faceSubMesh = theMesh->MeshElements( theFace );
01593 if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
01594 return false;
01595
01596
01597
01598 map< double, const SMDS_MeshNode*> sortedBaseNodes;
01599 if ( !SMESH_Algo::GetSortedNodesOnEdge( theMesh, theBaseEdge,true, sortedBaseNodes)
01600 || sortedBaseNodes.size() < 2 )
01601 return false;
01602
01603 int nbRows = faceSubMesh->NbElements() / ( sortedBaseNodes.size()-1 ) + 1;
01604 map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
01605 double f = u_n->first, range = sortedBaseNodes.rbegin()->first - f;
01606 for ( ; u_n != sortedBaseNodes.end(); u_n++ )
01607 {
01608 double par = ( u_n->first - f ) / range;
01609 vector<const SMDS_MeshNode*>& nCol = theParam2ColumnMap[ par ];
01610 nCol.resize( nbRows );
01611 nCol[0] = u_n->second;
01612 }
01613 TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin();
01614 if ( theProxyMesh )
01615 {
01616 for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 )
01617 {
01618 const SMDS_MeshNode* & n = par_nVec_1->second[0];
01619 n = theProxyMesh->GetProxyNode( n );
01620 }
01621 }
01622
01623
01624
01625
01626 par_nVec_2 = theParam2ColumnMap.begin();
01627 par_nVec_1 = par_nVec_2++;
01628 TIDSortedElemSet emptySet, avoidSet;
01629 for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
01630 {
01631 vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
01632 vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
01633
01634 int i1, i2, iRow = 0;
01635 const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0];
01636
01637 while ( const SMDS_MeshElement* face =
01638 SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2))
01639 {
01640 if ( faceSubMesh->Contains( face ))
01641 {
01642 int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
01643 if ( nbNodes != 4 )
01644 return false;
01645 n1 = face->GetNode( (i2+2) % 4 );
01646 n2 = face->GetNode( (i1+2) % 4 );
01647 if ( ++iRow >= nbRows )
01648 return false;
01649 nCol1[ iRow ] = n1;
01650 nCol2[ iRow ] = n2;
01651 avoidSet.clear();
01652 }
01653 avoidSet.insert( face );
01654 }
01655 if ( iRow + 1 < nbRows )
01656 nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 );
01657 }
01658 return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
01659 }
01660
01661
01662
01663
01664
01665
01666 int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
01667 const SMESH_Mesh& mesh,
01668 TopAbs_ShapeEnum ancestorType)
01669 {
01670 TopTools_MapOfShape ancestors;
01671 TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
01672 for ( ; ansIt.More(); ansIt.Next() ) {
01673 if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
01674 ancestors.Add( ansIt.Value() );
01675 }
01676 return ancestors.Extent();
01677 }
01678
01679
01680
01681
01682
01683
01684 TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
01685 const TopoDS_Shape& subShape)
01686 {
01687 TopAbs_Orientation ori = TopAbs_Orientation(-1);
01688 if ( !shape.IsNull() && !subShape.IsNull() )
01689 {
01690 TopExp_Explorer e( shape, subShape.ShapeType() );
01691 if ( shape.Orientation() >= TopAbs_INTERNAL )
01692 e.Init( shape.Oriented(TopAbs_FORWARD), subShape.ShapeType() );
01693 for ( ; e.More(); e.Next())
01694 if ( subShape.IsSame( e.Current() ))
01695 break;
01696 if ( e.More() )
01697 ori = e.Current().Orientation();
01698 }
01699 return ori;
01700 }
01701
01702
01703
01704
01705
01706
01707 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
01708 const TopoDS_Shape& mainShape )
01709 {
01710 if ( !shape.IsNull() && !mainShape.IsNull() )
01711 {
01712 for ( TopExp_Explorer exp( mainShape, shape.ShapeType());
01713 exp.More();
01714 exp.Next() )
01715 if ( shape.IsSame( exp.Current() ))
01716 return true;
01717 }
01718 SCRUTE((shape.IsNull()));
01719 SCRUTE((mainShape.IsNull()));
01720 return false;
01721 }
01722
01723
01724
01725
01726
01727
01728 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh )
01729 {
01730 if ( shape.IsNull() || !aMesh )
01731 return false;
01732 return
01733 aMesh->GetMeshDS()->ShapeToIndex( shape ) ||
01734
01735 (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ));
01736 }
01737
01738
01742
01743
01744 double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
01745 {
01746 double tol = Precision::Confusion();
01747 TopExp_Explorer exp;
01748 for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() )
01749 tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Face( exp.Current())));
01750 for ( exp.Init( shape, TopAbs_EDGE ); exp.More(); exp.Next() )
01751 tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Edge( exp.Current())));
01752 for ( exp.Init( shape, TopAbs_VERTEX ); exp.More(); exp.Next() )
01753 tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current())));
01754
01755 return tol;
01756 }
01757
01758
01764
01765
01766 bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge )
01767 {
01768 if ( anEdge.Orientation() >= TopAbs_INTERNAL )
01769 return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD )));
01770 return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge ));
01771 }
01772
01773
01778
01779
01780 TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool is2nd,
01781 TopoDS_Edge anEdge,
01782 const bool CumOri )
01783 {
01784 if ( anEdge.Orientation() >= TopAbs_INTERNAL )
01785 anEdge.Orientation( TopAbs_FORWARD );
01786
01787 const TopAbs_Orientation tgtOri = is2nd ? TopAbs_REVERSED : TopAbs_FORWARD;
01788 TopoDS_Iterator vIt( anEdge, CumOri );
01789 while ( vIt.More() && vIt.Value().Orientation() != tgtOri )
01790 vIt.Next();
01791
01792 return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() );
01793 }
01794
01795
01796
01797
01798
01799
01800
01801
01802 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
01803 {
01804 int NbAllEdgsAndFaces=0;
01805 int NbQuadFacesAndEdgs=0;
01806 int NbFacesAndEdges=0;
01807
01808 NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
01809
01810
01811 NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
01812
01813
01814 NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
01815
01816 if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
01817
01818 return SMESH_MesherHelper::QUADRATIC;
01819 }
01820 else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
01821
01822 return SMESH_MesherHelper::LINEAR;
01823 }
01824 else
01825
01826 return SMESH_MesherHelper::COMP;
01827 }
01828
01829
01830
01831
01832
01833
01834 double SMESH_MesherHelper::GetOtherParam(const double param) const
01835 {
01836 int i = myParIndex & U_periodic ? 0 : 1;
01837 return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
01838 }
01839
01840 namespace {
01841
01842
01846
01847
01848 struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
01849 {
01850 TopTools_ListIteratorOfListOfShape _ancIter;
01851 TopAbs_ShapeEnum _type;
01852 TopTools_MapOfShape _encountered;
01853 TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
01854 : _ancIter( ancestors ), _type( type )
01855 {
01856 if ( _ancIter.More() ) {
01857 if ( _ancIter.Value().ShapeType() != _type ) next();
01858 else _encountered.Add( _ancIter.Value() );
01859 }
01860 }
01861 virtual bool more()
01862 {
01863 return _ancIter.More();
01864 }
01865 virtual const TopoDS_Shape* next()
01866 {
01867 const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
01868 if ( _ancIter.More() )
01869 for ( _ancIter.Next(); _ancIter.More(); _ancIter.Next())
01870 if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
01871 break;
01872 return s;
01873 }
01874 };
01875
01876 }
01877
01878
01882
01883
01884 PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
01885 const SMESH_Mesh& mesh,
01886 TopAbs_ShapeEnum ancestorType)
01887 {
01888 return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
01889 }
01890
01891
01892
01893
01894 namespace {
01895
01896
01897 #define __DMP__(txt) \
01898 //cout << txt
01899 #define MSG(txt) __DMP__(txt<<endl)
01900 #define MSGBEG(txt) __DMP__(txt)
01901
01902
01903 bool isStraightLink(double linkLen2, double middleNodeMove2)
01904 {
01905
01906 return middleNodeMove2 < 1/15./15. * linkLen2;
01907 }
01908
01909 struct QFace;
01910
01914 struct QLink: public SMESH_TLink
01915 {
01916 const SMDS_MeshNode* _mediumNode;
01917 mutable vector<const QFace* > _faces;
01918 mutable gp_Vec _nodeMove;
01919 mutable int _nbMoves;
01920
01921 QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm):
01922 SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) {
01923 _faces.reserve(4);
01924
01925 _nodeMove = MediumPnt() - MiddlePnt();
01926 }
01927 void SetContinuesFaces() const;
01928 const QFace* GetContinuesFace( const QFace* face ) const;
01929 bool OnBoundary() const;
01930 gp_XYZ MiddlePnt() const { return ( XYZ( node1() ) + XYZ( node2() )) / 2.; }
01931 gp_XYZ MediumPnt() const { return XYZ( _mediumNode ); }
01932
01933 SMDS_TypeOfPosition MediumPos() const
01934 { return _mediumNode->GetPosition()->GetTypeOfPosition(); }
01935 SMDS_TypeOfPosition EndPos(bool isSecond) const
01936 { return (isSecond ? node2() : node1())->GetPosition()->GetTypeOfPosition(); }
01937 const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const
01938 { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; }
01939
01940 void Move(const gp_Vec& move, bool sum=false) const
01941 { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
01942 gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
01943 bool IsMoved() const { return (_nbMoves > 0 ); }
01944 bool IsStraight() const
01945 { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
01946 _nodeMove.SquareMagnitude());
01947 }
01948 bool operator<(const QLink& other) const {
01949 return (node1()->GetID() == other.node1()->GetID() ?
01950 node2()->GetID() < other.node2()->GetID() :
01951 node1()->GetID() < other.node1()->GetID());
01952 }
01953
01954
01955
01956 };
01957
01961 struct TChainLink
01962 {
01963 const QLink* _qlink;
01964 mutable const QFace* _qfaces[2];
01965
01966 TChainLink(const QLink* qlink=0):_qlink(qlink) {
01967 _qfaces[0] = _qfaces[1] = 0;
01968 }
01969 void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
01970
01971 bool IsBoundary() const { return !_qfaces[1]; }
01972
01973 void RemoveFace( const QFace* face ) const
01974 { _qfaces[(face == _qfaces[1])] = 0; if (!_qfaces[0]) std::swap(_qfaces[0],_qfaces[1]); }
01975
01976 const QFace* NextFace( const QFace* f ) const
01977 { return _qfaces[0]==f ? _qfaces[1] : _qfaces[0]; }
01978
01979 const SMDS_MeshNode* NextNode( const SMDS_MeshNode* n ) const
01980 { return n == _qlink->node1() ? _qlink->node2() : _qlink->node1(); }
01981
01982 bool operator<(const TChainLink& other) const { return *_qlink < *other._qlink; }
01983
01984 operator bool() const { return (_qlink); }
01985
01986 const QLink* operator->() const { return _qlink; }
01987
01988 gp_Vec Normal() const;
01989
01990 bool IsStraight() const;
01991 };
01992
01993 typedef list< TChainLink > TChain;
01994 typedef set < TChainLink > TLinkSet;
01995 typedef TLinkSet::const_iterator TLinkInSet;
01996
01997 const int theFirstStep = 5;
01998
01999 enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN };
02000
02004 struct QFace: public TIDSortedNodeSet
02005 {
02006 mutable const SMDS_MeshElement* _volumes[2];
02007 mutable vector< const QLink* > _sides;
02008 mutable bool _sideIsAdded[4];
02009 gp_Vec _normal;
02010 #ifdef _DEBUG_
02011 mutable const SMDS_MeshElement* _face;
02012 #endif
02013
02014 QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face=0 );
02015
02016 void SetVolume(const SMDS_MeshElement* v) const { _volumes[ _volumes[0] ? 1 : 0 ] = v; }
02017
02018 int NbVolumes() const { return !_volumes[0] ? 0 : !_volumes[1] ? 1 : 2; }
02019
02020 void AddSelfToLinks() const {
02021 for ( int i = 0; i < _sides.size(); ++i )
02022 _sides[i]->_faces.push_back( this );
02023 }
02024 int LinkIndex( const QLink* side ) const {
02025 for (int i=0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i;
02026 return -1;
02027 }
02028 bool GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& err) const;
02029
02030 bool GetLinkChain( TChainLink& link, TChain& chain, SMDS_TypeOfPosition pos, int& err) const
02031 {
02032 int i = LinkIndex( link._qlink );
02033 if ( i < 0 ) return true;
02034 _sideIsAdded[i] = true;
02035 link.SetFace( this );
02036
02037 return GetLinkChain( (i+2)%_sides.size(), chain, pos, err );
02038 }
02039 bool IsBoundary() const { return !_volumes[1]; }
02040
02041 bool Contains( const SMDS_MeshNode* node ) const { return count(node); }
02042
02043 bool IsSpoiled(const QLink* bentLink ) const;
02044
02045 TLinkInSet GetBoundaryLink( const TLinkSet& links,
02046 const TChainLink& avoidLink,
02047 TLinkInSet * notBoundaryLink = 0,
02048 const SMDS_MeshNode* nodeToContain = 0,
02049 bool * isAdjacentUsed = 0,
02050 int nbRecursionsLeft = -1) const;
02051
02052 TLinkInSet GetLinkByNode( const TLinkSet& links,
02053 const TChainLink& avoidLink,
02054 const SMDS_MeshNode* nodeToContain) const;
02055
02056 const SMDS_MeshNode* GetNodeInFace() const {
02057 for ( int iL = 0; iL < _sides.size(); ++iL )
02058 if ( _sides[iL]->MediumPos() == SMDS_TOP_FACE ) return _sides[iL]->_mediumNode;
02059 return 0;
02060 }
02061
02062 gp_Vec LinkNorm(const int i, SMESH_MesherHelper* theFaceHelper=0) const;
02063
02064 double MoveByBoundary( const TChainLink& theLink,
02065 const gp_Vec& theRefVec,
02066 const TLinkSet& theLinks,
02067 SMESH_MesherHelper* theFaceHelper=0,
02068 const double thePrevLen=0,
02069 const int theStep=theFirstStep,
02070 gp_Vec* theLinkNorm=0,
02071 double theSign=1.0) const;
02072 };
02073
02074
02078 ostream& operator << (ostream& out, const QLink& l)
02079 {
02080 out <<"QLink nodes: "
02081 << l.node1()->GetID() << " - "
02082 << l._mediumNode->GetID() << " - "
02083 << l.node2()->GetID() << endl;
02084 return out;
02085 }
02086 ostream& operator << (ostream& out, const QFace& f)
02087 {
02088 out <<"QFace nodes: ";
02089 for ( TIDSortedNodeSet::const_iterator n = f.begin(); n != f.end(); ++n )
02090 out << (*n)->GetID() << " ";
02091 out << " \tvolumes: "
02092 << (f._volumes[0] ? f._volumes[0]->GetID() : 0) << " "
02093 << (f._volumes[1] ? f._volumes[1]->GetID() : 0);
02094 out << " \tNormal: "<< f._normal.X() <<", "<<f._normal.Y() <<", "<<f._normal.Z() << endl;
02095 return out;
02096 }
02097
02098
02102
02103
02104 QFace::QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face )
02105 {
02106 _volumes[0] = _volumes[1] = 0;
02107 _sides = links;
02108 _sideIsAdded[0]=_sideIsAdded[1]=_sideIsAdded[2]=_sideIsAdded[3]=false;
02109 _normal.SetCoord(0,0,0);
02110 for ( int i = 1; i < _sides.size(); ++i ) {
02111 const QLink *l1 = _sides[i-1], *l2 = _sides[i];
02112 insert( l1->node1() ); insert( l1->node2() );
02113
02114 gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1()));
02115 gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2()));
02116 if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() )
02117 v1.Reverse();
02118 _normal += v1 ^ v2;
02119 }
02120 double normSqSize = _normal.SquareMagnitude();
02121 if ( normSqSize > numeric_limits<double>::min() )
02122 _normal /= sqrt( normSqSize );
02123 else
02124 _normal.SetCoord(1e-33,0,0);
02125
02126 #ifdef _DEBUG_
02127 _face = face;
02128 #endif
02129 }
02130
02140
02141
02142 bool QFace::GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& error) const
02143 {
02144 if ( iSide >= _sides.size() )
02145 return false;
02146 if ( _sideIsAdded[ iSide ])
02147 return true;
02148
02149 if ( _sides.size() != 4 ) {
02150 MSGBEG( *this );
02151 TLinkSet links;
02152 list< const QFace* > faces( 1, this );
02153 while ( !faces.empty() ) {
02154 const QFace* face = faces.front();
02155 for ( int i = 0; i < face->_sides.size(); ++i ) {
02156 if ( !face->_sideIsAdded[i] && face->_sides[i] ) {
02157 face->_sideIsAdded[i] = true;
02158
02159 TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first;
02160
02161
02162
02163
02164
02165
02166
02167 chLink->SetFace( face );
02168 if ( face->_sides[i]->MediumPos() == pos )
02169 if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face ))
02170 if ( contFace->_sides.size() == 3 )
02171 faces.push_back( contFace );
02172 }
02173 }
02174 faces.pop_front();
02175 }
02176 if ( error < ERR_TRI )
02177 error = ERR_TRI;
02178 chain.insert( chain.end(), links.begin(),links.end() );
02179 return false;
02180 }
02181 _sideIsAdded[iSide] = true;
02182 const QLink* link = _sides[iSide];
02183 if ( !link)
02184 return true;
02185
02186
02187 TChain::iterator chLink = chain.insert( chain.begin(), TChainLink(link));
02188 chLink->SetFace( this );
02189 MSGBEG( *this );
02190
02191
02192 if ( link->MediumPos() >= pos ) {
02193 int nbLinkFaces = link->_faces.size();
02194 if ( nbLinkFaces == 4 || (link->OnBoundary())) {
02195
02196 if ( const QFace* f = link->GetContinuesFace( this ))
02197 if ( f->_sides.size() == 4 )
02198 return f->GetLinkChain( *chLink, chain, pos, error );
02199 }
02200 else {
02201 TChainLink chLink(link);
02202 for ( int i = 0; i < nbLinkFaces; ++i )
02203 if ( link->_faces[i] )
02204 link->_faces[i]->GetLinkChain( chLink, chain, pos, error );
02205 if ( error < ERR_PRISM )
02206 error = ERR_PRISM;
02207 return false;
02208 }
02209 }
02210 return true;
02211 }
02212
02213
02224
02225
02226 TLinkInSet QFace::GetBoundaryLink( const TLinkSet& links,
02227 const TChainLink& avoidLink,
02228 TLinkInSet * notBoundaryLink,
02229 const SMDS_MeshNode* nodeToContain,
02230 bool * isAdjacentUsed,
02231 int nbRecursionsLeft) const
02232 {
02233 TLinkInSet linksEnd = links.end(), boundaryLink = linksEnd;
02234
02235 typedef list< pair< const QFace*, TLinkInSet > > TFaceLinkList;
02236 TFaceLinkList adjacentFaces;
02237
02238 for ( int iL = 0; iL < _sides.size(); ++iL )
02239 {
02240 if ( avoidLink._qlink == _sides[iL] )
02241 continue;
02242 TLinkInSet link = links.find( _sides[iL] );
02243 if ( link == linksEnd ) continue;
02244 if ( (*link)->MediumPos() > SMDS_TOP_FACE )
02245 continue;
02246
02247
02248 if ( link->IsBoundary() ) {
02249 if ( !nodeToContain ||
02250 (*link)->node1() == nodeToContain ||
02251 (*link)->node2() == nodeToContain )
02252 {
02253 boundaryLink = link;
02254 if ( !notBoundaryLink ) break;
02255 }
02256 }
02257 else if ( notBoundaryLink ) {
02258 *notBoundaryLink = link;
02259 if ( boundaryLink != linksEnd ) break;
02260 }
02261
02262 if ( boundaryLink == linksEnd && nodeToContain )
02263 if ( const QFace* adj = link->NextFace( this ))
02264 if ( adj->Contains( nodeToContain ))
02265 adjacentFaces.push_back( make_pair( adj, link ));
02266 }
02267
02268 if ( isAdjacentUsed ) *isAdjacentUsed = false;
02269 if ( boundaryLink == linksEnd && nodeToContain && nbRecursionsLeft)
02270 {
02271 if ( nbRecursionsLeft < 0 )
02272 nbRecursionsLeft = nodeToContain->NbInverseElements();
02273 TFaceLinkList::iterator adj = adjacentFaces.begin();
02274 for ( ; boundaryLink == linksEnd && adj != adjacentFaces.end(); ++adj )
02275 boundaryLink = adj->first->GetBoundaryLink( links, *(adj->second), 0, nodeToContain,
02276 isAdjacentUsed, nbRecursionsLeft-1);
02277 if ( isAdjacentUsed ) *isAdjacentUsed = true;
02278 }
02279 return boundaryLink;
02280 }
02281
02285
02286
02287 TLinkInSet QFace::GetLinkByNode( const TLinkSet& links,
02288 const TChainLink& avoidLink,
02289 const SMDS_MeshNode* nodeToContain) const
02290 {
02291 for ( int i = 0; i < _sides.size(); ++i )
02292 if ( avoidLink._qlink != _sides[i] &&
02293 (_sides[i]->node1() == nodeToContain || _sides[i]->node2() == nodeToContain ))
02294 return links.find( _sides[ i ]);
02295 return links.end();
02296 }
02297
02298
02302
02303
02304 gp_Vec QFace::LinkNorm(const int i, SMESH_MesherHelper* ) const
02305 {
02306 gp_Vec norm, vecOut;
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321 norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
02322 gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) +
02323 XYZ( _sides[0]->node2() ) +
02324 XYZ( _sides[1]->node1() )) / 3.;
02325 vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn );
02326
02327 if ( norm * vecOut < 0 )
02328 norm.Reverse();
02329 double mag2 = norm.SquareMagnitude();
02330 if ( mag2 > numeric_limits<double>::min() )
02331 norm /= sqrt( mag2 );
02332 return norm;
02333 }
02334
02347
02348
02349 double QFace::MoveByBoundary( const TChainLink& theLink,
02350 const gp_Vec& theRefVec,
02351 const TLinkSet& theLinks,
02352 SMESH_MesherHelper* theFaceHelper,
02353 const double thePrevLen,
02354 const int theStep,
02355 gp_Vec* theLinkNorm,
02356 double theSign) const
02357 {
02358 if ( !theStep )
02359 return thePrevLen;
02360
02361 int iL;
02362 for ( iL = 0; iL < _sides.size(); ++iL )
02363 if ( theLink._qlink == _sides[ iL ])
02364 break;
02365
02366 MSG(string(theStep,'.')<<" Ref( "<<theRefVec.X()<<","<<theRefVec.Y()<<","<<theRefVec.Z()<<" )"
02367 <<" thePrevLen " << thePrevLen);
02368 MSG(string(theStep,'.')<<" "<<*theLink._qlink);
02369
02370 gp_Vec linkNorm = -LinkNorm( iL );
02371 double refProj = theRefVec * linkNorm;
02372 if ( theStep == theFirstStep )
02373 theSign = refProj < 0. ? -1. : 1.;
02374 else if ( theSign * refProj < 0.4 * theRefVec.Magnitude())
02375 return thePrevLen;
02376
02377 int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3;
02378 TLinkInSet link1 = theLinks.find( _sides[iL1] );
02379 TLinkInSet link2 = theLinks.find( _sides[iL2] );
02380 if ( link1 == theLinks.end() || link2 == theLinks.end() )
02381 return thePrevLen;
02382 const QFace* f1 = link1->NextFace( this );
02383 const QFace* f2 = link2->NextFace( this );
02384
02385
02386 double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
02387 double len2 = thePrevLen + (theLink->MiddlePnt() - _sides[iL2]->MiddlePnt()).Modulus();
02388 gp_Vec linkDir1(0,0,0);
02389 gp_Vec linkDir2(0,0,0);
02390 try {
02391 OCC_CATCH_SIGNALS;
02392 if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() )
02393 len1 = f1->MoveByBoundary
02394 ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
02395 else
02396 linkDir1 = LinkNorm( iL1 );
02397 } catch (...) {
02398 MSG( " --------------- EXCEPTION");
02399 return thePrevLen;
02400 }
02401 try {
02402 OCC_CATCH_SIGNALS;
02403 if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() )
02404 len2 = f2->MoveByBoundary
02405 ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
02406 else
02407 linkDir2 = LinkNorm( iL2 );
02408 } catch (...) {
02409 MSG( " --------------- EXCEPTION");
02410 return thePrevLen;
02411 }
02412
02413 double fullLen = 0;
02414 if ( theStep != theFirstStep )
02415 {
02416
02417 bool choose1 = ( theRefVec * linkDir1 * theSign > theRefVec * linkDir2 * theSign );
02418 fullLen = choose1 ? len1 : len2;
02419 double r = thePrevLen / fullLen;
02420
02421 gp_Vec move = linkNorm * refProj * ( 1 - r );
02422 theLink->Move( move, true );
02423
02424 MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<<
02425 " by " << refProj * ( 1 - r ) << " following " <<
02426 (choose1 ? *link1->_qlink : *link2->_qlink));
02427
02428 if ( theLinkNorm ) *theLinkNorm = linkNorm;
02429 }
02430 return fullLen;
02431 }
02432
02433
02437
02438
02439 bool QFace::IsSpoiled(const QLink* bentLink ) const
02440 {
02441
02442 gp_XYZ gc(0,0,0);
02443 for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n)
02444 gc += XYZ( *n ) / size();
02445 for (unsigned i = 0; i < _sides.size(); ++i )
02446 {
02447 if ( _sides[i] == bentLink ) continue;
02448 gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
02449 gp_Vec vecOut( gc, _sides[i]->MiddlePnt() );
02450 if ( linkNorm * vecOut < 0 )
02451 linkNorm.Reverse();
02452 double mag2 = linkNorm.SquareMagnitude();
02453 if ( mag2 > numeric_limits<double>::min() )
02454 linkNorm /= sqrt( mag2 );
02455 gp_Vec vecBent ( _sides[i]->MiddlePnt(), bentLink->MediumPnt());
02456 gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt());
02457 if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() )
02458 return true;
02459 }
02460 return false;
02461
02462 }
02463
02464
02468
02469
02470 void QLink::SetContinuesFaces() const
02471 {
02472
02473
02474
02475
02476
02477
02478
02479
02480 if ( _faces.empty() )
02481 return;
02482 int iFaceCont = -1, nbBoundary = 0, iBoundary[2]={-1,-1};
02483 if ( _faces[0]->IsBoundary() )
02484 iBoundary[ nbBoundary++ ] = 0;
02485 for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF )
02486 {
02487
02488 bool sameVol = false;
02489 int nbVol = _faces[iF]->NbVolumes();
02490 for ( int iV = 0; !sameVol && iV < nbVol; ++iV )
02491 sameVol = ( _faces[iF]->_volumes[iV] == _faces[0]->_volumes[0] ||
02492 _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]);
02493 if ( !sameVol )
02494 iFaceCont = iF;
02495 if ( _faces[iF]->IsBoundary() )
02496 iBoundary[ nbBoundary++ ] = iF;
02497 }
02498
02499
02500
02501 if ( nbBoundary == 2 )
02502 {
02503 if (( iBoundary[0] < 2 ) != ( iBoundary[1] < 2 ))
02504 {
02505 int iNear0 = iBoundary[0] < 2 ? 1-iBoundary[0] : 5-iBoundary[0];
02506 std::swap( _faces[ iBoundary[1] ], _faces[iNear0] );
02507 }
02508 }
02509 else if ( iFaceCont > 0 )
02510 {
02511 if ( iFaceCont != 1 )
02512 std::swap( _faces[1], _faces[iFaceCont] );
02513 }
02514 else if ( _faces.size() > 1 )
02515 {
02516 _faces.insert( ++_faces.begin(), 0 );
02517 }
02518 }
02519
02523
02524
02525 const QFace* QLink::GetContinuesFace( const QFace* face ) const
02526 {
02527 for ( int i = 0; i < _faces.size(); ++i ) {
02528 if ( _faces[i] == face ) {
02529 int iF = i < 2 ? 1-i : 5-i;
02530 return iF < _faces.size() ? _faces[iF] : 0;
02531 }
02532 }
02533 return 0;
02534 }
02535
02539
02540
02541 bool QLink::OnBoundary() const
02542 {
02543 for ( int i = 0; i < _faces.size(); ++i )
02544 if (_faces[i] && _faces[i]->IsBoundary()) return true;
02545 return false;
02546 }
02547
02551
02552
02553 gp_Vec TChainLink::Normal() const {
02554 gp_Vec norm;
02555 if (_qfaces[0]) norm = _qfaces[0]->_normal;
02556 if (_qfaces[1]) norm += _qfaces[1]->_normal;
02557 return norm;
02558 }
02559
02563
02564
02565 bool TChainLink::IsStraight() const
02566 {
02567 bool isStraight = _qlink->IsStraight();
02568 if ( isStraight && _qfaces[0] && !_qfaces[1] )
02569 {
02570 int i = _qfaces[0]->LinkIndex( _qlink );
02571 int iOpp = ( i + 2 ) % _qfaces[0]->_sides.size();
02572 gp_XYZ mid1 = _qlink->MiddlePnt();
02573 gp_XYZ mid2 = _qfaces[0]->_sides[ iOpp ]->MiddlePnt();
02574 double faceSize2 = (mid1-mid2).SquareModulus();
02575 isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/3./3. * faceSize2;
02576 }
02577 return isStraight;
02578 }
02579
02580
02584
02585
02586 void fixPrism( TChain& allLinks )
02587 {
02588
02589 typedef set<const QLink*> QLinkSet;
02590 QLinkSet interLinks, bndLinks1, bndLink2;
02591
02592 bool isCurved = false;
02593 for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
02594 if ( (*lnk)->OnBoundary() )
02595 bndLinks1.insert( lnk->_qlink );
02596 else
02597 interLinks.insert( lnk->_qlink );
02598 isCurved = isCurved || !lnk->IsStraight();
02599 }
02600 if ( !isCurved )
02601 return;
02602
02603 QLinkSet *curBndLinks = &bndLinks1, *newBndLinks = &bndLink2;
02604
02605 while ( !interLinks.empty() && !curBndLinks->empty() )
02606 {
02607
02608 QLinkSet::iterator bnd = curBndLinks->begin(), bndEnd = curBndLinks->end();
02609 for ( ; bnd != bndEnd; ++bnd )
02610 {
02611 const QLink* bndLink = *bnd;
02612 for ( int i = 0; i < bndLink->_faces.size(); ++i )
02613 {
02614 const QFace* face = bndLink->_faces[i];
02615 if ( !face ) continue;
02616
02617 int interInd = ( face->LinkIndex( bndLink ) + 2 ) % face->_sides.size();
02618 const QLink* interLink = face->_sides[ interInd ];
02619 QLinkSet::iterator pInterLink = interLinks.find( interLink );
02620 if ( pInterLink == interLinks.end() ) continue;
02621 interLink->Move( bndLink->_nodeMove );
02622
02623 interLinks. erase( pInterLink );
02624 newBndLinks->insert( interLink );
02625 }
02626 }
02627 curBndLinks->clear();
02628 std::swap( curBndLinks, newBndLinks );
02629 }
02630 }
02631
02632
02636
02637
02638 void fixTriaNearBoundary( TChain & allLinks, SMESH_MesherHelper& )
02639 {
02640 if ( allLinks.empty() ) return;
02641
02642 TLinkSet linkSet( allLinks.begin(), allLinks.end());
02643 TLinkInSet linkIt = linkSet.begin(), linksEnd = linkSet.end();
02644
02645 for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt)
02646 {
02647 if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0])
02648 {
02649
02650 const QFace* face = linkIt->_qfaces[0];
02651 gp_XYZ pIn = ( face->_sides[0]->MiddlePnt() +
02652 face->_sides[1]->MiddlePnt() +
02653 face->_sides[2]->MiddlePnt() ) / 3.;
02654 gp_XYZ insideDir( pIn - (*linkIt)->MiddlePnt());
02655 bool linkBentInside = ((*linkIt)->_nodeMove.Dot( insideDir ) > 0 );
02656
02657 if ( linkBentInside )
02658 face->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet );
02659 }
02660 }
02661 }
02662
02663
02667
02668
02669 enum TSplitTriaResult {
02670 _OK, _NO_CORNERS, _FEW_ROWS, _MANY_ROWS, _NO_SIDELINK, _BAD_MIDQUAD, _NOT_RECT,
02671 _NO_MIDQUAD, _NO_UPTRIA, _BAD_SET_SIZE, _BAD_CORNER, _BAD_START, _NO_BOTLINK, _TWISTED_CHAIN };
02672
02673 TSplitTriaResult splitTrianglesIntoChains( TChain & allLinks,
02674 vector< TChain> & resultChains,
02675 SMDS_TypeOfPosition pos )
02676 {
02677
02678 TLinkSet linkSet;
02679 int nbBndLinks = 0;
02680 for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
02681 linkSet.insert( *lnk );
02682 nbBndLinks += lnk->IsBoundary();
02683 }
02684 resultChains.clear();
02685 resultChains.reserve( nbBndLinks / 2 );
02686
02687 TLinkInSet linkIt, linksEnd = linkSet.end();
02688
02689
02690
02691
02692 SMDS_TypeOfPosition cornerPos = SMDS_TypeOfPosition(pos-2);
02693 const SMDS_MeshNode* corner = 0;
02694 for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt )
02695 if ( linkIt->IsBoundary() && (corner = (*linkIt)->EndPosNode(cornerPos)))
02696 break;
02697 if ( !corner)
02698 return _NO_CORNERS;
02699
02700 TLinkInSet startLink = linkIt;
02701 const SMDS_MeshNode* startCorner = corner;
02702 vector< TChain* > rowChains;
02703 int iCol = 0;
02704
02705 while ( startLink != linksEnd)
02706 {
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721 if ( resultChains.size() == nbBndLinks / 2 )
02722 return _NOT_RECT;
02723 resultChains.push_back( TChain() );
02724 TChain& columnChain = resultChains.back();
02725
02726 TLinkInSet botLink = startLink;
02727 corner = startCorner;
02728 int iRow = 0;
02729 while ( botLink != linksEnd )
02730 {
02731
02732 columnChain.push_back( *botLink );
02733
02734 const QFace* botTria = botLink->_qfaces[0];
02735 if ( !botTria )
02736 {
02737 if ( botLink == startLink )
02738 return _TWISTED_CHAIN;
02739 linkSet.erase( botLink );
02740 if ( iRow != rowChains.size() )
02741 return _FEW_ROWS;
02742 break;
02743 }
02744
02745
02746
02747
02748
02749 bool isCase2;
02750 TLinkInSet midQuadLink = linksEnd;
02751 TLinkInSet sideLink = botTria->GetBoundaryLink( linkSet, *botLink, &midQuadLink,
02752 corner, &isCase2 );
02753 if ( isCase2 ) {
02754 midQuadLink = botTria->GetLinkByNode( linkSet, *botLink, corner );
02755 if ( midQuadLink->IsBoundary() )
02756 return _BAD_MIDQUAD;
02757 }
02758 if ( sideLink == linksEnd || midQuadLink == linksEnd || sideLink == midQuadLink )
02759 return sideLink == linksEnd ? _NO_SIDELINK : _NO_MIDQUAD;
02760
02761
02762 columnChain.push_back( *midQuadLink );
02763 if ( iRow >= rowChains.size() ) {
02764 if ( iCol > 0 )
02765 return _MANY_ROWS;
02766 if ( resultChains.size() == nbBndLinks / 2 )
02767 return _NOT_RECT;
02768 resultChains.push_back( TChain() );
02769 rowChains.push_back( & resultChains.back() );
02770 }
02771 rowChains[iRow]->push_back( *sideLink );
02772 rowChains[iRow]->push_back( *midQuadLink );
02773
02774 const QFace* upTria = midQuadLink->NextFace( botTria );
02775 if ( !upTria)
02776 return _NO_UPTRIA;
02777 if ( iRow == 0 ) {
02778
02779 startCorner = startLink->NextNode( startCorner );
02780 if (isCase2)
02781 startLink = botTria->GetBoundaryLink( linkSet, *botLink, 0, startCorner );
02782 else
02783 startLink = upTria->GetBoundaryLink( linkSet, *midQuadLink, 0, startCorner );
02784
02785 if ( startLink != linksEnd ) {
02786 const SMDS_MeshNode* botNode = startLink->NextNode( startCorner );
02787 if ( (isCase2 ? botTria : upTria)->Contains( botNode ))
02788 startLink = linksEnd;
02789 else if ( startLink == botLink || startLink == midQuadLink || startLink == sideLink )
02790 return _BAD_START;
02791 }
02792 }
02793
02794 corner = sideLink->NextNode( corner );
02795
02796 linkSet.erase( botLink );
02797 botLink = upTria->GetLinkByNode( linkSet, (isCase2 ? *sideLink : *midQuadLink), corner );
02798 if ( botLink == linksEnd || botLink == midQuadLink || botLink == sideLink)
02799 return _NO_BOTLINK;
02800 if ( midQuadLink == startLink || sideLink == startLink )
02801 return _TWISTED_CHAIN;
02802 linkSet.erase( midQuadLink );
02803 linkSet.erase( sideLink );
02804
02805
02806 if ( startLink != linksEnd ) {
02807 const QFace* tria = isCase2 ? botTria : upTria;
02808 for ( int iL = 0; iL < 3; ++iL ) {
02809 linkIt = linkSet.find( tria->_sides[iL] );
02810 if ( linkIt != linksEnd )
02811 linkIt->RemoveFace( tria );
02812 }
02813 }
02814 if ( botLink->_qfaces[0] == upTria || botLink->_qfaces[1] == upTria )
02815 botLink->RemoveFace( upTria );
02816
02817 iRow++;
02818 }
02819
02820 iCol++;
02821 }
02822
02823 if ( linkSet.size() != rowChains.size() )
02824 return _BAD_SET_SIZE;
02825 for ( int iRow = 0; iRow < rowChains.size(); ++iRow ) {
02826
02827 corner = 0;
02828 for ( startLink = linkSet.begin(); startLink != linksEnd; ++startLink ) {
02829 if ( (*startLink)->node1() == startCorner ) {
02830 corner = (*startLink)->node2(); break;
02831 }
02832 else if ( (*startLink)->node2() == startCorner) {
02833 corner = (*startLink)->node1(); break;
02834 }
02835 }
02836 if ( startLink == linksEnd )
02837 return _BAD_CORNER;
02838 rowChains[ iRow ]->push_back( *startLink );
02839 linkSet.erase( startLink );
02840 startCorner = corner;
02841 }
02842
02843 return _OK;
02844 }
02845 }
02846
02847
02854
02855
02856 void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
02857 {
02858
02859 if ( getenv("NO_FixQuadraticElements") )
02860 return;
02861
02862
02863
02864 if ( myShape.IsNull() ) {
02865 if ( !myMesh->HasShapeToMesh() ) return;
02866 SetSubShape( myMesh->GetShapeToMesh() );
02867
02868 #ifdef _DEBUG_
02869 int nbSolids = 0;
02870 TopTools_IndexedMapOfShape solids;
02871 TopExp::MapShapes(myShape,TopAbs_SOLID,solids);
02872 nbSolids = solids.Extent();
02873 #endif
02874 TopTools_MapOfShape faces;
02875 for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) {
02876 faces.Add( f.Current() );
02877 }
02878 for ( TopExp_Explorer s(myShape,TopAbs_SOLID); s.More(); s.Next() ) {
02879 if ( myMesh->GetSubMesh( s.Current() )->IsEmpty() ) {
02880 for ( TopExp_Explorer f( s.Current(), TopAbs_FACE); f.More(); f.Next() )
02881 faces.Add( f.Current() );
02882 }
02883 else {
02884 #ifdef _DEBUG_
02885 MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
02886 #endif
02887 SMESH_MesherHelper h(*myMesh);
02888 h.SetSubShape( s.Current() );
02889 h.FixQuadraticElements(false);
02890 }
02891 }
02892
02893 #ifdef _DEBUG_
02894 int nbfaces = faces.Extent(); nbfaces++, nbfaces--;
02895 #endif
02896 for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
02897 MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
02898 SMESH_MesherHelper h(*myMesh);
02899 h.SetSubShape( fIt.Key() );
02900 h.FixQuadraticElements(true);
02901 }
02902
02903 return;
02904 }
02905
02906
02907
02908
02909 SMDS_ElemIteratorPtr elemIt;
02910 SMDSAbs_ElementType elemType = SMDSAbs_All;
02911
02912 SMESH_subMesh* submesh = myMesh->GetSubMeshContaining( myShapeID );
02913 if ( !submesh )
02914 return;
02915 if ( SMESHDS_SubMesh* smDS = submesh->GetSubMeshDS() ) {
02916 elemIt = smDS->GetElements();
02917 if ( elemIt->more() ) {
02918 elemType = elemIt->next()->GetType();
02919 elemIt = smDS->GetElements();
02920 }
02921 }
02922 if ( !elemIt || !elemIt->more() || elemType < SMDSAbs_Face )
02923 return;
02924
02925
02926
02927
02928 set< QLink > links;
02929 set< QFace > faces;
02930 set< QLink >::iterator pLink;
02931 set< QFace >::iterator pFace;
02932
02933 bool isCurved = false;
02934
02935
02936 SMDS_VolumeTool volTool;
02937
02938 TIDSortedNodeSet apexOfPyramid;
02939 const int apexIndex = 4;
02940
02941 if ( elemType == SMDSAbs_Volume )
02942 {
02943 while ( elemIt->more() )
02944 {
02945 const SMDS_MeshElement* vol = elemIt->next();
02946 if ( !vol->IsQuadratic() || !volTool.Set( vol ))
02947 return;
02948 for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
02949 {
02950 int nbN = volTool.NbFaceNodes( iF );
02951
02952 const SMDS_MeshNode** faceNodes = volTool.GetFaceNodes( iF );
02953 vector< const QLink* > faceLinks( nbN/2 );
02954 for ( int iN = 0; iN < nbN; iN += 2 )
02955 {
02956
02957 QLink link( faceNodes[iN], faceNodes[iN+2], faceNodes[iN+1] );
02958 pLink = links.insert( link ).first;
02959 faceLinks[ iN/2 ] = & *pLink;
02960 if ( !isCurved )
02961 isCurved = !link.IsStraight();
02962 if ( link.MediumPos() == SMDS_TOP_3DSPACE && !link.IsStraight() )
02963 return;
02964 }
02965
02966 pFace = faces.insert( QFace( faceLinks )).first;
02967 if ( pFace->NbVolumes() == 0 )
02968 pFace->AddSelfToLinks();
02969 pFace->SetVolume( vol );
02970
02971
02972
02973 #ifdef _DEBUG_
02974 if ( nbN == 6 )
02975 pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
02976 else
02977 pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],
02978 faceNodes[4],faceNodes[6] );
02979 #endif
02980 }
02981
02982 if ( vol->NbCornerNodes() == 5 )
02983 apexOfPyramid.insert( vol->GetNode( apexIndex ));
02984 }
02985 set< QLink >::iterator pLink = links.begin();
02986 for ( ; pLink != links.end(); ++pLink )
02987 pLink->SetContinuesFaces();
02988 }
02989 else
02990 {
02991 while ( elemIt->more() )
02992 {
02993 const SMDS_MeshElement* face = elemIt->next();
02994 if ( !face->IsQuadratic() )
02995 continue;
02996
02997 int nbN = face->NbNodes()/2;
02998 vector< const QLink* > faceLinks( nbN );
02999 for ( int iN = 0; iN < nbN; ++iN )
03000 {
03001
03002 QLink link( face->GetNode(iN), face->GetNode((iN+1)%nbN), face->GetNode(iN+nbN) );
03003 pLink = links.insert( link ).first;
03004 faceLinks[ iN ] = & *pLink;
03005 if ( !isCurved )
03006 isCurved = !link.IsStraight();
03007 }
03008
03009 pFace = faces.insert( QFace( faceLinks )).first;
03010 pFace->AddSelfToLinks();
03011
03012 }
03013 }
03014 if ( !isCurved )
03015 return;
03016
03017
03018
03019
03020
03021 TopLoc_Location loc;
03022
03023 int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0;
03024 for ( ; isInside < 2; ++isInside ) {
03025 MSG( "--------------- LOOP (inside=" << isInside << ") ------------------");
03026 SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE;
03027 SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE;
03028
03029 for ( pFace = faces.begin(); pFace != faces.end(); ++pFace ) {
03030 if ( bool(isInside) == pFace->IsBoundary() )
03031 continue;
03032 for ( int dir = 0; dir < 2; ++dir )
03033 {
03034 MSG( "CHAIN");
03035
03036 int error = ERR_OK;
03037 TChain rawChain;
03038 if ( !pFace->GetLinkChain( dir, rawChain, pos, error) && error ==ERR_UNKNOWN ) continue;
03039 rawChain.reverse();
03040 if ( !pFace->GetLinkChain( dir+2, rawChain, pos, error ) && error ==ERR_UNKNOWN ) continue;
03041
03042 vector< TChain > chains;
03043 if ( error == ERR_OK ) {
03044 chains.resize(1);
03045 chains[0].splice( chains[0].begin(), rawChain );
03046 }
03047 else if ( error == ERR_TRI ) {
03048 TSplitTriaResult res = splitTrianglesIntoChains( rawChain, chains, pos );
03049 if ( res != _OK ) {
03050 fixTriaNearBoundary( rawChain, *this );
03051 break;
03052 }
03053 }
03054 else if ( error == ERR_PRISM ) {
03055 fixPrism( rawChain );
03056 break;
03057 }
03058 else {
03059 continue;
03060 }
03061 for ( int iC = 0; iC < chains.size(); ++iC )
03062 {
03063 TChain& chain = chains[iC];
03064 if ( chain.empty() ) continue;
03065 if ( chain.front().IsStraight() && chain.back().IsStraight() ) {
03066 MSG("3D straight - ignore");
03067 continue;
03068 }
03069 if ( chain.front()->MediumPos() > bndPos ||
03070 chain.back() ->MediumPos() > bndPos ) {
03071 MSG("Internal chain - ignore");
03072 continue;
03073 }
03074
03075 double chainLen = 0;
03076 vector< double > linkPos;
03077 MSGBEG( "Link medium nodes: ");
03078 TChain::iterator link0 = chain.begin(), link1 = chain.begin(), link2;
03079 for ( ++link1; link1 != chain.end(); ++link1, ++link0 ) {
03080 MSGBEG( (*link0)->_mediumNode->GetID() << "-" <<(*link1)->_mediumNode->GetID()<<" ");
03081 double len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
03082 while ( len < numeric_limits<double>::min() ) {
03083 link1 = chain.erase( link1 );
03084 if ( link1 == chain.end() )
03085 break;
03086 len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
03087 }
03088 chainLen += len;
03089 linkPos.push_back( chainLen );
03090 }
03091 MSG("");
03092 if ( linkPos.size() < 2 )
03093 continue;
03094
03095 gp_Vec move0 = chain.front()->_nodeMove;
03096 gp_Vec move1 = chain.back ()->_nodeMove;
03097
03098 TopoDS_Face face;
03099 bool checkUV = true;
03100 if ( !isInside )
03101 {
03102
03103 TChainLink& linkOnFace = *(++chain.begin());
03104 const SMDS_MeshNode* nodeOnFace = linkOnFace->_mediumNode;
03105 TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() );
03106 if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE )
03107 {
03108 face = TopoDS::Face( f );
03109 Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
03110 bool isStraight[2];
03111 for ( int is1 = 0; is1 < 2; ++is1 )
03112 {
03113 TChainLink& link = is1 ? chain.back() : chain.front();
03114 gp_XY uvm = GetNodeUV( face, link->_mediumNode, nodeOnFace, &checkUV);
03115 gp_XY uv1 = GetNodeUV( face, link->node1(), nodeOnFace, &checkUV);
03116 gp_XY uv2 = GetNodeUV( face, link->node2(), nodeOnFace, &checkUV);
03117 gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
03118
03119 gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, false);
03120 ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 );
03121 if ( !is1 )
03122 nodeOnFace = (*(++chain.rbegin()))->_mediumNode;
03123 isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(),
03124 10 * uvMove.SquareModulus());
03125 }
03126 if ( isStraight[0] && isStraight[1] ) {
03127 MSG("2D straight - ignore");
03128 continue;
03129 }
03130
03131
03132 gp_XY uvm = GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV);
03133 gp_XY uv1 = GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV);
03134 gp_XY uv2 = GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV);
03135 gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
03136 if (( uvm - uv12 ).SquareModulus() > 1e-10 )
03137 {
03138 MSG("Already fixed - ignore");
03139 continue;
03140 }
03141 }
03142 }
03143 gp_Trsf trsf;
03144 if ( isInside || face.IsNull() )
03145 {
03146
03147 {
03148 TChainLink& ln0 = chain.front(), ln1 = *(++chain.begin());
03149 trsf.SetTransformation( gp_Ax3( gp::Origin(), ln0.Normal(),
03150 gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
03151 move0.Transform(trsf);
03152 }
03153 {
03154 TChainLink& ln0 = *(++chain.rbegin()), ln1 = chain.back();
03155 trsf.SetTransformation( gp_Ax3( gp::Origin(), ln1.Normal(),
03156 gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
03157 move1.Transform(trsf);
03158 }
03159 }
03160
03161 link2 = chain.begin();
03162 link0 = link2++;
03163 link1 = link2++;
03164 for ( int i = 0; link2 != chain.end(); ++link0, ++link1, ++link2, ++i )
03165 {
03166 double r = linkPos[i] / chainLen;
03167
03168 gp_Vec move = (1. - r) * move0 + r * move1;
03169 if ( isInside || face.IsNull()) {
03170
03171 gp_Vec x01( (*link0)->MiddlePnt(), (*link1)->MiddlePnt() );
03172 gp_Vec x12( (*link1)->MiddlePnt(), (*link2)->MiddlePnt() );
03173 gp_Vec x = x01.Normalized() + x12.Normalized();
03174 trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() );
03175 move.Transform(trsf);
03176 }
03177 else {
03178
03179 Handle(Geom_Surface) s = BRep_Tool::Surface(face,loc);
03180 gp_XY oldUV = GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV);
03181 gp_XY newUV = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added);
03182 gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
03183 move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
03184 #ifdef _DEBUG_
03185 if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
03186 move.SquareMagnitude())
03187 {
03188 gp_XY uv0 = GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV);
03189 gp_XY uv2 = GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV);
03190 MSG( "TOO LONG MOVE \t" <<
03191 "uv0: "<<uv0.X()<<", "<<uv0.Y()<<" \t" <<
03192 "uv2: "<<uv2.X()<<", "<<uv2.Y()<<" \t" <<
03193 "uvOld: "<<oldUV.X()<<", "<<oldUV.Y()<<" \t" <<
03194 "newUV: "<<newUV.X()<<", "<<newUV.Y()<<" \t");
03195 }
03196 #endif
03197 }
03198 (*link1)->Move( move );
03199 MSG( "Move " << (*link1)->_mediumNode->GetID() << " following "
03200 << chain.front()->_mediumNode->GetID() <<"-"
03201 << chain.back ()->_mediumNode->GetID() <<
03202 " by " << move.Magnitude());
03203 }
03204 }
03205 }
03206 }
03207 }
03208
03209
03210
03211
03212
03213
03214
03215
03216 for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {
03217 if ( pLink->IsMoved() ) {
03218 gp_Pnt p = pLink->MiddlePnt() + pLink->Move();
03219 GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z());
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249 }
03250 }
03251
03252
03253
03254
03255
03256
03257 TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin();
03258 for ( ; apexIt != apexOfPyramid.end(); ++apexIt )
03259 {
03260 SMESH_TNodeXYZ apex = *apexIt;
03261
03262 gp_Vec maxMove( 0,0,0 );
03263 double maxMoveSize2 = 0;
03264
03265
03266 const int base2MediumShift = 5;
03267
03268
03269 SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume );
03270 vector< const SMDS_MeshElement* > pyramids;
03271 while ( volIt->more() )
03272 {
03273 const SMDS_MeshElement* pyram = volIt->next();
03274 if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue;
03275 pyramids.push_back( pyram );
03276
03277 for ( int iBase = 0; iBase < apexIndex; ++iBase )
03278 {
03279 SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift );
03280 if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
03281 {
03282 SMESH_TNodeXYZ n1 = pyram->GetNode( iBase );
03283 SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 );
03284 gp_Pnt middle = 0.5 * ( n1 + n2 );
03285 gp_Vec move( middle, medium );
03286 double moveSize2 = move.SquareMagnitude();
03287 if ( moveSize2 > maxMoveSize2 )
03288 maxMove = move, maxMoveSize2 = moveSize2;
03289 }
03290 }
03291 }
03292
03293
03294 if ( maxMoveSize2 > 1e-20 )
03295 {
03296 apex += maxMove.XYZ();
03297 GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z());
03298
03299
03300 const int base2MediumShift_2 = 9;
03301 for ( unsigned i = 0; i < pyramids.size(); ++i )
03302 for ( int iBase = 0; iBase < apexIndex; ++iBase )
03303 {
03304 SMESH_TNodeXYZ base = pyramids[i]->GetNode( iBase );
03305 const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 );
03306 gp_XYZ middle = 0.5 * ( apex + base );
03307 GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z());
03308 }
03309 }
03310 }
03311 }
03312