myTrianglePreference(false),
myTriaVertexID(-1),
myNeedSmooth(false),
+ myCheckOri(false),
myParams( NULL ),
myQuadType(QUAD_STANDARD),
myHelper( NULL )
_quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
myNeedSmooth = false;
+ myCheckOri = false;
FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true );
if (!quad)
if ( res == COMPUTE_OK && myNeedSmooth )
smooth( quad );
+ if ( res == COMPUTE_OK )
+ res = check();
+
return ( res == COMPUTE_OK );
}
FaceQuadStruct::Ptr newQuad = myQuadList.back();
if ( quad != newQuad ) // split done
{
- {
+ { // update left side limit till where to make triangles
FaceQuadStruct::Ptr botQuad = // a bottom part
( quad->side[ QUAD_LEFT_SIDE ].from == 0 ) ? quad : newQuad;
if ( botQuad->nbNodeOut( QUAD_LEFT_SIDE ) > 0 )
{
splitQuad( quad, quad->iSize-2, 0 );
}
- if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
+ if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
{
splitQuad( quad, 1, 0 );
+
+ if ( quad->nbNodeOut( QUAD_TOP_SIDE ))
+ {
+ newQuad = myQuadList.back();
+ if ( newQuad == quad ) // too narrow to split
+ {
+ // update left side limit till where to make triangles
+ quad->side[ QUAD_LEFT_SIDE ].to--;
+ }
+ else
+ {
+ FaceQuadStruct::Ptr leftQuad =
+ ( quad->side[ QUAD_BOTTOM_SIDE ].from == 0 ) ? quad : newQuad;
+ leftQuad->nbNodeOut( QUAD_TOP_SIDE ) = 0;
+ }
+ }
}
- return computeQuadDominant( aMesh, aFace );
+ if ( ! computeQuadDominant( aMesh, aFace ))
+ return false;
+
+ // try to fix zero-area triangles near straight-angle corners
+
+ return true;
}
//================================================================================
iup = nbhoriz - 1;
int stop = 0;
- // if left edge is out, we will stop at a second node
- //if (quad->nbNodeOut(3)) stop++;
- if ( quad->nbNodeOut( QUAD_RIGHT_SIDE ))
- quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node;
- if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
- quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node;
+ if ( quad->side[3].grid->Edge(0).IsNull() ) // left side is simulated one
+ {
+ // quad divided at I but not at J, as nbvertic==nbright==2
+ stop++; // we stop at a second node
+ }
+ else
+ {
+ if ( quad->nbNodeOut( QUAD_RIGHT_SIDE ))
+ quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node;
+ if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
+ quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node;
+ if ( nbright > 2 ) // there was a split at J
+ quad->nbNodeOut( QUAD_LEFT_SIDE ) = 0;
+ }
+ const SMDS_MeshNode *a, *b, *c, *d;
+ i = nbup - 1;
+ // avoid creating zero-area triangles near a straight-angle corner
+ {
+ a = uv_e2[i].node;
+ b = uv_e2[i-1].node;
+ c = uv_e1[nbright-2].node;
+ SMESH_TNodeXYZ pa( a ), pb( b ), pc( c );
+ double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus();
+ if ( Abs( area ) < 1e-20 )
+ {
+ --g;
+ d = quad->UVPt( g, nbvertic-2 ).node;
+ if ( myTrianglePreference )
+ {
+ if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c))
+ meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else
+ {
+ if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c))
+ meshDS->SetMeshElementOnShape(face, geomFaceID);
+ --i;
+ }
+ }
+ }
// for each node of the up edge find nearest node
// in the first row of the regular grid and link them
- for (i = nbup - 1; i > stop; i--) {
- const SMDS_MeshNode *a, *b, *c, *d;
+ for ( ; i > stop; i--) {
a = uv_e2[i].node;
b = uv_e2[i - 1].node;
gp_Pnt pb (b->X(), b->Y(), b->Z());
int g = nbvertic - 1; // last processed node in the grid
int stop = 0;
i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1;
- for (; i > stop; i--) {
- const SMDS_MeshNode *a, *b, *c, *d;
+
+ const SMDS_MeshNode *a, *b, *c, *d;
+ // avoid creating zero-area triangles near a straight-angle corner
+ {
+ a = uv_e3[i].node;
+ b = uv_e3[i-1].node;
+ c = quad->UVPt( 1, g ).node;
+ SMESH_TNodeXYZ pa( a ), pb( b ), pc( c );
+ double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus();
+ if ( Abs( area ) < 1e-20 )
+ {
+ --g;
+ d = quad->UVPt( 1, g ).node;
+ if ( myTrianglePreference )
+ {
+ if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c))
+ meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else
+ {
+ if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c))
+ meshDS->SetMeshElementOnShape(face, geomFaceID);
+ --i;
+ }
+ }
+ }
+ for (; i > stop; i--) // loop on nodes on the left side
+ {
a = uv_e3[i].node;
b = uv_e3[i - 1].node;
gp_Pnt pb (b->X(), b->Y(), b->Z());
if (i == stop + 1) { // down bondary reached
c = quad->uv_grid[nbhoriz*jlow + 1].node;
near = jlow;
- } else {
+ }
+ else {
double mind = RealLast();
for (int k = g; k >= jlow; k--) {
const SMDS_MeshNode *nk;
if (k > jup)
- nk = uv_e2[1].node;
+ nk = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
else
nk = quad->uv_grid[nbhoriz*k + 1].node;
gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
}
else { // make quadrangle
if (near + 1 > jup)
- d = uv_e2[1].node;
+ d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
else
d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
- //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
- if (!myTrianglePreference){
+ if (!myTrianglePreference) {
SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
for (int k = near + 1; k < g; k++) {
c = quad->uv_grid[nbhoriz*k + 1].node;
if (k + 1 > jup)
- d = uv_e2[1].node;
+ d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
else
d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
return true;
}
+//================================================================================
+/*!
+ * \brief Return true if the algorithm can mesh this shape
+ * \param [in] aShape - shape to check
+ * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+ * else, returns OK if at least one shape is OK
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
+{
+ int nbFoundFaces = 0;
+ for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
+ {
+ const TopoDS_Shape& aFace = exp.Current();
+ int nbWire = SMESH_MesherHelper::Count( aFace, TopAbs_WIRE, false );
+ if ( nbWire != 1 ) {
+ if ( toCheckAll ) return false;
+ continue;
+ }
+
+ int nbNoDegenEdges = 0;
+ TopExp_Explorer eExp( aFace, TopAbs_EDGE );
+ for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next() ) {
+ if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() )))
+ ++nbNoDegenEdges;
+ }
+ if ( toCheckAll && nbNoDegenEdges < 3 ) return false;
+ if ( !toCheckAll && nbNoDegenEdges >= 3 ) return true;
+ }
+ return ( toCheckAll && nbFoundFaces != 0 );
+}
//================================================================================
/*!
if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
const bool ignoreMediumNodes = _quadraticMesh;
- // verify 1 wire only, with 4 edges
+ // verify 1 wire only
list< TopoDS_Edge > edges;
list< int > nbEdgesInWire;
int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
double d = v1 ^ v2;
return d > 1e-100;
}
+ //================================================================================
+ /*!
+ * \brief Returns area of a triangle
+ */
+ //================================================================================
+
+ double getArea( const gp_UV uv1, const gp_UV uv2, const gp_UV uv3 )
+ {
+ gp_XY v1 = uv1 - uv2, v2 = uv3 - uv2;
+ double a = v2 ^ v1;
+ return a;
+ }
}
//================================================================================
}
}
+//================================================================================
+/*!
+ * \brief Checks validity of generated faces
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::check()
+{
+ const bool isOK = true;
+ if ( !myCheckOri || myQuadList.empty() || !myQuadList.front() || !myHelper )
+ return isOK;
+
+ TopoDS_Face geomFace = TopoDS::Face( myHelper->GetSubShape() );
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
+ bool toCheckUV;
+ if ( geomFace.Orientation() >= TopAbs_INTERNAL ) geomFace.Orientation( TopAbs_FORWARD );
+
+ // Get a reference orientation sign
+
+ double okSign;
+ {
+ TError err;
+ TSideVector wireVec =
+ StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err );
+ StdMeshers_FaceSidePtr wire = wireVec[0];
+
+ // find a right angle VERTEX
+ int iVertex;
+ double maxAngle = -1e100;
+ for ( int i = 0; i < wire->NbEdges(); ++i )
+ {
+ int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() );
+ const TopoDS_Edge& e1 = wire->Edge( iPrev );
+ const TopoDS_Edge& e2 = wire->Edge( i );
+ double angle = myHelper->GetAngle( e1, e2, geomFace );
+ if ( maxAngle < angle && angle < 0.9 * M_PI )
+ {
+ maxAngle = angle;
+ iVertex = i;
+ }
+ }
+ if ( maxAngle < -2*M_PI ) return isOK;
+
+ // get a sign of 2D area of a corner face
+
+ int iPrev = myHelper->WrapIndex( iVertex-1, wire->NbEdges() );
+ const TopoDS_Edge& e1 = wire->Edge( iPrev );
+ const TopoDS_Edge& e2 = wire->Edge( iVertex );
+
+ gp_Vec2d v1, v2; gp_Pnt2d p;
+ double u[2];
+ {
+ bool rev = ( e1.Orientation() == TopAbs_REVERSED );
+ Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e1, geomFace, u[0], u[1] );
+ c->D1( u[ !rev ], p, v1 );
+ if ( !rev )
+ v1.Reverse();
+ }
+ {
+ bool rev = ( e2.Orientation() == TopAbs_REVERSED );
+ Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e2, geomFace, u[0], u[1] );
+ c->D1( u[ rev ], p, v2 );
+ if ( rev )
+ v2.Reverse();
+ }
+
+ okSign = v2 ^ v1;
+
+ if ( maxAngle < 0 )
+ okSign *= -1;
+ }
+
+ // Look for incorrectly oriented faces
+
+ std::list<const SMDS_MeshElement*> badFaces;
+
+ const SMDS_MeshNode* nn [ 8 ]; // 8 is just for safety
+ gp_UV uv [ 8 ];
+ SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements();
+ while ( fIt->more() ) // loop on faces bound to a FACE
+ {
+ const SMDS_MeshElement* f = fIt->next();
+
+ const int nbN = f->NbCornerNodes();
+ for ( int i = 0; i < nbN; ++i )
+ nn[ i ] = f->GetNode( i );
+
+ const SMDS_MeshNode* nInFace = 0;
+ if ( myHelper->HasSeam() )
+ for ( int i = 0; i < nbN && !nInFace; ++i )
+ if ( !myHelper->IsSeamShape( nn[i]->getshapeId() ))
+ nInFace = nn[i];
+
+ for ( int i = 0; i < nbN; ++i )
+ uv[ i ] = myHelper->GetNodeUV( geomFace, nn[i], nInFace, &toCheckUV );
+
+ switch ( nbN ) {
+ case 4:
+ {
+ double sign1 = getArea( uv[0], uv[1], uv[2] );
+ double sign2 = getArea( uv[0], uv[2], uv[3] );
+ if ( sign1 * sign2 < 0 )
+ {
+ sign2 = getArea( uv[1], uv[2], uv[3] );
+ sign1 = getArea( uv[1], uv[3], uv[0] );
+ if ( sign1 * sign2 < 0 )
+ continue; // this should not happen
+ }
+ if ( sign1 * okSign < 0 )
+ badFaces.push_back ( f );
+ break;
+ }
+ case 3:
+ {
+ double sign = getArea( uv[0], uv[1], uv[2] );
+ if ( sign * okSign < 0 )
+ badFaces.push_back ( f );
+ break;
+ }
+ default:;
+ }
+ }
+
+ if ( !badFaces.empty() )
+ {
+ SMESH_subMesh* fSM = myHelper->GetMesh()->GetSubMesh( geomFace );
+ SMESH_ComputeErrorPtr& err = fSM->GetComputeError();
+ err.reset ( new SMESH_ComputeError( COMPERR_ALGO_FAILED,
+ "Inverted elements generated"));
+ err->myBadElements.swap( badFaces );
+
+ return !isOK;
+ }
+
+ return isOK;
+}
+
/*//================================================================================
/*!
* \brief Finds vertices at the most sharp face corners
isThereVariants = ( lostAngle * 1.1 >= lastAngle );
}
+ myCheckOri = ( vertexByAngle.size() > nbCorners ||
+ vertexByAngle.begin()->first < 5.* M_PI/180 );
+
// make theWire begin from a corner vertex or triaVertex
if ( nbCorners == 3 )
while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) ||
myQuadList.push_back( FaceQuadStruct::Ptr( newQuad ));
vector<UVPtStruct> points;
- if ( I > 0 )
+ if ( I > 0 && I <= quad->iSize-2 )
{
points.reserve( quad->jSize );
for ( int jP = 0; jP < quad->jSize; ++jP )
return QUAD_LEFT_SIDE;
}
- else if ( J > 0 ) //// split horizontally, a new quad is below an old one
+ else if ( J > 0 && J <= quad->jSize-2 ) //// split horizontally, a new quad is below an old one
{
points.reserve( quad->iSize );
for ( int iP = 0; iP < quad->iSize; ++iP )
return QUAD_TOP_SIDE;
}
+
+ myQuadList.pop_back();
+ return -1;
}
//================================================================================
oSide->contacts[iOC].other_side = this;
}
}
+ return *this;
}
//================================================================================
grid->Reverse();
}
}
+ return (bool)grid;
}
//================================================================================