# HG changeset patch # User Lisandro Dalcin # Date 1366813519 -10800 # Node ID 57ce076ce887613fde3f5849568f65996c5d63df # Parent 174fd488f719ac5c13f097d5e397c3518d0008db Add cad.refine() diff --git a/demo/refine.py b/demo/refine.py new file mode 100644 --- /dev/null +++ b/demo/refine.py @@ -0,0 +1,38 @@ +from igakit.cad import * + +# create two quarter circles +R0 = 1.0 +R1 = 2.0 +C0 = circle(radius=R0, angle=Pi/2) +C1 = circle(radius=R1, angle=Pi/2) + +# make a ruled surface out of the two arcs +srf = ruled(C0, C1) + +# make the radial direction first +srf.transpose() + +# refine the surface to have: +# * 3x6 nonempty knot spans +# * degree two in both directions +# * maximum allowable continuity +srf1 = refine(srf, factor=[3,6], degree=2) + +# refine the surface to have: +# * 5x10 nonempty knot spans +# * degree three in both directions +# * C^0 continuity +srf2 = refine(srf, factor=[5,10], degree=3, continuity=0) + +import sys +from igakit.plot import plt +try: + backend = sys.argv[1] + plt.use(backend) +except IndexError: + pass +for nrb in (srf, srf1, srf2): + plt.figure() + plt.cplot(nrb) + plt.kplot(nrb) +plt.show() diff --git a/makefile b/makefile --- a/makefile +++ b/makefile @@ -17,6 +17,7 @@ \$(PYTHON) demo/sweep.py \$(BACKEND) \$(PYTHON) demo/ruled.py \$(BACKEND) \$(PYTHON) demo/revolve.py \$(BACKEND) + \$(PYTHON) demo/refine.py \$(BACKEND) .PHONY: doctest doctest: \$(PYTHON) -m doctest src/igakit/transform.py @@ -41,7 +42,7 @@ \$(PYTHON) setup.py install --home=\${HOME} .PHONY: uninstall uninstall_home uninstall_user -uninstall: uninstall_user uninstall_home +uninstall: uninstall_user uninstall_user: -\$(RM) -r `\$(PYTHON) -m site --user-site`/\$(PACKAGE)* uninstall_home: diff --git a/src/igakit/cad.py b/src/igakit/cad.py --- a/src/igakit/cad.py +++ b/src/igakit/cad.py @@ -664,3 +664,97 @@ return nrb # ----- + +def refine(nrb, factor=None, degree=None, continuity=None): + """ + Refine a NURBS object by degree elevation and knot insertion. + + Parameters + ---------- + nrb : NURBS + factor : int or sequence of int, optional + degree : int or sequence of int, optional + continuity : int or sequence of int, optional + + """ + try: np_unique = np.lib.arraysetops.unique + except AttributeError: np_unique = np.unique1d + def Arg(arg, defval): + defval = tuple(defval) + dim = len(defval) + if arg is None: + return defval + arg = np.asarray(arg, dtype=object) + if arg.ndim == 0: + return arg.repeat(dim) + assert arg.shape == (dim,) + for i, val in enumerate(arg): + if val is None: + arg[i] = defval[i] + return arg + # + nrb = nrb.clone() + dim = nrb.dim + # clamping + for axis in range(dim): + nrb.clamp(axis) + # degree elevation + degree = Arg(degree, nrb.degree) + degree = np.asarray(degree, dtype='i') + assert np.all(degree > 0) + for axis in range(dim): + p = degree[axis] + t = p - nrb.degree[axis] + if t > 0: + nrb.elevate(axis, t) + p = nrb.degree[axis] + degree[axis] = p + # knot refinement + factor = Arg(factor, np.ones(dim, 'i')) + for i, fact in enumerate(factor): + factor[i] = np.asarray(fact, dtype='i') + assert np.all(factor[i] > 0) + continuity = Arg(continuity, [-1]*dim) + continuity = np.asarray(continuity, dtype='i') + continuity[continuity<0] += degree[continuity<0] + assert np.all(continuity >= 0) + assert np.all(continuity < degree) + for axis in range(dim): + N = factor[axis] + C = continuity[axis] + p = nrb.degree[axis] + # breaks & multiplicities + u, s = nrb.breaks(axis, mults=True) + # compute knots to insert + assert N.ndim == 0 or N.size == u.size-1 + delta_u = np.ediff1d(u)/N + u = u[:-1]; s = s[:-1]; + if N.ndim == 0: + U = u[:,np.newaxis].repeat(N, axis=1) + step = np.arange(1, N, dtype='d') + U[:,1:] += np.outer(delta_u, step) + S = np.empty((s.size, N), dtype='i') + S[:,0] = (p-C-s).clip(0, None) + S[:,1:] = p-C + u = U.ravel()[1:] + s = S.ravel()[1:] + else: + pos = np.empty(N.size+1, dtype='i') + pos[0] = 0; np.cumsum(N, out=pos[1:]) + U = np.empty(pos[-1], dtype='d') + S = np.empty(pos[-1], dtype='i') + for i, n in enumerate(N): + uu = np.arange(n, dtype='d') + uu *= delta_u[i]; uu += u[i] + a, b = pos[i], pos[i+1] + U[a:b] = uu + S[a] = max(p-C-s[i], 0) + S[a+1:b] = p-C + u = U[1:] + s = S[1:] + # insert knots + nrb.refine(axis, u.repeat(s)) + # + return nrb + +# -----