# igakit

committed 9b66773

Fix clamping on the right and add tests

# src/igakit/igalib.f90

`   real   (kind=8), intent(inout) :: U(0:n+p+1)`
`   real   (kind=8), intent(inout) :: Pw(d,0:n)`
`   logical(kind=4), intent(in)    :: l, r`
`+  real   (kind=8) :: uu`
`   integer(kind=4) :: k, s`
`   if (l) then ! Clamp at left end`
`+     uu = U(p)`
`+     s = FindMult(p,uu,p,U)`
`      k = p`
`-     s = FindMult(p,U(p),p,U)`
`-     if (s < p) call KntIns(d,n,p,U,Pw,k,s)`
`+     if (s < p) call KntIns(d,n,p,U,Pw,uu,k,p-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)`
`-     if (s < p) call KntIns(d,n,p,U,Pw,k,s)`
`+     uu = U(n+1)`
`+     s = FindMult(n,uu,p,U)`
`+     k = n+s`
`+     if (s < p) call KntIns(d,n,p,U,Pw,uu,k,p-s)`
`      U(n+2:n+p+1) = U(n+1)`
`   end if`
` contains`
`-  subroutine KntIns(d,n,p,U,Pw,k,s)`
`+  subroutine KntIns(d,n,p,U,Pw,uu,k,r)`
`       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(inout) :: Pw(d,0:n)`
`-      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)`
`-      uu = U(k)`
`-      r = p-s`
`-      Qw(:,0) = Pw(:,k-p)`
`-      Rw(:,0:p-s) = Pw(:,k-p:k-s)`
`+      real   (kind=8), intent(in)    :: uu`
`+      integer(kind=4), intent(in)    :: k, r`
`+      integer(kind=4) :: i, j, idx`
`+      real   (kind=8) :: alpha, Rw(d,0:r), Qw(d,0:r+r-1)`
`+      Qw(:,0)   = Pw(:,k-p)`
`+      Rw(:,0:r) = Pw(:,k-p:k-p+r)`
`       do j = 1, r`
`          idx = k-p+j`
`-         do i = 0, p-j-s`
`-            alpha = (uu-U(idx+i))/(U(i+k+1)-U(idx+i))`
`-            Rw(:,i) = alpha*Rw(:,i+1)+(1-alpha)*Rw(:,i)`
`+         do i = 0, r-j`
`+            alpha = (uu-U(idx+i))/(U(k+i+1)-U(idx+i))`
`+            Rw(:,i) = alpha*Rw(:,i+1)+(1.0-alpha)*Rw(:,i)`
`          end do`
`-         Qw(:,j) = Rw(:,0)`
`-         Qw(:,p-j-s+r) = Rw(:,p-j-s)`
`+         Qw(:,j)     = Rw(:,0)`
`+         Qw(:,r+r-j) = Rw(:,r-j)`
`       end do`
`       if (k == p) then ! left end`
`-         Pw(:,0:r-1) = Qw(:,r:r+r-1)`
`+         Pw(:,0:r-1)   = Qw(:,r:r+r-1)`
`       else             ! right end`
`-         Pw(:,n-r+1:n) = Qw(:,p-r:p-1)`
`+         Pw(:,n-r+1:n) = Qw(:,1:r)`
`       end if`
`     end subroutine KntIns`
` end subroutine ClampKnot`

# test/test_clamp.py

`+import numpy as np`
`+from igakit.cad import circle, Pi`
`+`
`+def make_crv(p,u):`
`+    c = circle(radius=1, angle=Pi/2)`
`+    c.rotate(Pi/4)`
`+    c.elevate(0,p-2)`
`+    c.refine(0,u)`
`+    return c`
`+`
`+def check_crv(c):`
`+    u0, u1 = c.breaks(0)[[0,-1]]`
`+    u = np.linspace(u0,u1,100)`
`+    x, y, z = c.evaluate(u).T`
`+    r = np.hypot(x,y)`
`+    return np.allclose(r, 1)`
`+`
`+def test_clamp():`
`+    for p in range(2,6):`
`+        for u in ([],[0.5],[1/3.0,2/3.0],[0.1,0.9]):`
`+            c = make_crv(p,u)`
`+            check_crv(c)`
`+            for continuity in range(c.degree[0]):`
`+                for side in (0, 1, None):`
`+                    cc = c.copy()`
`+                    cc.unclamp(0, continuity=continuity, side=side)`
`+                    check_crv(cc)`
`+                    cc.clamp(0, side=side)`
`+                    check_crv(cc)`
`+                    cc.clamp(0)`
`+                    check_crv(cc)`
`+                    assert np.allclose(cc.knots[0], c.knots[0])`
`+                    assert np.allclose(cc.array,    c.array)`
`+`
`+if __name__ == '__main__':`
`+    test_clamp()`
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.