+ multimap< double, ForcedPoint >::iterator d2uv = sortedFP.begin();
+ for ( ; d2uv != sortedFP.end(); ++d2uv )
+ myForcedPnts.push_back( (*d2uv).second );
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Splits quads by adding points of enforced nodes and create nodes on
+ * the sides shared by quads
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::addEnforcedNodes()
+{
+ // if ( myForcedPnts.empty() )
+ // return true;
+
+ // make a map of quads sharing a side
+ map< StdMeshers_FaceSidePtr, vector< FaceQuadStruct::Ptr > > quadsBySide;
+ list< FaceQuadStruct::Ptr >::iterator quadIt = myQuadList.begin();
+ for ( ; quadIt != myQuadList.end(); ++quadIt )
+ for ( size_t iSide = 0; iSide < (*quadIt)->side.size(); ++iSide )
+ {
+ if ( !setNormalizedGrid( *quadIt ))
+ return false;
+ quadsBySide[ (*quadIt)->side[iSide] ].push_back( *quadIt );
+ }
+
+ SMESH_Mesh* mesh = myHelper->GetMesh();
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ const TopoDS_Face& face = TopoDS::Face( myHelper->GetSubShape() );
+ Handle(Geom_Surface) surf = BRep_Tool::Surface( face );
+
+ for ( size_t iFP = 0; iFP < myForcedPnts.size(); ++iFP )
+ {
+ bool isNodeEnforced = false;
+
+ // look for a quad enclosing a enforced point
+ for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt )
+ {
+ FaceQuadStruct::Ptr quad = *quadIt;
+ if ( !setNormalizedGrid( *quadIt ))
+ return false;
+ int i,j;
+ if ( !quad->findCell( myForcedPnts[ iFP ], i, j ))
+ continue;
+
+ // a grid cell is found, select a node of the cell to move
+ // to the enforced point to and to split the quad at
+ multimap< double, pair< int, int > > ijByDist;
+ for ( int di = 0; di < 2; ++di )
+ for ( int dj = 0; dj < 2; ++dj )
+ {
+ double dist2 = ( myForcedPnts[ iFP ].uv - quad->UVPt( i+di,j+dj ).UV() ).SquareModulus();
+ ijByDist.insert( make_pair( dist2, make_pair( di,dj )));
+ }
+ // try all nodes starting from the closest one
+ set< FaceQuadStruct::Ptr > changedQuads;
+ multimap< double, pair< int, int > >::iterator d2ij = ijByDist.begin();
+ for ( ; !isNodeEnforced && d2ij != ijByDist.end(); ++d2ij )
+ {
+ int di = d2ij->second.first;
+ int dj = d2ij->second.second;
+
+ // check if a node is at a side
+ int iSide = -1;
+ if ( dj== 0 && j == 0 )
+ iSide = QUAD_BOTTOM_SIDE;
+ else if ( dj == 1 && j+2 == quad->jSize )
+ iSide = QUAD_TOP_SIDE;
+ else if ( di == 0 && i == 0 )
+ iSide = QUAD_LEFT_SIDE;
+ else if ( di == 1 && i+2 == quad->iSize )
+ iSide = QUAD_RIGHT_SIDE;
+
+ if ( iSide > -1 ) // ----- node is at a side
+ {
+ FaceQuadStruct::Side& side = quad->side[ iSide ];
+ // check if this node can be moved
+ if ( quadsBySide[ side ].size() < 2 )
+ continue; // its a face boundary -> can't move the node
+
+ int quadNodeIndex = ( iSide % 2 ) ? j : i;
+ int sideNodeIndex = side.ToSideIndex( quadNodeIndex );
+ if ( side.IsForced( sideNodeIndex ))
+ {
+ // the node is already moved to another enforced point
+ isNodeEnforced = quad->isEqual( myForcedPnts[ iFP ], i, j );
+ continue;
+ }
+ // make a node of a side forced
+ vector<UVPtStruct>& points = (vector<UVPtStruct>&) side.GetUVPtStruct();
+ points[ sideNodeIndex ].u = myForcedPnts[ iFP ].U();
+ points[ sideNodeIndex ].v = myForcedPnts[ iFP ].V();
+
+ updateSideUV( side, sideNodeIndex, quadsBySide );
+
+ // update adjacent sides
+ set< StdMeshers_FaceSidePtr > updatedSides;
+ updatedSides.insert( side );
+ for ( size_t i = 0; i < side.contacts.size(); ++i )
+ if ( side.contacts[i].point == sideNodeIndex )
+ {
+ const vector< FaceQuadStruct::Ptr >& adjQuads =
+ quadsBySide[ *side.contacts[i].other_side ];
+ if ( adjQuads.size() > 1 &&
+ updatedSides.insert( * side.contacts[i].other_side ).second )
+ {
+ updateSideUV( *side.contacts[i].other_side,
+ side.contacts[i].other_point,
+ quadsBySide );
+ }
+ changedQuads.insert( adjQuads.begin(), adjQuads.end() );
+ }
+ const vector< FaceQuadStruct::Ptr >& adjQuads = quadsBySide[ side ];
+ changedQuads.insert( adjQuads.begin(), adjQuads.end() );
+
+ isNodeEnforced = true;
+ }
+ else // ------------------ node is inside the quad
+ {
+ i += di;
+ j += dj;
+ // make a new side passing through IJ node and split the quad
+ int indForced, iNewSide;
+ if ( quad->iSize < quad->jSize ) // split vertically
+ {
+ quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/true );
+ indForced = j;
+ iNewSide = splitQuad( quad, i, 0 );
+ }
+ else
+ {
+ quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/false );
+ indForced = i;
+ iNewSide = splitQuad( quad, 0, j );
+ }
+ FaceQuadStruct::Ptr newQuad = myQuadList.back();
+ FaceQuadStruct::Side& newSide = newQuad->side[ iNewSide ];
+
+ newSide.forced_nodes.insert( indForced );
+ quad->side[( iNewSide+2 ) % 4 ].forced_nodes.insert( indForced );
+
+ quadsBySide[ newSide ].push_back( quad );
+ quadsBySide[ newQuad->side[0] ].push_back( newQuad );
+ quadsBySide[ newQuad->side[1] ].push_back( newQuad );
+ quadsBySide[ newQuad->side[2] ].push_back( newQuad );
+ quadsBySide[ newQuad->side[3] ].push_back( newQuad );
+
+ isNodeEnforced = true;
+
+ } // end of "node is inside the quad"
+
+ } // loop on nodes of the cell
+
+ // remove out-of-date uv grid of changedQuads
+ set< FaceQuadStruct::Ptr >::iterator qIt = changedQuads.begin();
+ for ( ; qIt != changedQuads.end(); ++qIt )
+ (*qIt)->uv_grid.clear();
+
+ if ( isNodeEnforced )
+ break;
+
+ } // loop on quads
+
+ if ( !isNodeEnforced )
+ {
+ if ( !myForcedPnts[ iFP ].vertex.IsNull() )
+ return error(TComm("Unable to move any node to vertex #")
+ <<myHelper->GetMeshDS()->ShapeToIndex( myForcedPnts[ iFP ].vertex ));
+ else
+ return error(TComm("Unable to move any node to point ( ")
+ << myForcedPnts[iFP].xyz.X() << ", "
+ << myForcedPnts[iFP].xyz.Y() << ", "
+ << myForcedPnts[iFP].xyz.Z() << " )");
+ }
+
+ } // loop on enforced points
+
+ // Compute nodes on all sides, where not yet present
+
+ for ( quadIt = myQuadList.begin(); quadIt != myQuadList.end(); ++quadIt )
+ {
+ FaceQuadStruct::Ptr quad = *quadIt;
+ for ( int iSide = 0; iSide < 4; ++iSide )
+ {
+ FaceQuadStruct::Side & side = quad->side[ iSide ];
+ if ( side.nbNodeOut > 0 )
+ continue; // emulated side
+ vector< FaceQuadStruct::Ptr >& quadVec = quadsBySide[ side ];
+ if ( quadVec.size() <= 1 )
+ continue; // outer side
+
+ bool missedNodesOnSide = false;
+ const vector<UVPtStruct>& points = side.grid->GetUVPtStruct();
+ for ( size_t iC = 0; iC < side.contacts.size(); ++iC )
+ {
+ const vector<UVPtStruct>& oGrid = side.contacts[iC].other_side->grid->GetUVPtStruct();
+ const UVPtStruct& uvPt = points[ side.contacts[iC].point ];
+ if ( side.contacts[iC].other_point >= oGrid.size() ||
+ side.contacts[iC].point >= points.size() )
+ throw SALOME_Exception( "StdMeshers_Quadrangle_2D::addEnforcedNodes(): wrong contact" );
+ if ( oGrid[ side.contacts[iC].other_point ].node )
+ (( UVPtStruct& ) uvPt).node = oGrid[ side.contacts[iC].other_point ].node;
+ }
+ for ( size_t iP = 0; iP < points.size(); ++iP )
+ if ( !points[ iP ].node )
+ {
+ UVPtStruct& uvPnt = ( UVPtStruct& ) points[ iP ];
+ gp_Pnt P = surf->Value( uvPnt.u, uvPnt.v );
+ uvPnt.node = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace( uvPnt.node, myHelper->GetSubShapeID(), uvPnt.u, uvPnt.v );
+ missedNodesOnSide = true;
+ }
+ if ( missedNodesOnSide )
+ {
+ // clear uv_grid where nodes are missing
+ for ( size_t iQ = 0; iQ < quadVec.size(); ++iQ )
+ quadVec[ iQ ]->uv_grid.clear();
+ }
+ }
+ }
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Splits a quad at I or J. Returns an index of a new side in the new quad
+ */
+//================================================================================
+
+int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J)
+{
+ FaceQuadStruct* newQuad = new FaceQuadStruct( quad->face );
+ myQuadList.push_back( FaceQuadStruct::Ptr( newQuad ));
+
+ vector<UVPtStruct> points;
+ if ( I > 0 && I <= quad->iSize-2 )
+ {
+ points.reserve( quad->jSize );
+ for ( int jP = 0; jP < quad->jSize; ++jP )
+ points.push_back( quad->UVPt( I, jP ));
+
+ newQuad->side.resize( 4 );
+ newQuad->side[ QUAD_BOTTOM_SIDE ] = quad->side[ QUAD_BOTTOM_SIDE ];
+ newQuad->side[ QUAD_RIGHT_SIDE ] = quad->side[ QUAD_RIGHT_SIDE ];
+ newQuad->side[ QUAD_TOP_SIDE ] = quad->side[ QUAD_TOP_SIDE ];
+ newQuad->side[ QUAD_LEFT_SIDE ] = StdMeshers_FaceSide::New( points, quad->face );
+
+ FaceQuadStruct::Side& newSide = newQuad->side[ QUAD_LEFT_SIDE ];
+ FaceQuadStruct::Side& newSide2 = quad->side [ QUAD_RIGHT_SIDE ];
+
+ quad->side[ QUAD_RIGHT_SIDE ] = newSide;
+
+ int iBot = quad->side[ QUAD_BOTTOM_SIDE ].ToSideIndex( I );
+ int iTop = quad->side[ QUAD_TOP_SIDE ].ToSideIndex( I );
+
+ newSide.AddContact ( 0, & quad->side[ QUAD_BOTTOM_SIDE ], iBot );
+ newSide2.AddContact( 0, & quad->side[ QUAD_BOTTOM_SIDE ], iBot );
+ newSide.AddContact ( quad->jSize - 1, & quad->side[ QUAD_TOP_SIDE ], iTop );
+ newSide2.AddContact( quad->jSize - 1, & quad->side[ QUAD_TOP_SIDE ], iTop );
+ // cout << "Contact: L " << &newSide << " "<< newSide.NbPoints()
+ // << " R " << &newSide2 << " "<< newSide2.NbPoints()
+ // << " B " << &quad->side[ QUAD_BOTTOM_SIDE ] << " "<< quad->side[ QUAD_BOTTOM_SIDE].NbPoints()
+ // << " T " << &quad->side[ QUAD_TOP_SIDE ] << " "<< quad->side[ QUAD_TOP_SIDE].NbPoints()<< endl;
+
+ newQuad->side[ QUAD_BOTTOM_SIDE ].from = iBot;
+ newQuad->side[ QUAD_TOP_SIDE ].from = iTop;
+ newQuad->name = ( TComm("Right of I=") << I );
+
+ quad->side[ QUAD_BOTTOM_SIDE ].to = iBot + 1;
+ quad->side[ QUAD_TOP_SIDE ].to = iTop + 1;
+ quad->uv_grid.clear();
+
+ return QUAD_LEFT_SIDE;
+ }
+ else if ( J > 0 && J <= quad->jSize-2 ) //// split horizontally, a new quad is below an old one
+ {
+ points.reserve( quad->iSize );
+ for ( int iP = 0; iP < quad->iSize; ++iP )
+ points.push_back( quad->UVPt( iP, J ));
+
+ newQuad->side.resize( 4 );
+ newQuad->side[ QUAD_BOTTOM_SIDE ] = quad->side[ QUAD_BOTTOM_SIDE ];
+ newQuad->side[ QUAD_RIGHT_SIDE ] = quad->side[ QUAD_RIGHT_SIDE ];
+ newQuad->side[ QUAD_TOP_SIDE ] = StdMeshers_FaceSide::New( points, quad->face );
+ newQuad->side[ QUAD_LEFT_SIDE ] = quad->side[ QUAD_LEFT_SIDE ];
+
+ FaceQuadStruct::Side& newSide = newQuad->side[ QUAD_TOP_SIDE ];
+ FaceQuadStruct::Side& newSide2 = quad->side [ QUAD_BOTTOM_SIDE ];
+
+ quad->side[ QUAD_BOTTOM_SIDE ] = newSide;
+
+ int iLft = quad->side[ QUAD_LEFT_SIDE ].ToSideIndex( J );
+ int iRgt = quad->side[ QUAD_RIGHT_SIDE ].ToSideIndex( J );
+
+ newSide.AddContact ( 0, & quad->side[ QUAD_LEFT_SIDE ], iLft );
+ newSide2.AddContact( 0, & quad->side[ QUAD_LEFT_SIDE ], iLft );
+ newSide.AddContact ( quad->iSize - 1, & quad->side[ QUAD_RIGHT_SIDE ], iRgt );
+ newSide2.AddContact( quad->iSize - 1, & quad->side[ QUAD_RIGHT_SIDE ], iRgt );
+ // cout << "Contact: T " << &newSide << " "<< newSide.NbPoints()
+ // << " B " << &newSide2 << " "<< newSide2.NbPoints()
+ // << " L " << &quad->side[ QUAD_LEFT_SIDE ] << " "<< quad->side[ QUAD_LEFT_SIDE].NbPoints()
+ // << " R " << &quad->side[ QUAD_RIGHT_SIDE ] << " "<< quad->side[ QUAD_RIGHT_SIDE].NbPoints()<< endl;
+
+ newQuad->side[ QUAD_RIGHT_SIDE ].to = iRgt+1;
+ newQuad->side[ QUAD_LEFT_SIDE ].to = iLft+1;
+ newQuad->name = ( TComm("Below J=") << J );
+
+ quad->side[ QUAD_RIGHT_SIDE ].from = iRgt;
+ quad->side[ QUAD_LEFT_SIDE ].from = iLft;
+ quad->uv_grid.clear();
+
+ return QUAD_TOP_SIDE;
+ }
+
+ myQuadList.pop_back();
+ return -1;
+}
+
+//================================================================================
+/*!
+ * \brief Updates UV of a side after moving its node
+ */
+//================================================================================
+
+void StdMeshers_Quadrangle_2D::updateSideUV( FaceQuadStruct::Side& side,
+ int iForced,
+ const TQuadsBySide& quadsBySide,
+ int * iNext)
+{
+ if ( !iNext )
+ {
+ side.forced_nodes.insert( iForced );
+
+ // update parts of the side before and after iForced
+
+ set<int>::iterator iIt = side.forced_nodes.upper_bound( iForced );
+ int iEnd = Min( side.NbPoints()-1, ( iIt == side.forced_nodes.end() ) ? int(1e7) : *iIt );
+ if ( iForced + 1 < iEnd )
+ updateSideUV( side, iForced, quadsBySide, &iEnd );
+
+ iIt = side.forced_nodes.lower_bound( iForced );
+ int iBeg = Max( 0, ( iIt == side.forced_nodes.begin() ) ? 0 : *--iIt );
+ if ( iForced - 1 > iBeg )
+ updateSideUV( side, iForced, quadsBySide, &iBeg );
+
+ return;
+ }
+
+ const int iFrom = Min ( iForced, *iNext );
+ const int iTo = Max ( iForced, *iNext ) + 1;
+ const int sideSize = iTo - iFrom;
+
+ vector<UVPtStruct> points[4]; // side points of a temporary quad
+
+ // from the quads get grid points adjacent to the side
+ // to make two sides of a temporary quad
+ vector< FaceQuadStruct::Ptr > quads = quadsBySide.find( side )->second; // copy!
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ points[ is2nd ].reserve( sideSize );
+ int nbLoops = 0;
+ while ( points[is2nd].size() < sideSize )
+ {
+ int iCur = iFrom + points[is2nd].size() - int( !points[is2nd].empty() );
+
+ // look for a quad adjacent to iCur-th point of the side
+ for ( size_t iQ = 0; iQ < quads.size(); ++iQ )
+ {
+ FaceQuadStruct::Ptr q = quads[ iQ ];
+ if ( !q )
+ continue;
+ size_t iS;
+ for ( iS = 0; iS < q->side.size(); ++iS )
+ if ( side.grid == q->side[ iS ].grid )
+ break;
+ bool isOut;
+ if ( !q->side[ iS ].IsReversed() )
+ isOut = ( q->side[ iS ].from > iCur || q->side[ iS ].to-1 <= iCur );
+ else
+ isOut = ( q->side[ iS ].to >= iCur || q->side[ iS ].from <= iCur );
+ if ( isOut )
+ continue;
+ if ( !setNormalizedGrid( q ))
+ continue;
+
+ // found - copy points
+ int i,j,di,dj,nb;
+ if ( iS % 2 ) // right or left
+ {
+ i = ( iS == QUAD_LEFT_SIDE ) ? 1 : q->iSize-2;
+ j = q->side[ iS ].ToQuadIndex( iCur );
+ di = 0;
+ dj = ( q->side[ iS ].IsReversed() ) ? -1 : +1;
+ nb = ( q->side[ iS ].IsReversed() ) ? j+1 : q->jSize-j;
+ }
+ else // bottom or top
+ {
+ i = q->side[ iS ].ToQuadIndex( iCur );
+ j = ( iS == QUAD_BOTTOM_SIDE ) ? 1 : q->jSize-2;
+ di = ( q->side[ iS ].IsReversed() ) ? -1 : +1;
+ dj = 0;
+ nb = ( q->side[ iS ].IsReversed() ) ? i+1 : q->iSize-i;
+ }
+ if ( !points[is2nd].empty() )
+ {
+ gp_UV lastUV = points[is2nd].back().UV();
+ gp_UV quadUV = q->UVPt( i, j ).UV();
+ if ( ( lastUV - quadUV ).SquareModulus() > 1e-10 )
+ continue; // quad is on the other side of the side
+ i += di; j += dj; --nb;
+ }
+ for ( ; nb > 0 ; --nb )
+ {
+ points[ is2nd ].push_back( q->UVPt( i, j ));
+ if ( points[is2nd].size() >= sideSize )
+ break;
+ i += di; j += dj;
+ }
+ quads[ iQ ].reset(); // not to use this quad anymore
+
+ if ( points[is2nd].size() >= sideSize )
+ break;
+ } // loop on quads
+
+ if ( nbLoops++ > quads.size() )
+ throw SALOME_Exception( "StdMeshers_Quadrangle_2D::updateSideUV() bug: infinite loop" );
+
+ } // while ( points[is2nd].size() < sideSize )
+ } // two loops to fill points[0] and points[1]
+
+ // points for other pair of opposite sides of the temporary quad
+
+ enum { L,R,B,T }; // side index of points[]
+
+ points[B].push_back( points[L].front() );
+ points[B].push_back( side.GetUVPtStruct()[ iFrom ]);
+ points[B].push_back( points[R].front() );
+
+ points[T].push_back( points[L].back() );
+ points[T].push_back( side.GetUVPtStruct()[ iTo-1 ]);
+ points[T].push_back( points[R].back() );
+
+ // make the temporary quad
+ FaceQuadStruct::Ptr tmpQuad
+ ( new FaceQuadStruct( TopoDS::Face( myHelper->GetSubShape() ), "tmpQuad"));
+ tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[B] )); // bottom
+ tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[R] )); // right
+ tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[T] ));
+ tmpQuad->side.push_back( StdMeshers_FaceSide::New( points[L] ));
+
+ // compute new UV of the side
+ setNormalizedGrid( tmpQuad );
+ gp_UV uv = tmpQuad->UVPt(1,0).UV();
+ tmpQuad->updateUV( uv, 1,0, /*isVertical=*/true );
+
+ // update UV of the side
+ vector<UVPtStruct>& sidePoints = (vector<UVPtStruct>&) side.GetUVPtStruct();
+ for ( int i = iFrom; i < iTo; ++i )
+ {
+ const uvPtStruct& uvPt = tmpQuad->UVPt( 1, i-iFrom );
+ sidePoints[ i ].u = uvPt.u;
+ sidePoints[ i ].v = uvPt.v;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Finds indices of a grid quad enclosing the given enforced UV
+ */
+//================================================================================
+
+bool FaceQuadStruct::findCell( const gp_XY& UV, int & I, int & J )
+{
+ // setNormalizedGrid() must be called before!
+ if ( uv_box.IsOut( UV ))
+ return false;
+
+ // find an approximate position
+ double x = 0.5, y = 0.5;
+ gp_XY t0 = UVPt( iSize - 1, 0 ).UV();
+ gp_XY t1 = UVPt( 0, jSize - 1 ).UV();
+ gp_XY t2 = UVPt( 0, 0 ).UV();
+ SMESH_MeshAlgos::GetBarycentricCoords( UV, t0, t1, t2, x, y );
+ x = Min( 1., Max( 0., x ));
+ y = Min( 1., Max( 0., y ));
+
+ // precise the position
+ normPa2IJ( x,y, I,J );
+ if ( !isNear( UV, I,J ))
+ {
+ // look for the most close IJ by traversing uv_grid in the middle
+ double dist2, minDist2 = ( UV - UVPt( I,J ).UV() ).SquareModulus();
+ for ( int isU = 0; isU < 2; ++isU )
+ {
+ int ind1 = isU ? 0 : iSize / 2;
+ int ind2 = isU ? jSize / 2 : 0;
+ int di1 = isU ? Max( 2, iSize / 20 ) : 0;
+ int di2 = isU ? 0 : Max( 2, jSize / 20 );
+ int i,nb = isU ? iSize / di1 : jSize / di2;
+ for ( i = 0; i < nb; ++i, ind1 += di1, ind2 += di2 )
+ if (( dist2 = ( UV - UVPt( ind1,ind2 ).UV() ).SquareModulus() ) < minDist2 )
+ {
+ I = ind1;
+ J = ind2;
+ if ( isNear( UV, I,J ))
+ return true;
+ minDist2 = ( UV - UVPt( I,J ).UV() ).SquareModulus();
+ }
+ }
+ if ( !isNear( UV, I,J, Max( iSize, jSize ) /2 ))
+ return false;
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Find indices (i,j) of a point in uv_grid by normalized parameters (x,y)
+ */
+//================================================================================
+
+void FaceQuadStruct::normPa2IJ(double X, double Y, int & I, int & J )
+{
+
+ I = Min( int ( iSize * X ), iSize - 2 );
+ J = Min( int ( jSize * Y ), jSize - 2 );
+
+ int oldI, oldJ;
+ do
+ {
+ oldI = I, oldJ = J;
+ while ( X <= UVPt( I,J ).x && I != 0 )
+ --I;
+ while ( X > UVPt( I+1,J ).x && I+2 < iSize )
+ ++I;
+ while ( Y <= UVPt( I,J ).y && J != 0 )
+ --J;
+ while ( Y > UVPt( I,J+1 ).y && J+2 < jSize )
+ ++J;
+ } while ( oldI != I || oldJ != J );
+}
+
+//================================================================================
+/*!
+ * \brief Looks for UV in quads around a given (I,J) and precise (I,J)
+ */
+//================================================================================
+
+bool FaceQuadStruct::isNear( const gp_XY& UV, int & I, int & J, int nbLoops )
+{
+ if ( I+1 >= iSize ) I = iSize - 2;
+ if ( J+1 >= jSize ) J = jSize - 2;
+
+ double bcI, bcJ;
+ gp_XY uvI, uvJ, uv0, uv1;
+ for ( int iLoop = 0; iLoop < nbLoops; ++iLoop )
+ {
+ int oldI = I, oldJ = J;
+
+ uvI = UVPt( I+1, J ).UV();
+ uvJ = UVPt( I, J+1 ).UV();
+ uv0 = UVPt( I, J ).UV();
+ SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv0, bcI, bcJ );
+ if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.)
+ return true;
+
+ if ( I > 0 && bcI < 0. ) --I;
+ if ( I+2 < iSize && bcI > 1. ) ++I;
+ if ( J > 0 && bcJ < 0. ) --J;
+ if ( J+2 < jSize && bcJ > 1. ) ++J;
+
+ uv1 = UVPt( I+1,J+1).UV();
+ if ( I != oldI || J != oldJ )
+ {
+ uvI = UVPt( I+1, J ).UV();
+ uvJ = UVPt( I, J+1 ).UV();
+ }
+ SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv1, bcI, bcJ );
+ if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.)
+ return true;
+
+ if ( I > 0 && bcI > 1. ) --I;
+ if ( I+2 < iSize && bcI < 0. ) ++I;
+ if ( J > 0 && bcJ > 1. ) --J;
+ if ( J+2 < jSize && bcJ < 0. ) ++J;
+
+ if ( I == oldI && J == oldJ )
+ return false;
+
+ if ( iLoop+1 == nbLoops )
+ {
+ uvI = UVPt( I+1, J ).UV();
+ uvJ = UVPt( I, J+1 ).UV();
+ uv0 = UVPt( I, J ).UV();
+ SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv0, bcI, bcJ );
+ if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.)
+ return true;
+
+ uv1 = UVPt( I+1,J+1).UV();
+ SMESH_MeshAlgos::GetBarycentricCoords( UV, uvI, uvJ, uv1, bcI, bcJ );
+ if ( bcI >= 0. && bcJ >= 0. && bcI + bcJ <= 1.)
+ return true;
+ }
+ }
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Checks if a given UV is equal to a given grid point
+ */
+//================================================================================
+
+bool FaceQuadStruct::isEqual( const gp_XY& UV, int I, int J )
+{
+ TopLoc_Location loc;
+ Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc );
+ gp_Pnt p1 = surf->Value( UV.X(), UV.Y() );
+ gp_Pnt p2 = surf->Value( UVPt( I,J ).u, UVPt( I,J ).v );
+
+ double dist2 = 1e100;
+ for ( int di = -1; di < 2; di += 2 )
+ {
+ int i = I + di;
+ if ( i < 0 || i+1 >= iSize ) continue;
+ for ( int dj = -1; dj < 2; dj += 2 )
+ {
+ int j = J + dj;
+ if ( j < 0 || j+1 >= jSize ) continue;
+
+ dist2 = Min( dist2,
+ p2.SquareDistance( surf->Value( UVPt( i,j ).u, UVPt( i,j ).v )));
+ }
+ }
+ double tol2 = dist2 / 1000.;
+ return p1.SquareDistance( p2 ) < tol2;
+}
+
+//================================================================================
+/*!
+ * \brief Recompute UV of grid points around a moved point in one direction
+ */
+//================================================================================
+
+void FaceQuadStruct::updateUV( const gp_XY& UV, int I, int J, bool isVertical )
+{
+ UVPt( I, J ).u = UV.X();
+ UVPt( I, J ).v = UV.Y();
+
+ if ( isVertical )
+ {
+ // above J
+ if ( J+1 < jSize-1 )
+ {
+ gp_UV a0 = UVPt( 0, J ).UV();
+ gp_UV a1 = UVPt( iSize-1, J ).UV();
+ gp_UV a2 = UVPt( iSize-1, jSize-1 ).UV();
+ gp_UV a3 = UVPt( 0, jSize-1 ).UV();
+
+ gp_UV p0 = UVPt( I, J ).UV();
+ gp_UV p2 = UVPt( I, jSize-1 ).UV();
+ const double y0 = UVPt( I, J ).y, dy = 1. - y0;
+ for (int j = J+1; j < jSize-1; j++)
+ {
+ gp_UV p1 = UVPt( iSize-1, j ).UV();
+ gp_UV p3 = UVPt( 0, j ).UV();
+
+ UVPtStruct& uvPt = UVPt( I, j );
+ gp_UV uv = calcUV( uvPt.x, ( uvPt.y - y0 ) / dy, a0,a1,a2,a3, p0,p1,p2,p3);
+ uvPt.u = uv.X();
+ uvPt.v = uv.Y();
+ }
+ }
+ // under J
+ if ( J-1 > 0 )
+ {
+ gp_UV a0 = UVPt( 0, 0 ).UV();
+ gp_UV a1 = UVPt( iSize-1, 0 ).UV();
+ gp_UV a2 = UVPt( iSize-1, J ).UV();
+ gp_UV a3 = UVPt( 0, J ).UV();
+
+ gp_UV p0 = UVPt( I, 0 ).UV();
+ gp_UV p2 = UVPt( I, J ).UV();
+ const double y0 = 0., dy = UVPt( I, J ).y - y0;
+ for (int j = 1; j < J; j++)
+ {
+ gp_UV p1 = UVPt( iSize-1, j ).UV();
+ gp_UV p3 = UVPt( 0, j ).UV();
+
+ UVPtStruct& uvPt = UVPt( I, j );
+ gp_UV uv = calcUV( uvPt.x, ( uvPt.y - y0 ) / dy, a0,a1,a2,a3, p0,p1,p2,p3);
+ uvPt.u = uv.X();
+ uvPt.v = uv.Y();
+ }
+ }
+ }
+ else // horizontally
+ {
+ // before I
+ if ( I-1 > 0 )
+ {
+ gp_UV a0 = UVPt( 0, 0 ).UV();
+ gp_UV a1 = UVPt( I, 0 ).UV();
+ gp_UV a2 = UVPt( I, jSize-1 ).UV();
+ gp_UV a3 = UVPt( 0, jSize-1 ).UV();
+
+ gp_UV p1 = UVPt( I, J ).UV();
+ gp_UV p3 = UVPt( 0, J ).UV();
+ const double x0 = 0., dx = UVPt( I, J ).x - x0;
+ for (int i = 1; i < I; i++)
+ {
+ gp_UV p0 = UVPt( i, 0 ).UV();
+ gp_UV p2 = UVPt( i, jSize-1 ).UV();
+
+ UVPtStruct& uvPt = UVPt( i, J );
+ gp_UV uv = calcUV(( uvPt.x - x0 ) / dx , uvPt.y, a0,a1,a2,a3, p0,p1,p2,p3);
+ uvPt.u = uv.X();
+ uvPt.v = uv.Y();
+ }
+ }
+ // after I
+ if ( I+1 < iSize-1 )
+ {
+ gp_UV a0 = UVPt( I, 0 ).UV();
+ gp_UV a1 = UVPt( iSize-1, 0 ).UV();
+ gp_UV a2 = UVPt( iSize-1, jSize-1 ).UV();
+ gp_UV a3 = UVPt( I, jSize-1 ).UV();
+
+ gp_UV p1 = UVPt( iSize-1, J ).UV();
+ gp_UV p3 = UVPt( I, J ).UV();
+ const double x0 = UVPt( I, J ).x, dx = 1. - x0;
+ for (int i = I+1; i < iSize-1; i++)
+ {
+ gp_UV p0 = UVPt( i, 0 ).UV();
+ gp_UV p2 = UVPt( i, jSize-1 ).UV();
+
+ UVPtStruct& uvPt = UVPt( i, J );
+ gp_UV uv = calcUV(( uvPt.x - x0 ) / dx , uvPt.y, a0,a1,a2,a3, p0,p1,p2,p3);
+ uvPt.u = uv.X();
+ uvPt.v = uv.Y();
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Side copying
+ */
+//================================================================================
+
+FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide)
+{
+ grid = otherSide.grid;
+ from = otherSide.from;
+ to = otherSide.to;
+ di = otherSide.di;
+ forced_nodes = otherSide.forced_nodes;
+ contacts = otherSide.contacts;
+ nbNodeOut = otherSide.nbNodeOut;
+
+ for ( size_t iC = 0; iC < contacts.size(); ++iC )
+ {
+ FaceQuadStruct::Side* oSide = contacts[iC].other_side;
+ for ( size_t iOC = 0; iOC < oSide->contacts.size(); ++iOC )
+ if ( oSide->contacts[iOC].other_side == & otherSide )
+ {
+ // cout << "SHIFT old " << &otherSide << " " << otherSide.NbPoints()
+ // << " -> new " << this << " " << this->NbPoints() << endl;
+ oSide->contacts[iOC].other_side = this;
+ }
+ }
+ return *this;
+}
+
+//================================================================================
+/*!
+ * \brief Converts node index of a quad to node index of this side
+ */
+//================================================================================
+
+int FaceQuadStruct::Side::ToSideIndex( int quadNodeIndex ) const
+{
+ return from + di * quadNodeIndex;
+}
+
+//================================================================================
+/*!
+ * \brief Converts node index of this side to node index of a quad
+ */
+//================================================================================
+
+int FaceQuadStruct::Side::ToQuadIndex( int sideNodeIndex ) const
+{
+ return ( sideNodeIndex - from ) * di;
+}
+
+//================================================================================
+/*!
+ * \brief Reverse the side
+ */
+//================================================================================
+
+bool FaceQuadStruct::Side::Reverse(bool keepGrid)
+{
+ if ( grid )
+ {
+ if ( keepGrid )
+ {
+ from -= di;
+ to -= di;
+ std::swap( from, to );
+ di *= -1;
+ }
+ else
+ {
+ grid->Reverse();
+ }
+ }
+ return (bool)grid;
+}
+
+//================================================================================
+/*!
+ * \brief Checks if a node is enforced
+ * \param [in] nodeIndex - an index of a node in a size
+ * \return bool - \c true if the node is forced
+ */
+//================================================================================
+
+bool FaceQuadStruct::Side::IsForced( int nodeIndex ) const
+{
+ if ( nodeIndex < 0 || nodeIndex >= grid->NbPoints() )
+ throw SALOME_Exception( " FaceQuadStruct::Side::IsForced(): wrong index" );
+
+ if ( forced_nodes.count( nodeIndex ) )
+ return true;
+
+ for ( size_t i = 0; i < this->contacts.size(); ++i )
+ if ( contacts[ i ].point == nodeIndex &&
+ contacts[ i ].other_side->forced_nodes.count( contacts[ i ].other_point ))
+ return true;
+
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Sets up a contact between this and another side
+ */
+//================================================================================
+
+void FaceQuadStruct::Side::AddContact( int ip, Side* side, int iop )
+{
+ if ( ip >= GetUVPtStruct().size() ||
+ iop >= side->GetUVPtStruct().size() )
+ throw SALOME_Exception( "FaceQuadStruct::Side::AddContact(): wrong point" );
+ {
+ contacts.resize( contacts.size() + 1 );
+ Contact& c = contacts.back();
+ c.point = ip;
+ c.other_side = side;
+ c.other_point = iop;
+ }
+ {
+ side->contacts.resize( side->contacts.size() + 1 );
+ Contact& c = side->contacts.back();
+ c.point = iop;
+ c.other_side = this;
+ c.other_point = ip;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Returns a normalized parameter of a point indexed within a quadrangle
+ */
+//================================================================================
+
+double FaceQuadStruct::Side::Param( int i ) const
+{
+ const vector<UVPtStruct>& points = GetUVPtStruct();
+ return (( points[ from + i * di ].normParam - points[ from ].normParam ) /
+ ( points[ to - 1 * di ].normParam - points[ from ].normParam ));
+}
+
+//================================================================================
+/*!
+ * \brief Returns UV by a parameter normalized within a quadrangle
+ */
+//================================================================================
+
+gp_XY FaceQuadStruct::Side::Value2d( double x ) const
+{
+ const vector<UVPtStruct>& points = GetUVPtStruct();
+ double u = ( points[ from ].normParam +
+ x * ( points[ to-di ].normParam - points[ from ].normParam ));
+ return grid->Value2d( u ).XY();
+}
+
+//================================================================================
+/*!
+ * \brief Returns side length
+ */
+//================================================================================
+
+double FaceQuadStruct::Side::Length(int theFrom, int theTo) const
+{
+ if ( IsReversed() != ( theTo < theFrom ))
+ std::swap( theTo, theFrom );
+
+ const vector<UVPtStruct>& points = GetUVPtStruct();
+ double r;
+ if ( theFrom == theTo && theTo == -1 )
+ r = Abs( First().normParam -
+ Last ().normParam );
+ else if ( IsReversed() )
+ r = Abs( points[ Max( to, theTo+1 ) ].normParam -
+ points[ Min( from, theFrom ) ].normParam );
+ else
+ r = Abs( points[ Min( to, theTo-1 ) ].normParam -
+ points[ Max( from, theFrom ) ].normParam );
+ return r * grid->Length();