#include "StdMeshers_QuadrangleParams.hxx"
#include "StdMeshers_ViscousLayers2D.hxx"
+#include <BRepBndLib.hxx>
#include <BRepClass_FaceClassifier.hxx>
#include <BRep_Tool.hxx>
+#include <Bnd_Box.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <Geom_Surface.hxx>
#include <NCollection_DefineArray2.hxx>
const TopoDS_Face& F = TopoDS::Face(aShape);
aMesh.GetSubMesh( F );
+ // do not initialize my fields before this as StdMeshers_ViscousLayers2D
+ // can call Compute() recursively
+ SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+ if ( !proxyMesh )
+ return false;
+
+ myProxyMesh = proxyMesh;
+
SMESH_MesherHelper helper (aMesh);
myHelper = &helper;
- myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
- if ( !myProxyMesh )
- return false;
-
_quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
myNeedSmooth = false;
if ( !setNormalizedGrid( quad ))
return false;
- if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE ))
- {
- splitQuad( quad, 0, 1 );
- }
if ( quad->nbNodeOut( QUAD_TOP_SIDE ))
{
splitQuad( quad, 0, quad->jSize-2 );
}
+ if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE )) // this should not happen
+ {
+ splitQuad( quad, 0, 1 );
+ }
FaceQuadStruct::Ptr newQuad = myQuadList.back();
if ( quad != newQuad ) // split done
{
+ {
+ FaceQuadStruct::Ptr botQuad = // a bottom part
+ ( quad->side[ QUAD_LEFT_SIDE ].from == 0 ) ? quad : newQuad;
+ if ( botQuad->nbNodeOut( QUAD_LEFT_SIDE ) > 0 )
+ botQuad->side[ QUAD_LEFT_SIDE ].to += botQuad->nbNodeOut( QUAD_LEFT_SIDE );
+ else if ( botQuad->nbNodeOut( QUAD_RIGHT_SIDE ) > 0 )
+ botQuad->side[ QUAD_RIGHT_SIDE ].to += botQuad->nbNodeOut( QUAD_RIGHT_SIDE );
+ }
// make quad be a greatest one
if ( quad->side[ QUAD_LEFT_SIDE ].NbPoints() == 2 ||
quad->side[ QUAD_RIGHT_SIDE ].NbPoints() == 2 )
int nbright = (int) uv_e1.size();
int nbleft = (int) uv_e3.size();
- if (quad->nbNodeOut(0) && nbvertic == 2)
+ if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occure
{
// Down edge is out
//
}
// right or left boundary quadrangles
- if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2)
+ if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) // this should not occure
{
int g = 0; // last processed node in the grid
int stop = nbright - 1;
// MESSAGE("left edge is out");
int g = nbvertic - 1; // last processed node in the grid
int stop = 0;
- i = nbleft - 1;
- if (quad->side[3].from != stop ) stop++;
- if (quad->side[3].to != i ) i--;
+ i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1;
for (; i > stop; i--) {
const SMDS_MeshNode *a, *b, *c, *d;
a = uv_e3[i].node;
//purpose : auxilary function for computeQuadPref
//=======================================================================
-static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num)
+void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int num )
{
- quad->shift( num, /*ori=*/true );
+ quad->shift( num, /*ori=*/true, /*keepGrid=*/myQuadList.size() > 1 );
}
//================================================================================
* \brief Rotate sides of a quad by given nb of quartes
* \param nb - number of rotation quartes
* \param ori - to keep orientation of sides as in an unit quad or not
+ * \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to
+ * are altered instead
*/
//================================================================================
-void FaceQuadStruct::shift( size_t nb, bool ori )
+void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid )
{
if ( nb == 0 ) return;
bool wasForward = (i < QUAD_TOP_SIDE);
bool newForward = (id < QUAD_TOP_SIDE);
if ( wasForward != newForward )
- side[ i ].Reverse();
+ side[ i ].Reverse( keepGrid );
}
newSides[ id ] = side[ i ];
sidePtrs[ i ] = & side[ i ];
{
// rotate the quad to have nt > nb [and nr > nl]
if ( nb > nt )
- quad->shift( nr > nl ? 1 : 2, true );
+ shiftQuad ( quad, nr > nl ? 1 : 2 );
else if ( nr > nl )
- quad->shift( nb == nt ? 1 : 0, true );
+ shiftQuad( quad, nb == nt ? 1 : 0 );
else if ( nl > nr )
- quad->shift( 3, true );
+ shiftQuad( quad, 3 );
}
nb = quad->side[0].NbPoints();
// 0------------0
// 0 bottom 1
- const vector<UVPtStruct>& uv_eb_vec = quad->side[0].GetUVPtStruct(true,0);
- const vector<UVPtStruct>& uv_er_vec = quad->side[1].GetUVPtStruct(false,1);
- const vector<UVPtStruct>& uv_et_vec = quad->side[2].GetUVPtStruct(true,1);
- const vector<UVPtStruct>& uv_el_vec = quad->side[3].GetUVPtStruct(false,0);
+
const int bfrom = quad->side[0].from;
const int rfrom = quad->side[1].from;
const int tfrom = quad->side[2].from;
const int lfrom = quad->side[3].from;
- if (uv_eb_vec.size() < nb + bfrom ||
- uv_er_vec.size() < nr + rfrom ||
- uv_et_vec.size() < nt + tfrom ||
- uv_el_vec.size() < nl + lfrom)
- return error(COMPERR_BAD_INPUT_MESH);
-
- const UVPtStruct * uv_eb = & uv_eb_vec[0] + bfrom;
- const UVPtStruct * uv_er = & uv_er_vec[0] + rfrom;
- const UVPtStruct * uv_et = & uv_et_vec[0] + tfrom;
- const UVPtStruct * uv_el = & uv_el_vec[0] + lfrom;
+ {
+ const vector<UVPtStruct>& uv_eb_vec = quad->side[0].GetUVPtStruct(true,0);
+ const vector<UVPtStruct>& uv_er_vec = quad->side[1].GetUVPtStruct(false,1);
+ const vector<UVPtStruct>& uv_et_vec = quad->side[2].GetUVPtStruct(true,1);
+ const vector<UVPtStruct>& uv_el_vec = quad->side[3].GetUVPtStruct(false,0);
+ if (uv_eb_vec.empty() ||
+ uv_er_vec.empty() ||
+ uv_et_vec.empty() ||
+ uv_el_vec.empty())
+ return error(COMPERR_BAD_INPUT_MESH);
+ }
+ FaceQuadStruct::SideIterator uv_eb, uv_er, uv_et, uv_el;
+ uv_eb.Init( quad->side[0] );
+ uv_er.Init( quad->side[1] );
+ uv_et.Init( quad->side[2] );
+ uv_el.Init( quad->side[3] );
gp_UV a0,a1,a2,a3, p0,p1,p2,p3, uv;
double x,y;
if ( !myForcedPnts.empty() )
{
- if ( dv != 0 && dh != 0 )
+ if ( dv != 0 && dh != 0 ) // here myQuadList.size() == 1
{
const int dmin = Min( dv, dh );
x = uv_et[ dmin ].normParam;
p0 = quad->side[0].grid->Value2d( x ).XY();
p2 = uv_et[ dmin ].UV();
- for ( int i = 1; i < nl; ++i )
+ double y0 = uv_er[ dmin ].normParam;
+ for ( int i = 1; i < nl-1; ++i )
{
- y = uv_er[ i + dmin ].normParam;
- p1 = uv_er[ i + dmin ].UV();
+ y = y0 + i / ( nl-1. ) * ( 1. - y0 );
+ p1 = quad->side[1].grid->Value2d( y ).XY();
p3 = quad->side[3].grid->Value2d( y ).XY();
uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 );
pointsLCt[ i ].u = uv.X();
p1 = uv_er[ dmin ].UV();
p3 = quad->side[3].grid->Value2d( y ).XY();
double x0 = uv_et[ dmin ].normParam;
- double xn = uv_et[ dmin+nb-1 ].normParam;
- double kx = ( 1. - x0 ) / ( xn - x0 );
for ( int i = 1; i < nb-1; ++i )
{
- x = x0 + ( uv_et[ i + dmin ].normParam - x0 ) * kx;
+ x = x0 + i / ( nb-1. ) * ( 1. - x0 );
p2 = quad->side[2].grid->Value2d( x ).XY();
p0 = quad->side[0].grid->Value2d( x ).XY();
uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 );
} // if ( dv != 0 && dh != 0 )
- // Case dv == 0
+ const int db = quad->side[0].IsReversed() ? -1 : +1;
+ const int dr = quad->side[1].IsReversed() ? -1 : +1;
+ const int dt = quad->side[2].IsReversed() ? -1 : +1;
+ const int dl = quad->side[3].IsReversed() ? -1 : +1;
+
+ // Case dv == 0, here possibly myQuadList.size() > 1
//
// lw nb lw = dh/2
// +------------+
// +------------+
const int lw = dh/2; // lateral width
- const double lL = quad->side[3].Length();
- const double lLwL = quad->side[2].Length( tfrom, tfrom + lw );
- const double yCbL = lLwL / ( lLwL + lL );
-
- const double lR = quad->side[1].Length();
- const double lLwR = quad->side[2].Length( tfrom + lw + nb-1,
- tfrom + lw + nb-1 + lw );
- const double yCbR = lLwR / ( lLwR + lR );
-
+ double yCbL, yCbR;
+ {
+ double lL = quad->side[3].Length();
+ double lLwL = quad->side[2].Length( tfrom,
+ tfrom + ( lw ) * dt );
+ yCbL = lLwL / ( lLwL + lL );
+
+ double lR = quad->side[1].Length();
+ double lLwR = quad->side[2].Length( tfrom + ( lw + nb-1 ) * dt,
+ tfrom + ( lw + nb-1 + lw ) * dt);
+ yCbR = lLwR / ( lLwR + lR );
+ }
// Make sides separating domains Cb and L and R
StdMeshers_FaceSidePtr sideLCb, sideRCb;
UVPtStruct pTBL, pTBR; // points where 3 domains meat
x = quad->side[2].Param( i2 );
y = yCbR * i / lw;
+ p1 = quad->side[1].Value2d( y );
p0 = quad->side[0].Value2d( x );
p2 = uv_et[ i2 ].UV();
+ p3 = quad->side[3].Value2d( y );
uv = calcUV( x,y, a0,a1,a2,a3, p0,p1,p2,p3 );
pointsRCb[ i ].u = uv.X();
pointsRCb[ i ].v = uv.Y();
qL->side[1] = sideLCt;
qL->side[2] = quad->side[2];
qL->side[3] = quad->side[3];
- qL->side[2].to = lw + 1 + tfrom;
+ qL->side[2].to = ( lw + 1 ) * dt + tfrom;
// Make R quad
FaceQuadStruct* qR = new FaceQuadStruct( quad->face, "R" );
myQuadList.push_back( FaceQuadStruct::Ptr( qR ));
qR->side[0] = sideRCb;
qR->side[0].from = lw;
qR->side[0].to = -1;
+ qR->side[0].di = -1;
qR->side[1] = quad->side[1];
qR->side[2] = quad->side[2];
- qR->side[2].from = lw + nb-1 + tfrom;
+ qR->side[2].from = ( lw + nb-1 ) * dt + tfrom;
qR->side[3] = sideRCt;
// Make Ct from the main quad
FaceQuadStruct::Ptr qCt = quad;
qCt->side[0] = sideCbCt;
qCt->side[1] = sideRCt;
- qCt->side[2].from = lw + tfrom;
- qCt->side[2].to = lw + nb + tfrom;
+ qCt->side[2].from = ( lw ) * dt + tfrom;
+ qCt->side[2].to = ( lw + nb ) * dt + tfrom;
qCt->side[3] = sideLCt;
qCt->uv_grid.clear();
qCt->name = "Ct";
// Get nodes to smooth
+ // TODO: do not smooth fixed nodes
+
typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
TNo2SmooNoMap smooNoMap;
TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID );
if ( !triaVertex.IsNull() &&
triaVertex.ShapeType() == TopAbs_VERTEX &&
- helper.IsSubShape( triaVertex, theFace ))
+ helper.IsSubShape( triaVertex, theFace ) &&
+ ( vertexByAngle.size() != 4 || vertexByAngle.begin()->first < 5 * M_PI/180. ))
nbCorners = 3;
else
triaVertex.Nullify();
//================================================================================
FaceQuadStruct::Side::Side(StdMeshers_FaceSidePtr theGrid)
- : grid(theGrid), nbNodeOut(0), from(0), to(theGrid ? theGrid->NbPoints() : 0 )
+ : grid(theGrid), nbNodeOut(0), from(0), to(theGrid ? theGrid->NbPoints() : 0 ), di(1)
{
}
surf->Bounds( u1,u2,v1,v2 );
GeomAPI_ProjectPointOnSurf project;
project.Init(surf, u1,u2, v1,v2, tol );
+ Bnd_Box bbox;
+ BRepBndLib::Add( face, bbox );
+ double farTol = 0.01 * sqrt( bbox.SquareExtent() );
for ( size_t iP = 0; iP < points.size(); ++iP )
{
<< points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )");
continue;
}
- if ( project.LowerDistance() > tol*1000 )
+ if ( project.LowerDistance() > farTol )
{
if ( isStrictCheck && iP < nbPoints )
return error
for ( int dj = 0; dj < 2; ++dj )
{
double dist2 = ( myForcedPnts[ iFP ].uv - quad->UVPt( i+di,j+dj ).UV() ).SquareModulus();
- ijByDist.insert( make_pair( dist2, make_pair( i+di,j+dj )));
+ ijByDist.insert( make_pair( dist2, make_pair( di,dj )));
}
// try all nodes starting from the closest one
set< FaceQuadStruct::Ptr > changedQuads;
multimap< double, pair< int, int > >::iterator d2ij = ijByDist.begin();
for ( ; !isNodeEnforced && d2ij != ijByDist.end(); ++d2ij )
{
- i = d2ij->second.first;
- j = d2ij->second.second;
+ int di = d2ij->second.first;
+ int dj = d2ij->second.second;
// check if a node is at a side
int iSide = -1;
- if ( j == 0 )
+ if ( dj== 0 && j == 0 )
iSide = QUAD_BOTTOM_SIDE;
- else if ( j+1 == quad->jSize )
+ else if ( dj == 1 && j+2 == quad->jSize )
iSide = QUAD_TOP_SIDE;
- else if ( i == 0 )
+ else if ( di == 0 && i == 0 )
iSide = QUAD_LEFT_SIDE;
- else if ( i+1 == quad->iSize )
+ else if ( di == 1 && i+2 == quad->iSize )
iSide = QUAD_RIGHT_SIDE;
if ( iSide > -1 ) // ----- node is at a side
}
else // ------------------ node is inside the quad
{
+ i += di;
+ j += dj;
// make a new side passing through IJ node and split the quad
int indForced, iNewSide;
if ( quad->iSize < quad->jSize ) // split vertically
{
quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/true );
- indForced = i;
+ indForced = j;
iNewSide = splitQuad( quad, i, 0 );
}
else
{
quad->updateUV( myForcedPnts[ iFP ].uv, i, j, /*isVert=*/false );
- indForced = j;
+ indForced = i;
iNewSide = splitQuad( quad, 0, j );
}
FaceQuadStruct::Ptr newQuad = myQuadList.back();
// update UV of the side
vector<UVPtStruct>& sidePoints = (vector<UVPtStruct>&) side.GetUVPtStruct();
for ( int i = iFrom; i < iTo; ++i )
- sidePoints[ i ] = tmpQuad->UVPt( 1, i-iFrom );
+ {
+ const uvPtStruct& uvPt = tmpQuad->UVPt( 1, i-iFrom );
+ sidePoints[ i ].u = uvPt.u;
+ sidePoints[ i ].v = uvPt.v;
+ }
}
//================================================================================
grid = otherSide.grid;
from = otherSide.from;
to = otherSide.to;
+ di = otherSide.di;
forced_nodes = otherSide.forced_nodes;
contacts = otherSide.contacts;
nbNodeOut = otherSide.nbNodeOut;
int FaceQuadStruct::Side::ToSideIndex( int quadNodeIndex ) const
{
- return ( from > to ) ? ( from - quadNodeIndex ) : ( quadNodeIndex + from );
+ return from + di * quadNodeIndex;
}
//================================================================================
int FaceQuadStruct::Side::ToQuadIndex( int sideNodeIndex ) const
{
- return ( from > to ) ? ( from - sideNodeIndex ) : ( sideNodeIndex - from );
+ return ( sideNodeIndex - from ) * di;
}
//================================================================================
*/
//================================================================================
-bool FaceQuadStruct::Side::Reverse()
+bool FaceQuadStruct::Side::Reverse(bool keepGrid)
{
if ( grid )
{
- // if ( nbNodeOut == 0 )
- // {
- // if ( from > to )
- // {
- // from++;
- // to++;
- // }
- // else
- // {
- // from--;
- // to--;
- // }
- // std::swap( from, to );
- // }
- // else
+ if ( keepGrid )
+ {
+ from -= di;
+ to -= di;
+ std::swap( from, to );
+ di *= -1;
+ }
+ else
{
grid->Reverse();
}
double FaceQuadStruct::Side::Param( int i ) const
{
const vector<UVPtStruct>& points = GetUVPtStruct();
- return (( points[ from + i ].normParam - points[ from ].normParam ) /
- ( points[ to - 1 ].normParam - points[ from ].normParam ));
+ return (( points[ from + i * di ].normParam - points[ from ].normParam ) /
+ ( points[ to - 1 * di ].normParam - points[ from ].normParam ));
}
//================================================================================
{
const vector<UVPtStruct>& points = GetUVPtStruct();
double u = ( points[ from ].normParam +
- x * ( points[ to-1 ].normParam - points[ from ].normParam ));
+ x * ( points[ to-di ].normParam - points[ from ].normParam ));
return grid->Value2d( u ).XY();
}
