1function skl = nrbsurfderiveval (srf, uv, d) 2 3% 4% NRBSURFDERIVEVAL: Evaluate n-th order derivatives of a NURBS surface 5% 6% usage: skl = nrbsurfderiveval (srf, [u; v], d) 7% 8% INPUT: 9% 10% srf : NURBS surface structure, see nrbmak 11% 12% u, v : parametric coordinates of the point where we compute the 13% derivatives 14% 15% d : number of partial derivatives to compute 16% 17% 18% OUTPUT: 19% 20% skl (i, j, k, l) = i-th component derived j-1,k-1 times at the 21% l-th point. 22% 23% Adaptation of algorithm A4.4 from the NURBS book, pg137 24% 25% Copyright (C) 2009 Carlo de Falco 26% 27% This program is free software: you can redistribute it and/or modify 28% it under the terms of the GNU General Public License as published by 29% the Free Software Foundation, either version 3 of the License, or 30% (at your option) any later version. 31 32% This program is distributed in the hope that it will be useful, 33% but WITHOUT ANY WARRANTY; without even the implied warranty of 34% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 35% GNU General Public License for more details. 36% 37% You should have received a copy of the GNU General Public License 38% along with this program. If not, see <http://www.gnu.org/licenses/>. 39 40 skl = zeros (3, d+1, d+1, size (uv, 2)); 41 42 for iu = 1:size(uv, 2); 43 wders = squeeze (surfderiveval (srf.number(1)-1, srf.order(1)-1, ... 44 srf.knots{1}, srf.number(2)-1, ... 45 srf.order(2)-1, srf.knots{2}, ... 46 squeeze (srf.coefs(4, :, :)), uv(1,iu), ... 47 uv(2,iu), d)); 48 49 for idim = 1:3 50 Aders = squeeze (surfderiveval (srf.number(1)-1, srf.order(1)-1, ... 51 srf.knots{1}, srf.number(2)-1, ... 52 srf.order(2)-1, srf.knots{2}, ... 53 squeeze (srf.coefs(idim, :, :)), uv(1,iu),... 54 uv(2,iu), d)); 55 56 for k=0:d 57 for l=0:d-k 58 v = Aders(k+1, l+1); 59 for j=1:l 60 v = v - nchoosek(l,j)*wders(1,j+1)*skl(idim, k+1, l-j+1,iu); 61 end 62 for i=1:k 63 v = v - nchoosek(k,i)*wders(i+1,1)*skl(idim, k-i+1, l+1, iu); 64 v2 =0; 65 for j=1:l 66 v2 = v2 + nchoosek(l,j)*wders(i+1,j+1)*skl(idim, k-i+1, l-j+1, iu); 67 end 68 v = v - nchoosek(k,i)*v2; 69 end 70 skl(idim, k+1, l+1, iu) = v/wders(1,1); 71 end 72 end 73 end 74 75 end 76 77end 78 79%!test 80%! k = [0 0 1 1]; 81%! c = [0 1]; 82%! [coef(2,:,:), coef(1,:,:)] = meshgrid (c, c); 83%! coef(3,:,:) = coef(1,:,:); 84%! srf = nrbmak (coef, {k, k}); 85%! [u, v] = meshgrid (linspace(0,1,11)); 86%! uv = [u(:)';v(:)']; 87%! skl = nrbsurfderiveval (srf, uv, 0); 88%! aux = nrbeval(srf,uv); 89%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 90 91 92%!test 93%! k = [0 0 1 1]; 94%! c = [0 1]; 95%! [coef(2,:,:), coef(1,:,:)] = meshgrid (c, c); 96%! coef(3,:,:) = coef(1,:,:); 97%! srf = nrbmak (coef, {k, k}); 98%! srf = nrbkntins (srf, {[], rand(2,1)}); 99%! [u, v] = meshgrid (linspace(0,1,11)); 100%! uv = [u(:)';v(:)']; 101%! skl = nrbsurfderiveval (srf, uv, 0); 102%! aux = nrbeval(srf,uv); 103%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 104 105%!shared srf, uv 106%!test 107%! k = [0 0 0 1 1 1]; 108%! c = [0 1/2 1]; 109%! [coef(1,:,:), coef(2,:,:)] = meshgrid (c, c); 110%! coef(3,:,:) = coef(1,:,:); 111%! srf = nrbmak (coef, {k, k}); 112%! ders= nrbderiv (srf); 113%! [u, v] = meshgrid (linspace(0,1,11)); 114%! uv = [u(:)';v(:)']; 115%! skl = nrbsurfderiveval (srf, uv, 1); 116%! [fun, der] = nrbdeval (srf, ders, uv); 117%! assert (squeeze (skl (1:2,1,1,:)), fun(1:2,:), 1e3*eps) 118%! assert (squeeze (skl (1:2,2,1,:)), der{1}(1:2,:), 1e3*eps) 119%! assert (squeeze (skl (1:2,1,2,:)), der{2}(1:2,:), 1e3*eps) 120%! 121%!test 122%! srf = nrbdegelev (srf, [3, 1]); 123%! ders= nrbderiv (srf); 124%! [fun, der] = nrbdeval (srf, ders, uv); 125%! skl = nrbsurfderiveval (srf, uv, 1); 126%! assert (squeeze (skl (1:2,1,1,:)), fun(1:2,:), 1e3*eps) 127%! assert (squeeze (skl (1:2,2,1,:)), der{1}(1:2,:), 1e3*eps) 128%! assert (squeeze (skl (1:2,1,2,:)), der{2}(1:2,:), 1e3*eps) 129 130%!shared uv 131%!test 132%! k = [0 0 0 1 1 1]; 133%! c = [0 1/2 1]; 134%! [coef(2,:,:), coef(1,:,:)] = meshgrid (c, c); 135%! coef(3,:,:) = coef(1,:,:); 136%! srf = nrbmak (coef, {k, k}); 137%! ders= nrbderiv (srf); 138%! [u, v] = meshgrid (linspace(0,1,11)); 139%! uv = [u(:)';v(:)']; 140%! skl = nrbsurfderiveval (srf, uv, 1); 141%! [fun, der] = nrbdeval (srf, ders, uv); 142%! assert (squeeze (skl (1:2,1,1,:)), fun(1:2,:), 1e3*eps) 143%! assert (squeeze (skl (1:2,2,1,:)), der{1}(1:2,:), 1e3*eps) 144%! assert (squeeze (skl (1:2,1,2,:)), der{2}(1:2,:), 1e3*eps) 145%! 146%!test 147%! p = 3; q = 3; 148%! mcp = 5; ncp = 5; 149%! Lx = 10*rand(1); Ly = Lx; 150%! srf = nrbdegelev (nrb4surf ([0 0], [Lx, 0], [0 Ly], [Lx Ly]), [p-1, q-1]); 151%! %%srf = nrbkntins (srf, {linspace(0,1,mcp-p+2)(2:end-1), linspace(0,1,ncp-q+2)(2:end-1)}); 152%! %%srf.coefs = permute (srf.coefs, [1 3 2]); 153%! ders= nrbderiv (srf); 154%! [fun, der] = nrbdeval (srf, ders, uv); 155%! skl = nrbsurfderiveval (srf, uv, 1); 156%! assert (squeeze (skl (1:2,1,1,:)), fun(1:2,:), 1e3*eps) 157%! assert (squeeze (skl (1:2,2,1,:)), der{1}(1:2,:), 1e3*eps) 158%! assert (squeeze (skl (1:2,1,2,:)), der{2}(1:2,:), 1e3*eps) 159 160%!shared srf, uv, P, dPdx, d2Pdx2, c1, c2 161%!test 162%! [u, v] = meshgrid (linspace(0,1,10)); 163%! uv = [u(:)';v(:)']; 164%! c1 = nrbmak([0 1/2 1; 0 1 0],[0 0 0 1 1 1]); 165%! c1 = nrbtform (c1, vecrotx (pi/2)); 166%! c2 = nrbtform(c1, vectrans([0 1 0])); 167%! srf = nrbdegelev (nrbruled (c1, c2), [3, 1]); 168%! skl = nrbsurfderiveval (srf, uv, 2); 169%! P = squeeze(skl(:,1,1,:)); 170%! dPdx = squeeze(skl(:,2,1,:)); 171%! d2Pdx2 = squeeze(skl(:,3,1,:)); 172%!assert(P(3,:), 2*(P(1,:)-P(1,:).^2),100*eps) 173%!assert(dPdx(3,:), 2-4*P(1,:), 100*eps) 174%!assert(d2Pdx2(3,:), -4+0*P(1,:), 100*eps) 175%! srf = nrbdegelev (nrbruled (c1, c2), [5, 6]); 176%! skl = nrbsurfderiveval (srf, uv, 2); 177%! P = squeeze(skl(:,1,1,:)); 178%! dPdx = squeeze(skl(:,2,1,:)); 179%! d2Pdx2 = squeeze(skl(:,3,1,:)); 180%! aux = nrbeval(srf,uv); 181%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 182%!assert(P(3,:), 2*(P(1,:)-P(1,:).^2),100*eps) 183%!assert(dPdx(3,:), 2-4*P(1,:), 100*eps) 184%!assert(d2Pdx2(3,:), -4+0*P(1,:), 100*eps) 185%! 186%!test 187%! skl = nrbsurfderiveval (srf, uv, 0); 188%! aux = nrbeval (srf, uv); 189%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 190 191%!shared dPdu, d2Pdu2, P, srf, uv 192%!test 193%! [u, v] = meshgrid (linspace(0,1,10)); 194%! uv = [u(:)';v(:)']; 195%! c1 = nrbmak([0 1/2 1; 0.1 1.6 1.1; 0 0 0],[0 0 0 1 1 1]); 196%! c2 = nrbmak([0 1/2 1; 0.1 1.6 1.1; 1 1 1],[0 0 0 1 1 1]); 197%! srf = nrbdegelev (nrbruled (c1, c2), [0, 1]); 198%! skl = nrbsurfderiveval (srf, uv, 2); 199%! P = squeeze(skl(:,1,1,:)); 200%! dPdu = squeeze(skl(:,2,1,:)); 201%! dPdv = squeeze(skl(:,1,2,:)); 202%! d2Pdu2 = squeeze(skl(:,3,1,:)); 203%! aux = nrbeval(srf,uv); 204%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 205%!assert(dPdu(2,:), 3-4*P(1,:),100*eps) 206%!assert(d2Pdu2(2,:), -4+0*P(1,:),100*eps) 207%! 208%!test 209%! skl = nrbsurfderiveval (srf, uv, 0); 210%! aux = nrbeval(srf,uv); 211%! assert (squeeze (skl (1:2,1,1,:)), aux(1:2,:), 1e3*eps) 212 213%!test 214%! srf = nrb4surf([0 0], [1 0], [0 1], [1 1]); 215%! geo = nrbdegelev (srf, [3 3]); 216%1 geo.coefs (4, 2:end-1, 2:end-1) = geo.coefs (4, 2:end-1, 2:end-1) + .1 * rand (1, geo.number(1)-2, geo.number(2)-2); 217%! geo = nrbkntins (geo, {[.1:.1:.9], [.2:.2:.8]}); 218%! [u, v] = meshgrid (linspace(0,1,10)); 219%! uv = [u(:)';v(:)']; 220%! skl = nrbsurfderiveval (geo, uv, 2); 221%! dgeo = nrbderiv (geo); 222%! [pnts, ders] = nrbdeval (geo, dgeo, uv); 223%! assert (ders{1}, squeeze(skl(:,2,1,:)), 1e-9) 224%! assert (ders{2}, squeeze(skl(:,1,2,:)), 1e-9) 225 226%!test 227%! crv = nrbline ([1 0], [2 0]); 228%! srf = nrbrevolve (crv, [0 0 0], [0 0 1], pi/2); 229%! srf = nrbtransp (srf); 230%! [v, u] = meshgrid (linspace (0, 1, 11)); 231%! uv = [u(:)'; v(:)']; 232%! skl = nrbsurfderiveval (srf, uv, 2); 233%! c = sqrt(2); 234%! w = @(x, y) (2 - c)*y.^2 + (c-2)*y + 1; 235%! dwdy = @(x, y) 2*(2-c)*y + c - 2; 236%! d2wdy2 = @(x, y) 2*(2-c); 237%! F1 = @(x, y) (x+1) .* ((1-y).^2 + c*y.*(1-y)) ./ w(x,y); 238%! F2 = @(x, y) (x+1) .* (y.^2 + c*y.*(1-y)) ./ w(x,y); 239%! dF1dx = @(x, y) ((1-y).^2 + c*y.*(1-y)) ./ w(x,y); 240%! dF2dx = @(x, y) (y.^2 + c*y.*(1-y)) ./ w(x,y); 241%! dF1dy = @(x, y) (x+1) .* ((2 - 2*c)*y + c - 2) ./ w(x,y) - (x+1) .* ((1-y).^2 + c*y.*(1-y)) .* dwdy(x,y) ./ w(x,y).^2; 242%! dF2dy = @(x, y) (x+1) .* ((2 - 2*c)*y + c) ./ w(x,y) - (x+1) .* (y.^2 + c*y.*(1-y)) .* dwdy(x,y) ./ w(x,y).^2; 243%! d2F1dx2 = @(x, y) zeros (size (x)); 244%! d2F2dx2 = @(x, y) zeros (size (x)); 245%! d2F1dxdy = @(x, y) ((2 - 2*c)*y + c - 2) ./ w(x,y) - ((1-y).^2 + c*y.*(1-y)) .* dwdy(x,y) ./ w(x,y).^2; 246%! d2F2dxdy = @(x, y) ((2 - 2*c)*y + c) ./ w(x,y) - (y.^2 + c*y.*(1-y)) .* dwdy(x,y) ./ w(x,y).^2; 247%! d2F1dy2 = @(x, y) (x+1)*(2 - 2*c) ./ w(x,y) - 2*(x+1) .* ((2 - 2*c)*y + c - 2) .* dwdy(x,y) ./ w(x,y).^2 - ... 248%! (x+1) .* ((1-y).^2 + c*y.*(1-y)) * d2wdy2(x,y) ./ w(x,y).^2 + ... 249%! 2 * (x+1) .* ((1-y).^2 + c*y.*(1-y)) .* w(x,y) .*dwdy(x,y).^2 ./ w(x,y).^4; 250%! d2F2dy2 = @(x, y) (x+1)*(2 - 2*c) ./ w(x,y) - 2*(x+1) .* ((2 - 2*c)*y + c) .* dwdy(x,y) ./ w(x,y).^2 - ... 251%! (x+1) .* (y.^2 + c*y.*(1-y)) * d2wdy2(x,y) ./ w(x,y).^2 + ... 252%! 2 * (x+1) .* (y.^2 + c*y.*(1-y)) .* w(x,y) .*dwdy(x,y).^2 ./ w(x,y).^4; 253%! assert ([F1(u(:),v(:)), F2(u(:),v(:))], squeeze(skl(1:2,1,1,:))', 1e2*eps); 254%! assert ([dF1dx(u(:),v(:)), dF2dx(u(:),v(:))], squeeze(skl(1:2,2,1,:))', 1e2*eps); 255%! assert ([dF1dy(u(:),v(:)), dF2dy(u(:),v(:))], squeeze(skl(1:2,1,2,:))', 1e2*eps); 256%! assert ([d2F1dx2(u(:),v(:)), d2F2dx2(u(:),v(:))], squeeze(skl(1:2,3,1,:))', 1e2*eps); 257%! assert ([d2F1dxdy(u(:),v(:)), d2F2dxdy(u(:),v(:))], squeeze(skl(1:2,2,2,:))', 1e2*eps); 258%! assert ([d2F1dy2(u(:),v(:)), d2F2dy2(u(:),v(:))], squeeze(skl(1:2,1,3,:))', 1e2*eps); 259 260%!test 261%! knots = {[0 0 1 1] [0 0 1 1]}; 262%! coefs(:,1,1) = [0;0;0;1]; 263%! coefs(:,2,1) = [1;0;0;1]; 264%! coefs(:,1,2) = [0;1;0;1]; 265%! coefs(:,2,2) = [1;1;1;2]; 266%! srf = nrbmak (coefs, knots); 267%! [v, u] = meshgrid (linspace (0, 1, 3)); 268%! uv = [u(:)'; v(:)']; 269%! skl = nrbsurfderiveval (srf, uv, 2); 270%! w = @(x, y) x.*y + 1; 271%! F1 = @(x, y) x ./ w(x,y); 272%! F2 = @(x, y) y ./ w(x,y); 273%! F3 = @(x, y) x .* y ./ w(x,y); 274%! dF1dx = @(x, y) 1./w(x,y) - x.*y./w(x,y).^2; 275%! dF1dy = @(x, y) - x.^2./w(x,y).^2; 276%! dF2dx = @(x, y) - y.^2./w(x,y).^2; 277%! dF2dy = @(x, y) 1./w(x,y) - x.*y./w(x,y).^2; 278%! dF3dx = @(x, y) y./w(x,y) - x.*(y./w(x,y)).^2; 279%! dF3dy = @(x, y) x./w(x,y) - y.*(x./w(x,y)).^2; 280%! d2F1dx2 = @(x, y) -2*y./w(x,y).^2 + 2*x.*y.^2./w(x,y).^3; 281%! d2F1dy2 = @(x, y) 2*x.^3./w(x,y).^3; 282%! d2F1dxdy = @(x, y) -x./w(x,y).^2 - x./w(x,y).^2 + 2*x.^2.*y./w(x,y).^3; 283%! d2F2dx2 = @(x, y) 2*y.^3./w(x,y).^3; 284%! d2F2dy2 = @(x, y) -2*x./w(x,y).^2 + 2*y.*x.^2./w(x,y).^3; 285%! d2F2dxdy = @(x, y) -y./w(x,y).^2 - y./w(x,y).^2 + 2*y.^2.*x./w(x,y).^3; 286%! d2F3dx2 = @(x, y) -2*y.^2./w(x,y).^2 + 2*x.*y.^3./w(x,y).^3; 287%! d2F3dy2 = @(x, y) -2*x.^2./w(x,y).^2 + 2*y.*x.^3./w(x,y).^3; 288%! d2F3dxdy = @(x, y) 1./w(x,y) - 3*x.*y./w(x,y).^2 + 2*(x.*y).^2./w(x,y).^3; 289%! assert ([F1(u(:),v(:)), F2(u(:),v(:)), F3(u(:),v(:))], squeeze(skl(1:3,1,1,:))', 1e2*eps); 290%! assert ([dF1dx(u(:),v(:)), dF2dx(u(:),v(:)), dF3dx(u(:),v(:))], squeeze(skl(1:3,2,1,:))', 1e2*eps); 291%! assert ([dF1dy(u(:),v(:)), dF2dy(u(:),v(:)), dF3dy(u(:),v(:))], squeeze(skl(1:3,1,2,:))', 1e2*eps); 292%! assert ([d2F1dx2(u(:),v(:)), d2F2dx2(u(:),v(:)), d2F3dx2(u(:),v(:))], squeeze(skl(1:3,3,1,:))', 1e2*eps); 293%! assert ([d2F1dy2(u(:),v(:)), d2F2dy2(u(:),v(:)), d2F3dy2(u(:),v(:))], squeeze(skl(1:3,1,3,:))', 1e2*eps); 294%! assert ([d2F1dxdy(u(:),v(:)), d2F2dxdy(u(:),v(:)), d2F3dxdy(u(:),v(:))], squeeze(skl(1:3,2,2,:))', 1e2*eps); 295