-// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
#include "StdMeshers_Quadrangle_2D.hxx"
-#include "StdMeshers_FaceSide.hxx"
-
-#include "StdMeshers_QuadrangleParams.hxx"
-
+#include "SMDS_EdgePosition.hxx"
+#include "SMDS_FacePosition.hxx"
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMESH_Block.hxx"
+#include "SMESH_Comment.hxx"
#include "SMESH_Gen.hxx"
#include "SMESH_Mesh.hxx"
-#include "SMESH_subMesh.hxx"
#include "SMESH_MesherHelper.hxx"
-#include "SMESH_Block.hxx"
-#include "SMESH_Comment.hxx"
-
-#include "SMDS_MeshElement.hxx"
-#include "SMDS_MeshNode.hxx"
-#include "SMDS_EdgePosition.hxx"
-#include "SMDS_FacePosition.hxx"
+#include "SMESH_subMesh.hxx"
+#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_QuadrangleParams.hxx"
+#include "StdMeshers_ViscousLayers2D.hxx"
#include <BRep_Tool.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <Geom_Surface.hxx>
#include <NCollection_DefineArray2.hxx>
#include <Precision.hxx>
-#include <TColStd_SequenceOfReal.hxx>
+#include <Quantity_Parameter.hxx>
#include <TColStd_SequenceOfInteger.hxx>
+#include <TColStd_SequenceOfReal.hxx>
#include <TColgp_SequenceOfXY.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
+#include <TopTools_DataMapOfShapeReal.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopTools_MapOfShape.hxx>
#include <TopoDS.hxx>
//=============================================================================
/*!
- *
+ *
*/
//=============================================================================
StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
SMESH_Gen* gen)
- : SMESH_2D_Algo(hypId, studyId, gen)
+ : SMESH_2D_Algo(hypId, studyId, gen),
+ myQuadranglePreference(false),
+ myTrianglePreference(false),
+ myTriaVertexID(-1),
+ myNeedSmooth(false),
+ myQuadType(QUAD_STANDARD),
+ myHelper( 0 )
{
MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
_name = "Quadrangle_2D";
_compatibleHypothesis.push_back("QuadrangleParams");
_compatibleHypothesis.push_back("QuadranglePreference");
_compatibleHypothesis.push_back("TrianglePreference");
- myHelper = 0;
+ _compatibleHypothesis.push_back("ViscousLayers2D");
}
//=============================================================================
/*!
- *
+ *
*/
//=============================================================================
myQuadType = QUAD_STANDARD;
myQuadranglePreference = false;
myTrianglePreference = false;
+ myQuadStruct.reset();
bool isFirstParams = true;
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
- const TopoDS_Shape& aShape)// throw (SALOME_Exception)
+bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape)
{
- // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
- //Unexpect aCatchSalomeException);
+ const TopoDS_Face& F = TopoDS::Face(aShape);
+ Handle(Geom_Surface) S = BRep_Tool::Surface(F);
SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
aMesh.GetSubMesh(aShape);
SMESH_MesherHelper helper (aMesh);
myHelper = &helper;
+ myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+ if ( !myProxyMesh )
+ return false;
+
_quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
myNeedSmooth = false;
- FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
- std::auto_ptr<FaceQuadStruct> quadDeleter (quad); // to delete quad at exit from Compute()
+ FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape);
if (!quad)
return false;
+ myQuadStruct = quad;
if (myQuadranglePreference) {
int n1 = quad->side[0]->NbPoints();
Smooth( quad );
return ok;
}
+ if ( n1 != n3 && n2 != n4 )
+ error( COMPERR_WARNING,
+ "To use 'Reduced' transition, "
+ "two opposite sides should have same number of segments, "
+ "but actual number of segments is different on all sides. "
+ "'Standard' transion has been used.");
+ else
+ error( COMPERR_WARNING,
+ "To use 'Reduced' transition, "
+ "two opposite sides should have an even difference in number of segments. "
+ "'Standard' transion has been used.");
}
// set normalized grid on unit square in parametric domain
int nbhoriz = Min(nbdown, nbup);
int nbvertic = Min(nbright, nbleft);
- const TopoDS_Face& F = TopoDS::Face(aShape);
- Handle(Geom_Surface) S = BRep_Tool::Surface(F);
-
// internal mesh nodes
int i, j, geomFaceID = meshDS->ShapeToIndex(F);
for (i = 1; i < nbhoriz - 1; i++) {
for (i = ilow; i < iup; i++) {
for (j = jlow; j < jup; j++) {
const SMDS_MeshNode *a, *b, *c, *d;
- a = quad->uv_grid[j * nbhoriz + i].node;
- b = quad->uv_grid[j * nbhoriz + i + 1].node;
+ a = quad->uv_grid[j * nbhoriz + i ].node;
+ b = quad->uv_grid[j * nbhoriz + i + 1].node;
c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
- d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
+ d = quad->uv_grid[(j + 1) * nbhoriz + i ].node;
SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
if (face) {
meshDS->SetMeshElementOnShape(face, geomFaceID);
*/
//=============================================================================
-FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
- const TopoDS_Shape & aShape)
- //throw(SALOME_Exception)
+FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
+ const TopoDS_Shape & aShape)
{
+ if ( myQuadStruct && myQuadStruct->face.IsSame( aShape ))
+ return myQuadStruct;
+
TopoDS_Face F = TopoDS::Face(aShape);
if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
const bool ignoreMediumNodes = _quadraticMesh;
// verify 1 wire only, with 4 edges
- TopoDS_Vertex V;
list< TopoDS_Edge > edges;
list< int > nbEdgesInWire;
- int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
if (nbWire != 1) {
error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
- return 0;
+ return FaceQuadStruct::Ptr();
}
- FaceQuadStruct* quad = new FaceQuadStruct;
+
+ // find corner vertices of the quad
+ vector<TopoDS_Vertex> corners;
+ int nbDegenEdges, nbSides = GetCorners( F, aMesh, edges, corners, nbDegenEdges );
+ if ( nbSides == 0 )
+ {
+ return FaceQuadStruct::Ptr();
+ }
+ FaceQuadStruct::Ptr quad( new FaceQuadStruct );
quad->uv_grid = 0;
quad->side.reserve(nbEdgesInWire.front());
quad->face = F;
- int nbSides = 0;
list< TopoDS_Edge >::iterator edgeIt = edges.begin();
- if (nbEdgesInWire.front() == 3) // exactly 3 edges
+ if ( nbSides == 3 ) // 3 sides and corners[0] is a vertex with myTriaVertexID
{
- SMESH_Comment comment;
- SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
- if (myTriaVertexID == -1)
+ for ( int iSide = 0; iSide < 3; ++iSide )
{
- comment << "No Base vertex parameter provided for a trilateral geometrical face";
- }
- else
+ list< TopoDS_Edge > sideEdges;
+ TopoDS_Vertex nextSideV = corners[( iSide + 1 ) % 3 ];
+ while ( edgeIt != edges.end() &&
+ !nextSideV.IsSame( SMESH_MesherHelper::IthVertex( 0, *edgeIt )))
+ if ( SMESH_Algo::isDegenerated( *edgeIt ))
+ ++edgeIt;
+ else
+ sideEdges.push_back( *edgeIt++ );
+ if ( !sideEdges.empty() )
+ quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE,
+ ignoreMediumNodes, myProxyMesh));
+ else
+ --iSide;
+ }
+ const vector<UVPtStruct>& UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
+ /* vector<UVPtStruct>& UVPStop = */quad->side[1]->GetUVPtStruct(false,1);
+ /* vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
+ const SMDS_MeshNode* aNode = UVPSleft[0].node;
+ gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v);
+ quad->side.push_back(new StdMeshers_FaceSide(quad->side[1], aNode, &aPnt2d));
+ myNeedSmooth = ( nbDegenEdges > 0 );
+ return quad;
+ }
+ else // 4 sides
+ {
+ myNeedSmooth = ( corners.size() == 4 && nbDegenEdges > 0 );
+ int iSide = 0, nbUsedDegen = 0, nbLoops = 0;
+ for ( ; edgeIt != edges.end(); ++nbLoops )
{
- TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
- if (!V.IsNull()) {
- TopoDS_Edge E1,E2,E3;
- for (; edgeIt != edges.end(); ++edgeIt) {
- TopoDS_Edge E = *edgeIt;
- TopoDS_Vertex VF, VL;
- TopExp::Vertices(E, VF, VL, true);
- if (VF.IsSame(V))
- E1 = E;
- else if (VL.IsSame(V))
- E3 = E;
+ list< TopoDS_Edge > sideEdges;
+ TopoDS_Vertex nextSideV = corners[( iSide + 1 - nbUsedDegen ) % corners.size() ];
+ while ( edgeIt != edges.end() &&
+ !nextSideV.IsSame( myHelper->IthVertex( 0, *edgeIt )))
+ {
+ if ( SMESH_Algo::isDegenerated( *edgeIt ))
+ {
+ if ( myNeedSmooth )
+ {
+ ++edgeIt; // no side on the degenerated EDGE
+ }
else
- E2 = E;
+ {
+ if ( sideEdges.empty() )
+ {
+ ++nbUsedDegen;
+ sideEdges.push_back( *edgeIt++ ); // a degenerated side
+ break;
+ }
+ else
+ {
+ break; // do not append a degenerated EDGE to a regular side
+ }
+ }
}
- if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull())
+ else
{
- quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes));
- quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes));
- quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes));
- const vector<UVPtStruct>& UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
- /* vector<UVPtStruct>& UVPStop = */quad->side[1]->GetUVPtStruct(false,1);
- /* vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
- const SMDS_MeshNode* aNode = UVPSleft[0].node;
- gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v);
- quad->side.push_back(new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]));
- return quad;
- }
- }
- comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among [";
- TopTools_MapOfShape vMap;
- for (TopExp_Explorer v(aShape, TopAbs_VERTEX); v.More(); v.Next())
- if (vMap.Add(v.Current()))
- comment << meshDS->ShapeToIndex(v.Current()) << (vMap.Extent()==3 ? "]" : ", ");
- }
- error(comment);
- delete quad;
- return quad = 0;
- }
- else if (nbEdgesInWire.front() == 4) // exactly 4 edges
- {
- for (; edgeIt != edges.end(); ++edgeIt, nbSides++)
- quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
- nbSides<TOP_SIDE, ignoreMediumNodes));
- }
- else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some
- {
- list< TopoDS_Edge > sideEdges;
- vector< int > degenSides;
- while (!edges.empty()) {
- sideEdges.clear();
- sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
- bool sameSide = true;
- while (!edges.empty() && sameSide) {
- sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
- if (sameSide)
- sideEdges.splice(sideEdges.end(), edges, edges.begin());
- }
- if (nbSides == 0) { // go backward from the first edge
- sameSide = true;
- while (!edges.empty() && sameSide) {
- sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
- if (sameSide)
- sideEdges.splice(sideEdges.begin(), edges, --edges.end());
+ sideEdges.push_back( *edgeIt++ );
}
}
- if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() ))
- degenSides.push_back( nbSides );
-
- quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
- nbSides<TOP_SIDE, ignoreMediumNodes));
- ++nbSides;
- }
- if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
- {
- myNeedSmooth = true;
- for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
- quad->side[i]->Reverse();
-
- for ( int i = degenSides.size()-1; i > -1; --i )
+ if ( !sideEdges.empty() )
{
- StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]];
- delete degenSide;
- quad->side.erase( quad->side.begin() + degenSides[ i ] );
+ quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE,
+ ignoreMediumNodes, myProxyMesh));
+ ++iSide;
}
- for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
- quad->side[i]->Reverse();
-
- nbSides -= degenSides.size();
- }
- // issue 20222. Try to unite only edges shared by two same faces
- if (nbSides < 4) {
- // delete found sides
- { FaceQuadStruct cleaner(*quad); }
- quad->side.clear();
- quad->side.reserve(nbEdgesInWire.front());
- nbSides = 0;
-
- SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
- while (!edges.empty()) {
- sideEdges.clear();
- sideEdges.splice(sideEdges.end(), edges, edges.begin());
- bool sameSide = true;
- while (!edges.empty() && sameSide) {
- sameSide =
- SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
- twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
- if (sameSide)
- sideEdges.splice(sideEdges.end(), edges, edges.begin());
- }
- if (nbSides == 0) { // go backward from the first edge
- sameSide = true;
- while (!edges.empty() && sameSide) {
- sameSide =
- SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
- twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
- if (sameSide)
- sideEdges.splice(sideEdges.begin(), edges, --edges.end());
- }
- }
- quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
- nbSides<TOP_SIDE, ignoreMediumNodes));
- ++nbSides;
+ if ( nbLoops > 8 )
+ {
+ error(TComm("Bug: infinite loop in StdMeshers_Quadrangle_2D::CheckNbEdges()"));
+ quad.reset();
+ break;
}
}
- }
- if (nbSides != 4) {
-#ifdef _DEBUG_
- MESSAGE ("StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n");
- for (int i = 0; i < nbSides; ++i) {
- MESSAGE (" (");
- for (int e = 0; e < quad->side[i]->NbEdges(); ++e)
- MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " ");
- MESSAGE (")\n");
- }
- //cout << endl;
-#endif
- if (!nbSides)
- nbSides = nbEdgesInWire.front();
- error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
- delete quad;
- quad = 0;
+ if ( quad->side.size() != 4 )
+ {
+ error(TComm("Bug: ") << quad->side.size() << " sides found instead of 4");
+ quad.reset();
+ }
}
return quad;
const TopoDS_Face & F = TopoDS::Face(aShape);
// verify 1 wire only, with 4 edges
- TopoDS_Vertex V;
list< TopoDS_Edge > edges;
list< int > nbEdgesInWire;
- int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
if (nbWire != 1) {
return false;
}
// issue 20222. Try to unite only edges shared by two same faces
if (nbSides < 4) {
nbSides = 0;
- SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
while (!edges.empty()) {
sideEdges.clear();
sideEdges.splice(sideEdges.end(), edges, edges.begin());
*/
//=============================================================================
-FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
- (SMESH_Mesh & aMesh,
- const TopoDS_Shape & aShape,
- const bool CreateQuadratic) //throw(SALOME_Exception)
+FaceQuadStruct::Ptr
+StdMeshers_Quadrangle_2D::CheckAnd2Dcompute (SMESH_Mesh & aMesh,
+ const TopoDS_Shape & aShape,
+ const bool CreateQuadratic)
{
_quadraticMesh = CreateQuadratic;
- FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
-
- if (!quad) return 0;
-
- // set normalized grid on unit square in parametric domain
- bool stat = SetNormalizedGrid(aMesh, aShape, quad);
- if (!stat) {
- if (quad) delete quad;
- quad = 0;
+ FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape);
+ if ( quad )
+ {
+ // set normalized grid on unit square in parametric domain
+ if (!SetNormalizedGrid(aMesh, aShape, quad))
+ quad.reset();
}
-
return quad;
}
//=============================================================================
/*!
- *
+ *
*/
//=============================================================================
faceQuadStruct::~faceQuadStruct()
{
- for (int i = 0; i < side.size(); i++) {
- if (side[i]) delete side[i];
+ for (size_t i = 0; i < side.size(); i++) {
+ if (side[i]) {
+ delete side[i];
+ for (size_t j = i+1; j < side.size(); j++)
+ if ( side[i] == side[j] )
+ side[j] = 0;
+ }
+ }
+ side.clear();
+
+ if (uv_grid) {
+ delete [] uv_grid;
+ uv_grid = 0;
}
- if (uv_grid) delete [] uv_grid;
}
-namespace {
- inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
+namespace
+{
+ inline const vector<UVPtStruct>& getUVPtStructIn(FaceQuadStruct::Ptr& quad, int i, int nbSeg)
{
- bool isXConst = (i == BOTTOM_SIDE || i == TOP_SIDE);
- double constValue = (i == BOTTOM_SIDE || i == LEFT_SIDE) ? 0 : 1;
+ bool isXConst = (i == QUAD_BOTTOM_SIDE || i == QUAD_TOP_SIDE);
+ double constValue = (i == QUAD_BOTTOM_SIDE || i == QUAD_LEFT_SIDE) ? 0 : 1;
return
quad->isEdgeOut[i] ?
quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
quad->side[i]->GetUVPtStruct(isXConst,constValue);
}
- inline gp_UV CalcUV(double x, double y,
+ inline gp_UV calcUV(double x, double y,
const gp_UV& a0,const gp_UV& a1,const gp_UV& a2,const gp_UV& a3,
const gp_UV& p0,const gp_UV& p1,const gp_UV& p2,const gp_UV& p3)
{
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
- const TopoDS_Shape& aShape,
- FaceQuadStruct* & quad) //throw (SALOME_Exception)
+bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
+ const TopoDS_Shape& aShape,
+ FaceQuadStruct::Ptr & quad)
{
// Algorithme décrit dans "Génération automatique de maillages"
// P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
// 3 --- 2D normalized values on unit square [0..1][0..1]
+ UpdateDegenUV( quad );
+
int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
- const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn(quad, 0, nbhoriz - 1);
- const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn(quad, 1, nbvertic - 1);
- const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn(quad, 2, nbhoriz - 1);
- const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn(quad, 3, nbvertic - 1);
+ const vector<UVPtStruct>& uv_e0 = getUVPtStructIn(quad, 0, nbhoriz - 1);
+ const vector<UVPtStruct>& uv_e1 = getUVPtStructIn(quad, 1, nbvertic - 1);
+ const vector<UVPtStruct>& uv_e2 = getUVPtStructIn(quad, 2, nbhoriz - 1);
+ const vector<UVPtStruct>& uv_e3 = getUVPtStructIn(quad, 3, nbvertic - 1);
if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
//return error("Can't find nodes on sides");
return error(COMPERR_BAD_INPUT_MESH);
- if ( myNeedSmooth )
- UpdateDegenUV( quad );
-
- // nodes Id on "in" edges
- if (! quad->isEdgeOut[0]) {
- int j = 0;
- for (int i = 0; i < nbhoriz; i++) { // down
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e0[i].node;
- }
+ // copy data of face boundary
+ /*if (! quad->isEdgeOut[0])*/ {
+ const int j = 0;
+ for (int i = 0; i < nbhoriz; i++) // down
+ uv_grid[ j * nbhoriz + i ] = uv_e0[i];
}
- if (! quad->isEdgeOut[1]) {
- int i = nbhoriz - 1;
- for (int j = 0; j < nbvertic; j++) { // right
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e1[j].node;
- }
+ /*if (! quad->isEdgeOut[1])*/ {
+ const int i = nbhoriz - 1;
+ for (int j = 0; j < nbvertic; j++) // right
+ uv_grid[ j * nbhoriz + i ] = uv_e1[j];
}
- if (! quad->isEdgeOut[2]) {
- int j = nbvertic - 1;
- for (int i = 0; i < nbhoriz; i++) { // up
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e2[i].node;
- }
+ /*if (! quad->isEdgeOut[2])*/ {
+ const int j = nbvertic - 1;
+ for (int i = 0; i < nbhoriz; i++) // up
+ uv_grid[ j * nbhoriz + i ] = uv_e2[i];
}
- if (! quad->isEdgeOut[3]) {
+ /*if (! quad->isEdgeOut[3])*/ {
int i = 0;
- for (int j = 0; j < nbvertic; j++) { // left
- int ij = j * nbhoriz + i;
- uv_grid[ij].node = uv_e3[j].node;
- }
+ for (int j = 0; j < nbvertic; j++) // left
+ uv_grid[ j * nbhoriz + i ] = uv_e3[j];
}
- // normalized 2d values on grid
+ // normalized 2d parameters on grid
for (int i = 0; i < nbhoriz; i++) {
for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
}
// 4 --- projection on 2d domain (u,v)
- gp_UV a0(uv_e0.front().u, uv_e0.front().v);
- gp_UV a1(uv_e0.back().u, uv_e0.back().v);
- gp_UV a2(uv_e2.back().u, uv_e2.back().v);
- gp_UV a3(uv_e2.front().u, uv_e2.front().v);
+ gp_UV a0 (uv_e0.front().u, uv_e0.front().v);
+ gp_UV a1 (uv_e0.back().u, uv_e0.back().v );
+ gp_UV a2 (uv_e2.back().u, uv_e2.back().v );
+ gp_UV a3 (uv_e2.front().u, uv_e2.front().v);
+
+ for (int i = 0; i < nbhoriz; i++)
+ {
+ gp_UV p0( uv_e0[i].u, uv_e0[i].v );
+ gp_UV p2( uv_e2[i].u, uv_e2[i].v );
+ for (int j = 0; j < nbvertic; j++)
+ {
+ gp_UV p1( uv_e1[j].u, uv_e1[j].v );
+ gp_UV p3( uv_e3[j].u, uv_e3[j].v );
- for (int i = 0; i < nbhoriz; i++) {
- for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
double x = uv_grid[ij].x;
double y = uv_grid[ij].y;
- double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam); // sud
- double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam); // nord
- double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam); // est
- double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam); // ouest
-
- //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
- gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
- gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
- gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
- gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
- gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
+ gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
uv_grid[ij].u = uv.X();
uv_grid[ij].v = uv.Y();
//purpose : auxilary function for ComputeQuadPref
//=======================================================================
-static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
+static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num, bool)
{
- StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
- for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i) {
- int id = (i + num) % NB_SIDES;
- bool wasForward = (i < TOP_SIDE);
- bool newForward = (id < TOP_SIDE);
- if (wasForward != newForward)
- side[ i ]->Reverse();
- quad->side[ id ] = side[ i ];
+ quad->shift( num, /*ori=*/true );
+}
+
+//================================================================================
+/*!
+ * \brief Rotate sides of a quad by nb
+ * \param nb - number of rotation quartes
+ * \param ori - to keep orientation of sides as in an unit quad or not
+ */
+//================================================================================
+
+void FaceQuadStruct::shift( size_t nb, bool ori )
+{
+ if ( nb == 0 ) return;
+ StdMeshers_FaceSide* sideArr[4] = { side[0], side[1], side[2], side[3] };
+ for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) {
+ int id = (i + nb) % NB_QUAD_SIDES;
+ bool wasForward = (i < QUAD_TOP_SIDE);
+ bool newForward = (id < QUAD_TOP_SIDE);
+ if (ori && wasForward != newForward)
+ sideArr[ i ]->Reverse();
+ side[ id ] = sideArr[ i ];
}
}
//=======================================================================
-//function : CalcUV
+//function : calcUV
//purpose : auxilary function for ComputeQuadPref
//=======================================================================
-static gp_UV CalcUV(double x0, double x1, double y0, double y1,
- FaceQuadStruct* quad,
+static gp_UV calcUV(double x0, double x1, double y0, double y1,
+ FaceQuadStruct::Ptr& quad,
const gp_UV& a0, const gp_UV& a1,
const gp_UV& a2, const gp_UV& a3)
{
- const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
- const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
- const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
- const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
-
double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
double y = y0 + x * (y1 - y0);
- double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
- double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
- double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
- double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
+ gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY();
+ gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY();
+ gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY();
+ gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY();
- gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
- gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
- gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(param_t).XY();
- gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(param_l).XY();
-
- gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
+ gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
return uv;
}
//=======================================================================
-//function : CalcUV2
+//function : calcUV2
//purpose : auxilary function for ComputeQuadPref
//=======================================================================
-static gp_UV CalcUV2(double x, double y,
- FaceQuadStruct* quad,
+static gp_UV calcUV2(double x, double y,
+ FaceQuadStruct::Ptr& quad,
const gp_UV& a0, const gp_UV& a1,
const gp_UV& a2, const gp_UV& a3)
{
- gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(x).XY();
- gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(y).XY();
- gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(x).XY();
- gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(y).XY();
+ gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY();
+ gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY();
+ gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY();
+ gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY();
- gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
+ gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
return uv;
}
bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh,
const TopoDS_Shape& aShape,
- FaceQuadStruct* quad)
+ FaceQuadStruct::Ptr quad)
{
// Auxilary key in order to keep old variant
// of meshing after implementation new variant
bool WisF = true;
int i,j,geomFaceID = meshDS->ShapeToIndex(F);
+ UpdateDegenUV( quad );
+
int nb = quad->side[0]->NbPoints();
int nr = quad->side[1]->NbPoints();
int nt = quad->side[2]->NbPoints();
if (dh>=dv) {
if (nt>nb) {
// it is a base case => not shift quad but me be replacement is need
- ShiftQuad(quad,0,WisF);
+ shiftQuad(quad,0,WisF);
}
else {
// we have to shift quad on 2
- ShiftQuad(quad,2,WisF);
+ shiftQuad(quad,2,WisF);
}
}
else {
if (nr>nl) {
// we have to shift quad on 1
- ShiftQuad(quad,1,WisF);
+ shiftQuad(quad,1,WisF);
}
else {
// we have to shift quad on 3
- ShiftQuad(quad,3,WisF);
+ shiftQuad(quad,3,WisF);
}
}
int addh = 0;
int addv = 0;
+ // Orientation of face and 3 main domain for future faces
// ----------- Old version ---------------
- // orientation of face and 3 main domain for future faces
// 0 top 1
// 1------------1
// | | | |
- // | | | |
+ // | |C | |
// | L | | R |
- // left | | | | rigth
+ // left | |__| | rigth
// | / \ |
// | / C \ |
// |/ \|
// 0 bottom 1
// ----------- New version ---------------
- // orientation of face and 3 main domain for future faces
// 0 top 1
// 1------------1
- // | |____| |
+ // | |_C__| |
// | / \ |
// | / C \ |
// left |/________\| rigth
// | |
- // | |
+ // | C |
// | |
// 0------------0
// 0 bottom 1
if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
return error(COMPERR_BAD_INPUT_MESH);
- if ( myNeedSmooth )
- UpdateDegenUV( quad );
-
// arrays for normalized params
- //cout<<"Dump B:"<<endl;
TColStd_SequenceOfReal npb, npr, npt, npl;
for (i=0; i<nb; i++) {
npb.Append(uv_eb[i].normParam);
- //cout<<"i="<<i<<" par="<<uv_eb[i].normParam<<" npar="<<uv_eb[i].normParam;
- //const SMDS_MeshNode* N = uv_eb[i].node;
- //cout<<" node("<<N->X()<<","<<N->Y()<<","<<N->Z()<<")"<<endl;
}
for (i=0; i<nr; i++) {
npr.Append(uv_er[i].normParam);
npl.InsertAfter(1,npl.Value(2)-dpr);
}
}
- //cout<<"npb:";
- //for (i=1; i<=npb.Length(); i++) {
- // cout<<" "<<npb.Value(i);
- //}
- //cout<<endl;
gp_XY a0(uv_eb.front().u, uv_eb.front().v);
gp_XY a1(uv_eb.back().u, uv_eb.back().v);
gp_XY a2(uv_et.back().u, uv_et.back().v);
gp_XY a3(uv_et.front().u, uv_et.front().v);
- //cout<<" a0("<<a0.X()<<","<<a0.Y()<<")"<<" a1("<<a1.X()<<","<<a1.Y()<<")"
- // <<" a2("<<a2.X()<<","<<a2.Y()<<")"<<" a3("<<a3.X()<<","<<a3.Y()<<")"<<endl;
int nnn = Min(nr,nl);
// auxilary sequence of XY for creation nodes
// diagonal node
double y0 = npl.Value(i+1);
double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=2; j<nl; j++) {
double y0 = npl.Value(dl+j);
double y1 = npr.Value(dl+j);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn; i++) {
UVL.Append(UVtmp.Value(i));
}
- //cout<<"Dump NodesL:"<<endl;
- //for (i=1; i<=dl+1; i++) {
- // cout<<"i="<<i;
- // for (j=1; j<=nl; j++) {
- // cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
- // }
- // cout<<endl;
- //}
// create faces
for (i=1; i<=dl; i++) {
for (j=1; j<nl; j++) {
// diagonal node
double y0 = npl.Value(i+1);
double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=2; j<nr; j++) {
double y0 = npl.Value(nbv-j+1);
double y1 = npr.Value(nbv-j+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=1; j<nnn; j++) {
double y0 = npl.Value(nbv-nnn+j);
double y1 = npr.Value(nbv-nnn+j);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
}
}
// add diagonal layers
- //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
- //cout<<"Dump UVL:"<<endl;
- //for (i=1; i<=UVL.Length(); i++) {
- // cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
- //}
- //cout<<endl;
gp_UV A2 = UVR.Value(nbv-nnn);
gp_UV A3 = UVL.Value(nbv-nnn);
for (i=1; i<nbv-nnn; i++) {
double x = npb.Value(j);
gp_UV p0( uv_eb[j-1].u, uv_eb[j-1].v );
gp_UV p2 = UVT.Value( j-1 );
- gp_UV UV = CalcUV(x, y, a0, a1, A2, A3, p0,p1,p2,p3 );
+ gp_UV UV = calcUV(x, y, a0, a1, A2, A3, p0,p1,p2,p3 );
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
for (j=2; j<nb; j++) {
double x = npb.Value(j);
double y = (1-x) * npl.Value(i+1) + x * npr.Value(i+1);
- gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
for (j=1; j<nb; j++) {
double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
double y = y0 + dy*x;
- gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
double x = npt.Value(i+1+drl) +
npb.Value(j) * (npt.Value(nt-i) - npt.Value(i+1+drl));
double y = yy0 + dyy*x;
- gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=2; j<=nb; j++) {
double x = npb.Value(j)*npt.Value(nt-i);
double y = y0 + dy*x;
- gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
double x = npt.Value(i+1) +
npb.Value(j) * (npt.Value(nt-i-drl) - npt.Value(i+1));
double y = yy0 + dyy*x;
- gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
inline SMDS_MeshNode* makeNode( UVPtStruct & uvPt,
const double y,
- FaceQuadStruct* quad,
+ FaceQuadStruct::Ptr& quad,
const gp_UV* UVs,
SMESH_MesherHelper* helper,
Handle(Geom_Surface) S)
{
- const vector<UVPtStruct>& uv_eb = quad->side[BOTTOM_SIDE]->GetUVPtStruct();
- const vector<UVPtStruct>& uv_et = quad->side[TOP_SIDE ]->GetUVPtStruct();
+ const vector<UVPtStruct>& uv_eb = quad->side[QUAD_BOTTOM_SIDE]->GetUVPtStruct();
+ const vector<UVPtStruct>& uv_et = quad->side[QUAD_TOP_SIDE ]->GetUVPtStruct();
double rBot = ( uv_eb.size() - 1 ) * uvPt.normParam;
double rTop = ( uv_et.size() - 1 ) * uvPt.normParam;
int iBot = int( rBot );
double xTop = uv_et[ iTop ].normParam + ( rTop - iTop ) * ( uv_et[ iTop+1 ].normParam - uv_et[ iTop ].normParam );
double x = xBot + y * ( xTop - xBot );
- gp_UV uv = CalcUV(/*x,y=*/x, y,
+ gp_UV uv = calcUV(/*x,y=*/x, y,
/*a0,...=*/UVs[UV_A0], UVs[UV_A1], UVs[UV_A2], UVs[UV_A3],
- /*p0=*/quad->side[BOTTOM_SIDE]->Value2d( x ).XY(),
+ /*p0=*/quad->side[QUAD_BOTTOM_SIDE]->Value2d( x ).XY(),
/*p1=*/UVs[ UV_R ],
- /*p2=*/quad->side[TOP_SIDE ]->Value2d( x ).XY(),
+ /*p2=*/quad->side[QUAD_TOP_SIDE ]->Value2d( x ).XY(),
/*p3=*/UVs[ UV_L ]);
gp_Pnt P = S->Value( uv.X(), uv.Y() );
uvPt.u = uv.X();
vector<UVPtStruct>& next_base,
const int j,
int & next_base_len,
- FaceQuadStruct* quad,
+ FaceQuadStruct::Ptr& quad,
gp_UV* UVs,
const double y,
SMESH_MesherHelper* helper,
vector<UVPtStruct>& next_base,
const int j,
int & next_base_len,
- FaceQuadStruct* quad,
+ FaceQuadStruct::Ptr& quad,
gp_UV* UVs,
const double y,
SMESH_MesherHelper* helper,
vector<UVPtStruct>& next_base,
const int j,
int & next_base_len,
- FaceQuadStruct* quad,
+ FaceQuadStruct::Ptr & quad,
gp_UV* UVs,
const double y,
SMESH_MesherHelper* helper,
* Implementation of Reduced algorithm (meshing with quadrangles only)
*/
//=======================================================================
+
bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh,
const TopoDS_Shape& aShape,
- FaceQuadStruct* quad)
+ FaceQuadStruct::Ptr quad)
{
- SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
- const TopoDS_Face& F = TopoDS::Face(aShape);
+ SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
+ const TopoDS_Face& F = TopoDS::Face(aShape);
Handle(Geom_Surface) S = BRep_Tool::Surface(F);
- int i,j,geomFaceID = meshDS->ShapeToIndex(F);
+ int i,j,geomFaceID = meshDS->ShapeToIndex(F);
- int nb = quad->side[0]->NbPoints();
- int nr = quad->side[1]->NbPoints();
- int nt = quad->side[2]->NbPoints();
- int nl = quad->side[3]->NbPoints();
+ int nb = quad->side[0]->NbPoints(); // bottom
+ int nr = quad->side[1]->NbPoints(); // right
+ int nt = quad->side[2]->NbPoints(); // top
+ int nl = quad->side[3]->NbPoints(); // left
// Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step)
//
}
// number of rows and columns
- int nrows = nr1 - 1;
+ int nrows = nr1 - 1;
int ncol_top = nt1 - 1;
int ncol_bot = nb1 - 1;
// number of rows needed to reduce ncol_bot to ncol_top using simple 3->1 "tree" (see below)
- int nrows_tree31 = int( log( (double)(ncol_bot / ncol_top) ) / log((double) 3 )); // = log x base 3
+ int nrows_tree31 =
+ int( ceil( log( double(ncol_bot) / ncol_top) / log( 3.))); // = log x base 3
if ( nrows < nrows_tree31 )
+ {
MultipleReduce = true;
+ error( COMPERR_WARNING,
+ SMESH_Comment("To use 'Reduced' transition, "
+ "number of face rows should be at least ")
+ << nrows_tree31 << ". Actual number of face rows is " << nrows << ". "
+ "'Quadrangle preference (reversed)' transion has been used.");
+ }
}
if (MultipleReduce) { // == ComputeQuadPref QUAD_QUADRANGLE_PREF_REVERSED
if (dh >= dv) {
if (nt > nb) {
// it is a base case => not shift quad but may be replacement is need
- ShiftQuad(quad,0,true);
+ shiftQuad(quad,0,true);
}
else {
// we have to shift quad on 2
- ShiftQuad(quad,2,true);
+ shiftQuad(quad,2,true);
}
}
else {
if (nr > nl) {
// we have to shift quad on 1
- ShiftQuad(quad,1,true);
+ shiftQuad(quad,1,true);
}
else {
// we have to shift quad on 3
- ShiftQuad(quad,3,true);
+ shiftQuad(quad,3,true);
}
}
if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
return error(COMPERR_BAD_INPUT_MESH);
- if ( myNeedSmooth )
- UpdateDegenUV( quad );
+ UpdateDegenUV( quad );
// arrays for normalized params
TColStd_SequenceOfReal npb, npr, npt, npl;
// diagonal node
double y0 = npl.Value(i+1);
double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=2; j<nl; j++) {
double y0 = npl.Value(dl+j);
double y1 = npr.Value(dl+j);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
// diagonal node
double y0 = npl.Value(i+1);
double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
for (j=2; j<nr; j++) {
double y0 = npl.Value(nbv-j+1);
double y1 = npr.Value(nbv-j+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
// add bottom nodes (first columns)
for (i=2; i<nb; i++)
NodesC.SetValue(i,1,uv_eb[i-1].node);
-
+
// create and add needed nodes
// add linear layers
for (i=2; i<nb; i++) {
for (j=1; j<nnn; j++) {
double y0 = npl.Value(nbv-nnn+j);
double y1 = npr.Value(nbv-nnn+j);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
// create faces
for (i=1; i<nb; i++) {
for (j=1; j<nbv; j++) {
- SMDS_MeshFace* F =
- myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
+ SMDS_MeshFace* F =
+ myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
- if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
}
- // TODO ???
} // end Multiple Reduce implementation
else { // Simple Reduce (!MultipleReduce)
//=========================================================
if (nr == nl) {
if (nt < nb) {
// it is a base case => not shift quad
- //ShiftQuad(quad,0,true);
+ //shiftQuad(quad,0,true);
}
else {
// we have to shift quad on 2
- ShiftQuad(quad,2,true);
+ shiftQuad(quad,2,true);
}
}
else {
if (nl > nr) {
// we have to shift quad on 1
- ShiftQuad(quad,1,true);
+ shiftQuad(quad,1,true);
}
else {
// we have to shift quad on 3
- ShiftQuad(quad,3,true);
+ shiftQuad(quad,3,true);
}
}
vector<int> nb_col_by_row;
- int delta_all = nb - nt;
+ int delta_all = nb - nt;
int delta_one_col = nrows * 2;
- int nb_col = delta_all / delta_one_col;
- int remainder = delta_all - nb_col * delta_one_col;
+ int nb_col = delta_all / delta_one_col;
+ int remainder = delta_all - nb_col * delta_one_col;
if (remainder > 0) {
nb_col++;
}
nb_col = ( nt - 1 ) / col_top_size;
nb_col_by_row.resize( nrows, nb_col );
int nbrows_not_full = nrows - 1;
- int cur_top_size = nt - 1;
+ int cur_top_size = nt - 1;
remainder = delta_all - nb_col * delta_one_col;
while ( remainder > 0 )
{
- delta_one_col = nbrows_not_full * 2;
- int nb_col_add = remainder / delta_one_col;
- cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ];
+ delta_one_col = nbrows_not_full * 2;
+ int nb_col_add = remainder / delta_one_col;
+ cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ];
int nb_col_free = cur_top_size / col_top_size - nb_col_by_row[ nbrows_not_full-1 ];
if ( nb_col_add > nb_col_free )
nb_col_add = nb_col_free;
*/
struct TSmoothNode
{
- gp_XY _uv;
+ gp_XY _uv;
+ gp_XYZ _xyz;
vector< TTriangle > _triangles; // if empty, then node is not movable
};
// --------------------------------------------------------------------------------
*/
//================================================================================
-void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad)
+void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct::Ptr quad)
{
- for ( unsigned i = 0; i < quad->side.size(); ++i )
- {
- StdMeshers_FaceSide* side = quad->side[i];
- const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
-
- // find which end of the side is on degenerated shape
- int degenInd = -1;
- if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
- degenInd = 0;
- else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
- degenInd = uvVec.size() - 1;
- else
- continue;
+ if ( myNeedSmooth )
- // find another side sharing the degenerated shape
- bool isPrev = ( degenInd == 0 );
- if ( i >= TOP_SIDE )
- isPrev = !isPrev;
- int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
- StdMeshers_FaceSide* side2 = quad->side[ i2 ];
- const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
- int degenInd2 = -1;
- if ( uvVec[ degenInd ].node == uvVec2[0].node )
- degenInd2 = 0;
- else if ( uvVec[ degenInd ].node == uvVec2.back().node )
- degenInd2 = uvVec2.size() - 1;
- else
- throw SALOME_Exception( LOCALIZED( "Logical error" ));
+ // Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE
+ // --------------------------------------------------------------------------
+ for ( unsigned i = 0; i < quad->side.size(); ++i )
+ {
+ StdMeshers_FaceSide* side = quad->side[i];
+ const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
+
+ // find which end of the side is on degenerated shape
+ int degenInd = -1;
+ if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
+ degenInd = 0;
+ else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
+ degenInd = uvVec.size() - 1;
+ else
+ continue;
+
+ // find another side sharing the degenerated shape
+ bool isPrev = ( degenInd == 0 );
+ if ( i >= QUAD_TOP_SIDE )
+ isPrev = !isPrev;
+ int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
+ StdMeshers_FaceSide* side2 = quad->side[ i2 ];
+ const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
+ int degenInd2 = -1;
+ if ( uvVec[ degenInd ].node == uvVec2[0].node )
+ degenInd2 = 0;
+ else if ( uvVec[ degenInd ].node == uvVec2.back().node )
+ degenInd2 = uvVec2.size() - 1;
+ else
+ throw SALOME_Exception( LOCALIZED( "Logical error" ));
- // move UV in the middle
- uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd ]);
- uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
- uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
- uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
- }
+ // move UV in the middle
+ uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd ]);
+ uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
+ uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
+ uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
+ }
+
+ else if ( quad->side.size() == 4 )
+
+ // Set number of nodes on a degenerated side to be same as on an opposite side
+ // ----------------------------------------------------------------------------
+ for ( unsigned i = 0; i < quad->side.size(); ++i )
+ {
+ StdMeshers_FaceSide* degSide = quad->side[i];
+ if ( !myHelper->IsDegenShape( degSide->EdgeID(0) ))
+ continue;
+ StdMeshers_FaceSide* oppSide = quad->side[( i+2 ) % quad->side.size() ];
+ if ( degSide->NbSegments() == oppSide->NbSegments() )
+ continue;
+
+ // make new side data
+ const vector<UVPtStruct>& uvVecDegOld = degSide->GetUVPtStruct();
+ const SMDS_MeshNode* n = uvVecDegOld[0].node;
+ Handle(Geom2d_Curve) c2d = degSide->Curve2d(0);
+ double f = degSide->FirstU(0), l = degSide->LastU(0);
+ gp_Pnt2d p1( uvVecDegOld.front().u, uvVecDegOld.front().v );
+ gp_Pnt2d p2( uvVecDegOld.back().u, uvVecDegOld.back().v );
+
+ delete degSide;
+ quad->side[i] = new StdMeshers_FaceSide( oppSide, n, &p1, &p2, c2d, f, l );
+ }
}
//================================================================================
*/
//================================================================================
-void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad)
+void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct::Ptr quad)
{
if ( !myNeedSmooth ) return;
typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
TNo2SmooNoMap smooNoMap;
- const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
- SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
- SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
- SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
+ const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
+ Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
+ double U1, U2, V1, V2;
+ surface->Bounds(U1, U2, V1, V2);
+ GeomAPI_ProjectPointOnSurf proj;
+ proj.Init( surface, U1, U2, V1, V2, BRep_Tool::Tolerance( geomFace ) );
+
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
+ SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
while ( nIt->more() ) // loop on nodes bound to a FACE
{
const SMDS_MeshNode* node = nIt->next();
TSmoothNode & sNode = smooNoMap[ node ];
- sNode._uv = myHelper->GetNodeUV( geomFace, node );
+ sNode._uv = myHelper->GetNodeUV( geomFace, node );
+ sNode._xyz = SMESH_TNodeXYZ( node );
// set sNode._triangles
SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
{
TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
+ sNode._xyz = SMESH_TNodeXYZ( uvVec[j].node );
}
}
if ( sNode._triangles.empty() )
continue; // not movable node
- // compute a new UV
- gp_XY newUV (0,0);
+ // compute a new XYZ
+ gp_XYZ newXYZ (0,0,0);
for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
- newUV += sNode._triangles[i]._n1->_uv;
- newUV /= sNode._triangles.size();
+ newXYZ += sNode._triangles[i]._n1->_xyz;
+ newXYZ /= sNode._triangles.size();
- // check validity of the newUV
- bool isValid = true;
- for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
- isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
+ // compute a new UV by projection
+ gp_XY newUV;
+ proj.Perform( newXYZ );
+ bool isValid = ( proj.IsDone() && proj.NbPoints() > 0 );
+ if ( isValid )
+ {
+ // check validity of the newUV
+ Quantity_Parameter u,v;
+ proj.LowerDistanceParameters( u, v );
+ newUV.SetCoord( u, v );
+ for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
+ isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
+ }
+ if ( !isValid )
+ {
+ // compute a new UV by averaging
+ newUV.SetCoord(0.,0.);
+ for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
+ newUV += sNode._triangles[i]._n1->_uv;
+ newUV /= sNode._triangles.size();
+ // check validity of the newUV
+ isValid = true;
+ for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
+ isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
+ }
if ( isValid )
+ {
sNode._uv = newUV;
+ sNode._xyz = surface->Value( newUV.X(), newUV.Y() ).XYZ();
+ }
}
}
// Set new XYZ to the smoothed nodes
- Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
-
for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
{
TSmoothNode& sNode = n2sn->second;
}
}
}
+
+/*//================================================================================
+/*!
+ * \brief Finds vertices at the most sharp face corners
+ * \param [in] theFace - the FACE
+ * \param [in,out] theWire - the ordered edges of the face. It can be modified to
+ * have the first VERTEX of the first EDGE in \a vertices
+ * \param [out] theVertices - the found corner vertices in the order corresponding to
+ * the order of EDGEs in \a theWire
+ * \param [out] theNbDegenEdges - nb of degenerated EDGEs in theFace
+ * \return int - number of quad sides found: 0, 3 or 4
+ */
+//================================================================================
+
+int StdMeshers_Quadrangle_2D::GetCorners(const TopoDS_Face& theFace,
+ SMESH_Mesh & theMesh,
+ std::list<TopoDS_Edge>& theWire,
+ std::vector<TopoDS_Vertex>& theVertices,
+ int & theNbDegenEdges)
+{
+ theNbDegenEdges = 0;
+
+ SMESH_MesherHelper helper( theMesh );
+
+ // sort theVertices by angle
+ multimap<double, TopoDS_Vertex> vertexByAngle;
+ TopTools_DataMapOfShapeReal angleByVertex;
+ TopoDS_Edge prevE = theWire.back();
+ if ( SMESH_Algo::isDegenerated( prevE ))
+ {
+ list<TopoDS_Edge>::reverse_iterator edge = ++theWire.rbegin();
+ while ( SMESH_Algo::isDegenerated( *edge ))
+ ++edge;
+ if ( edge == theWire.rend() )
+ return false;
+ prevE = *edge;
+ }
+ list<TopoDS_Edge>::iterator edge = theWire.begin();
+ for ( ; edge != theWire.end(); ++edge )
+ {
+ if ( SMESH_Algo::isDegenerated( *edge ))
+ {
+ ++theNbDegenEdges;
+ continue;
+ }
+ TopoDS_Vertex v = helper.IthVertex( 0, *edge );
+ if ( SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
+ {
+ double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace );
+ vertexByAngle.insert( make_pair( angle, v ));
+ angleByVertex.Bind( v, angle );
+ }
+ prevE = *edge;
+ }
+
+ // find out required nb of corners (3 or 4)
+ int nbCorners = 4;
+ TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID );
+ if ( !triaVertex.IsNull() &&
+ triaVertex.ShapeType() == TopAbs_VERTEX &&
+ helper.IsSubShape( triaVertex, theFace ))
+ nbCorners = 3;
+ else
+ triaVertex.Nullify();
+
+ // check nb of available corners
+ if ( nbCorners == 3 )
+ {
+ if ( vertexByAngle.size() < 3 )
+ return error(COMPERR_BAD_SHAPE,
+ TComm("Face must have 3 sides but not ") << vertexByAngle.size() );
+ }
+ else
+ {
+ if ( vertexByAngle.size() == 3 && theNbDegenEdges == 0 )
+ {
+ if ( myTriaVertexID < 1 )
+ return error(COMPERR_BAD_PARMETERS,
+ "No Base vertex provided for a trilateral geometrical face");
+
+ TComm comment("Invalid Base vertex: ");
+ comment << myTriaVertexID << " its ID is not among [ ";
+ multimap<double, TopoDS_Vertex>::iterator a2v = vertexByAngle.begin();
+ comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++;
+ comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << ", "; a2v++;
+ comment << helper.GetMeshDS()->ShapeToIndex( a2v->second ) << " ]";
+ return error(COMPERR_BAD_PARMETERS, comment );
+ }
+ if ( vertexByAngle.size() + ( theNbDegenEdges > 0 ) < 4 &&
+ vertexByAngle.size() + theNbDegenEdges != 4 )
+ return error(COMPERR_BAD_SHAPE,
+ TComm("Face must have 4 sides but not ") << vertexByAngle.size() );
+ }
+
+ // put all corner vertices in a map
+ TopTools_MapOfShape vMap;
+ if ( nbCorners == 3 )
+ vMap.Add( triaVertex );
+ multimap<double, TopoDS_Vertex>::reverse_iterator a2v = vertexByAngle.rbegin();
+ for ( ; a2v != vertexByAngle.rend() && vMap.Extent() < nbCorners; ++a2v )
+ vMap.Add( (*a2v).second );
+
+ // check if there are possible variations in choosing corners
+ bool isThereVariants = false;
+ if ( vertexByAngle.size() > nbCorners )
+ {
+ double lostAngle = a2v->first;
+ double lastAngle = ( --a2v, a2v->first );
+ isThereVariants = ( lostAngle * 1.1 >= lastAngle );
+ }
+
+ // make theWire begin from a corner vertex or triaVertex
+ if ( nbCorners == 3 )
+ while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) ||
+ SMESH_Algo::isDegenerated( theWire.front() ))
+ theWire.splice( theWire.end(), theWire, theWire.begin() );
+ else
+ while ( !vMap.Contains( helper.IthVertex( 0, theWire.front() )) ||
+ SMESH_Algo::isDegenerated( theWire.front() ))
+ theWire.splice( theWire.end(), theWire, theWire.begin() );
+
+ // fill the result vector and prepare for its refinement
+ theVertices.clear();
+ vector< double > angles;
+ vector< TopoDS_Edge > edgeVec;
+ vector< int > cornerInd;
+ angles.reserve( vertexByAngle.size() );
+ edgeVec.reserve( vertexByAngle.size() );
+ cornerInd.reserve( nbCorners );
+ for ( edge = theWire.begin(); edge != theWire.end(); ++edge )
+ {
+ if ( SMESH_Algo::isDegenerated( *edge ))
+ continue;
+ TopoDS_Vertex v = helper.IthVertex( 0, *edge );
+ bool isCorner = vMap.Contains( v );
+ if ( isCorner )
+ {
+ theVertices.push_back( v );
+ cornerInd.push_back( angles.size() );
+ }
+ angles.push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI );
+ edgeVec.push_back( *edge );
+ }
+
+ // refine the result vector - make sides elual by length if
+ // there are several equal angles
+ if ( isThereVariants )
+ {
+ if ( nbCorners == 3 )
+ angles[0] = 2 * M_PI; // not to move the base triangle VERTEX
+
+ set< int > refinedCorners;
+ for ( size_t iC = 0; iC < cornerInd.size(); ++iC )
+ {
+ int iV = cornerInd[iC];
+ if ( !refinedCorners.insert( iV ).second )
+ continue;
+ list< int > equalVertices;
+ equalVertices.push_back( iV );
+ int nbC[2] = { 0, 0 };
+ // find equal angles backward and forward from the iV-th corner vertex
+ for ( int isFwd = 0; isFwd < 2; ++isFwd )
+ {
+ int dV = isFwd ? +1 : -1;
+ int iCNext = helper.WrapIndex( iC + dV, cornerInd.size() );
+ int iVNext = helper.WrapIndex( iV + dV, angles.size() );
+ while ( iVNext != iV )
+ {
+ bool equal = Abs( angles[iV] - angles[iVNext] ) < 0.1 * angles[iV];
+ if ( equal )
+ equalVertices.insert( isFwd ? equalVertices.end() : equalVertices.begin(), iVNext );
+ if ( iVNext == cornerInd[ iCNext ])
+ {
+ if ( !equal )
+ break;
+ nbC[ isFwd ]++;
+ refinedCorners.insert( cornerInd[ iCNext ] );
+ iCNext = helper.WrapIndex( iCNext + dV, cornerInd.size() );
+ }
+ iVNext = helper.WrapIndex( iVNext + dV, angles.size() );
+ }
+ }
+ // move corners to make sides equal by length
+ int nbEqualV = equalVertices.size();
+ int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] );
+ if ( nbExcessV > 0 )
+ {
+ // calculate normalized length of each side enclosed between neighbor equalVertices
+ vector< double > curLengths;
+ double totalLen = 0;
+ vector< int > evVec( equalVertices.begin(), equalVertices.end() );
+ int iEV = 0;
+ int iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )];
+ int iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )];
+ while ( curLengths.size() < nbEqualV + 1 )
+ {
+ curLengths.push_back( totalLen );
+ do {
+ curLengths.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]);
+ iE = helper.WrapIndex( iE + 1, edgeVec.size());
+ if ( iEV < evVec.size() && iE == evVec[ iEV++ ] )
+ break;
+ }
+ while( iE != iEEnd );
+ totalLen = curLengths.back();
+ }
+ curLengths.resize( equalVertices.size() );
+ for ( size_t iS = 0; iS < curLengths.size(); ++iS )
+ curLengths[ iS ] /= totalLen;
+
+ // find equalVertices most close to the ideal sub-division of all sides
+ int iBestEV = 0;
+ int iCorner = helper.WrapIndex( iC - nbC[0], cornerInd.size() );
+ int nbSides = 2 + nbC[0] + nbC[1];
+ for ( int iS = 1; iS < nbSides; ++iS, ++iBestEV )
+ {
+ double idealLen = iS / double( nbSides );
+ double d, bestDist = 1.;
+ for ( iEV = iBestEV; iEV < curLengths.size(); ++iEV )
+ if (( d = Abs( idealLen - curLengths[ iEV ])) < bestDist )
+ {
+ bestDist = d;
+ iBestEV = iEV;
+ }
+ if ( iBestEV > iS-1 + nbExcessV )
+ iBestEV = iS-1 + nbExcessV;
+ theVertices[ iCorner ] = helper.IthVertex( 0, edgeVec[ evVec[ iBestEV ]]);
+ iCorner = helper.WrapIndex( iCorner + 1, cornerInd.size() );
+ }
+ }
+ }
+ }
+
+ return nbCorners;
+}
