Mercurial > forge
annotate extra/nurbs/src/tbasisfun.cc @ 9403:2108d03d50de octave-forge
Fixed bug in 1D, added 3D case, added derivatives
author | rafavzqz |
---|---|
date | Mon, 06 Feb 2012 13:48:17 +0000 |
parents | 937278c100fd |
children | 869d0f838b63 |
rev | line source |
---|---|
6958 | 1 /* Copyright (C) 2009 Carlo de Falco |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
2 Copyright (C) 2012 Rafael Vazquez |
6958 | 3 |
4 This program is free software: you can redistribute it and/or modify | |
5 it under the terms of the GNU General Public License as published by | |
6 the Free Software Foundation, either version 2 of the License, or | |
7 (at your option) any later version. | |
6052 | 8 |
6958 | 9 This program is distributed in the hope that it will be useful, |
10 but WITHOUT ANY WARRANTY; without even the implied warranty of | |
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
12 GNU General Public License for more details. | |
6052 | 13 |
6958 | 14 You should have received a copy of the GNU General Public License |
15 along with this program. If not, see <http://www.gnu.org/licenses/>. | |
6052 | 16 */ |
17 | |
18 #include <octave/oct.h> | |
19 #include <iostream> | |
20 | |
21 double onebasisfun__ (double u, octave_idx_type p, RowVector U) | |
22 { | |
23 | |
24 //std::cout << "u=" << u << " " << "p=" << p << " \n" << "U=" << U; | |
25 | |
26 double N = 0.0; | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
27 if ((u <= U.min ()) || ( u > U.max ())) |
6052 | 28 return (N); |
29 else if (p == 0) | |
30 return (1.0); | |
31 | |
32 double ln = u - U(0); | |
33 double ld = U(U.length () - 2) - U(0); | |
34 if (ld != 0) | |
35 N += ln * onebasisfun__ (u, p-1, U.extract (0, U.length () - 2))/ ld; | |
36 | |
37 double dn = U(U.length () - 1) - u; | |
38 double dd = U(U.length () - 1) - U(1); | |
39 if (dd != 0) | |
40 N += dn * onebasisfun__ (u, p-1, U.extract (1, U.length () - 1))/ dd; | |
41 | |
42 return (N); | |
43 } | |
44 | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
45 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
46 double onebasisfunder__ (double u, octave_idx_type p, RowVector U) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
47 { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
48 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
49 //std::cout << "u=" << u << " " << "p=" << p << " \n" << "U=" << U; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
50 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
51 double N = 0.0; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
52 if ((u <= U.min ()) || ( u > U.max ())) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
53 return (N); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
54 else if (p == 0) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
55 return (0.0); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
56 else { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
57 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
58 double ld = U(U.length () - 2) - U(0); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
59 if (ld != 0) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
60 N += p * onebasisfun__ (u, p-1, U.extract (0, U.length () - 2))/ ld; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
61 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
62 double dd = U(U.length () - 1) - U(1); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
63 if (dd != 0) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
64 N -= p * onebasisfun__ (u, p-1, U.extract (1, U.length () - 1))/ dd; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
65 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
66 return (N); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
67 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
68 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
69 |
6052 | 70 |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
71 DEFUN_DLD(tbasisfunder, args, nargout,"\ |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
72 TBASISFUN: Compute a B- or T-Spline basis function, and its derivatives, from its local knot vector.\n\ |
6464 | 73 \n\ |
6052 | 74 usage:\n\ |
75 \n\ | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
76 [N, Nder] = tbasisfun (u, p, U)\n\ |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
77 [N, Nder] = tbasisfun ([u; v], [p q], {U, V})\n\ |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
78 [N, Nder] = tbasisfun ([u; v; w], [p q r], {U, V, W})\n\ |
6052 | 79 \n\ |
80 INPUT:\n\ | |
81 u or [u; v] : points in parameter space where the basis function is to be\n\ | |
82 evaluated \n\ | |
83 \n\ | |
84 U or {U, V} : local knot vector\n\ | |
85 \n\ | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
86 p or [p q] : polynomial order of the basis function\n\ |
6052 | 87 \n\ |
88 OUTPUT:\n\ | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
89 N : basis function evaluated at the given parametric points\n\ |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
90 Nder : gradient of the basis function evaluated at the given points\n") |
6052 | 91 |
92 { | |
93 | |
7204 | 94 octave_value_list retval; |
6052 | 95 Matrix u = args(0).matrix_value (); |
96 | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
97 RowVector N(u.cols ()); |
6052 | 98 if (! args(2).is_cell ()) |
99 { | |
100 | |
101 double p = args(1).idx_type_value (); | |
102 RowVector U = args(2).row_vector_value (true, true); | |
103 assert (U.numel () == p+2); | |
104 | |
105 for (octave_idx_type ii=0; ii<u.numel (); ii++) | |
106 N(ii) = onebasisfun__ (u(ii), p, U); | |
107 | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
108 if (nargout == 2) { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
109 RowVector Nder(u.cols ()); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
110 for (octave_idx_type ii=0; ii<u.numel (); ii++) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
111 Nder(ii) = onebasisfunder__ (u(ii), p, U); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
112 retval(1) = Nder; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
113 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
114 |
6052 | 115 } else { |
116 RowVector p = args(1).row_vector_value (); | |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
117 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
118 if (p.length() == 2) { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
119 Cell C = args(2).cell_value (); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
120 RowVector U = C(0).row_vector_value (true, true); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
121 RowVector V = C(1).row_vector_value (true, true); |
6052 | 122 |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
123 for (octave_idx_type ii=0; ii<u.cols (); ii++) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
124 { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
125 N(ii) = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
126 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
127 //std::cout << "N=" << N(ii) << "\n\n\n"; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
128 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
129 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
130 if (nargout == 2) { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
131 Matrix Nder (2, u.cols()); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
132 for (octave_idx_type ii=0; ii<u.cols (); ii++) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
133 { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
134 Nder(0,ii) = onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
135 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
136 Nder(1,ii) = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
137 onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
138 //std::cout << "N=" << N(ii) << "\n\n\n"; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
139 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
140 retval(1) = Nder; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
141 } |
6052 | 142 |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
143 } else if (p.length() == 3) { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
144 Cell C = args(2).cell_value (); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
145 RowVector U = C(0).row_vector_value (true, true); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
146 RowVector V = C(1).row_vector_value (true, true); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
147 RowVector W = C(2).row_vector_value (true, true); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
148 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
149 for (octave_idx_type ii=0; ii<u.cols (); ii++) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
150 { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
151 N(ii) = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
152 onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
153 onebasisfun__ (u(2, ii), octave_idx_type(p(2)), W); |
6052 | 154 //std::cout << "N=" << N(ii) << "\n\n\n"; |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
155 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
156 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
157 if (nargout == 2) { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
158 Matrix Nder (3, u.cols()); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
159 double Nu, Nv, Nw; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
160 for (octave_idx_type ii=0; ii<u.cols (); ii++) |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
161 { |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
162 Nu = onebasisfun__ (u(0, ii), octave_idx_type(p(0)), U); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
163 Nv = onebasisfun__ (u(1, ii), octave_idx_type(p(1)), V); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
164 Nw = onebasisfun__ (u(2, ii), octave_idx_type(p(2)), W); |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
165 Nder(0,ii) = onebasisfunder__ (u(0, ii), octave_idx_type(p(0)), U) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
166 Nv * Nw; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
167 Nder(1,ii) = onebasisfunder__ (u(1, ii), octave_idx_type(p(1)), V) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
168 Nu * Nw; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
169 Nder(2,ii) = onebasisfunder__ (u(2, ii), octave_idx_type(p(2)), W) * |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
170 Nu * Nv; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
171 //std::cout << "N=" << N(ii) << "\n\n\n"; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
172 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
173 retval(1) = Nder; |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
174 } |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
175 |
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
176 } |
6052 | 177 } |
9403
2108d03d50de
Fixed bug in 1D, added 3D case, added derivatives
rafavzqz
parents:
7204
diff
changeset
|
178 retval(0) = octave_value (N); |
7204 | 179 return retval; |
6052 | 180 } |