G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMvMonomialBasis.h
Go to the documentation of this file.
1
15#pragma once
16
17#include <gsCore/gsLinearAlgebra.h>
18#include <gsCore/gsBasisFun.h>
19
20#include <gsCore/gsDebug.h>
21
22#include <gsCore/gsBoundary.h>
23
24#include <gsCore/gsBasis.h>
25
26#include <gsUtils/gsCombinatorics.h>
27
28namespace gismo
29{
30
31template<short_t d,class T>
32class gsMvMonomialBasis : public gsBasis<T>
33{
34public:
35#define Eigen gsEigen
36 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
37#undef Eigen
38
39public:
41 typedef memory::shared_ptr< gsMvMonomialBasis > Ptr;
42
44 typedef memory::unique_ptr< gsMvMonomialBasis > uPtr;
45
46 typedef T Scalar_t;
47
48 static const bool IsRational = false;
49
50 typedef memory::unique_ptr< gsDomainIterator<T> > domainIter;
51
53 static const short_t Dim = d;
54
55private:
56
57 typedef gsEigen::internal::variable_if_dynamic<unsigned,d> dimType;
58 dimType m_d;
59 short_t m_degree;
60 std::vector<gsVector<index_t> > m_compositions;
61
62public:
63 gsMvMonomialBasis() : m_d(-1), m_degree(-1) { }
64
65 gsMvMonomialBasis(unsigned _d, unsigned p) : m_d(_d), m_degree(p)
66 {
67 getCompositions(m_compositions);
68 }
69
70 // enable_if<d"=-1>
71 gsMvMonomialBasis(unsigned p) : m_d(d), m_degree(p)
72 {
73 getCompositions(m_compositions);
74 }
75
77 ~gsMvMonomialBasis() { };
78
79public:
80
81 // Look at gsBasis class for a description
82 short_t domainDim() const;
83
84 // Look at gsBasis class for a description
85 void active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const
86 {
87 const int sz = size();
88 result.resize( sz, u.cols() );
89 for ( int i = 0; i< sz; ++i )
90 result.row(i).setConstant(i);// globally active
91 }
92
93 // Look at gsBasis class for a description
94 gsMatrix<T> support() const;
95
96 // Look at gsBasis class for a description
97 gsMatrix<T> support(const index_t & i) const;
98
99 // Look at gsBasis class for a description
100 void eval_into(const gsMatrix<T> & u, gsMatrix<T>& result) const;
101
102 // Look at gsBasis class for a description
103 void evalSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result) const;
104
105 // Look at gsBasis class for a description
106 void deriv_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const;
107
108 // Look at gsBasis class for a description
109 void derivSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
110
111 // Look at gsBasis class for a description
112 void deriv2_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const;
113
114 // Look at gsBasis class for a description
115 void deriv2Single_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const;
116
117 GISMO_CLONE_FUNCTION(gsMvMonomialBasis)
118
119 // Look at gsBasis class for a description
120 memory::unique_ptr<gsGeometry<T> > makeGeometry(gsMatrix<T> coefs ) const
121 { GISMO_NO_IMPLEMENTATION }
122
123 // Look at gsBasis class for a description
124 domainIter makeDomainIterator() const
125 { GISMO_NO_IMPLEMENTATION }
126
127 // Look at gsBasis class for a description
128 std::ostream &print(std::ostream &os) const
129 {
130 os <<"Multivariate monomial basis basis of degree "
131 << m_degree <<" and "<<m_d.value()<<" variables\n";
132 return os;
133 }
134
136 std::string detail() const
137 {
138 // By default just uses print(..)
139 std::ostringstream os;
140 print(os);
141 return os.str();
142 }
143
144 /*
145 Member functions that may be implemented or not in the derived class
146 */
147
149 index_t size() const;
150
151 short_t totalDegree() const
152 { return m_degree; }
153
154 short_t degree(short_t i = 0) const
155 { return m_degree; }
156
157 int degreeOf(index_t k) const
158 { return m_compositions[k].sum(); }
159
160private:
161 //create compositions for basis functions for m_degree
162 void getCompositions( std::vector<gsVector<index_t> > & compos)const;
163
164 // create compositions for basis functions for arbitrary degree
165 void getCompositions( std::vector<gsVector<index_t> > & compos, index_t degree)const;
166
167 //rt derivative
168 //void rThDerivSingle(index_t r,index_t i, const gsMatrix<T> & dir, const gsMatrix<T> & u, gsMatrix<T>& result)const;
169
171 int findIndex(gsVector<index_t> const p, std::vector<gsVector<index_t> >const compos )const;
172
173}; // class gsMvMonomialBasis
174
175
176template<short_t d,class T>
177short_t gsMvMonomialBasis<d,T>::domainDim() const { return m_d.value(); }
178
179template<short_t d,class T>
180index_t gsMvMonomialBasis<d,T>::size() const{
181 return m_compositions.size();
182}
183
184template<short_t d,class T>
185gsMatrix<T> gsMvMonomialBasis<d,T>::support(const index_t & i) const
186{return gsMatrix<T>();}
187
188template<short_t d,class T>
189gsMatrix<T> gsMvMonomialBasis<d,T>:: support() const
190{return gsMatrix<T>();}
191
192template<short_t d,class T>
193void gsMvMonomialBasis<d,T>::evalSingle_into(index_t i,
194 const gsMatrix<T> & u,
195 gsMatrix<T>& result) const
196{
197 gsVector<T> point;
198 result.setZero(1,u.cols() );
199 for(int j = 0; j < u.cols(); j++)
200 {
201 point = u.col(j);
202 T val = 1;
203 for(int k = 0; k < m_compositions[i].size(); k++)
204 {
205 val = val * math::pow(point[k], static_cast<int>(m_compositions[i][k]));
206 }
207 result.at(j) = val;
208 }
209}
210
211template<short_t d,class T>
212void gsMvMonomialBasis<d,T>::eval_into(const gsMatrix<T> & u, gsMatrix<T>& result) const
213{
214 result.setZero(m_compositions.size(), u.cols());
215 gsMatrix<T> single_result;
216 for(unsigned i = 0; i < m_compositions.size(); i++)
217 {
218 evalSingle_into(i, u, single_result);
219 result.row(i) = single_result.row(0);
220 }
221}
222
223template<short_t d,class T>
224void gsMvMonomialBasis<d,T>::derivSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const
225{
226GISMO_NO_IMPLEMENTATION
227}
228
229template<short_t d,class T>
230void gsMvMonomialBasis<d,T>::deriv_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const
231{
232 result.setZero(m_compositions.size()*this->dim(), u.cols());
233 gsMatrix<T> single_result;
234 for(unsigned i = 0; i < m_compositions.size(); i++)
235 {
236 derivSingle_into(i, u, single_result);
237 for(int j = 0; j < single_result.rows();j++)
238 {
239 result.row(m_d.value()*i+j) = single_result.row(j); // Map
240 }
241 }
242}
243
244template<short_t d,class T>
245void gsMvMonomialBasis<d,T>::deriv2Single_into(index_t i, const gsMatrix<T> & u,
246 gsMatrix<T>& result ) const
247{
248 result.setZero(m_d.value(), u.cols());
249 gsVector<T> point;
250
251 for(index_t j = 0; j < u.cols(); ++j)
252 {
253 point = u.col(j);
254 for(index_t l = 0; l < m_compositions[i].size(); ++l)
255 {
256 T val = 1;
257 for(index_t k = 0; k < m_compositions[i].size(); ++k)
258 {
259 const unsigned ex = m_compositions[i][k];
260 if (ex<1 ) { val=0; break; }
261 if (ex==1) { continue; }
262 val = val * ex * math::pow(point[k], static_cast<int>(ex-1));
263 }
264 result(l,j) = val;
265 }
266 }
267}
268
269template<short_t d,class T>
270void gsMvMonomialBasis<d,T>::deriv2_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const
271{
272 result.setZero(m_compositions.size()*((d * (d + 1)) / 2), u.cols());
273 gsMatrix<T> single_result;
274 for(unsigned i = 0; i < m_compositions.size(); i++)
275 {
276 deriv2Single_into(i, u, single_result);
277 for(int j = 0; j < single_result.rows();j++){
278 result.row(single_result.rows()*i+j) = single_result.row(j);
279 }
280 }
281}
282
283template<short_t d,class T>
284void gsMvMonomialBasis<d,T>::getCompositions(std::vector<gsVector<index_t> > & compos)const{
285 this->getCompositions(compos,m_degree);
286}
287
288
289namespace
290{
291bool comp_deg(const gsVector<index_t> & a, const gsVector<index_t>& b)
292{ return a.sum()<b.sum(); }
293}
294
295template<short_t d,class T>
296void gsMvMonomialBasis<d,T>::getCompositions(std::vector<gsVector<index_t> > & compos, index_t degree)const
297{
298 compos.reserve(numCompositions(degree,m_d.value()+1));
299 gsVector<index_t> RR;
300 firstComposition(degree,m_d.value()+1,RR);
301 do
302 {
303 compos.push_back(RR.tail(m_d.value()));
304 } while ( nextComposition(RR) );
305
306 std::sort(compos.begin(), compos.end(), comp_deg );
307}
308
309template<short_t d,class T>
310int gsMvMonomialBasis<d,T>::findIndex(gsVector<index_t> const p,
311 std::vector<gsVector<index_t> >const compos )const
312{
313 //std::find
314 for(index_t i = 0; i < compos.size();i++)
315 {
316 if(p==compos[i])return i;
317 }
318 return -1;
319}
320
321}; // namespace gismo
322
323
324// #ifndef GISMO_BUILD_LIB
325// #include GISMO_HPP_HEADER(gsMvMonomialBasis.hpp)
326// #endif
gismo::gsBasis::makeDomainIterator
virtual domainIter makeDomainIterator(const boxSide &s) const
Create a boundary domain iterator for the computational mesh this basis, that points to the first ele...
Definition gsBasis.hpp:547
gismo::nextComposition
bool nextComposition(Vec &v)
Returns (inplace) the next composition in lexicographic order.
Definition gsCombinatorics.h:667
gismo::firstComposition
void firstComposition(typename Vec::Scalar sum, index_t dim, Vec &res)
Construct first composition of sum into dim integers.
Definition gsCombinatorics.h:657
gismo::numCompositions
unsigned numCompositions(int sum, short_t dim)
Number of compositions of sum into dim integers.
Definition gsCombinatorics.h:691
gsBasisFun.h
Provides definition of the BasisFun class.
gsBasis.h
Provides declaration of Basis abstract interface.
gsBoundary.h
Provides structs and classes related to interfaces and boundaries.
gsCombinatorics.h
Provides combinatorial unitilies.
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
gsDebug.h
This file contains the debugging and messaging system of G+Smo.
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo
The G+Smo namespace, containing all definitions for the library.