Commits

Lisandro Dalcin committed db9ad98

Fixes and enhancements to clamping/unclamping

  • Participants
  • Parent commits 77731d1

Comments (0)

Files changed (4)

File demo/unclamp.py

 from igakit.cad import *
 
 c1 = circle(radius=1, angle=Pi/2)
-c2 = c1.copy().unclamp()
+c2 = c1.copy().unclamp(0)
 
 from igakit.plot import plt
 import sys

File src/igakit/igalib.f90

   if (l) then ! Clamp at left end
      k = p
      s = FindMult(p,U(p),p,U)
-     call KntIns(d,n,p,U,Pw,k,s)
+     if (s < p) call KntIns(d,n,p,U,Pw,k,s)
      U(0:p-1) = U(p)
   end if
   if (r) then ! Clamp at right end
      k = n+1
      s = FindMult(n,U(n+1),p,U)
-     call KntIns(d,n,p,U,Pw,k,s)
+     if (s < p) call KntIns(d,n,p,U,Pw,k,s)
      U(n+2:n+p+1) = U(n+1)
   end if
 contains
       integer(kind=4), intent(in)    :: k, s
       integer(kind=4) :: r, i, j, idx
       real   (kind=8) :: uu, alpha, Rw(d,0:p), Qw(d,0:2*p)
-      if (s >= p) return
       uu = U(k)
       r = p-s
       Qw(:,0) = Pw(:,k-p)
     end subroutine KntIns
 end subroutine ClampKnot
 
-subroutine UnclampKnot(d,n,p,U,Pw,l,r)
+subroutine UnclampKnot(d,n,p,U,Pw,C,l,r)
   implicit none
   integer(kind=4), intent(in)    :: d
   integer(kind=4), intent(in)    :: n, p
   real   (kind=8), intent(inout) :: U(0:n+p+1)
   real   (kind=8), intent(inout) :: Pw(d,0:n)
+  integer(kind=4), intent(in)    :: C
   logical(kind=4), intent(in)    :: l, r
-  integer(kind=4) :: i, j, k
+  integer(kind=4) :: m, i, j
   real   (kind=8) :: alpha
   if (l) then ! Unclamp at left end
-     do i = 0, p-2
-        U(p-i-1) = U(p-i) - (U(n-i+1)-U(n-i))
-        k = p-1
+     do i = 0, C
+        U(C-i) = U(p) - U(n+1) + U(n-i)
+     end do
+     do i = p-C-1, p-2
         do j = i, 0, -1
-           alpha = (U(p)-U(k))/(U(p+j+1)-U(k))
+           alpha = (U(p)-U(p+j-i-1))/(U(p+j+1)-U(p+j-i-1))
            Pw(:,j) = (Pw(:,j)-alpha*Pw(:,j+1))/(1-alpha)
-           k = k-1
         end do
      end do
-     U(0) = U(1) - (U(n-p+2)-U(n-p+1)) ! Set first knot
   end if
   if (r) then ! Unclamp at right end
-     do i = 0, p-2
-        U(n+i+2) = U(n+i+1) + (U(p+i+1)-U(p+i))
+     m = n+p+1
+     do i = 0, C
+        U(m-C+i) = U(n+1) - U(p) + U(p+i+1)
+     end do
+     do i = p-C-1, p-2
         do j = i, 0, -1
            alpha = (U(n+1)-U(n-j))/(U(n-j+i+2)-U(n-j))
            Pw(:,n-j) = (Pw(:,n-j)-(1-alpha)*Pw(:,n-j-1))/alpha
         end do
      end do
-     U(n+p+1) = U(n+p) + (U(2*p)-U(2*p-1)) ! Set last knot
   end if
 end subroutine UnclampKnot
 
   call ClampKnot(d,n,p,V,Qw,l,r)
 end subroutine Clamp
 
-subroutine Unclamp(d,n,p,U,Pw,l,r,V,Qw)
+subroutine Unclamp(d,n,p,U,Pw,C,l,r,V,Qw)
   use bspline
   implicit none
   integer(kind=4), intent(in)  :: d
   integer(kind=4), intent(in)  :: n, p
   real   (kind=8), intent(in)  :: U(0:n+p+1)
   real   (kind=8), intent(in)  :: Pw(d,0:n)
