From: eap Date: Fri, 30 Dec 2005 07:16:14 +0000 (+0000) Subject: PAL10467 Meshing in quadrangles even if the number of nodes on opposite edges is... X-Git-Tag: T_Before_Join_BR-D5-38-2003~14 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=584d5262ff0c36b1431cfc36b1f6591c467423df;p=modules%2Fsmesh.git PAL10467 Meshing in quadrangles even if the number of nodes on opposite edges is not the same --- diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index e40390c6e..9ffa37d1f 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -51,10 +51,21 @@ using namespace std; #include #include #include +#include +#include #include "utilities.h" #include "Utils_ExceptHandlers.hxx" +#ifndef StdMeshers_Array2OfNode_HeaderFile +#define StdMeshers_Array2OfNode_HeaderFile +typedef const SMDS_MeshNode* SMDS_MeshNodePtr; +#include +DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +DEFINE_ARRAY2(StdMeshers_Array2OfNode, + StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr) +#endif + //============================================================================= /*! @@ -117,7 +128,28 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); aMesh.GetSubMesh(aShape); - FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); + //FaceQuadStruct *quad = CheckAnd2Dcompute(aMesh, aShape); + FaceQuadStruct* quad = CheckNbEdges(aMesh, aShape); + + if (!quad) + return false; + + if(myQuadranglePreference) { + int n1 = quad->nbPts[0]; + int n2 = quad->nbPts[1]; + int n3 = quad->nbPts[2]; + int n4 = quad->nbPts[3]; + int nfull = n1+n2+n3+n4; + int ntmp = nfull/2; + ntmp = ntmp*2; + if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) { + // special path for using only quandrangle faces + return ComputeQuadPref(aMesh, aShape, quad); + } + } + + // set normalized grid on unit square in parametric domain + SetNormalizedGrid(aMesh, aShape, quad); if (!quad) return false; @@ -490,14 +522,16 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh, return isOk; } + //============================================================================= /*! * */ //============================================================================= -FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute - (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception) +FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape) + throw(SALOME_Exception) { Unexpect aCatch(SalomeException); @@ -505,42 +539,56 @@ FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute // verify 1 wire only, with 4 edges - if (NumberOfWires(F) != 1) - { + if (NumberOfWires(F) != 1) { INFOS("only 1 wire by face (quadrangles)"); return 0; } const TopoDS_Wire& W = BRepTools::OuterWire(F); BRepTools_WireExplorer wexp (W, F); - FaceQuadStruct *quad = new FaceQuadStruct; + FaceQuadStruct* quad = new FaceQuadStruct; for (int i = 0; i < 4; i++) quad->uv_edges[i] = 0; quad->uv_grid = 0; int nbEdges = 0; - for (wexp.Init(W, F); wexp.More(); wexp.Next()) - { + for (wexp.Init(W, F); wexp.More(); wexp.Next()) { const TopoDS_Edge& E = wexp.Current(); int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); - if (nbEdges < 4) - { + if (nbEdges < 4) { quad->edge[nbEdges] = E; quad->nbPts[nbEdges] = nb + 2; // internal points + 2 extrema } nbEdges++; } - if (nbEdges != 4) - { + if (nbEdges != 4) { INFOS("face must have 4 edges /quadrangles"); QuadDelete(quad); return 0; } - // set normalized grid on unit square in parametric domain + return quad; +} + + +//============================================================================= +/*! + * + */ +//============================================================================= + +FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute + (SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) throw(SALOME_Exception) +{ + Unexpect aCatch(SalomeException); + + FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape); + + if(!quad) return 0; - SetNormalizedGrid(aMesh, F, quad); + // set normalized grid on unit square in parametric domain + SetNormalizedGrid(aMesh, aShape, quad); return quad; } @@ -789,6 +837,654 @@ void StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh, } } + +//======================================================================= +//function : ShiftQuad +//purpose : auxilary function for ComputeQuadPref +//======================================================================= +static void ShiftQuad(FaceQuadStruct* quad, const int num, bool WisF) +{ + if(num>3) return; + int i; + for(i=1; i<=num; i++) { + int nbPts3 = quad->nbPts[0]; + quad->nbPts[0] = quad->nbPts[1]; + quad->nbPts[1] = quad->nbPts[2]; + quad->nbPts[2] = quad->nbPts[3]; + quad->nbPts[3] = nbPts3; + TopoDS_Edge edge3 = quad->edge[0]; + quad->edge[0] = quad->edge[1]; + quad->edge[1] = quad->edge[2]; + quad->edge[2] = quad->edge[3]; + quad->edge[3] = edge3; + double first3 = quad->first[0]; + quad->first[0] = quad->first[1]; + quad->first[1] = quad->first[2]; + quad->first[2] = quad->first[3]; + quad->first[3] = first3; + double last3 = quad->last[0]; + quad->last[0] = quad->last[1]; + quad->last[1] = quad->last[2]; + quad->last[2] = quad->last[3]; + quad->last[3] = last3; + bool isEdgeForward3 = quad->isEdgeForward[0]; + quad->isEdgeForward[0] = quad->isEdgeForward[1]; + quad->isEdgeForward[1] = quad->isEdgeForward[2]; + quad->isEdgeForward[2] = quad->isEdgeForward[3]; + quad->isEdgeForward[3] = isEdgeForward3; + bool isEdgeOut3 = quad->isEdgeOut[0]; + quad->isEdgeOut[0] = quad->isEdgeOut[1]; + quad->isEdgeOut[1] = quad->isEdgeOut[2]; + quad->isEdgeOut[2] = quad->isEdgeOut[3]; + quad->isEdgeOut[3] = isEdgeOut3; + UVPtStruct* uv_edges3 = quad->uv_edges[0]; + quad->uv_edges[0] = quad->uv_edges[1]; + quad->uv_edges[1] = quad->uv_edges[2]; + quad->uv_edges[2] = quad->uv_edges[3]; + quad->uv_edges[3] = uv_edges3; + } + if(!WisF) { + // replacement left and right edges + int nbPts3 = quad->nbPts[1]; + quad->nbPts[1] = quad->nbPts[3]; + quad->nbPts[3] = nbPts3; + TopoDS_Edge edge3 = quad->edge[1]; + quad->edge[1] = quad->edge[3]; + quad->edge[3] = edge3; + double first3 = quad->first[1]; + quad->first[1] = quad->first[3]; + quad->first[3] = first3; + double last3 = quad->last[1]; + quad->last[1] = quad->last[2]; + quad->last[3] = last3; + bool isEdgeForward3 = quad->isEdgeForward[1]; + quad->isEdgeForward[1] = quad->isEdgeForward[3]; + quad->isEdgeForward[3] = isEdgeForward3; + bool isEdgeOut3 = quad->isEdgeOut[1]; + quad->isEdgeOut[1] = quad->isEdgeOut[3]; + quad->isEdgeOut[3] = isEdgeOut3; + UVPtStruct* uv_edges3 = quad->uv_edges[1]; + quad->uv_edges[1] = quad->uv_edges[3]; + quad->uv_edges[3] = uv_edges3; + } +} + + +//======================================================================= +//function : CalcUV +//purpose : auxilary function for ComputeQuadPref +//======================================================================= +static gp_XY CalcUV(double x0, double x1, double y0, double y1, + FaceQuadStruct* quad, + const gp_Pnt2d& a0, const gp_Pnt2d& a1, + const gp_Pnt2d& a2, const gp_Pnt2d& a3, + const Handle(Geom2d_Curve)& c2db, + const Handle(Geom2d_Curve)& c2dr, + const Handle(Geom2d_Curve)& c2dt, + const Handle(Geom2d_Curve)& c2dl) +{ + int nb = quad->nbPts[0]; + int nr = quad->nbPts[1]; + int nt = quad->nbPts[2]; + int nl = quad->nbPts[3]; + + UVPtStruct* uv_eb = quad->uv_edges[0]; + UVPtStruct* uv_er = quad->uv_edges[1]; + UVPtStruct* uv_et = quad->uv_edges[2]; + UVPtStruct* uv_el = quad->uv_edges[3]; + + double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + double y = y0 + x * (y1 - y0); + + double param_b = uv_eb[0].param + x * (uv_eb[nb-1].param - uv_eb[0].param); + double param_t = uv_et[0].param + x * (uv_et[nt-1].param - uv_et[0].param); + double param_r = uv_er[0].param + y * (uv_er[nr-1].param - uv_er[0].param); + double param_l = uv_el[0].param + y * (uv_el[nl-1].param - uv_el[0].param); + + gp_Pnt2d p0 = c2db->Value(param_b); + gp_Pnt2d p1 = c2dr->Value(param_r); + gp_Pnt2d p2 = c2dt->Value(param_t); + gp_Pnt2d p3 = c2dl->Value(param_l); + + double u = (1 - y) * p0.X() + x * p1.X() + y * p2.X() + (1 - x) * p3.X(); + double v = (1 - y) * p0.Y() + x * p1.Y() + y * p2.Y() + (1 - x) * p3.Y(); + + u -= (1 - x) * (1 - y) * a0.X() + x * (1 - y) * a1.X() + + x * y * a2.X() + (1 - x) * y * a3.X(); + v -= (1 - x) * (1 - y) * a0.Y() + x * (1 - y) * a1.Y() + + x * y * a2.Y() + (1 - x) * y * a3.Y(); + + //cout<<"x0="<ShapeToIndex( F ); + + int nb = quad->nbPts[0]; + int nr = quad->nbPts[1]; + int nt = quad->nbPts[2]; + int nl = quad->nbPts[3]; + 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 3 + ShiftQuad(quad,3,WisF); + } + else { + // we have to shift quad on 1 + ShiftQuad(quad,1,WisF); + } + } + + nb = quad->nbPts[0]; + nr = quad->nbPts[1]; + nt = quad->nbPts[2]; + nl = quad->nbPts[3]; + dh = abs(nb-nt); + dv = abs(nr-nl); + int nbh = Max(nb,nt); + int nbv = Max(nr,nl); + int addh = 0; + int addv = 0; + + // orientation of face and 3 main domain for future faces + // 0 top 1 + // 1------------1 + // | | | | + // | | | | + // | L | | R | + // left | | | | rigth + // | / \ | + // | / C \ | + // |/ \| + // 0------------0 + // 0 bottom 1 + + if(dh>dv) { + addv = (dh-dv)/2; + nbv = nbv + addv; + } + else { // dv>=dh + addh = (dv-dh)/2; + nbh = nbh + addh; + } + + Handle(Geom2d_Curve) c2d[4]; + for(i=0; i<4; i++) { + c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i], F, + quad->first[i], quad->last[i]); + } + + bool loadOk = true; + for(i=0; i<2; i++) { + quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], false); + if(!quad->uv_edges[i]) loadOk = false; + } + for(i=2; i<4; i++) { + quad->uv_edges[i] = LoadEdgePoints2(aMesh, F, quad->edge[i], true); + if (!quad->uv_edges[i]) loadOk = false; + } + if (!loadOk) { + INFOS("StdMeshers_Quadrangle_2D::ComputeQuadPref - LoadEdgePoints failed"); + QuadDelete( quad ); + quad = 0; + return false; + } + + UVPtStruct* uv_eb = quad->uv_edges[0]; + UVPtStruct* uv_er = quad->uv_edges[1]; + UVPtStruct* uv_et = quad->uv_edges[2]; + UVPtStruct* uv_el = quad->uv_edges[3]; + + // arrays for normalized params + //cout<<"Dump B:"<X()<<","<Y()<<","<Z()<<")"<D0(uv_eb[0].param,a[0]); + c2d[0]->D0(uv_eb[nb-1].param,a[1]); + c2d[2]->D0(uv_et[nt-1].param,a[2]); + c2d[2]->D0(uv_et[0].param,a[3]); + //cout<<" a[0]("<0) { + // add top nodes + for(i=1; i<=dl; i++) + NodesL.SetValue(i+1,nl,uv_et[i].node); + // create and add needed nodes + TColgp_SequenceOfXY UVtmp; + for(i=1; i<=dl; i++) { + double x0 = npt.Value(i+1); + double x1 = x0; + // diagonal node + double y0 = npl.Value(i+1); + double y1 = npr.Value(i+1); + gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3], + c2d[0], c2d[1], c2d[2], c2d[3]); + 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()); + NodesL.SetValue(i+1,1,N); + if(UVL.Length()Value(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesL.SetValue(i+1,j,N); + if( i==dl ) UVtmp.Append(UV); + } + } + for(i=1; i<=UVtmp.Length() && UVL.Length()X()<<","<Y()<<","<Z()<<")"; + // } + // cout<AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j), + NodesL.Value(i+1,j+1), NodesL.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + meshDS->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1), + NodesL.Value(i+1,j+1), NodesL.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVL using c2d + for(i=1; iD0(uv_el[i].param,p2d); + UVL.Append(p2d.XY()); + } + } + + // step2: create faces for right domain + StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr); + // add right nodes + for(j=1; j<=nr; j++) + NodesR.SetValue(1,j,uv_er[nr-j].node); + if(dr>0) { + // add top nodes + for(i=1; i<=dr; i++) + NodesR.SetValue(i+1,1,uv_et[nt-1-i].node); + // create and add needed nodes + TColgp_SequenceOfXY UVtmp; + for(i=1; i<=dr; i++) { + double x0 = npt.Value(nt-i); + double x1 = x0; + // diagonal node + double y0 = npl.Value(i+1); + double y1 = npr.Value(i+1); + gp_XY UV = CalcUV(x0, x1, y0, y1, quad, a[0], a[1], a[2], a[3], + c2d[0], c2d[1], c2d[2], c2d[3]); + 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()); + NodesR.SetValue(i+1,nr,N); + if(UVR.Length()Value(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesR.SetValue(i+1,j,N); + if( i==dr ) UVtmp.Prepend(UV); + } + } + for(i=1; i<=UVtmp.Length() && UVR.Length()AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j), + NodesR.Value(i+1,j+1), NodesR.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + meshDS->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1), + NodesR.Value(i+1,j+1), NodesR.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + } + else { + // fill UVR using c2d + for(i=1; iD0(uv_er[i].param,p2d); + UVR.Append(p2d.XY()); + } + } + + // step3: create faces for central domain + StdMeshers_Array2OfNode NodesC(1,nb,1,nbv); + // add first string using NodesL + for(i=1; i<=dl+1; i++) + NodesC.SetValue(1,i,NodesL(i,1)); + for(i=2; i<=nl; i++) + NodesC.SetValue(1,dl+i,NodesL(dl+1,i)); + // add last string using NodesR + for(i=1; i<=dr+1; i++) + NodesC.SetValue(nb,i,NodesR(i,nr)); + for(i=1; iD0(uv_eb[i-1].param,p2d); + } + // create and add needed nodes + // add linear layers + for(i=2; iValue(UV.X(),UV.Y()); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y()); + NodesC.SetValue(i,nbv-nnn+j,N); + } + } + // add diagonal layers + //cout<<"UVL.Length()="<Value(u,v); + SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace(N, geomFaceID, u, v); + NodesC.SetValue(j,i+1,N); + } + } + // create faces + for(i=1; iAddFace(NodesC.Value(i,j), NodesC.Value(i+1,j), + NodesC.Value(i+1,j+1), NodesC.Value(i,j+1)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + else { + SMDS_MeshFace* F = + meshDS->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1), + NodesC.Value(i+1,j+1), NodesC.Value(i+1,j)); + meshDS->SetMeshElementOnShape(F, geomFaceID); + } + } + } + + QuadDelete(quad); + bool isOk = true; + return isOk; +} + + +//============================================================================= +/*! + * LoadEdgePoints2 + */ +//============================================================================= +UVPtStruct* StdMeshers_Quadrangle_2D::LoadEdgePoints2 (SMESH_Mesh & aMesh, + const TopoDS_Face& F, + const TopoDS_Edge& E, + bool IsReverse) +{ + //MESSAGE("StdMeshers_Quadrangle_2D::LoadEdgePoints"); + // --- IDNodes of first and last Vertex + TopoDS_Vertex VFirst, VLast; + TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l + + ASSERT(!VFirst.IsNull()); + SMDS_NodeIteratorPtr lid = aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes(); + if (!lid->more()) { + MESSAGE ( "NO NODE BUILT ON VERTEX" ); + return 0; + } + const SMDS_MeshNode* idFirst = lid->next(); + + ASSERT(!VLast.IsNull()); + lid = aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes(); + if (!lid->more()) { + MESSAGE ( "NO NODE BUILT ON VERTEX" ); + return 0; + } + const SMDS_MeshNode* idLast = lid->next(); + + // --- edge internal IDNodes (relies on good order storage, not checked) + + map params; + SMDS_NodeIteratorPtr ite = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes(); + + while(ite->more()) { + const SMDS_MeshNode* node = ite->next(); + const SMDS_EdgePosition* epos = + static_cast(node->GetPosition().get()); + double param = epos->GetUParameter(); + params[param] = node; + } + + int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes(); + if (nbPoints != params.size()) { + MESSAGE( "BAD NODE ON EDGE POSITIONS" ); + return 0; + } + UVPtStruct* uvslf = new UVPtStruct[nbPoints + 2]; + + double f, l; + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l); + + const TopoDS_Wire& W = BRepTools::OuterWire(F); + bool FisF = (F.Orientation()==TopAbs_FORWARD); + bool WisF = (W.Orientation()==TopAbs_FORWARD); + bool isForward = (E.Orientation()==TopAbs_FORWARD); + //if(isForward) cout<<"E is FORWARD"<Value(f); // first point = Vertex Forward + uvslf[0].x = p.X(); + uvslf[0].y = p.Y(); + uvslf[0].param = f; + uvslf[0].node = idFirst; + //MESSAGE("__ f "<::iterator itp = params.begin(); + for (int i = 1; i <= nbPoints; i++) { // nbPoints internal + double param = (*itp).first; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].node = (*itp).second; + //MESSAGE("__ "<Value(l); // last point = Vertex Reversed + uvslf[nbPoints + 1].x = p.X(); + uvslf[nbPoints + 1].y = p.Y(); + uvslf[nbPoints + 1].param = l; + uvslf[nbPoints + 1].node = idLast; + //MESSAGE("__ l "<Value(l); // first point = Vertex Reversed + uvslf[0].x = p.X(); + uvslf[0].y = p.Y(); + uvslf[0].param = l; + uvslf[0].node = idLast; + //MESSAGE("__ l "<::reverse_iterator itp = params.rbegin(); + for (int j = nbPoints; j >= 1; j--) { // nbPoints internal + double param = (*itp).first; + int i = nbPoints + 1 - j; + gp_Pnt2d p = C2d->Value(param); + uvslf[i].x = p.X(); + uvslf[i].y = p.Y(); + uvslf[i].param = param; + uvslf[i].node = (*itp).second; + //MESSAGE("__ "<Value(f); // last point = Vertex Forward + uvslf[nbPoints + 1].x = p.X(); + uvslf[nbPoints + 1].y = p.Y(); + uvslf[nbPoints + 1].param = f; + uvslf[nbPoints + 1].node = idFirst; + //MESSAGE("__ f "<