G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsLagrangeBasis.hpp
1
2#pragma once
3
4namespace gismo
5{
6
7template<class T>
9{
10 result.resize(1,m_p+1);
11 gsMatrix<> res(1,1);
12 for(int j = 0; j<=m_p;++j)
13 {
14 anchor_into(j,res);
15 result(0,j)=res(0,0);
16 }
17}
18
19template<class T>
21{
22 result.resize(m_p+1,u.cols());
23 result.colwise() = gsVector<index_t>::LinSpaced(m_p+1, 0, m_p);
24}
25
26template<class T>
28{
29 return new gsLagrangeBasis<T>(this->m_breaks,this->m_start,this->m_end);
30}
31
32template<class T>
34{
35 gsMatrix<T> res(1,2);
36 if(dir==0)
37 {
38 (res)(0,0)=m_start;
39 (res)(0,1)=m_end;
40 }
41 else
42 {
43 (res)(0,1)=m_start;
44 (res)(0,0)=m_end;
45 }
46 return res;
47}
48
49template<class T>
51{
52 result.resize(m_p+1,u.cols());
53 for (int i=0;i<=m_p;++i)
54 {
55 gsMatrix<> res(1,u.cols());
56 evalSingle_into(i,u,res);
57 for(int j=0;j<u.cols();++j)
58 {
59 result(i,j)=res(0,j);
60 }
61 }
62}
63
64template<class T>
66{
67 result.resize(m_p+1,u.cols());
68 for (int i=0;i<=m_p;++i)
69 {
70 gsMatrix<> res(1,u.cols());
71 derivSingle_into(i,u,res);
72 for(int j=0;j<u.cols();++j)
73 {
74 result(i,j)=res(0,j);
75 }
76 }
77}
78
79template<class T>
81{
82 result.resize(m_p+1,u.cols());
83 for(int i=0;i<=m_p;i++)
84 {
85 gsMatrix<> res(1,u.cols());
86 deriv2Single_into(i,u,res);
87 for(int j=0;j<u.cols();j++)
88 {
89 result(i,j)=res(0,j);
90 }
91 }
92}
93
94template<class T>
96{
97 result.resize((n+1)*(m_p + 2),u.cols());
98 for (int i=0;i<=m_p+1;++i)
99 {
100 gsMatrix<> res(n+1,u.cols());
101 evalAllDersSingle_into(i,u,n,res);
102 for(int j=0;j<u.cols();++j)
103 {
104 for(int k=0;k<n+1;k++)
105 {
106 result(i*(n+1)+k,j)=res(1,j);
107 }
108 }
109 }
110}
111
112template<class T>
114{
115 result.resize(n+1, u.cols() );
116 result.setZero();
117 for(int order = 0;order<n+1;order++)
118 {
119 gsMatrix<> res(1,u.cols());
120 evalDerSingle_into(i,u,order,res);
121 for(int j=0;j<u.cols();j++)
122 result(order,j)=res(0,j);
123 }
124}
125
126template<class T>
128{
129 result.resize(1, u.cols() );
130 result.setZero();
131 if(n<this->degree())
132 {
133 T fact = _getFactor(i);
134 std::vector<int> vec(m_p-n);
135 for(int k = 0; k<=m_p-n-1;k++) vec[k]=k;
136
137 do
138 {
139 std::vector<real_t> mults(u.cols());
140 for(int j=0;j<u.cols();j++)
141 {
142 mults[j]=1;
143 for(int k = 0; k<=m_p-n-1;k++)
144 {
145 int index = vec[k];
146 if(index>=static_cast<int>(i))
147 index++;
148 mults[j]*=(u(0,j)-m_breaks[index]);
149 }
150 result(0,j)+=mults[j];
151 }
152 }
153 while(_nextPoint(vec, m_p));
154 for(int j=0;j<u.cols();j++)
155 {
156 if(n>1)
157 result(0,j)/=(fact/std::pow(2,n-1));
158 else
159 result(0,j)/=fact;
160 }
161 }
162 else if(n==this->degree())
163 {
164 T fact = _getFactor(i);
165 for(int j=0;j<u.cols();j++)
166 {
167 if(n>1)
168 result(0,j)=m_p*std::pow(2,n-1)/fact;
169 else
170 result(0,j)=m_p/fact;
171 }
172 }
173 else
174 {
175 for(int j=0;j<u.cols();j++)
176 result(0,j)=0;
177 }
178}
179
180template<class T>
182{
183 T old_length = m_end-m_start;
184 for(unsigned i=0;i<m_breaks.size();i++)
185 {
186 m_breaks[i]=(m_breaks[i]-m_start)/old_length;
187 }
188 m_start = 0;
189 m_end = 1;
190}
191
192template<class T>
194{
195 int basis_size = m_breaks.size();
196 result.resize(basis_size, basis_size);
197 result.setZero();
198 for(int i = 0; i < basis_size; i++)
199 {
200 T fact = _getFactor(i);
201 for(int j = 0; j < basis_size-1; j++)
202 {
203 std::vector<int> vec(m_p-j);
204 for(int k = 0; k<m_p-j;k++)
205 {
206 vec[k]=k;
207 }
208 real_t sum = 0;
209 do
210 {
211 real_t prod=1;
212 for(int k = 0; k<m_p-j;k++)
213 {
214 int index = vec[k];
215 if(index>=i)
216 index++;
217 prod*=m_breaks[index];
218 }
219 sum+=prod;
220 }
221 while(_nextPoint(vec,basis_size-1));
222 result(i,j)=sum;
223 }
224 result(i,basis_size-1)=1;
225 int sign = 1;
226 for(int j = basis_size-1; j >= 0; j--)
227 {
228 result(i,j) *= sign / fact;
229 sign *= -1;
230 }
231 }
232 result.transposeInPlace(); // transpose for control points
233}
234
235template<class T>
237{
238 std::vector<T> breaks = m_breaks;
239 int basis_size = breaks.size();
240 result.resize(basis_size, basis_size);
241 result.setZero();
242 int n=this->degree();
243 for(int j = 0; j < basis_size; j++)
244 {
245 for(int i = 0; i <= j; i++)
246 {
247 result(i,j)=1./binomial(n,j)*binomial(n-i,j-i);
248 }
249 }
250 result.transposeInPlace(); // transpose for control points
251}
252
253template<class T>
254bool gsLagrangeBasis<T>::_nextPoint(std::vector<int> & vec, int end) const
255{
256 int n=vec.size(), changed=-1, i = 0;
257 for(int j = n-1 ; j>=0 ; j--)
258 {
259 if(vec[j]+1<end-i){
260 vec[j]++;
261 changed=j;
262 break;
263 }
264 i++;
265 }
266 if (changed == -1)
267 return false;
268 for(int j=changed;j<n-1;++j)
269 {
270 vec[j+1]=vec[j]+1;
271 }
272 return true;
273}
274
275} // namespace gismo
276
Struct which represents a certain side of a box.
Definition gsBoundary.h:85
A basis represents a family of scalar basis functions defined over a common parameter domain.
Definition gsBasis.h:79
A univariate Lagrange basis.
Definition gsLagrangeBasis.h:46
void deriv2_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates the (partial) derivatives of the nonzero basis functions at points u into result.
Definition gsLagrangeBasis.hpp:80
void anchors_into(gsMatrix< T > &result) const
Returns the anchors (greville points) of the basis.
Definition gsLagrangeBasis.hpp:8
void evalAllDersSingle_into(index_t i, const gsMatrix< T > &u, int n, gsMatrix< T > &result) const
Evaluate the basis function i and its derivatives up to order n at points u into result.
Definition gsLagrangeBasis.hpp:113
void reparameterizeToZeroOne()
Changes the basis so the curve is defined in the interval [0,1].
Definition gsLagrangeBasis.hpp:181
void _getTransformationLagrangeMonomial(gsMatrix< T > &result) const
Definition gsLagrangeBasis.hpp:193
void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates the non-zero basis functions at value u.
Definition gsLagrangeBasis.hpp:50
void active_into(const gsMatrix< T > &u, gsMatrix< index_t > &result) const
Returns the indices of active (non zero) basis functions at points (columns of) u,...
Definition gsLagrangeBasis.hpp:20
bool _nextPoint(std::vector< int > &vec, int end) const
Definition gsLagrangeBasis.hpp:254
void evalAllDers_into(const gsMatrix< T > &u, int n, gsMatrix< T > &result) const
Evaluate the nonzero basis functions and their derivatives up to order n at points u into result.
Definition gsLagrangeBasis.hpp:95
void deriv_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates the (partial) derivatives of non-zero basis functions at (the columns of) u.
Definition gsLagrangeBasis.hpp:65
gsMatrix< T > supportInterval(index_t dir) const
Returns an interval that contains the parameter values in direction \ dir.
Definition gsLagrangeBasis.hpp:33
void _getTransformationMonomialBezier(gsMatrix< T > &result) const
Definition gsLagrangeBasis.hpp:236
void evalDerSingle_into(index_t i, const gsMatrix< T > &u, int n, gsMatrix< T > &result) const
Evaluate the (partial) derivative(s) of order n the i-th basis function at points u into result.
Definition gsLagrangeBasis.hpp:127
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
unsigned binomial()
Returns binomial(n,r) as a compile time constant.
Definition gsCombinatorics.h:125
#define index_t
Definition gsConfig.h:32
The G+Smo namespace, containing all definitions for the library.