+  integer(kind=4), intent(in)  :: C
   logical(kind=4), intent(in)  :: l, r
   real   (kind=8), intent(out) :: V(0:n+p+1)
   real   (kind=8), intent(out) :: Qw(d,0:n)
   V  = U
   Qw = Pw
-  call UnclampKnot(d,n,p,V,Qw,l,r)
+  call ClampKnot(d,n,p,V,Qw,l,r)
+  call UnclampKnot(d,n,p,V,Qw,C,l,r)
 end subroutine Unclamp
 
 subroutine RefineKnotVector(d,n,p,U,Pw,r,X,Ubar,Qw)

File src/igakit/igalib.pyf

        real   (kind=8), intent(c,out) :: Qw(0:n,d)
      end subroutine Unclamp
 
-     subroutine Unclamp(d,n,p,U,Pw,l,r,V,Qw)
+     subroutine Unclamp(d,n,p,U,Pw,C,l,r,V,Qw)
        !f2py threadsafe
        integer(kind=4), intent(hide), check(d>=1) :: d
        integer(kind=4), intent(hide), check(n>=1) :: n
        integer(kind=4), intent(in),   check(p>=1) :: p
        real   (kind=8), intent(c,in)  :: U(0:n+p+1)
        real   (kind=8), intent(c,in)  :: Pw(0:n,d)
+       integer(kind=4), optional,     :: C = p - 1
        logical(kind=4), optional      :: l = 1
        logical(kind=4), optional      :: r = 1
        real   (kind=8), intent(c,out) :: V(0:n+p+1)

File src/igakit/nurbs.py

         self._knots = tuple(knots)
         return self
 
-    def clamp(self, *axes, **kargs):
+    def clamp(self, axis, side=None):
         """
 
         Examples
 
         >>> C = np.random.rand(3,3)
         >>> U = [0,0,0,1,1,1]
-        >>> c1 = NURBS([U], C)
+        >>> c1 = NURBS([U], C).unclamp(0)
         >>> c1.knots[0].tolist()
-        [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
-        >>> c2 = c1.clone().unclamp()
+        [-2.0, -1.0, 0.0, 1.0, 2.0, 3.0]
+        >>> c2 = c1.clone().clamp(0)
         >>> c2.knots[0].tolist()
-        [-2.0, -1.0, 0.0, 1.0, 2.0, 3.0]
-        >>> c3 = c2.clone().clamp()
-        >>> c3.knots[0].tolist()
         [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
 
         >>> C = np.random.rand(4,4)
         >>> U = [0,0,0,0.5,1,1,1]
-        >>> c1 = NURBS([U], C)
+        >>> c1 = NURBS([U], C).unclamp(0)
         >>> c1.knots[0].tolist()
-        [0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]
-        >>> c2 = c1.clone().unclamp()
+        [-1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0]
+        >>> c2 = c1.clone().clamp(0, side=0)
         >>> c2.knots[0].tolist()
-        [-1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0]
-        >>> c3 = c2.clone().clamp(side=0)
+        [0.0, 0.0, 0.0, 0.5, 1.0, 1.5, 2.0]
+        >>> c3 = c1.clone().clamp(0, side=1)
         >>> c3.knots[0].tolist()
-        [0.0, 0.0, 0.0, 0.5, 1.0, 1.5, 2.0]
-        >>> c4 = c3.clone().clamp(axis=0, side=1)
-        >>> c4.knots[0].tolist()
-        [0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]
+        [-1.0, -0.5, 0.0, 0.5, 1.0, 1.0, 1.0]
 
         """
-        allaxes = range(self.dim)
-        left, right = True, True
-        if not axes:
-            axis = kargs.pop('axis', None)
-            if axis is not None:
-                axes = (axis,)
-            else:
-                axes = allaxes
-        side = kargs.pop('side', None)
+        axis = range(self.dim)[axis]
+        array = self.array
+        knots = list(self.knots)
+        p = self.degree[axis]
+        U = knots[axis]
+        #
+        l, r = True, True
         if side is not None:
-            assert len(axes) == 1
             assert side in (0, 1)
