octave_packages/nurbs-1.3.6/nrbsurfderiveval.m

   1 function 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 2 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
  77 end
  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);