From bb848edb77dee2cde85e47a5bff1eaf45ddbef90 Mon Sep 17 00:00:00 2001 From: benoit128 Date: Thu, 5 Sep 2024 15:56:11 +0200 Subject: [PATCH] CPlot: update getActivePointIndex, MapEdge: update local refinement for CAD --- Cassiopee/CPlot/CPlot/Plugins/mouseClick.cpp | 60 +++++--- Cassiopee/CPlot/CPlot/get.cpp | 34 ++++- Cassiopee/CPlot/apps/tkMapEdge.py | 150 ++++++++++++++++++- Cassiopee/Geom/Geom/MapEdge.py | 17 ++- Cassiopee/Intersector/test/closeCells.py | 5 +- Cassiopee/Intersector/test/closeCells_t1.py | 5 +- 6 files changed, 227 insertions(+), 44 deletions(-) diff --git a/Cassiopee/CPlot/CPlot/Plugins/mouseClick.cpp b/Cassiopee/CPlot/CPlot/Plugins/mouseClick.cpp index 390d6d237..93818e2b2 100644 --- a/Cassiopee/CPlot/CPlot/Plugins/mouseClick.cpp +++ b/Cassiopee/CPlot/CPlot/Plugins/mouseClick.cpp @@ -267,55 +267,71 @@ E_Int Data::findBlockContaining(double x, double y, double z, double xmi, ymi, zmi, xma, yma, zma; double d, dn, de; E_Int nzmin = -1; - double dmin = 1.e6; // distMin node/element - double dminNode = 1.e6; // distMin node - double dminElt = 1.e6; // distMin element + double dmin = K_CONST::E_MAX_FLOAT; // distMin node/element + double dminNode = K_CONST::E_MAX_FLOAT; // distMin node + double dminElt = K_CONST::E_MAX_FLOAT; // distMin element ind = 0; indE = 0; double eps = 0.05; d = MAX(xmax-xmin, ymax-ymin); d = MAX(d, zmax-zmin); eps = d*eps; - //printf("eps: %f\n", eps); ncon = 0; for (nz = 0; nz < _numberOfZones; nz++) { - Zone* zone = _zones[nz]; - if (zone->active == 1 || - (zone->active == 0 && ptrState->ghostifyDeactivatedZones == 1)) + Zone* zonep = _zones[nz]; + if (zonep->active == 1 || + (zonep->active == 0 && ptrState->ghostifyDeactivatedZones == 1)) { - xma = zone->xmax + eps; - xmi = zone->xmin - eps; - yma = zone->ymax + eps; - ymi = zone->ymin - eps; - zma = zone->zmax + eps; - zmi = zone->zmin - eps; - //printf("try: zone %f %f %f %f %f %f\n", xmi, xma, ymi, yma, zmi, zma); + xma = zonep->xmax + eps; + xmi = zonep->xmin - eps; + yma = zonep->ymax + eps; + ymi = zonep->ymin - eps; + zma = zonep->zmax + eps; + zmi = zonep->zmin - eps; if (x >= xmi && x <= xma && y >= ymi && y <= yma && z >= zmi && z <= zma) { if (nz < _numberOfStructZones) { - findNearestPoint(x, y, z, (StructZone*)zone, indl, dn); - inde = findElement(x, y, z, (StructZone*)zone, de); + findNearestPoint(x, y, z, (StructZone*)zonep, indl, dn); + inde = findElement(x, y, z, (StructZone*)zonep, de); } else { - findNearestPoint(x, y, z, (UnstructZone*)zone, indl, dn); - inde = findElement(x, y, z, (UnstructZone*)zone, de, ncon); + findNearestPoint(x, y, z, (UnstructZone*)zonep, indl, dn); + inde = findElement(x, y, z, (UnstructZone*)zonep, de, ncon); } d = MIN(dn, de); - if (zone->active == 0) d = d + eps*0.01; // malus + if (zonep->active == 0) d = d + eps*0.01; // malus - if (d < dmin-1.e-10) + if (d < dmin-1.e-10) // node ou centre plus proche { dmin = d; nzmin = nz; ind = indl; indE = inde; dminElt = de; dminNode = dn; } - else if (d <= dmin) // meme point mini: compare dn ou de + else if (d <= dmin+1.e-10) // meme point mini: compare dn ou de { - if (de < dminElt && _zones[nzmin]->dim != 1) + if (zonep->dim == 1 && _zones[nzmin]->dim == 1) + { + if (de < dminElt) + { + dmin = d; nzmin = nz; ind = indl; indE = inde; + dminElt = de; dminNode = dn; + } + else if (dn < dminNode) + { + dmin = d; nzmin = nz; ind = indl; indE = inde; + dminElt = de; dminNode = dn; + } + } + else if (zonep->dim == 1 && _zones[nzmin]->dim != 1) + { + dmin = d; nzmin = nz; ind = indl; indE = inde; + dminElt = de; dminNode = dn; + } + else if (de < dminElt && _zones[nzmin]->dim != 1) { dmin = d; nzmin = nz; ind = indl; indE = inde; dminElt = de; dminNode = dn; diff --git a/Cassiopee/CPlot/CPlot/get.cpp b/Cassiopee/CPlot/CPlot/get.cpp index 267e61310..cc8dc9df5 100644 --- a/Cassiopee/CPlot/CPlot/get.cpp +++ b/Cassiopee/CPlot/CPlot/get.cpp @@ -19,8 +19,16 @@ #include "cplot.h" #include "Data.h" -int findFace(double xp, double yp, double zp, E_Int elt, - UnstructZone* zone, double& dist); +E_Int findNearestPoint(double xp, double yp, double zp, + StructZone* zone, E_Int& ind, double& dist); +E_Int findNearestPoint(double xp, double yp, double zp, + UnstructZone* zone, E_Int& ind, double& dist); +E_Int findElement(double xp, double yp, double zp, + UnstructZone* zone, double& dist, E_Int& ncon); +E_Int findElement(double xp, double yp, double zp, + StructZone* zone, double& dist); +E_Int findFace(double xp, double yp, double zp, E_Int elt, + UnstructZone* zone, double& dist); //====================================================== // Get mode from PyObject, return an int (0: mesh, 1: solid, @@ -390,15 +398,31 @@ PyObject* K_CPLOT::getActivePointIndex(PyObject* self, PyObject* args) if (nz == 0) return l; else { - // Re-check (if zone has been replaced) double posX, posY, posZ; E_Int zone, ind, indE, ncon; double dist; posX = d->ptrState->activePointX; posY = d->ptrState->activePointY; posZ = d->ptrState->activePointZ; - d->findBlockContaining(posX, posY, posZ, - zone, ind, indE, dist, ncon); + + // recompute zone and index + //d->findBlockContaining(posX, posY, posZ, + // zone, ind, indE, dist, ncon); + + // only recompute index + double dn, de; + zone = d->ptrState->selectedZone-1; Zone* z = d->_zones[zone]; + if (zone < d->_numberOfStructZones) + { + findNearestPoint(posX, posY, posZ, (StructZone*)z, ind, dn); + indE = findElement(posX, posY, posZ, (StructZone*)z, de); + } + else + { + findNearestPoint(posX, posY, posZ, (UnstructZone*)z, ind, dn); + indE = findElement(posX, posY, posZ, (UnstructZone*)z, de, ncon); + } + if (zone < d->_numberOfStructZones) { StructZone* zz = (StructZone*)z; diff --git a/Cassiopee/CPlot/apps/tkMapEdge.py b/Cassiopee/CPlot/apps/tkMapEdge.py index 00cd8505c..ee8ea175e 100644 --- a/Cassiopee/CPlot/apps/tkMapEdge.py +++ b/Cassiopee/CPlot/apps/tkMapEdge.py @@ -117,11 +117,11 @@ def setEnforce(event=None): nob = CTK.Nb[nz]+1 noz = CTK.Nz[nz] z = CTK.t[2][nob][2][noz] - setEnforceZ(z) + h = CTK.varsFromWidget(VARS[1].get(), 1) + setEnforceZ(z, h, VARS[8].get()) # Perform a setH on z with clicked point -def setEnforceZ(z): - h = CTK.varsFromWidget(VARS[1].get(), 1) +def setEnforceZ(z, h, mode): if len(h) != 1: CTK.TXT.insert('START', 'Invalid spacing.\n') CTK.TXT.insert('START', 'Error: ', 'Error') @@ -131,7 +131,7 @@ def setEnforceZ(z): ind = CPlot.getActivePointIndex() if ind == []: return True ind = ind[0] - if VARS[8].get() == 'HFactor': + if mode == 'HFactor': P0 = C.getValue(z, 'GridCoordinates', ind) if ind == npts-1: P1 = C.getValue(z, 'GridCoordinates', ind-1) else: P1 = C.getValue(z, 'GridCoordinates', ind+1) @@ -166,7 +166,7 @@ def setEnforceMode(event=None): nob = CTK.Nb[nz]+1 noz = CTK.Nz[nz] z = CTK.t[2][nob][2][noz] - setEnforceZ(z) + setEnforceZ(z, CTK.varsFromWidget(VARS[1].get(), 1), VARS[8].get()) CTK.__BUSY__ = False TTK.raiseButton(W) @@ -1061,6 +1061,7 @@ def copyDistrib(): CTK.setCursor(0, WIDGETS['frame']) +#============================================================================== # get selection from CTK.t def getSelection(nzs): zones = [] @@ -1070,7 +1071,9 @@ def getSelection(nzs): zones.append(CTK.t[2][nob][2][noz]) return zones +#============================================================================== # remesh CAD when an edge is modified +#============================================================================== def remeshCAD(edges): try: import OCC.PyTree as OCC except: return @@ -1088,7 +1091,125 @@ def remeshCAD(edges): [h,hmax,hausd] = CTK.CADHOOK OCC._remeshTreeFromEdges(h, CTK.t, valids) CTK.display(CTK.t) + +#============================================================================== +# enforce h in edge locally +#============================================================================== +def enforceLocal(event=None): + v = CTK.varsFromWidget(VARS[12].get(), 1) + if len(v) != 1: + CTK.TXT.insert('START', 'Invalid h or hfactor.\n') + CTK.TXT.insert('START', 'Error: ', 'Error') + return + v = v[0] + width = WIDGETS['widthScale'].get() / 100. + width = max(width, 0.1) + mode = VARS[11].get() + + # get edge + nzs = CPlot.getSelectedZones() + nz = nzs[0] + nob = CTK.Nb[nz]+1 + noz = CTK.Nz[nz] + z = CTK.t[2][nob][2][noz] + dim = Internal.getZoneDim(z) + if dim[4] != 1: + CTK.TXT.insert('START', 'Zone must be an edge.\n') + CTK.TXT.insert('START', 'Error: ', 'Error') + return + npts = C.getNPts(z) + + CTK.saveTree() + CPlot.setState(cursor=1) + # split zone of width + ind = CPlot.getActivePointIndex() + if ind == []: return + ind = ind[0] + + P0 = C.getValue(z, 'GridCoordinates', ind) + if ind == npts-1: P1 = C.getValue(z, 'GridCoordinates', ind-1) + else: P1 = C.getValue(z, 'GridCoordinates', ind+1) + hloc = Vector.norm(Vector.sub(P1,P0)) + if mode == 'HFactor': + h = v*hloc; factor = v + else: h = v; factor = h / hloc + + delta = int(width*npts)+1 + imin = ind+1-delta + imax = ind+1+delta + + CAD = Internal.getNodeFromName1(z, 'CAD') + + if imin > 1: z0 = T.subzone(z, (1,1,1), (imin,-1,-1)) + else: z0 = None + if ind > 2: zp1 = T.subzone(z, (max(imin,1),1,1), (ind+1,-1,-1)) + else: zp1 = None + if ind < npts-1: zp2 = T.subzone(z, (ind+1,1,1), (min(imax, npts),-1,-1)) + else: zp2 = None + if imax < npts: z1 = T.subzone(z, (imax,1,1), (npts,-1,-1)) + else: z1 = None + + if zp1 is not None: + P0 = C.getValue(zp1, 'GridCoordinates', 0) + P1 = C.getValue(zp1, 'GridCoordinates', 1) + h1 = Vector.norm(Vector.sub(P1,P0)) + L = D.getLength(zp1) + if h+h1 > L: h1 = h + + if zp2 is not None: + P0 = C.getValue(zp2, 'GridCoordinates', -1) + P1 = C.getValue(zp2, 'GridCoordinates', -2) + h2 = Vector.norm(Vector.sub(P1,P0)) + L = D.getLength(zp2) + if h+h2 > L: h2 = h + + if zp1 is None: h1 = h2 + if zp2 is None: h2 = h1 + + #if z0 is None: h1 = (h1+h)*0.5 + #if z1 is None: h2 = (h2+h)*0.5 + + # D._setH(zp, ind-imin+1, h) + # # guess a cool number of points + # if factor == 1.: + # N = npts + # elif factor < 1.: + # N = npts + (1./factor)/100.*npts + # N = int(N)+1 + # else: + # N = npts - (1./factor)/100.*npts + # N = int(N)+1 + # print("nbre de points=", N) + # D._enforceh(zp, N=N) + + if zp1 is not None: + d2 = D.distrib2(zp1, h1, h, algo=1) + zp1 = G.map(zp1, d2) + if zp2 is not None: + d2 = D.distrib2(zp2, h, h2, algo=1) + zp2 = G.map(zp2, d2) + + zo = None + if zp1 is not None: zo = zp1 + if zp2 is not None: + if zo is not None: zo = T.join(zo, zp2) + else: zo = zp2 + if z0 is not None: zo = T.join(z0, zo) + if z1 is not None: zo = T.join(zo, z1) + zo[0] = z[0] # keep orig name and CAD + zo[2].append(CAD) + + CTK.replace(CTK.t, nob, noz, zo) + (CTK.Nb, CTK.Nz) = CPlot.updateCPlotNumbering(CTK.t) + CTK.TKTREE.updateApp() + CPlot.render() + CTK.TXT.insert('START', 'Local spacing enforced.\n') + + # add CAD remesh if possible + edges = getSelection(nzs) + remeshCAD(edges) + CPlot.setState(cursor=0) #============================================================================== # Create app widgets @@ -1145,6 +1266,10 @@ def createApp(win): V = TK.StringVar(win); V.set('NFactor'); VARS.append(V) # -10- Nbre de pts pour h enforce V = TK.StringVar(win); V.set('1.'); VARS.append(V) + # -11- Type of local enforce + V = TK.StringVar(win); V.set('HFactor'); VARS.append(V) + # -12- Factor for local enforce + V = TK.StringVar(win); V.set('1.'); VARS.append(V) # - Uniformize - B = TTK.Button(Frame, text="Uniformize", command=uniformize) @@ -1217,6 +1342,21 @@ def createApp(win): B.grid(row=5, column=2, columnspan=2, sticky=TK.EW) BB = CTK.infoBulle(parent=B, text='Enforced number of points.') B.bind('', enforceH) + + # - Enforce local - + B = TTK.Button(Frame, text="Enforce", command=enforceLocal) + B.grid(row=6, column=0, columnspan=1, sticky=TK.EW) + BB = CTK.infoBulle(parent=B, text='Enforce local spacing.') + B = TTK.OptionMenu(Frame, VARS[11], 'HFactor', 'H') + B.grid(row=6, column=1, sticky=TK.EW) + B = TTK.Entry(Frame, textvariable=VARS[12], background='White', width=7) + B.grid(row=6, column=2, columnspan=1, sticky=TK.EW) + BB = CTK.infoBulle(parent=B, text='Enforced variable.') + B.bind('', enforceLocal) + B = TTK.Scale(Frame, from_=0, to=100, orient=TK.HORIZONTAL, + showvalue=0, borderwidth=1, value=50) + WIDGETS['widthScale'] = B + B.grid(row=6, column=3, sticky=TK.EW) #============================================================================== # Called to display widgets diff --git a/Cassiopee/Geom/Geom/MapEdge.py b/Cassiopee/Geom/Geom/MapEdge.py index 753c406f3..2132adbb9 100644 --- a/Cassiopee/Geom/Geom/MapEdge.py +++ b/Cassiopee/Geom/Geom/MapEdge.py @@ -432,12 +432,17 @@ def distrib2(a, h1, h2, add=20, forceAdd=False, normalized=True, algo=0, verbose if normalized: a = D.getDistribution(a) return a else: # geometrique - q = (L-h1)/(L-h2) - if verbose > 0: print("Info: distrib2: geometric progression: %g"%q) - N = numpy.log(h2/h1) / numpy.log(q) - N = int(T.kround(N))+2 - if N <= 2: print('Error: distrib2: not enough point to remesh.') - a = G.cartr1((0,0,0), (h1,1,1), (q,1,1), (N,1,1)) + if abs(h2-h1) < 1.e-6: # constant step + N = int(T.kround(L / h1)) + h = L/N + a = G.cart((0,0,0), (h,1,1), (N,1,1)) + else: + q = (L-h1)/(L-h2) + if verbose > 0: print("Info: distrib2: geometric progression: %g"%q) + N = numpy.log(h2/h1) / numpy.log(q) + N = int(T.kround(N))+2 + if N <= 2: print('Error: distrib2: not enough point to remesh.') + a = G.cartr1((0,0,0), (h1,1,1), (q,1,1), (N,1,1)) if normalized: a = D.getDistribution(a) return a diff --git a/Cassiopee/Intersector/test/closeCells.py b/Cassiopee/Intersector/test/closeCells.py index 61fa2db42..2d381556b 100644 --- a/Cassiopee/Intersector/test/closeCells.py +++ b/Cassiopee/Intersector/test/closeCells.py @@ -1,4 +1,4 @@ -# - triangulateExteriorFaces (array) - +# - closeCells (array) - import Intersector as XOR import Converter as C @@ -6,5 +6,4 @@ m = C.convertArray2NGon(m[0]) m = XOR.closeCells(m) -C.convertArrays2File([m], 'out.plt') - +C.convertArrays2File([m], 'out.plt') \ No newline at end of file diff --git a/Cassiopee/Intersector/test/closeCells_t1.py b/Cassiopee/Intersector/test/closeCells_t1.py index 974f7003a..9550f3aec 100644 --- a/Cassiopee/Intersector/test/closeCells_t1.py +++ b/Cassiopee/Intersector/test/closeCells_t1.py @@ -1,4 +1,4 @@ -# - triangulateExteriorFaces (array) - +# - closeCells (array) - import Intersector as XOR import Converter as C import KCore.test as test @@ -8,5 +8,4 @@ m = XOR.closeCells(m) -test.testA([m],1) - +test.testA([m],1) \ No newline at end of file