diff --git a/Cassiopee/CPlot/apps/tkCADFix.py b/Cassiopee/CPlot/apps/tkCADFix.py index ac64e171f..c49f0f596 100644 --- a/Cassiopee/CPlot/apps/tkCADFix.py +++ b/Cassiopee/CPlot/apps/tkCADFix.py @@ -30,7 +30,7 @@ def readCAD(event=None): [hmin, hmax, hausd] = OCC.getCADcontainer(CTK.t) if hmax is None or (hmax < 0 and hausd < 0): (hmax,hmin,hausd) = OCC.occ.analyseEdges(hook) - OCC._setCADcontainer(CTK.t, fileName, fileFmt, hmax, hausd) + OCC._setCADcontainer(CTK.t, fileName, fileFmt, hmin, hmax, hausd) CTK.CADHOOK = hook # remesh and redisplay CTK.setCursor(2, WIDGETS['frame']) @@ -410,15 +410,16 @@ def checkWatertight(event=None): if CAD is not None: hmax = Internal.getNodeFromName1(CAD, 'hmax') hmax = Internal.getValue(hmax) - tol = hmax/10. + tol = hmax/100. CTK.t = C.addBase2PyTree(CTK.t, 'LEAKS', 1) p = Internal.getNodeFromName1(CTK.t, 'LEAKS') gnob = C.getNobOfBase(p, CTK.t) f = Internal.getZones(b) + f = G.zip(f, tol) f = T.join(f) - f = G.close(f, tol) + #f = G.close(f, tol) ef = T.splitConnexity(f) VARS[6].set('Components: %d'%(len(ef))) diff --git a/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp b/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp index 38ef999ff..689bb25e7 100644 --- a/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp +++ b/Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp @@ -253,170 +253,11 @@ void geom4(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u printf("h1=%f real=%f\n", h1, ue[N-1]-ue[N-2]); } -// ============================================================================ -// Return the nbPoints and ue for meshing E with best of deflection and hmax -// ============================================================================ -E_Int __getParamHmaxHausd(const TopoDS_Edge& E, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue) -{ - // First call param hausd - E_Int ret = __getParamHausd(E, hausd, nbPoints, ue); - - BRepAdaptor_Curve C0(E); - GeomAdaptor_Curve geomAdap(C0.Curve()); - Standard_Real u0 = geomAdap.FirstParameter(); - Standard_Real u1 = geomAdap.LastParameter(); - E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1); - - // Then split in region h > hmax and h < hmax - E_Float delta; - E_Int state = -1; // if 0, we are in a h < hmax zone, if 1 in a h >= hmax zone - std::vector index; - for (E_Int i = 1; i < nbPoints; i++) - { - delta = (ue[i]-ue[i-1])/(u1-u0)*L; - //printf("%f %f\n", (ue[i]-u0)/(u1-u0), delta); - if (state == -1) - { - if (delta < hmax) state = 0; - else state = 1; - } - else if (state == 0 && delta >= hmax) - { - state = 1; - index.push_back(i); - } - else if (state == 1 && delta < hmax) - { - state = 0; - index.push_back(i); - } - } - for (size_t i = 0; i < index.size(); i++) printf("split %zu\n", i); - - E_Int size = index.size(); - if (size == 0 && state == 0) - { - // One zone, all well refined, nothing to do - printf("already fine\n"); - } - else if (size == 0 && state == 1) - { - // One zone but too coarse, we regenerate a full hmax distribution - E_Int np = E_Int(L/hmax)+1; - if (np == 1) np = 2; - printf("remesh 1 zone with hmax (%d)\n", np); - E_Float* ue2 = new E_Float [np]; - //E_Float hreg = L/(np-1); - for (E_Int i = 0; i < np; i++) ue2[i] = (i*1.)/(np-1)*(u1-u0)+u0; - delete [] ue; - ue = ue2; - nbPoints = np; - } - else if (size == 1 && state == 0) - { - E_Int is = index[0]; - // Two zones, the last is already refined - printf("2 zones, last refined, split in %d\n", is); - E_Float hr = ue[is+1]-ue[is]; - E_Float h0 = ue[1]-ue[0]; - E_Float* ur; E_Int N; - geom2(ue[0], ue[is], h0, hr, N, ur); - // merge distrib - E_Int Nf = nbPoints-is+N; - E_Float* uf = new E_Float [Nf]; - for (E_Int i = 0; i < N; i++) uf[i] = ur[i]; - for (E_Int i = N; i < Nf; i++) uf[i] = ue[is-N+i]; - // switch - delete [] ue; delete [] ur; - ue = uf; - nbPoints = Nf; - } - else - { - // general case - printf("general case =========================\n"); - std::vector ul(size+1); - std::vector Nl(size+1); - E_Int isp = 0; - - if (state%2 != 0) - { - if (state == 1) state = 0; - else state = 1; - } - printf("starting state=%d\n", state); - - for (E_Int l = 0; l < size; l++) - { - E_Int is = index[l]; - if (state == 1) // remesh - { - E_Float u0 = ue[isp]; - E_Float u1 = ue[is]; - E_Float h0 = ue[isp+1]-ue[isp]; - E_Float h1 = ue[is+1]-ue[is]; - geom2(u0, u1, h0, h1, Nl[l], ul[l]); - state = 0; - } - else // copy - { - Nl[l] = is-isp+1; - ul[l] = new E_Float [Nl[l]]; - E_Float* ull = ul[l]; - for (E_Int i = isp; i <= is; i++) ull[i] = ue[i]; - state = 1; - } - isp = is; - } - E_Int l = size; - E_Int is = index[l-1]; - if (state == 1) - { - E_Float u0 = ue[is]; - E_Float u1 = ue[nbPoints-1]; - E_Float h0 = ue[is+1]-ue[is]; - E_Float h1 = ue[nbPoints-1]-ue[nbPoints-2]; - geom1(u0, u1, h0, h1, Nl[l], ul[l]); - } - else - { - Nl[l] = nbPoints-isp+1; - ul[l] = new E_Float [Nl[l]]; - E_Float* ull = ul[l]; - for (E_Int i = is; i < nbPoints; i++) ull[i] = ue[i]; - } - - // Merge - E_Int Nf = 0; - for (E_Int l = 0; l <= size; l++) - { - Nf += Nl[l]; - } - E_Float* uf = new E_Float [Nf]; - - E_Int b = 0; - for (E_Int l = 0; l <= size; l++) - { - E_Float* ull = ul[l]; - for (E_Int i = 0; i < Nl[l]; i++) uf[b] = ull[i]; - b += Nl[l]-1; - } - // switch - for (E_Int l = 0; l <= size; l++) delete [] ul[l]; - delete [] ue; - ue = uf; - nbPoints = Nf; - } - - return ret; -} - // ============================================================================ // Return the nbPoints and ue for meshing E with hmin/hmax/hausd evaluated -// only on edge -// using geometric progression between points +// only on edge using geometric progression between points // ============================================================================ -E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue) +E_Int __getParamHminHmaxHausdE(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue) { BRepAdaptor_Curve C0(E); GeomAdaptor_Curve geomAdap(C0.Curve()); @@ -524,9 +365,9 @@ E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, // ============================================================================ // Return the nbPoints and ue for meshing E with hmin/hmax/hausd evaluated -// only on edge using linear evaluation between points +// only on edge using progressive walk // ============================================================================ -E_Int __getParamHminHmaxHausd4(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue) +E_Int __getParamHminHmaxHausdE4(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue) { BRepAdaptor_Curve C0(E); GeomAdaptor_Curve geomAdap(C0.Curve()); @@ -541,31 +382,30 @@ E_Int __getParamHminHmaxHausd4(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, } E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1); L = K_FUNC::E_max(L, 1.e-12); - - // Compute local h + + gp_Pnt point; + gp_Vec dE, d2E; + E_Float ues, dE1, dE2, dE3, d2E1, d2E2, d2E3; + E_Float h, t1, t2, t3, a, b, rho; + Handle(Geom_Curve) curve; Standard_Real first, last; curve = BRep_Tool::Curve(E, first, last); - gp_Pnt point; - gp_Vec dE, d2E; - E_Float U, dE1, dE2, dE3, d2E1, d2E2, d2E3; - E_Float t1, t2, t3, a, b, rho; - E_Int npts = 5; - std::vector Us(npts); - for (E_Int i = 0; i < npts; i++) - { - Us[i] = u0+(i/(npts-1.))*(u1-u0); - //GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, i*L/(npts-1.), u0); - //Us[i] = Pt.Parameter(); - } - std::vector h(npts); + E_Float Ltot = 0.; + E_Int N = 0; + E_Int sizeCont = 5000; + E_Float* hi = new E_Float [sizeCont]; // a dimensionner dynamiquement + ues = u0; h = 0.; - for (E_Int i = 0; i < npts; i++) + while (Ltot < L) { - U = Us[i]; + GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, h, ues); + ues = Pt.Parameter(); + + // Compute local h // Compute the derivatives - curve->D2(U, point, dE, d2E); + curve->D2(ues, point, dE, d2E); // compute curvature dE1 = dE.X(); @@ -586,115 +426,68 @@ E_Int __getParamHminHmaxHausd4(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, else rho = 1.e12; // Compute local h - h[i] = std::sqrt(8.*hausd*rho); - h[i] = K_FUNC::E_min(h[i], hmax); - h[i] = K_FUNC::E_max(h[i], hmin); - //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(8.*hausd*rho), h[i]); - } - - printf("hi dimensioned= "); - for (E_Int i = 0; i < npts; i++) - { printf("%g ", h[i]); } - printf("\n"); - - for (E_Int i = 0; i < npts; i++) - { - // passage en parametres - h[i] = h[i]/L; - } - - //printf("h parameter= "); - //for (E_Int i = 0; i < npts; i++) - //{ printf("%g ", h[i]); } - //printf("\n"); - - // compute integral of h dxsi - E_Float h1, h2, Nf, res, xi, r; - E_Int N, i0, i1; - E_Float dxi = 1./(npts-1); - - E_Float integ = 0.; - for (E_Int i = 0; i < npts-1; i++) - { - h1 = h[i]; h2 = h[i+1]; - integ += dxi*(h1+h2)*0.5; // mid - } - //printf("integ=%g\n", integ); - - Nf = 1./integ; - N = std::rint(Nf); - N = std::max(N, E_Int(2)); - - printf("N= %d, Nf = %g\n", N, Nf); - - // scale h to match N - res = Nf/N; - if (std::abs(res) > 1.e-10) - { - for (E_Int i = 0; i < npts; i++) h[i] = h[i]*res; + h = std::sqrt(8.*hausd*rho); + h = K_FUNC::E_min(h, hmax); + h = K_FUNC::E_max(h, hmin); + hi[N] = h; + Ltot += h; N += 1; + if (N >= sizeCont) + { + sizeCont += 5000; + E_Float* hi2 = new E_Float [sizeCont]; + for (E_Int i = 0; i < N; i++) hi2[i] = hi[i]; + delete [] hi; + hi = hi2; + } } - // DBX + // smooth hi /* - integ = 0.; - for (E_Int i = 0; i < npts-1; i++) + for (E_Int i = 0; i < N-1; i++) { - h1 = h[i]; h2 = h[i+1]; - integ += dxi*(h1+h2)*0.5; // mid + hi[i] = 0.5*(hi[i+1]+hi[i]); } - Nf = 1./integ+1.; - printf("Nf after rescale=%g\n", Nf); + hi[N-1] = 0.5*(hi[N-2]+hi[N-1]); */ - // END DBX - // get linear interpolation of h on new discretization - std::vector hn(N); - dxi = 1./(N-1); - for (E_Int j = 0; j < N; j++) - { - xi = j*dxi; - i0 = xi*(npts-1); - i1 = std::min(i0+1, npts-1); - r = xi-i0/(npts-1); - hn[j] = (1.-r)*h[i0]+r*h[i1]; - } - - printf("hn remesh\n"); - for (E_Int i = 0; i < N; i++) - { printf("%g ", hn[i]); } - printf("\n"); + E_Float r = Ltot-L; + if (r > 0.5*hi[N-1] && N >= 3) { N = N-1; Ltot = Ltot - hi[N]; } + + //printf("hi = "); + //for (E_Int i = 0; i < N; i++) printf("%g ", hi[i]); + //printf("\n"); - // for each point, get x - ue = new E_Float [N]; + // scale hi to match full L + E_Float s = L-Ltot; + s = s/N; + for (E_Int i = 0; i < N; i++) hi[i] = hi[i]+s; - ue[0] = 0.; - for (E_Int j = 1; j < N; j++) + // DBX + //Ltot = 0.; + //for (E_Int i = 0; i < N; i++) Ltot += hi[i]; + //printf("Ltot=%g %g\n", Ltot, L); + // ENDDBX + + // Get ue + ue = new E_Float [N+1]; + ues = u0; + ue[0] = u0; + for (E_Int i = 0; i < N; i++) { - ue[j] = ue[j-1] + N*hn[j]*dxi; + GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, hi[i], ues); + ues = Pt.Parameter(); + ue[i+1] = ues; } + ue[N] = u1; // forced - if (std::abs(ue[N-1]-1.) > 1.e-10) - { - // rescale - res = 1./ue[N-1]; - for (E_Int j = 1; j < N; j++) - { - ue[j] = ue[j]*res; - } - } + //printf("N= %d\n", N); + //printf("ue5= "); + //for (E_Int i = 0; i < N; i++) printf("%g ", ue[i]); + //printf("\n"); - printf("ue end\n"); - for (E_Int i = 0; i < N; i++) - { printf("%g ", ue[i]); } - printf("\n"); + delete [] hi; + nbPoints = N+1; - for (E_Int j = 0; j < N; j++) - { - //ue[j] = ue[j]*(u1-u0)+ u0; - GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, ue[j]*L, u0); - ue[j] = Pt.Parameter(); - } - nbPoints = N; return 1; } @@ -702,7 +495,7 @@ E_Int __getParamHminHmaxHausd4(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, // Return the nbPoints and ue for meshing E with hmin/hmax/hausd evaluated // as max of faces using geometric progression between points // =============================================================================== -E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue, +E_Int __getParamHminHmaxHausdF(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue, TopoDS_Shape* shape) { TopTools_IndexedDataMapOfShapeListOfShape edgeFaceMap; @@ -727,7 +520,7 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, gp_Pnt2d Puv; E_Float u,v; E_Float ues; - Standard_Real aFirst = C0.FirstParameter(), aEnd=C0.LastParameter(); + Standard_Real aFirst=C0.FirstParameter(), aEnd=C0.LastParameter(); //printf("afirst=%g aend=%g\n", aFirst, aEnd); Standard_Real pFirst = aFirst, pEnd=aEnd; @@ -738,10 +531,7 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Int npts = 2; std::vector Us(npts); for (E_Int i = 0; i < npts; i++) Us[i] = i/(npts-1.); - - //Us[0] = Us[0]+0.001*(Us[1]-Us[0]); - //Us[npts-1] = Us[npts-1]-0.001*(Us[npts-1]-Us[npts-2]); - + std::vector h(npts); for (E_Int i = 0; i < npts; i++) h[i] = K_CONST::E_MAX_FLOAT; @@ -779,7 +569,6 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float M11, M12, M21, M22, detM; E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS; - //printf("surf: u=%g, v=%g\n", u, v); // estimer les derivees sur la surface surface->D2(u, v, point, DU1, DV1, DU2, DV2, DUV); //printf("point=%g %g %g\n", point.X(), point.Y(), point.Z()); @@ -937,9 +726,7 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, //printf("final: "); //for (E_Int i = 0; i < Ntot; i++) - //{ - // printf("%g ", ue[i]); - //} + //{ printf("%g ", ue[i]); } //printf("\n"); for (E_Int i = 0; i < npts; i++) @@ -951,311 +738,11 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, } -// =============================================================================== -// Return the nbPoints and ue for meshing E with hmin/hmax/hausd evaluated -// as max of faces using linear between points evaluation -// =============================================================================== -E_Int __getParamHminHmaxHausd3(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue, - TopoDS_Shape* shape) -{ - TopTools_IndexedDataMapOfShapeListOfShape edgeFaceMap; - TopExp::MapShapesAndAncestors(*shape, TopAbs_EDGE, TopAbs_FACE, edgeFaceMap); - - BRepAdaptor_Curve C0(E); - GeomAdaptor_Curve geomAdap(C0.Curve()); - Standard_Real u0 = geomAdap.FirstParameter(); - Standard_Real u1 = geomAdap.LastParameter(); - //printf("===============edge =========================\n"); - //printf("u0=%g u1=%g\n", u0, u1); - if (BRep_Tool::Degenerated(E)) - { - nbPoints = 2; - ue = new E_Float [nbPoints]; - for (E_Int i = 0; i < nbPoints; i++) ue[i] = u0; - return 1; - } - E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1); - L = K_FUNC::E_max(L, 1.e-12); - - gp_Pnt2d Puv; - E_Float u,v, ues; - Standard_Real aFirst=C0.FirstParameter(), aEnd=C0.LastParameter(); - //printf("afirst=%g aend=%g\n", aFirst, aEnd); - Standard_Real pFirst=aFirst, pEnd=aEnd; - - gp_Pnt point; - gp_Vec DU1, DV1, DU2, DV2, DUV; - E_Float rho, hh, K1, K2; - - E_Int npts = 20; // number of evaluation points - std::vector Us(npts); - for (E_Int i = 0; i < npts; i++) { Us[i] = i/(npts-1.); } - //Us[0] = Us[0]+0.001*(Us[1]-Us[0]); - //Us[npts-1] = Us[npts-1]-0.001*(Us[npts-1]-Us[npts-2]); - - std::vector h(npts); - for (E_Int i = 0; i < npts; i++) h[i] = K_CONST::E_MAX_FLOAT; - - // Find faces of E - const TopTools_ListOfShape& connectedFaces = edgeFaceMap.FindFromKey(E); - for (TopTools_ListIteratorOfListOfShape it(connectedFaces); it.More(); it.Next()) - { - const TopoDS_Face& F = TopoDS::Face(it.Value()); - Handle(Geom_Surface) surface = BRep_Tool::Surface(F); - - Standard_Real BU1, BU2, BV1, BV2; surface->Bounds(BU1, BU2, BV1, BV2); - //printf("face bound = %g %g and %g %g\n", BU1, BU2, BV1, BV2); - - // p curve on F - Handle(Geom2d_Curve) pCurve = BRep_Tool::CurveOnSurface(E, F, pFirst, pEnd); - //printf("pfirst=%g pend=%g\n", pFirst, pEnd); - //printf("pCurve = %p\n", (void*)pCurve); - // estimate at some points - for (E_Int i = 0; i < npts; i++) - { - //ues = Us[i]*(pEnd-pFirst)+pFirst; - GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, Us[i]*L, u0); - ues = Pt.Parameter(); - C0.D0(ues, point); - printf("param=%g, %g %g %g\n", ues, point.X(), point.Y(), point.Z()); - - pCurve->D0(ues, Puv); - u = Puv.X(); v = Puv.Y(); // u,v edge on surface - if (u <= BU1) u = BU1+1.e-1; - if (u >= BU2) u = BU2-1.e-1; - if (v <= BV1) v = BV1+1.e-1; - if (v >= BV2) v = BV2-1.e-1; - -// must be MINE (0) or OCC (1) -#define CURVATURE 1 -#if CURVATURE == 0 - E_Float nx, ny, nz, n2; - E_Float G11, G12, G21, G22, detG; - E_Float M11, M12, M21, M22, detM; - E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS; - - //printf("surf: u=%g, v=%g\n", u, v); - // estimer les derivees sur la surface - surface->D2(u, v, point, DU1, DV1, DU2, DV2, DUV); - //printf("point=%g %g %g\n", point.X(), point.Y(), point.Z()); - //printf("dU1=%g %g %g\n", DU1.X(), DU1.Y(), DU1.Z()); - //printf("dV1=%g %g %g\n", DV1.X(), DV1.Y(), DV1.Z()); - - // get normal - nx = DU1.Y()*DV1.Z()-DU1.Z()*DV1.Y(); - ny = DU1.Z()*DV1.X()-DU1.X()*DV1.Z(); - nz = DU1.X()*DV1.Y()-DU1.Y()*DV1.X(); - n2 = nx*nx+ny*ny+nz*nz; - n2 = std::sqrt(n2); - if (n2 < 1.e-24) n2 = 1.e12; - else n2 = 1./n2; - nx = nx*n2; - ny = ny*n2; - nz = nz*n2; - //printf("normale= %g %g %g norm=%g\n", nx, ny, nz, nx*nx+ny*ny+nz*nz); - - // Metric - M11 = DU1.X()*DU1.X()+DU1.Y()*DU1.Y()+DU1.Z()*DU1.Z(); - M22 = DV1.X()*DV1.X()+DV1.Y()*DV1.Y()+DV1.Z()*DV1.Z(); - M12 = M21 = DU1.X()*DV1.X()+DU1.Y()*DV1.Y()+DU1.Z()*DV1.Z(); - //printf("matrice M:\n"); - //printf(" %g %g\n", M11, M12); - //printf(" %g %g\n", M21, M22); - - detM = M11*M22-M12*M21; - if (std::abs(detM) < 1.e-12) det = 1.e12; // invalid metric - else det = 1./detM; - Mi11 = det*M22; - Mi21 = -det*M21; - Mi12 = -det*M12; - Mi22 = det*M11; - detM = std::abs(detM); - //printf("detM=%g\n", detM); - - // matrice de courbure - G11 = DU2.X()*nx+DU2.Y()*ny+DU2.Z()*nz; - G22 = DV2.X()*nx+DV2.Y()*ny+DV2.Z()*nz; - G12 = G21 = DUV.X()*nx+DUV.Y()*ny+DUV.Z()*nz; - detG = G11*G22-G12*G21; - detG = std::abs(detG); - - S11 = G11*Mi11 + G12*Mi21; - S12 = G11*Mi12 + G12*Mi22; - S21 = G21*Mi11 + G22*Mi21; - S22 = G21*Mi12 + G22*Mi22; - - detS = S11*S22-S12*S21; - //printf("dU2=%g %g %g\n", DU2.X(), DU2.Y(), DU2.Z()); - //printf("dV2=%g %g %g\n", DV2.X(), DV2.Y(), DV2.Z()); - //printf("dUV=%g %g %g\n", DUV.X(), DUV.Y(), DUV.Z()); - //printf("detG=%g\n", detG); - - // courbure de gauss - rayon K1*K2 - //if (detS > 1.e-12) rho = 1./detS; - //else rho = 1.e12; - //printf("gauss=%g ou bien %g\n", detG/detM, detS); - - // K1 et K2 - //printf("matrice G:\n"); - //printf(" %g %g\n", G11, G12); - //printf(" %g %g\n", G21, G22); - - K_LINEAR::DelaunayMath sl; - sl.eigen_values(S11, S22, S12, K2, K1); - //printf("K1=%g K2=%g K1K2=%g\n", K1, K2, K1*K2); - //printf("mean=%g\n", 0.5*(S11+S22)); - - rho = K_FUNC::E_max(std::abs(K1), std::abs(K2)); - if (rho > 1.e-12) rho = 1./rho; - else rho = 1.e12; - -#else - // with open cascade - GeomLProp_SLProps props(surface, u, v, 2, Precision::Confusion()); - - if (props.IsCurvatureDefined()) - { - Standard_Real gaussianCurvature = props.GaussianCurvature(); - Standard_Real meanCurvature = props.MeanCurvature(); - - //std::cout << "OCC: gauss: " << gaussianCurvature << std::endl; - //std::cout << "OCC: mean: " << meanCurvature << std::endl; - - // Principal curvatures - K1 = meanCurvature + sqrt(meanCurvature * meanCurvature - gaussianCurvature); - K2 = meanCurvature - sqrt(meanCurvature * meanCurvature - gaussianCurvature); - //std::cout << "OCC: k1: " << K1 << " k2: " < 1.e-12) rho = 1./rho; - else rho = 1.e12; - } - else - { - rho = 1.e12; - //std::cout << "OCC: Curvature is not defined at the given parameters." << std::endl; - } -#endif - - hh = std::sqrt(8.*hausd*rho); - hh = K_FUNC::E_min(hh, hmax); - hh = K_FUNC::E_max(hh, hmin); - //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(1.*hausd*rho), hh); - h[i] = K_FUNC::E_min(h[i], hh); - } - } - - printf("Length = %g\n", L); - printf("hi dimensioned= "); - for (E_Int i = 0; i < npts; i++) - { printf("%g ", h[i]); } - printf("\n"); - - for (E_Int i = 0; i < npts; i++) - { - // passage en parametres - h[i] = h[i]/L; - } - - //printf("h parameter= "); - //for (E_Int i = 0; i < npts; i++) - //{ printf("%g ", h[i]); } - //printf("\n"); - - // compute integral of h dxsi - E_Float h1, h2, Nf, res, xi, r; - E_Int N, i0, i1; - E_Float dxi = 1./(npts-1); - - E_Float integ = 0.; - for (E_Int i = 0; i < npts-1; i++) - { - h1 = h[i]; h2 = h[i+1]; - integ += dxi*(h1+h2)*0.5; // mid - } - printf("integ=%g\n", integ); - - Nf = 1./integ; - N = std::rint(Nf); - N = std::max(N, E_Int(2)); - - printf("N= %d, Nf = %g\n", N, Nf); - - // scale h to match N - res = Nf/N; - if (std::abs(res) > 1.e-10) - { - for (E_Int i = 0; i < npts; i++) h[i] = h[i]*res; - } - - // DBX - /* - integ = 0.; - for (E_Int i = 0; i < npts-1; i++) - { - h1 = h[i]; h2 = h[i+1]; - integ += dxi*(h1+h2)*0.5; // mid - } - Nf = 1./integ+1.; - printf("Nf after rescale=%g\n", Nf); - */ - // END DBX - - // get linear interpolation of h on new discretization - std::vector hn(N); - dxi = 1./(N-1); - for (E_Int j = 0; j < N; j++) - { - xi = j*dxi; - i0 = xi*(npts-1); - i1 = std::min(i0+1, npts-1); - r = xi-i0/(npts-1); - hn[j] = (1.-r)*h[i0]+r*h[i1]; - } - - printf("hn remesh\n"); - for (E_Int i = 0; i < N; i++) - { printf("%g ", hn[i]); } - printf("\n"); - - // for each point, get x - ue = new E_Float [N]; - - ue[0] = 0.; - for (E_Int j = 1; j < N; j++) - { - ue[j] = ue[j-1]+N*hn[j]*dxi; - } - - if (std::abs(ue[N-1]-1.) > 1.e-10) - { - // rescale - res = 1./ue[N-1]; - for (E_Int j = 1; j < N; j++) - { - ue[j] = ue[j]*res; - } - } - - //printf("ue end\n"); - //for (E_Int i = 0; i < N; i++) - //{ printf("%g ", ue[i]); } - //printf("\n"); - - for (E_Int j = 0; j < N; j++) - { - //ue[j] = ue[j]*(aEnd-aFirst)+ aFirst; - GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, ue[j]*L, u0); - ue[j] = Pt.Parameter(); - } - nbPoints = N; - return 1; -} - // =============================================================================== // Return the nbPoints and ue for meshing E with hmin/hmax/hausd evaluated // as max of faces using progressive walk // =============================================================================== -E_Int __getParamHminHmaxHausd5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue, +E_Int __getParamHminHmaxHausdF5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue, TopoDS_Shape* shape) { TopTools_IndexedDataMapOfShapeListOfShape edgeFaceMap; @@ -1266,7 +753,7 @@ E_Int __getParamHminHmaxHausd5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, Standard_Real u0 = geomAdap.FirstParameter(); Standard_Real u1 = geomAdap.LastParameter(); //printf("===============edge =========================\n"); - //printf("u0=%g u1=%g\n", u0, u1); + //printf("edge u0=%g u1=%g\n", u0, u1); if (BRep_Tool::Degenerated(E)) { nbPoints = 2; @@ -1274,6 +761,7 @@ E_Int __getParamHminHmaxHausd5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, for (E_Int i = 0; i < nbPoints; i++) ue[i] = u0; return 1; } + E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1); L = K_FUNC::E_max(L, 1.e-12); @@ -1286,19 +774,21 @@ E_Int __getParamHminHmaxHausd5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, gp_Pnt point; gp_Vec DU1, DV1, DU2, DV2, DUV; E_Float rho, hh, K1, K2, h; - + + E_Float Ltot = 0.; E_Int N = 0; - ue = new E_Float [1000]; + E_Int sizeCont = 5000; + E_Float* hi = new E_Float [sizeCont]; // a dimensionner dynamiquement + ues = u0; h = 0.; while (Ltot < L) { - h = K_CONST::E_MAX_FLOAT; - - GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, Ltot, u0); + GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, h, ues); ues = Pt.Parameter(); - ue[N] = ues; + h = K_CONST::E_MAX_FLOAT; + // Find faces of E const TopTools_ListOfShape& connectedFaces = edgeFaceMap.FindFromKey(E); for (TopTools_ListIteratorOfListOfShape it(connectedFaces); it.More(); it.Next()) @@ -1434,21 +924,68 @@ E_Int __getParamHminHmaxHausd5(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, hh = std::sqrt(8.*hausd*rho); hh = K_FUNC::E_min(hh, hmax); hh = K_FUNC::E_max(hh, hmin); - //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(1.*hausd*rho), hh); + //printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(1.*hausd*rho), hh); fflush(stdout); h = K_FUNC::E_min(h, hh); } + hi[N] = h; Ltot += h; N += 1; + if (N >= sizeCont) + { + sizeCont += 5000; + E_Float* hi2 = new E_Float [sizeCont]; + for (E_Int i = 0; i < N; i++) hi2[i] = hi[i]; + delete [] hi; + hi = hi2; + } } - printf("N= %d\n", N); - for (E_Int i = 0; i < N; i++) printf("%g ", ue[i]); - printf("\n"); - // scale ue to match u0, u1 exactly - + // smooth hi + /* + for (E_Int i = 0; i < N-1; i++) + { + hi[i] = 0.5*(hi[i+1]+hi[i]); + } + hi[N-1] = 0.5*(hi[N-2]+hi[N-1]); + */ + + E_Float r = Ltot-L; + if (r > 0.5*hi[N-1] && N >= 3) { N = N-1; Ltot = Ltot - hi[N]; } + //printf("hi = "); + //for (E_Int i = 0; i < N; i++) printf("%g ", hi[i]); + //printf("\n"); + + // scale hi to match full L + E_Float s = L-Ltot; + s = s/N; + for (E_Int i = 0; i < N; i++) hi[i] = hi[i]+s; - nbPoints = N; + // DBX + //Ltot = 0.; + //for (E_Int i = 0; i < N; i++) Ltot += hi[i]; + //printf("Ltot=%g %g\n", Ltot, L); + // ENDDBX + + // Get ue + ue = new E_Float [N+1]; + ues = u0; + ue[0] = u0; + for (E_Int i = 0; i < N; i++) + { + GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, hi[i], ues); + ues = Pt.Parameter(); + ue[i+1] = ues; + } + ue[N] = u1; // forced + + //printf("N= %d\n", N); + //printf("ue5= "); + //for (E_Int i = 0; i < N; i++) printf("%g ", ue[i]); + //printf("\n"); + + delete [] hi; + nbPoints = N+1; return 1; } @@ -1462,15 +999,21 @@ E_Int __getParamExt(const TopoDS_Edge& E, E_Int nbPoints, E_Float* uext, E_Float Standard_Real u0 = geomAdap.FirstParameter(); Standard_Real u1 = geomAdap.LastParameter(); ue = new E_Float [nbPoints]; - //for (E_Int i = 0; i < nbPoints; i++) ue[i] = uext[i]*(u1-u0)+u0; E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1); - E_Float abscissa; + E_Float abscissa = 0.; for (E_Int i = 0; i < nbPoints; i++) { abscissa = uext[i]*L; - GCPnts_AbscissaPoint Pt(1.e-10, geomAdap, abscissa, u0); + GCPnts_AbscissaPoint Pt(1.e-12, geomAdap, abscissa, u0); ue[i] = Pt.Parameter(); + + // Maybe faster but less accurate + //h = uext[i]*L - abscissa; + //GCPnts_AbscissaPoint Pt(1.e-12, geomAdap, h, us); + //us = Pt.Parameter(); + //ue[i] = us; + //abscissa = uext[i]*L; } //for (E_Int i = 0; i < nbPoints; i++) printf("ue: %d : %f %f\n", i, ue[i], uext[i]); @@ -1648,7 +1191,7 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args) else if (not useFaces && hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd { // courbure uniquement sur edge - __getParamHminHmaxHausd(E, hmin, hmax, hausd, nbPoints, ue); + __getParamHminHmaxHausdE4(E, hmin, hmax, hausd, nbPoints, ue); PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1); FldArrayF* f; K_ARRAY::getFromArray2(o, f); __meshEdge(E, nbPoints, ue, *f, false); @@ -1659,10 +1202,7 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args) else if (useFaces && hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd { // courbure sur le max des faces - //__getParamHminHmaxHausd2(E, hmin, hmax, hausd, nbPoints, ue, shape); - __getParamHminHmaxHausd3(E, hmin, hmax, hausd, nbPoints, ue, shape); - //__getParamHminHmaxHausd4(E, hmin, hmax, hausd, nbPoints, ue); - + __getParamHminHmaxHausdF5(E, hmin, hmax, hausd, nbPoints, ue, shape); PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1); FldArrayF* f; K_ARRAY::getFromArray2(o, f); __meshEdge(E, nbPoints, ue, *f, false);