[Pkg-octave-commit] [SCM] UNNAMED PROJECT branch, master, updated. 9e4ab9ff627998fc000541cb525ab88e146843f9
Rafael Laboissiere
rafael at laboissiere.net
Fri Mar 30 15:32:15 UTC 2012
The following commit has been merged in the master branch:
commit 21ae4a2a3b932f6ab08c907600626c1702bde45a
Author: Rafael Laboissiere <rafael at laboissiere.net>
Date: Sun Mar 18 12:24:13 2012 +0100
Imported Upstream version 1.3.6
diff --git a/DESCRIPTION b/DESCRIPTION
index ce4e260..36b13f8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Name: Nurbs
-Version: 1.3.5
-Date: 2011-08-26
+Version: 1.3.6
+Date: 2012-03-17
Author: Mark Spink, Daniel Claxton, Carlo de Falco, Rafael Vazquez
Maintainer: Carlo de Falco and Rafael Vazquez
Title: Nurbs.
@@ -10,5 +10,5 @@ Depends: octave (>= 3.2)
Autoload: no
License: GPL version 2 or later
Url: http://octave.sf.net
-SVNRelease: 8477
+SVNRelease: 9939
diff --git a/INDEX b/INDEX
index 2c736f9..a4f1138 100644
--- a/INDEX
+++ b/INDEX
@@ -62,6 +62,7 @@ Vector and Transformation Utilities
vecrotx
vecroty
vecrotz
+ vecrot
vecscale
vectrans
Misc Utilities
diff --git a/inst/basisfunder.m b/inst/basisfunder.m
index ba320f9..c6f0a89 100644
--- a/inst/basisfunder.m
+++ b/inst/basisfunder.m
@@ -20,7 +20,7 @@ function dersv = basisfunder (ii, pl, uu, u_knotl, nders)
%
% Adapted from Algorithm A2.3 from 'The NURBS BOOK' pg72.
%
-% Copyright (C) 2009 Rafael Vazquez
+% Copyright (C) 2009,2011 Rafael Vazquez
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
@@ -35,6 +35,8 @@ function dersv = basisfunder (ii, pl, uu, u_knotl, nders)
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
+ dersv = zeros(numel(uu), nders+1, pl+1);
+
for jj = 1:numel(uu)
i = ii(jj)+1; %% convert to base-1 numbering of knot spans
diff --git a/inst/kntuniform.m b/inst/kntuniform.m
index 6805ca4..2e6df46 100644
--- a/inst/kntuniform.m
+++ b/inst/kntuniform.m
@@ -14,6 +14,7 @@
% zeta: breaks = knots without repetitions
%
% Copyright (C) 2009, 2010 Carlo de Falco
+% Copyright (C) 2011 Rafael Vazquez
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
@@ -43,5 +44,9 @@ function [csi, zeta] = kntuniform (num, degree, regularity)
error ('kntuniform: regularity requested is too high')
end
end
+ if (numel(num) == 1)
+ csi = csi{1};
+ zeta = zeta{1};
+ end
end
end
diff --git a/inst/nrbexport.m b/inst/nrbexport.m
index 1644730..b53057c 100644
--- a/inst/nrbexport.m
+++ b/inst/nrbexport.m
@@ -1,79 +1,80 @@
-function nrbexport (nurbs, filename)
-
-%
-% NRBEXPORT: export NURBS geometries to a format compatible with the one used in GeoPDEs (version 0.6).
-%
-% Calling Sequence:
-%
-% nrbexport (nurbs, filename);
-%
-% INPUT:
-%
-% nurbs : NURBS curve, surface or volume, see nrbmak.
-% filename : name of the output file.
-%
-%
-% Description:
-%
-% The data of the nurbs structure is written in the file, in a format
-% that can be read by GeoPDEs.
-%
-% Copyright (C) 2011 Rafael Vazquez
-%
-% This program is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 2 of the License, or
-% (at your option) any later version.
-
-% This program is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-fid = fopen (filename, 'w');
-if (fid < 0)
- error ('nrbexport: cannot open file %s', filename);
-end
-
-ndim = numel (nurbs(1).order);
-npatch = numel (nurbs);
-fprintf (fid, '%s\n', '# nurbs mesh v.0.7');
-fprintf (fid, '%s\n', '#');
-fprintf (fid, '%s\n', ['# ' date]);
-fprintf (fid, '%s\n', '#');
-
-fprintf (fid, '%2i', ndim, 1);
-fprintf (fid, '\n');
-for iptc = 1:npatch
- fprintf (fid, '%s %i', 'PATCH', iptc);
- fprintf (fid, '\n');
- fprintf (fid, '%2i', nurbs(iptc).order-1);
- fprintf (fid, '\n');
- fprintf (fid, '%2i', nurbs(iptc).number);
- fprintf (fid, '\n');
- for ii = 1:ndim
- fprintf (fid, '%1.7f ', nurbs(iptc).knots{ii});
- fprintf (fid, '\n');
- end
-
- if (ndim == 2)
- for ii = 1:ndim
- fprintf (fid, '%1.15f ', nurbs(iptc).coefs(ii,:,:));
- fprintf (fid, '\n');
- end
- fprintf (fid, '%1.15f ', nurbs(iptc).coefs(4,:,:));
- fprintf (fid, '\n');
- elseif (ndim == 3)
- for ii = 1:ndim
- fprintf (fid, '%1.15f ', nurbs(iptc).coefs(ii,:,:,:));
- fprintf (fid, '\n');
- end
- fprintf (fid, '%1.15f ', nurbs(iptc).coefs(4,:,:,:));
- fprintf (fid, '\n');
- end
-end
-
-end
+function nrbexport (nurbs, filename)
+
+%
+% NRBEXPORT: export NURBS geometries to a format compatible with the one used in GeoPDEs (version 0.6).
+%
+% Calling Sequence:
+%
+% nrbexport (nurbs, filename);
+%
+% INPUT:
+%
+% nurbs : NURBS curve, surface or volume, see nrbmak.
+% filename : name of the output file.
+%
+%
+% Description:
+%
+% The data of the nurbs structure is written in the file, in a format
+% that can be read by GeoPDEs.
+%
+% Copyright (C) 2011 Rafael Vazquez
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 2 of the License, or
+% (at your option) any later version.
+
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+fid = fopen (filename, 'w');
+if (fid < 0)
+ error ('nrbexport: cannot open file %s', filename);
+end
+
+ndim = numel (nurbs(1).order);
+npatch = numel (nurbs);
+fprintf (fid, '%s\n', '# nurbs mesh v.0.7');
+fprintf (fid, '%s\n', '#');
+fprintf (fid, '%s\n', ['# ' date]);
+fprintf (fid, '%s\n', '#');
+
+fprintf (fid, '%2i', ndim, 1);
+fprintf (fid, '\n');
+for iptc = 1:npatch
+ fprintf (fid, '%s %i', 'PATCH', iptc);
+ fprintf (fid, '\n');
+ fprintf (fid, '%2i', nurbs(iptc).order-1);
+ fprintf (fid, '\n');
+ fprintf (fid, '%2i', nurbs(iptc).number);
+ fprintf (fid, '\n');
+ for ii = 1:ndim
+ fprintf (fid, '%1.7f ', nurbs(iptc).knots{ii});
+ fprintf (fid, '\n');
+ end
+
+ if (ndim == 2)
+ for ii = 1:ndim
+ fprintf (fid, '%1.15f ', nurbs(iptc).coefs(ii,:,:));
+ fprintf (fid, '\n');
+ end
+ fprintf (fid, '%1.15f ', nurbs(iptc).coefs(4,:,:));
+ fprintf (fid, '\n');
+ elseif (ndim == 3)
+ for ii = 1:ndim
+ fprintf (fid, '%1.15f ', nurbs(iptc).coefs(ii,:,:,:));
+ fprintf (fid, '\n');
+ end
+ fprintf (fid, '%1.15f ', nurbs(iptc).coefs(4,:,:,:));
+ fprintf (fid, '\n');
+ end
+end
+
+fclose (fid);
+end
diff --git a/inst/private/onebasisfun__.m b/inst/private/onebasisfun__.m
index 3fd962b..bbd1764 100644
--- a/inst/private/onebasisfun__.m
+++ b/inst/private/onebasisfun__.m
@@ -1,30 +1,56 @@
-function N = onebasisfun__ (u, p, U)
+function Nip = onebasisfun__ (u, p, U)
% __ONEBASISFUN__: Undocumented internal function
%
+% Adapted from Algorithm A2.4 from 'The NURBS BOOK' pg74.
+%
% Copyright (C) 2009 Carlo de Falco
+% Copyright (C) 2012 Rafael Vazquez
% This software comes with ABSOLUTELY NO WARRANTY; see the file
% COPYING for details. This is free software, and you are welcome
% to distribute it under the conditions laid out in COPYING.
- N = 0;
- if (~ any (U <= u)) || (~ any (U > u))
- return;
- elseif (p == 0)
- N = 1;
- return;
+ Nip = zeros (size (u));
+ N = zeros (p+1, 1);
+
+ for ii = 1:numel(u)
+ if ((u(ii) == U(1)) && (U(1) == U(end-1)) || ...
+ (u(ii) == U(end)) && (U(end) == U(2)))
+ Nip(ii) = 1;
+ continue
+ end
+ if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
+ continue;
+ end
+ for jj = 1:p+1 % Initialize zero-th degree functions
+ if (u(ii) > U(jj) && u(ii) < U(jj+1))
+ N(jj) = 1;
+ else
+ N(jj) = 0;
+ end
+ end
+ for k = 1:p
+ if (N(1) == 0)
+ saved = 0;
+ else
+ saved = (u(ii) - U(1))*N(1) / (U(k+1)-U(1));
+ end
+
+ for jj = 1:p-k+1
+ Uleft = U(1+jj);
+ Uright = U(1+jj+k);
+ if (N(jj+1) == 0)
+ N(jj) = saved;
+ saved = 0;
+ else
+ temp = N(jj+1)/(Uright-Uleft);
+ N(jj) = saved + (Uright - u(ii))*temp;
+ saved = (u(ii) - Uleft)*temp;
+ end
+ end
+ end
+ Nip(ii) = N(1);
end
- ln = u - U(1);
- ld = U(end-1) - U(1);
- if (ld ~= 0)
- N = N + ln * onebasisfun__ (u, p-1, U(1:end-1))/ ld;
- end
- dn = U(end) - u;
- dd = U(end) - U(2);
- if (dd ~= 0)
- N = N + dn * onebasisfun__ (u, p-1, U(2:end))/ dd;
- end
-
end
diff --git a/inst/private/onebasisfunder__.m b/inst/private/onebasisfunder__.m
new file mode 100644
index 0000000..b397d10
--- /dev/null
+++ b/inst/private/onebasisfunder__.m
@@ -0,0 +1,39 @@
+function [N, Nder] = onebasisfunder__ (u, p, U)
+
+% __ONEBASISFUNDER__: Undocumented internal function
+%
+% Copyright (C) 2012 Rafael Vazquez
+% This software comes with ABSOLUTELY NO WARRANTY; see the file
+% COPYING for details. This is free software, and you are welcome
+% to distribute it under the conditions laid out in COPYING.
+
+ N = zeros (size (u));
+ Nder = zeros (size (u));
+
+ for ii = 1:numel (u)
+ if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
+ continue;
+ elseif (p == 0)
+ N(ii) = 1;
+ Nder(ii) = 0;
+ continue;
+ else
+ ln = u(ii) - U(1);
+ ld = U(end-1) - U(1);
+ if (ld ~= 0)
+ aux = onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld;
+ N(ii) = N(ii) + ln * aux;
+ Nder(ii) = Nder(ii) + p * aux;
+ end
+
+ dn = U(end) - u(ii);
+ dd = U(end) - U(2);
+ if (dd ~= 0)
+ aux = onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
+ N(ii) = N(ii) + dn * aux;
+ Nder(ii) = Nder(ii) - p * aux;
+ end
+ end
+ end
+
+end
diff --git a/inst/tbasisfun.m b/inst/tbasisfun.m
index 408d64e..10d3d80 100644
--- a/inst/tbasisfun.m
+++ b/inst/tbasisfun.m
@@ -1,12 +1,12 @@
-function N = tbasisfun (u, p, U)
+function [N, Nder] = tbasisfun (u, p, U)
%
-% TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.
+% TBASISFUN: Compute a B- or T-Spline basis function, and its derivatives, from its local knot vector.
%
% usage:
%
-% N = tbasisfun (u, p, U)
-% N = tbasisfun ([u; v], [p q], {U, V})
-% N = tbasisfun ([u; v; w], [p q r], {U, V, W})
+% [N, Nder] = tbasisfun (u, p, U)
+% [N, Nder] = tbasisfun ([u; v], [p q], {U, V})
+% [N, Nder] = tbasisfun ([u; v; w], [p q r], {U, V, W})
%
% INPUT:
%
@@ -19,9 +19,11 @@ function N = tbasisfun (u, p, U)
%
% OUTPUT:
%
-% N : basis function evaluated at the given parametric points
+% N : basis function evaluated at the given parametric points
+% Nder : basis function gradient evaluated at the given parametric points
%
% Copyright (C) 2009 Carlo de Falco
+% Copyright (C) 2012 Rafael Vazquez
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
@@ -38,46 +40,58 @@ function N = tbasisfun (u, p, U)
if (~ iscell (U))
U = sort (U);
- assert (numel (U) == p+2)
-
- N = zeros(1,numel(u));
- for ii=1:numel(u)
- N(ii) = onebasisfun__ (u(ii), p, U);
+ if (numel (U) ~= p+2)
+ error ('tbasisfun: knot vector and degree do not correspond')
end
-
- elseif size(U,2) == 2
- U{1} = sort(U{1}); U{2} = sort(U{2});
- assert (numel(U{1}) == p(1)+2 && numel(U{2}) == p(2)+2)
- Nu = zeros(1,numel(u(1,:))); Nv = zeros(1,numel(u(1,:)));
- for ii=1:numel(u(1,:))
- Nu(ii) = onebasisfun__ (u(1,ii), p(1), U{1});
+ if (nargout == 1)
+ N = onebasisfun__ (u, p, U);
+ else
+ [N, Nder] = onebasisfunder__ (u, p, U);
end
-
- for ii=1:numel(u(1,:))
- Nv(ii) = onebasisfun__ (u(2,ii), p(2), U{2});
+
+ elseif (size(U,2) == 2)
+ U{1} = sort(U{1}); U{2} = sort(U{2});
+ if (numel(U{1}) ~= p(1)+2 || numel(U{2}) ~= p(2)+2)
+ error ('tbasisfun: knot vector and degree do not correspond')
end
+
+ if (nargout == 1)
+ Nu = onebasisfun__ (u(1,:), p(1), U{1});
+ Nv = onebasisfun__ (u(2,:), p(2), U{2});
- N = Nu.*Nv;
+ N = Nu.*Nv;
+ elseif (nargout == 2)
+ [Nu, Ndu] = onebasisfunder__ (u(1,:), p(1), U{1});
+ [Nv, Ndv] = onebasisfunder__ (u(2,:), p(2), U{2});
- elseif size(U,2) == 3
- U{1} = sort(U{1}); U{2} = sort(U{2}); U{3} = sort(U{3});
- assert (numel(U{1}) == p(1)+2 && numel(U{2}) == p(2)+2 && numel(U{3}) == p(3)+2)
+ N = Nu.*Nv;
+ Nder(1,:) = Ndu.*Nv;
+ Nder(2,:) = Nu.*Ndv;
+ end
- Nu = zeros(1,numel(u(1,:))); Nv = zeros(1,numel(u(1,:))); Nw = zeros(1,numel(u(1,:)));
- for ii=1:numel(u(1,:))
- Nu(ii) = onebasisfun__ (u(1,ii), p(1), U{1});
+ elseif (size(U,2) == 3)
+ U{1} = sort(U{1}); U{2} = sort(U{2}); U{3} = sort(U{3});
+ if (numel(U{1}) ~= p(1)+2 || numel(U{2}) ~= p(2)+2 || numel(U{3}) ~= p(3)+2)
+ error ('tbasisfun: knot vector and degree do not correspond')
end
- for ii=1:numel(u(1,:))
- Nv(ii) = onebasisfun__ (u(2,ii), p(2), U{2});
- end
+ if (nargout == 1)
+ Nu = onebasisfun__ (u(1,:), p(1), U{1});
+ Nv = onebasisfun__ (u(2,:), p(2), U{2});
+ Nw = onebasisfun__ (u(3,:), p(3), U{3});
- for ii=1:numel(u(1,:))
- Nw(ii) = onebasisfun__ (u(3,ii), p(3), U{3});
+ N = Nu.*Nv.*Nw;
+ else
+ [Nu, Ndu] = onebasisfunder__ (u(1,:), p(1), U{1});
+ [Nv, Ndv] = onebasisfunder__ (u(2,:), p(2), U{2});
+ [Nw, Ndw] = onebasisfunder__ (u(3,:), p(3), U{3});
+
+ N = Nu.*Nv.*Nw;
+ Nder(1,:) = Ndu.*Nv.*Nw;
+ Nder(2,:) = Nu.*Ndv.*Nw;
+ Nder(3,:) = Nu.*Nv.*Ndw;
end
-
- N = Nu.*Nv.*Nw;
end
end
@@ -91,3 +105,27 @@ end
%! surf (X, Y, reshape (N, size(X)))
%! title('Basis function associated to a local knot vector')
%! hold off
+
+%!test
+%! U = [0 1/2 1];
+%! p = 1;
+%! u = [0.3 0.4 0.6 0.7];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.6 0.8 0.8 0.6], 1e-12);
+%! assert (Nder, [2 2 -2 -2], 1e-12);
+
+%!test
+%! U = {[0 1/2 1] [0 1/2 1]};
+%! p = [1 1];
+%! u = [0.3 0.4 0.6 0.7; 0.3 0.4 0.6 0.7];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.36 0.64 0.64 0.36], 1e-12);
+%! assert (Nder, [1.2 1.6 -1.6 -1.2; 1.2 1.6 -1.6 -1.2], 1e-12);
+
+%!test
+%! U = {[0 1/2 1] [0 1/2 1] [0 1/2 1]};
+%! p = [1 1 1];
+%! u = [0.4 0.4 0.6 0.6; 0.4 0.4 0.6 0.6; 0.4 0.6 0.4 0.6];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.512 0.512 0.512 0.512], 1e-12);
+%! assert (Nder, [1.28 1.28 -1.28 -1.28; 1.28 1.28 -1.28 -1.28; 1.28 -1.28 1.28 -1.28], 1e-12);
diff --git a/inst/vecrotx.m b/inst/vecrot.m
similarity index 58%
copy from inst/vecrotx.m
copy to inst/vecrot.m
index 8b03721..d8ed20d 100644
--- a/inst/vecrotx.m
+++ b/inst/vecrot.m
@@ -1,14 +1,15 @@
-function rx = vecrotx(angle)
+function rx = vecrot(angle, vector)
%
-% VECROTX: Transformation matrix for a rotation around the x axis.
+% VECROT: Transformation matrix for a rotation around the axis given by a vector.
%
% Calling Sequence:
%
-% rx = vecrotx(angle);
+% rx = vecrot (angle, vector);
%
% INPUT:
%
% angle : rotation angle defined in radians
+% vector : vector defining the rotation axis
%
% OUTPUT:
%
@@ -20,27 +21,11 @@ function rx = vecrotx(angle)
% Return the (4x4) Transformation matrix for a rotation about the x axis
% by the defined angle.
%
-% The matrix is:
-%
-% [ 1 0 0 0]
-% [ 0 cos(angle) -sin(angle) 0]
-% [ 0 sin(angle) cos(angle) 0]
-% [ 0 0 0 1]
-%
-% Examples:
-%
-% Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees
-% around the x-axis
-%
-% line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]);
-% trans = vecrotx(%pi/4);
-% rline = nrbtform(line, trans);
-%
% See also:
%
% nrbtform
%
-% Copyright (C) 2000 Mark Spink
+% Copyright (C) 2011 Rafael Vazquez
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
@@ -55,8 +40,14 @@ function rx = vecrotx(angle)
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
-sn = sin(angle);
-cn = cos(angle);
-rx = [1 0 0 0; 0 cn -sn 0; 0 sn cn 0; 0 0 0 1];
+% Normalize the vector
+vec = vector / norm (vector);
+
+sn = sin (angle);
+cn = cos (angle);
+rx = [cn+vec(1)^2*(1-cn), vec(1)*vec(2)*(1-cn)-vec(3)*sn, vec(1)*vec(3)*(1-cn)+vec(2)*sn, 0;
+ vec(1)*vec(2)*(1-cn)+vec(3)*sn, cn+vec(2)^2*(1-cn), vec(2)*vec(3)*(1-cn)-vec(1)*sn, 0;
+ vec(1)*vec(3)*(1-cn)-vec(2)*sn, vec(2)*vec(3)*(1-cn)+vec(1)*sn, cn+vec(3)^2*(1-cn), 0;
+ 0 0 0 1];
end
diff --git a/inst/vectrans.m b/inst/vectrans.m
index 6faf64b..91a634c 100644
--- a/inst/vectrans.m
+++ b/inst/vectrans.m
@@ -22,8 +22,8 @@ function dd = vectrans(vector)
% The matrix is:
%
% [ 1 0 0 tx ]
-% [ 0 0 0 ty ]
-% [ 0 0 0 tz ]
+% [ 0 1 0 ty ]
+% [ 0 0 1 tz ]
% [ 0 0 0 1 ]
%
% Examples:
diff --git a/src/tbasisfun.cc b/src/tbasisfun.cc
index 5d33b17..a17e324 100644
--- a/src/tbasisfun.cc
+++ b/src/tbasisfun.cc
@@ -1,4 +1,5 @@
/* Copyright (C) 2009 Carlo de Falco
+ Copyright (C) 2012 Rafael Vazquez
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@@ -17,38 +18,123 @@
#include <octave/oct.h>
#include <iostream>
-double onebasisfun__ (double u, octave_idx_type p, RowVector U)
+void onebasisfun__ (double u, octave_idx_type p, RowVector U, double *N)
{
-
- //std::cout << "u=" << u << " " << "p=" << p << " \n" << "U=" << U;
-
- double N = 0.0;
- if ((u < U.min ()) || ( u > U.max ()))
- return (N);
+ *N = 0.0;
+ if ((u <= U.min ()) || ( u > U.max ()))
+ return;
else if (p == 0)
- return (1.0);
+ {
+ *N = 1.0;
+ return;
+ }
+ else if (p == 1)
+ {
+ if (u < U(1))
+ {
+ *N = (u - U(0)) / (U(1) - U(0));
+ return;
+ }
+ else
+ {
+ *N = (U(2) - u) / (U(2) - U(1));
+ return;
+ }
+ }
+ else if (p == 2)
+ {
+ double ln = u - U(0);
+ double dn = U(3) - u;
+ double ld = U(2) - U(0);
+ double dd = U(3) - U(1);
+ if (u < U(1))
+ {
+ *N = ln*ln / (ld * (U(1) - U(0)));
+ return;
+ }
+ else if (u > U(2))
+ {
+ *N = dn*dn / (dd * (U(3) - U(2)));
+ return;
+ }
+ else
+ {
+ if (ld != 0)
+ *N += ln * (U(2) - u) / ((U(2) - U(1)) * ld);
+
+ if (dd != 0)
+ *N += dn * (u - U(1)) / ((U(2) - U(1)) * dd);
+ return;
+ }
+ }
double ln = u - U(0);
double ld = U(U.length () - 2) - U(0);
if (ld != 0)
- N += ln * onebasisfun__ (u, p-1, U.extract (0, U.length () - 2))/ ld;
-
+ {
+ double tmp;
+ onebasisfun__ (u, p-1, U.extract (0, U.length () - 2), &tmp);
+ *N += ln * tmp / ld;
+ }
double dn = U(U.length () - 1) - u;
double dd = U(U.length () - 1) - U(1);
if (dd != 0)
- N += dn * onebasisfun__ (u, p-1, U.extract (1, U.length () - 1))/ dd;
-
- return (N);
+ {
+ double tmp;
+ onebasisfun__ (u, p-1, U.extract (1, U.length () - 1), &tmp);
+ *N += dn * tmp / dd;
+ }
+ return;
}
+
+void onebasisfunder__ (double u, octave_idx_type p, RowVector U, double *N, double *Nder)
+{
+ double aux;
+ *N = 0.0; *Nder = 0.0;
+ if ((u <= U.min ()) || ( u > U.max ()))
+ return;
+
+ else if (p == 0)
+ {
+ *N = 1.0;
+ *Nder = 0.0;
+ return;
+ }
+
+ else {
+
+ double ln = u - U(0);
+ double ld = U(U.length () - 2) - U(0);
+
+ if (ld != 0)
+ {
+ onebasisfun__ (u, p-1, U.extract (0, U.length () - 2), &aux);
+ aux = aux / ld;
+ *N += ln * aux;
+ *Nder += p * aux;
+ }
+
+ double dn = U(U.length () - 1) - u;
+ double dd = U(U.length () - 1) - U(1);
+ if (dd != 0)
+ {
+ onebasisfun__ (u, p-1, U.extract (1, U.length () - 1), &aux);
+ aux = aux / dd;
+ *N += dn *aux;
+ *Nder -= p * aux;
+ }
+ }
+}
DEFUN_DLD(tbasisfun, args, nargout,"\
-TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.\n\
+TBASISFUN: Compute a B- or T-Spline basis function, and its derivatives, from its local knot vector.\n\
\n\
usage:\n\
\n\
- N = tbasisfun (u, p, U)\n\
- N = tbasisfun ([u; v], [p q], {U, V})\n\
+ [N, Nder] = tbasisfun (u, p, U)\n\
+ [N, Nder] = tbasisfun ([u; v], [p q], {U, V})\n\
+ [N, Nder] = tbasisfun ([u; v; w], [p q r], {U, V, W})\n\
\n\
INPUT:\n\
u or [u; v] : points in parameter space where the basis function is to be\n\
@@ -56,16 +142,20 @@ TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.\n
\n\
U or {U, V} : local knot vector\n\
\n\
- p or [p q] : polynomial order of the basis function\n\
+ p or [p q] : polynomial order of the basis function\n\
\n\
OUTPUT:\n\
- N : basis function evaluated at the given parametric points\n")
+ N : basis function evaluated at the given parametric points\n\
+ Nder : gradient of the basis function evaluated at the given points\n")
{
octave_value_list retval;
Matrix u = args(0).matrix_value ();
+ RowVector N(u.cols ());
+ double *Nptr = N.fortran_vec ();
+
if (! args(2).is_cell ())
{
@@ -73,32 +163,102 @@ TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.\n
RowVector U = args(2).row_vector_value (true, true);
assert (U.numel () == p+2);
- RowVector N(u.cols ());
- for (octave_idx_type ii=0; ii<u.numel (); ii++)
- N(ii) = onebasisfun__ (u(ii), p, U);
+ if (nargout == 1)
+ for (octave_idx_type ii = 0; ii < u.numel (); ii++)
+ onebasisfun__ (u(ii), p, U, &(Nptr[ii]));
- } else {
+ if (nargout == 2)
+ {
+ RowVector Nder(u.cols ());
+ double *Nderptr = Nder.fortran_vec ();
- RowVector p = args(1).row_vector_value ();
- Cell C = args(2).cell_value ();
- RowVector U = C(0).row_vector_value (true, true);
- RowVector V = C(1).row_vector_value (true, true);
-
+ for (octave_idx_type ii=0; ii<u.numel (); ii++)
+ onebasisfunder__ (u(ii), p, U, &(Nptr[ii]), &(Nderptr[ii]));
- RowVector N(u.cols ());
- for (octave_idx_type ii=0; ii<u.cols (); ii++)
- {
- N(ii) = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U) *
- onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V);
- //std::cout << "N=" << N(ii) << "\n\n\n";
+ retval(1) = Nder;
+ }
+ }
+ else
+ {
+ RowVector p = args(1).row_vector_value ();
+ if (p.length() == 2)
+ {
+ Cell C = args(2).cell_value ();
+ RowVector U = C(0).row_vector_value (true, true);
+ RowVector V = C(1).row_vector_value (true, true);
+
+ if (nargout == 1)
+ {
+ for (octave_idx_type ii=0; ii<u.cols (); ii++)
+ {
+ double Nu, Nv;
+ onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U, &Nu);
+ onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V, &Nv);
+ Nptr[ii] = Nu * Nv;
+ }
+ }
+ else if (nargout == 2)
+ {
+ double Nu, Nv, Ndu, Ndv;
+ Matrix Nder (2, u.cols());
+ double *Nderptr = Nder.fortran_vec ();
+ for (octave_idx_type ii = 0; ii < u.cols (); ii++)
+ {
+ onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U, &Nu, &Ndu);
+ onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V, &Nv, &Ndv);
+ Nptr[ii] = Nu * Nv;
+
+ Nderptr[0 + (2 * ii)] = Ndu * Nv;
+ Nderptr[1 + (2 * ii)] = Ndv * Nu;
+ }
+ retval(1) = Nder;
+ }
+
+ }
+ else if (p.length() == 3)
+ {
+ Cell C = args(2).cell_value ();
+ RowVector U = C(0).row_vector_value (true, true);
+ RowVector V = C(1).row_vector_value (true, true);
+ RowVector W = C(2).row_vector_value (true, true);
+
+ if (nargout == 1)
+ {
+ for (octave_idx_type ii = 0; ii < u.cols (); ii++)
+ {
+ double Nu, Nv, Nw;
+ onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U, &Nu);
+ onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V, &Nv);
+ onebasisfun__ (u(2, ii), octave_idx_type(p(2)), W, &Nw);
+ Nptr[ii] = Nu * Nv * Nw;
+ }
+ }
+ else if (nargout == 2)
+ {
+ double Nu, Nv, Nw, Ndu, Ndv, Ndw;
+ Matrix Nder (3, u.cols());
+ double *Nderptr = Nder.fortran_vec ();
+ for (octave_idx_type ii=0; ii<u.cols (); ii++)
+ {
+ onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U, &Nu, &Ndu);
+ onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V, &Nv, &Ndv);
+ onebasisfunder__ (u(2, ii), octave_idx_type(p(2)), W, &Nw, &Ndw);
+ Nptr[ii] = Nu * Nv * Nw;
+ Nderptr[0 + (3 * ii)] = Ndu * Nv * Nw;
+ Nderptr[1 + (3 * ii)] = Ndv * Nu * Nw;
+ Nderptr[2 + (3 * ii)] = Ndw * Nu * Nv;
+ }
+ retval(1) = Nder;
+ }
+
}
- retval(0) = octave_value (N);
- }
+ }
+ retval(0) = octave_value (N);
return retval;
}
-
/*
+
%!demo
%! U = {[0 0 1/2 1 1], [0 0 0 1 1]};
%! p = [3, 3];
@@ -106,4 +266,31 @@ TBASISFUN: Compute a B- or T-Spline basis function from its local knot vector.\n
%! u = [X(:), Y(:)]';
%! N = tbasisfun (u, p, U);
%! surf (X, Y, reshape (N, size(X)))
+%! title('Basis function associated to a local knot vector')
+%! hold off
+
+%!test
+%! U = [0 1/2 1];
+%! p = 1;
+%! u = [0.3 0.4 0.6 0.7];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.6 0.8 0.8 0.6], 1e-12);
+%! assert (Nder, [2 2 -2 -2], 1e-12);
+
+%!test
+%! U = {[0 1/2 1] [0 1/2 1]};
+%! p = [1 1];
+%! u = [0.3 0.4 0.6 0.7; 0.3 0.4 0.6 0.7];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.36 0.64 0.64 0.36], 1e-12);
+%! assert (Nder, [1.2 1.6 -1.6 -1.2; 1.2 1.6 -1.6 -1.2], 1e-12);
+
+%!test
+%! U = {[0 1/2 1] [0 1/2 1] [0 1/2 1]};
+%! p = [1 1 1];
+%! u = [0.4 0.4 0.6 0.6; 0.4 0.4 0.6 0.6; 0.4 0.6 0.4 0.6];
+%! [N, Nder] = tbasisfun (u, p, U);
+%! assert (N, [0.512 0.512 0.512 0.512], 1e-12);
+%! assert (Nder, [1.28 1.28 -1.28 -1.28; 1.28 1.28 -1.28 -1.28; 1.28 -1.28 1.28 -1.28], 1e-12);
+
*/
--
UNNAMED PROJECT
More information about the Pkg-octave-commit
mailing list