-            if side == 0: right = False
-            else:         left  = False
-        assert not kargs
+            if side == 0: r = False
+            else:         l = False
         #
         Clamp = _bsp.Clamp
-        array = self.array
-        knots = list(self.knots)
-        degree = self.degree
-        #
-        for axis in axes:
-            axis = allaxes[axis]
-            p = degree[axis]
-            U = knots[axis]
-            Pw = np.rollaxis(array, axis, 0).copy()
-            shape = Pw.shape
-            Pw.shape = (shape[0], -1)
-            V, Qw = Clamp(p, U, Pw, left, right)
-            Qw.shape = shape
-            array = np.rollaxis(Qw, 0, axis+1)
-            knots[axis] = V
+        Pw = np.rollaxis(array, axis, 0).copy()
+        shape = Pw.shape
+        Pw.shape = (shape[0], -1)
+        V, Qw = Clamp(p, U, Pw, l, r)
+        Qw.shape = shape
+        array = np.rollaxis(Qw, 0, axis+1)
+        knots[axis] = V
         #
         self._array = np.ascontiguousarray(array)
         self._knots = tuple(knots)
         return self
 
-    def unclamp(self, *axes, **kargs):
+    def unclamp(self, axis, side=None, continuity=None):
         """
 
         Examples
         >>> c1.knots[0].tolist()
         [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
 
-        >>> c2 = c1.clone().unclamp()
+        >>> c2 = c1.clone().unclamp(0)
         >>> c2.knots[0].tolist()
         [-2.0, -1.0, 0.0, 1.0, 2.0, 3.0]
         >>> u = np.linspace(0,1,100)
         >>> np.allclose(xyz1, xyz2, rtol=0, atol=1e-15)
         True
 
-        >>> c3 = c1.clone().unclamp(side=0)
+        >>> c3 = c1.clone().unclamp(0,side=0)
         >>> c3.knots[0].tolist()
         [-2.0, -1.0, 0.0, 1.0, 1.0, 1.0]
         >>> u = np.linspace(0,1,100)
         >>> np.allclose(xyz1, xyz2, rtol=0, atol=1e-15)
         True
 
-        >>> c4 = c1.clone().unclamp(axis=0,side=1)
+        >>> c4 = c1.clone().unclamp(0, side=1)
         >>> c4.knots[0].tolist()
         [0.0, 0.0, 0.0, 1.0, 2.0, 3.0]
         >>> u = np.linspace(0,1,100)
         True
 
         """
-        allaxes = range(self.dim)
-        left, right = True, True
-        if not axes:
-            axis = kargs.pop('axis', None)
-            if axis is not None:
-                axes = (axis,)
-            else:
-                axes = allaxes
-        side = kargs.pop('side', None)
+        axis = range(self.dim)[axis]
+        array = self.array
+        knots = list(self.knots)
+        p = self.degree[axis]
+        U = knots[axis]
+        #
+        l, r = True, True
         if side is not None:
-            assert len(axes) == 1
             assert side in (0, 1)
-            if side == 0: right = False
-            else:         left  = False
-        assert not kargs
+            if side == 0: r = False
+            else:         l = False
+        C = p - 1
+        if continuity is not None:
+            assert 0 <= continuity <= p-1
+            C = continuity
         #
         Unclamp = _bsp.Unclamp
-        array = self.array
-        knots = list(self.knots)
-        degree = self.degree
-        #
-        for axis in axes:
-            axis = allaxes[axis]
-            p = self.degree[axis]
-            U = knots[axis]
-            Pw = np.rollaxis(array, axis, 0).copy()
-            shape = Pw.shape
-            Pw.shape = (shape[0], -1)
-            V, Qw = Unclamp(p, U, Pw, left, right)
-            Qw.shape = shape
-            array = np.rollaxis(Qw, 0, axis+1)
-            knots[axis] = V
+        Pw = np.rollaxis(array, axis, 0).copy()
+        shape = Pw.shape
+        Pw.shape = (shape[0], -1)
+        V, Qw = Unclamp(p, U, Pw, C, l, r)
+        Qw.shape = shape
+        array = np.rollaxis(Qw, 0, axis+1)
+        knots[axis] = V
         #
         self._array = np.ascontiguousarray(array)
         self._knots = tuple(knots)