#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>
//=============================================================================
/*!
- *
+ *
*/
//=============================================================================
//=============================================================================
/*!
- *
+ *
*/
//=============================================================================
const TopoDS_Shape& aShape,
SMESH_Hypothesis::Hypothesis_Status& aStatus)
{
+ myTriaVertexID = -1;
+ myQuadType = QUAD_STANDARD;
+ myQuadranglePreference = false;
+ myTrianglePreference = false;
+ myQuadStruct.reset();
+ myHelper = NULL;
+
bool isOk = true;
aStatus = SMESH_Hypothesis::HYP_OK;
GetUsedHypothesis(aMesh, aShape, false);
const SMESHDS_Hypothesis * aHyp = 0;
- myTriaVertexID = -1;
- myQuadType = QUAD_STANDARD;
- myQuadranglePreference = false;
- myTrianglePreference = false;
- myQuadStruct.reset();
-
bool isFirstParams = true;
// First assigned hypothesis (if any) is processed now
const TopoDS_Shape& aShape)
{
const TopoDS_Face& F = TopoDS::Face(aShape);
- Handle(Geom_Surface) S = BRep_Tool::Surface(F);
-
- SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
- aMesh.GetSubMesh(aShape);
+ aMesh.GetSubMesh( F );
SMESH_MesherHelper helper (aMesh);
myHelper = &helper;
_quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
myNeedSmooth = false;
- FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape);
+ FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true );
if (!quad)
return false;
myQuadStruct = quad;
- if (myQuadranglePreference) {
- int n1 = quad->side[0]->NbPoints();
- int n2 = quad->side[1]->NbPoints();
- int n3 = quad->side[2]->NbPoints();
- int n4 = quad->side[3]->NbPoints();
+ updateDegenUV( quad );
+
+ enum { NOT_COMPUTED = -1, COMPUTE_FAILED = 0, COMPUTE_OK = 1 };
+ int res = NOT_COMPUTED;
+ if (myQuadranglePreference)
+ {
+ int n1 = quad->side[0]->NbPoints();
+ int n2 = quad->side[1]->NbPoints();
+ int n3 = quad->side[2]->NbPoints();
+ int n4 = quad->side[3]->NbPoints();
int nfull = n1+n2+n3+n4;
- int ntmp = nfull/2;
+ int ntmp = nfull/2;
ntmp = ntmp*2;
- if (nfull == ntmp && ((n1 != n3) || (n2 != n4))) {
- // special path for using only quandrangle faces
- bool ok = ComputeQuadPref(aMesh, aShape, quad);
- if ( ok && myNeedSmooth )
- Smooth( quad );
- return ok;
- }
- }
- else if (myQuadType == QUAD_REDUCED) {
- int n1 = quad->side[0]->NbPoints();
- int n2 = quad->side[1]->NbPoints();
- int n3 = quad->side[2]->NbPoints();
- int n4 = quad->side[3]->NbPoints();
- int n13 = n1 - n3;
- int n24 = n2 - n4;
+ if (nfull == ntmp && ((n1 != n3) || (n2 != n4)))
+ {
+ // special path genarating only quandrangle faces
+ res = computeQuadPref( aMesh, F, quad );
+ }
+ }
+ else if (myQuadType == QUAD_REDUCED)
+ {
+ int n1 = quad->side[0]->NbPoints();
+ int n2 = quad->side[1]->NbPoints();
+ int n3 = quad->side[2]->NbPoints();
+ int n4 = quad->side[3]->NbPoints();
+ int n13 = n1 - n3;
+ int n24 = n2 - n4;
int n13tmp = n13/2; n13tmp = n13tmp*2;
int n24tmp = n24/2; n24tmp = n24tmp*2;
if ((n1 == n3 && n2 != n4 && n24tmp == n24) ||
- (n2 == n4 && n1 != n3 && n13tmp == n13)) {
- bool ok = ComputeReduced(aMesh, aShape, quad);
- if ( ok && myNeedSmooth )
- Smooth( quad );
- return ok;
+ (n2 == n4 && n1 != n3 && n13tmp == n13))
+ {
+ res = computeReduced( aMesh, F, quad );
}
+ else
+ {
+ 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.");
+ }
+ }
+
+ if ( res == NOT_COMPUTED )
+ {
+ res = computeQuadDominant( aMesh, F, quad );
}
+ if ( res == COMPUTE_OK && myNeedSmooth )
+ smooth( quad );
+
+ return ( res == COMPUTE_OK );
+}
+
+//================================================================================
+/*!
+ * \brief Compute quadrangles and possibly triangles
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh& aMesh,
+ const TopoDS_Face& aFace,
+ FaceQuadStruct::Ptr quad)
+{
// set normalized grid on unit square in parametric domain
-
- if (!SetNormalizedGrid(aMesh, aShape, quad))
+
+ if (!setNormalizedGrid(aMesh, aFace, quad))
return false;
// --- compute 3D values on points, store points & quadrangles
int nbvertic = Min(nbright, nbleft);
// internal mesh nodes
- int i, j, geomFaceID = meshDS->ShapeToIndex(F);
+ SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
+ Handle(Geom_Surface) S = BRep_Tool::Surface(aFace);
+ int i, j, geomFaceID = meshDS->ShapeToIndex(aFace);
for (i = 1; i < nbhoriz - 1; i++) {
for (j = 1; j < nbvertic - 1; j++) {
int ij = j * nbhoriz + i;
if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else {
- SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ splitQuad(meshDS, geomFaceID, a, b, c, d);
}
// if node d is not at position g - make additional triangles
if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else {
- SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ splitQuad(meshDS, geomFaceID, a, b, c, d);
}
if (near + 1 < g) { // if d not is at g - make additional triangles
if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else {
- SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ splitQuad(meshDS, geomFaceID, a, b, c, d);
}
if (near - 1 > g) { // if d not is at g - make additional triangles
if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else {
- SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ splitQuad(meshDS, geomFaceID, a, b, c, d);
}
if (near + 1 < g) { // if d not is at g - make additional triangles
}
}
- if ( myNeedSmooth )
- Smooth( quad );
-
bool isOk = true;
return isOk;
}
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
- const TopoDS_Shape& aShape,
- MapShapeNbElems& aResMap)
+bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aFace,
+ MapShapeNbElems& aResMap)
{
- aMesh.GetSubMesh(aShape);
+ aMesh.GetSubMesh(aFace);
std::vector<int> aNbNodes(4);
bool IsQuadratic = false;
- if (!CheckNbEdgesForEvaluate(aMesh, aShape, aResMap, aNbNodes, IsQuadratic)) {
+ if (!checkNbEdgesForEvaluate(aMesh, aFace, aResMap, aNbNodes, IsQuadratic)) {
std::vector<int> aResVec(SMDSEntity_Last);
for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
- SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aFace);
aResMap.insert(std::make_pair(sm,aResVec));
SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
ntmp = ntmp*2;
if (nfull==ntmp && ((n1!=n3) || (n2!=n4))) {
// special path for using only quandrangle faces
- return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
+ return evaluateQuadPref(aMesh, aFace, aNbNodes, aResMap, IsQuadratic);
//return true;
}
}
aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
}
}
- SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aFace);
aResMap.insert(std::make_pair(sm,aVec));
return true;
//=============================================================================
FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
- const TopoDS_Shape & aShape)
+ const TopoDS_Shape & aShape,
+ const bool considerMesh)
{
if ( myQuadStruct && myQuadStruct->face.IsSame( aShape ))
return myQuadStruct;
error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
return FaceQuadStruct::Ptr();
}
+
+ // find corner vertices of the quad
+ vector<TopoDS_Vertex> corners;
+ int nbDegenEdges, nbSides = getCorners( F, aMesh, edges, corners, nbDegenEdges, considerMesh );
+ 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, myProxyMesh));
- quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true,
- ignoreMediumNodes, myProxyMesh));
- quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,
- ignoreMediumNodes, myProxyMesh));
- 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);
- quad.reset();
- return quad;
- }
- 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 < QUAD_TOP_SIDE,
- ignoreMediumNodes, myProxyMesh));
- }
- 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 < QUAD_TOP_SIDE,
- ignoreMediumNodes, myProxyMesh));
- ++nbSides;
- }
- if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
- {
- myNeedSmooth = true;
- for ( unsigned i = QUAD_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 = QUAD_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)
- {
- quad.reset( new FaceQuadStruct );
- quad->side.reserve(nbEdgesInWire.front());
- nbSides = 0;
-
- SMESH_Block::GetOrderedEdges (F, 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 < QUAD_TOP_SIDE,
+ else if ( !SMESH_Algo::isDegenerated( *edgeIt ) && // closed EDGE
+ myHelper->IthVertex( 0, *edgeIt ).IsSame( myHelper->IthVertex( 1, *edgeIt )))
+ {
+ quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt++, &aMesh, iSide < QUAD_TOP_SIDE,
ignoreMediumNodes, myProxyMesh));
- ++nbSides;
+ ++iSide;
+ }
+ if ( quad->side.size() == 4 )
+ break;
+ 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");
+ if ( quad && quad->side.size() != 4 )
+ {
+ error(TComm("Bug: ") << quad->side.size() << " sides found instead of 4");
+ quad.reset();
}
-#endif
- if (!nbSides)
- nbSides = nbEdgesInWire.front();
- error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
- quad.reset();
}
return quad;
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
+bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh& aMesh,
const TopoDS_Shape & aShape,
- MapShapeNbElems& aResMap,
- std::vector<int>& aNbNodes,
- bool& IsQuadratic)
+ MapShapeNbElems& aResMap,
+ std::vector<int>& aNbNodes,
+ bool& IsQuadratic)
{
const TopoDS_Face & F = TopoDS::Face(aShape);
if ( quad )
{
// set normalized grid on unit square in parametric domain
- if (!SetNormalizedGrid(aMesh, aShape, quad))
+ if ( ! setNormalizedGrid( aMesh, TopoDS::Face( aShape ), quad))
quad.reset();
}
return quad;
*/
//=============================================================================
-bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
- const TopoDS_Shape& aShape,
+bool StdMeshers_Quadrangle_2D::setNormalizedGrid (SMESH_Mesh & aMesh,
+ const TopoDS_Face& aFace,
FaceQuadStruct::Ptr & quad)
{
// Algorithme décrit dans "Génération automatique de maillages"
// traitement dans le domaine paramétrique 2d u,v
// transport - projection sur le carré unité
-// MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
-// const TopoDS_Face& F = TopoDS::Face(aShape);
-
- // 1 --- find orientation of the 4 edges, by test on extrema
-
// max min 0 x1 1
// |<----north-2-------^ a3 -------------> a2
// | | ^1 1^
// =down
//
- // 3 --- 2D normalized values on unit square [0..1][0..1]
-
int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
//return error("Can't find nodes on sides");
return error(COMPERR_BAD_INPUT_MESH);
- if ( myNeedSmooth )
- UpdateDegenUV( quad );
-
// 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])*/ {
+ {
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])*/ {
+ {
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])*/ {
- int i = 0;
+ {
+ const int i = 0;
for (int j = 0; j < nbvertic; j++) // left
uv_grid[ j * nbhoriz + i ] = uv_e3[j];
}
// normalized 2d parameters on grid
+
for (int i = 0; i < nbhoriz; i++) {
for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
// --- droite i cste : x = x0 + y(x1-x0)
- double x0 = uv_e0[i].normParam; // bas - sud
+ double x0 = uv_e0[i].normParam; // bas - sud
double x1 = uv_e2[i].normParam; // haut - nord
// --- droite j cste : y = y0 + x(y1-y0)
- double y0 = uv_e3[j].normParam; // gauche-ouest
+ double y0 = uv_e3[j].normParam; // gauche - ouest
double y1 = uv_e1[j].normParam; // droite - est
// --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
double y = y0 + x * (y1 - y0);
uv_grid[ij].x = x;
uv_grid[ij].y = y;
- //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
- //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
}
}
- // 4 --- projection on 2d domain (u,v)
+ // 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 );
//=======================================================================
//function : ShiftQuad
-//purpose : auxilary function for ComputeQuadPref
+//purpose : auxilary function for computeQuadPref
//=======================================================================
-static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num, bool)
+static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num)
{
quad->shift( num, /*ori=*/true );
}
//================================================================================
/*!
* \brief Rotate sides of a quad by nb
- * \param nb - number of rotation quartes
+ * \param nb - number of rotation quartes
* \param ori - to keep orientation of sides as in an unit quad or not
*/
//================================================================================
//=======================================================================
//function : calcUV
-//purpose : auxilary function for ComputeQuadPref
+//purpose : auxilary function for computeQuadPref
//=======================================================================
static gp_UV calcUV(double x0, double x1, double y0, double y1,
//=======================================================================
//function : calcUV2
-//purpose : auxilary function for ComputeQuadPref
+//purpose : auxilary function for computeQuadPref
//=======================================================================
static gp_UV calcUV2(double x, double y,
*/
//=======================================================================
-bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh,
- const TopoDS_Shape& aShape,
+bool StdMeshers_Quadrangle_2D::computeQuadPref (SMESH_Mesh & aMesh,
+ const TopoDS_Face& aFace,
FaceQuadStruct::Ptr quad)
{
// Auxilary key in order to keep old variant
// of meshing after implementation new variant
// for bug 0016220 from Mantis.
- bool OldVersion = false;
- if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
- OldVersion = true;
+ bool OldVersion = (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED);
- SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
- const TopoDS_Face& F = TopoDS::Face(aShape);
- Handle(Geom_Surface) S = BRep_Tool::Surface(F);
+ SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
+ Handle(Geom_Surface) S = BRep_Tool::Surface(aFace);
bool WisF = true;
- int i,j,geomFaceID = meshDS->ShapeToIndex(F);
+ int i,j,geomFaceID = meshDS->ShapeToIndex(aFace);
int nb = quad->side[0]->NbPoints();
int nr = quad->side[1]->NbPoints();
int dh = abs(nb-nt);
int dv = abs(nr-nl);
- if (dh>=dv) {
- if (nt>nb) {
- // it is a base case => not shift quad but me be replacement is need
- shiftQuad(quad,0,WisF);
- }
- else {
- // we have to shift quad on 2
- shiftQuad(quad,2,WisF);
- }
- }
- else {
- if (nr>nl) {
- // we have to shift quad on 1
- shiftQuad(quad,1,WisF);
- }
- else {
- // we have to shift quad on 3
- shiftQuad(quad,3,WisF);
- }
- }
+ // rotate sides to be as in the picture below and to have
+ // dh >= dv and nt > nb
+ if ( dh >= dv )
+ shiftQuad( quad, ( nt > nb ) ? 0 : 2 );
+ else
+ shiftQuad( quad, ( nr > nl ) ? 1 : 3 );
nb = quad->side[0]->NbPoints();
nr = quad->side[1]->NbPoints();
dh = abs(nb-nt);
dv = abs(nr-nl);
int nbh = Max(nb,nt);
- int nbv = Max(nr,nl);
+ int nbv = Max(nr,nl);
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 \ |
// left |/________\| rigth
// | |
- // | |
+ // | C |
// | |
// 0------------0
// 0 bottom 1
- if (dh>dv) {
+ if ( dh > dv ) {
addv = (dh-dv)/2;
- nbv = nbv + addv;
+ nbv = nbv + addv;
}
- else { // dv>=dh
+ else { // dv >= dh
addh = (dv-dh)/2;
- nbh = nbh + addh;
+ nbh = nbh + addh;
}
const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
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 );
+ if ( !OldVersion )
+ {
+ // dh/2, Min(nb,nt), dh - dh/2, dv
+ }
// arrays for normalized params
TColStd_SequenceOfReal npb, npr, npt, npl;
*/
//=======================================================================
-bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh,
+bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh & aMesh,
const TopoDS_Shape& aShape,
- std::vector<int>& aNbNodes,
- MapShapeNbElems& aResMap,
- bool IsQuadratic)
+ std::vector<int>& aNbNodes,
+ MapShapeNbElems& aResMap,
+ bool IsQuadratic)
{
// Auxilary key in order to keep old variant
// of meshing after implementation new variant
return true;
}
-
//=============================================================================
/*! Split quadrangle in to 2 triangles by smallest diagonal
*
*/
//=============================================================================
-void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
- int theFaceID,
+
+void StdMeshers_Quadrangle_2D::splitQuad(SMESHDS_Mesh * theMeshDS,
+ int theFaceID,
const SMDS_MeshNode* theNode1,
const SMDS_MeshNode* theNode2,
const SMDS_MeshNode* theNode3,
const SMDS_MeshNode* theNode4)
{
- gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
- gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
- gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
- gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
SMDS_MeshFace* face;
- if (a.Distance(c) > b.Distance(d)){
+ if ( SMESH_TNodeXYZ( theNode1 ).SquareDistance( theNode3 ) >
+ SMESH_TNodeXYZ( theNode2 ).SquareDistance( theNode4 ) )
+ {
face = myHelper->AddFace(theNode2, theNode4 , theNode1);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
face = myHelper->AddFace(theNode2, theNode3, theNode4);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
-
}
- else{
+ else
+ {
face = myHelper->AddFace(theNode1, theNode2 ,theNode3);
if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
face = myHelper->AddFace(theNode1, theNode3, theNode4);
* Implementation of Reduced algorithm (meshing with quadrangles only)
*/
//=======================================================================
-bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh,
- const TopoDS_Shape& aShape,
+
+bool StdMeshers_Quadrangle_2D::computeReduced (SMESH_Mesh & aMesh,
+ const TopoDS_Face& aFace,
FaceQuadStruct::Ptr quad)
{
- 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);
+ SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
+ Handle(Geom_Surface) S = BRep_Tool::Surface(aFace);
+ int i,j,geomFaceID = meshDS->ShapeToIndex(aFace);
- 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 (MultipleReduce) { // == computeQuadPref QUAD_QUADRANGLE_PREF_REVERSED
//==================================================
int dh = abs(nb-nt);
int dv = abs(nr-nl);
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);
}
else {
// we have to shift quad on 2
- shiftQuad(quad,2,true);
+ shiftQuad(quad,2);
}
}
else {
if (nr > nl) {
// we have to shift quad on 1
- shiftQuad(quad,1,true);
+ shiftQuad(quad,1);
}
else {
// we have to shift quad on 3
- shiftQuad(quad,3,true);
+ shiftQuad(quad,3);
}
}
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
TColStd_SequenceOfReal npb, npr, npt, npl;
for (j = 0; j < nb; j++) {
// 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++) {
// 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)
//=========================================================
}
else {
// we have to shift quad on 2
- shiftQuad(quad,2,true);
+ shiftQuad(quad,2);
}
}
else {
if (nl > nr) {
// we have to shift quad on 1
- shiftQuad(quad,1,true);
+ shiftQuad(quad,1);
}
else {
// we have to shift quad on 3
- shiftQuad(quad,3,true);
+ shiftQuad(quad,3);
}
}
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::Ptr 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 >= 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" ));
+ // 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 && myQuadType == QUAD_STANDARD)
+
+ // 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::Ptr 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 );
while ( fIt->more() )
{
const SMDS_MeshElement* face = fIt->next();
- const int nbN = face->NbCornerNodes();
- const int nInd = face->GetNodeIndex( node );
+ const int nbN = face->NbCornerNodes();
+ const int nInd = face->GetNodeIndex( node );
const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
{
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);
- for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
- newUV += sNode._triangles[i]._n1->_uv;
- newUV /= sNode._triangles.size();
+ gp_XY newUV;
+ bool isValid = false;
+ bool use3D = ( iLoop > 2 ); // 3 loops in 2D and 2, in 3D
- // 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 );
+ if ( use3D )
+ {
+ // compute a new XYZ
+ gp_XYZ newXYZ (0,0,0);
+ for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
+ newXYZ += sNode._triangles[i]._n1->_xyz;
+ newXYZ /= sNode._triangles.size();
+
+ // compute a new UV by projection
+ proj.Perform( newXYZ );
+ 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
+ * \param [in] theConsiderMesh - if \c true, only meshed VERTEXes are considered
+ * as possible corners
+ * \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,
+ const bool theConsiderMesh)
+{
+ 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 ( !theConsiderMesh || 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, nbSeg;
+ angles.reserve( vertexByAngle.size() );
+ edgeVec.reserve( vertexByAngle.size() );
+ nbSeg.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 );
+ if ( theConsiderMesh && isThereVariants )
+ {
+ if ( SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( *edge ))
+ nbSeg.push_back( sm->NbNodes() + 1 );
+ else
+ nbSeg.push_back( 0 );
+ }
+ }
+
+ // 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;
+}
