29 #ifndef __OrthonormalPolynomialSet_H
30 #define __OrthonormalPolynomialSet_H 1
44 template<uint16_type Dim,
47 template<u
int16_type>
class PolySetType = Scalar,
49 uint16_type TheTAG = 0,
50 template<u
int16_type,u
int16_type,u
int16_type>
class Convex = Simplex>
51 class OrthonormalPolynomialSet
54 template<uint16_type Dim,
57 template<u
int16_type>
class PolySetType,
60 class OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Simplex>
62 public PolynomialSet<Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas>, PolySetType >
64 typedef PolynomialSet<Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas>, PolySetType > super;
67 static const uint16_type TAG = TheTAG;
68 static const uint16_type nDim = Dim;
69 static const uint16_type nOrder = Order;
70 static const uint16_type nRealDim = RealDim;
71 static const bool isTransformationEquivalent =
true;
72 typedef OrthonormalPolynomialSet<Dim, Order,RealDim, PolySetType, T, TheTAG, Simplex> self_type;
73 typedef self_type component_basis_type;
75 typedef typename super::polyset_type polyset_type;
76 static const bool is_tensor2 = polyset_type::is_tensor2;
77 static const bool is_vectorial = polyset_type::is_vectorial;
78 static const bool is_scalar = polyset_type::is_scalar;
79 static const bool is_continuous =
false;
80 static const bool is_modal =
true;
81 static const uint16_type nComponents = polyset_type::nComponents;
82 static const bool is_product =
true;
83 static const bool isContinuous =
false;
84 typedef Discontinuous continuity_type;
86 typedef typename super::component_type component_type;
89 typedef Dubiner<Dim, RealDim, Order, Normalized<true>, T, StorageUBlas> basis_type;
90 typedef Simplex<Dim, Order, Dim> convex_type;
94 typedef Simplex<Dim, O, Dim> type;
96 typedef Reference<convex_type, nDim, nOrder, nDim, value_type> reference_convex_type;
98 typedef typename super::polynomial_type polynomial_type;
101 static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
103 static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
105 static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
108 static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
110 static const uint16_type nLocalDof = convex_type::numPoints;
112 static const uint16_type nDof = nLocalDof;
113 static const uint16_type nNodes = nDof;
114 static const uint16_type nDofGrad = super::nDim*nDof;
115 static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
116 typedef typename matrix_node<value_type>::type points_type;
118 OrthonormalPolynomialSet()
120 super( basis_type() )
123 ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
126 std::cout <<
"[orthonormalpolynomialset] m = " << m <<
"\n";
128 if ( !( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
130 std::cout <<
"m1=" << m <<
"\n"
131 <<
"m2=" << polyset_type::toMatrix( polyset_type::toType( m ) ) <<
"\n"
132 << ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) - m ) <<
"\n";
134 FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
135 m ) < 1e-10 )( m ).warn (
"invalid transformation" );
136 this->setCoefficient( polyset_type::toType( m ),
true );
139 OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Simplex > toScalar()
const
141 return OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Simplex >();
147 std::string familyName()
const
153 points_type
points()
const
155 return points_type();
157 points_type
points(
int f )
const
159 return points_type();
163 template<uint16_type Dim,
166 template<u
int16_type>
class PolySetType,
169 const uint16_type OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType,T, TheTAG, Simplex>::nLocalDof;
172 template<uint16_type Dim,
175 template<u
int16_type>
class PolySetType,
178 class OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Hypercube>
180 public PolynomialSet<Legendre<Dim, RealDim, Order, Normalized<true>, T>, PolySetType >
182 typedef PolynomialSet<Legendre<Dim, RealDim, Order, Normalized<true>, T>, PolySetType > super;
185 static const uint16_type nDim = Dim;
186 static const uint16_type nOrder = Order;
187 static const uint16_type nRealDim = RealDim;
188 static const bool isTransformationEquivalent =
true;
190 typedef OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType, T, TheTAG, Hypercube> self_type;
191 typedef self_type component_basis_type;
193 typedef typename super::polyset_type polyset_type;
194 static const bool is_tensor2 = polyset_type::is_tensor2;
195 static const bool is_vectorial = polyset_type::is_vectorial;
196 static const bool is_scalar = polyset_type::is_scalar;
197 static const bool is_continuous =
false;
198 static const bool is_modal =
true;
199 static const uint16_type nComponents = polyset_type::nComponents;
200 static const bool is_product =
true;
201 static const bool isContinuous =
false;
202 typedef Discontinuous continuity_type;
204 typedef typename super::component_type component_type;
205 typedef T value_type;
206 typedef Legendre<Dim, RealDim, Order, Normalized<true>, T> basis_type;
207 typedef Hypercube<Dim, Order, Dim> convex_type;
208 typedef typename matrix_node<value_type>::type points_type;
213 typedef Hypercube<Dim, O, nDim> type;
215 typedef Reference<convex_type, nDim, nOrder, nDim, value_type> reference_convex_type;
217 typedef typename super::polynomial_type polynomial_type;
220 static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
222 static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
224 static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
227 static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
229 static const uint16_type nLocalDof = convex_type::numPoints;
231 static const uint16_type nDof = nLocalDof;
232 static const uint16_type nNodes = nDof;
233 static const uint16_type nDofGrad = super::nDim*nDof;
234 static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
236 OrthonormalPolynomialSet()
238 super( basis_type() )
241 ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
244 std::cout <<
"[orthonormalpolynomialset] m = " << m <<
"\n";
246 FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
247 m ) < 1e-10 )( m ).warn (
"invalid transformation" );
248 this->setCoefficient( polyset_type::toType( m ),
true );
251 OrthonormalPolynomialSet<Dim, Order, RealDim,Scalar,T, TheTAG, Hypercube > toScalar()
const
253 return OrthonormalPolynomialSet<Dim, Order, RealDim, Scalar,T, TheTAG, Hypercube >();
255 std::string familyName()
const
259 points_type
points()
const
261 return points_type();
263 points_type
points(
int f )
const
265 return points_type();
269 template<uint16_type Dim,
272 template<u
int16_type>
class PolySetType,
275 const uint16_type OrthonormalPolynomialSet<Dim, Order, RealDim, PolySetType,T, TheTAG, Hypercube>::nLocalDof;
279 template<uint16_type Order,
280 template<u
int16_type Dim>
class PolySetType = Scalar,
281 uint16_type TheTAG=0 >
282 class OrthonormalPolynomialSet
285 template<uint16_type N,
288 typename Convex = Simplex<N> >
291 typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
292 mpl::identity<detail::OrthonormalPolynomialSet<N,Order,RealDim,PolySetType,T,TheTAG,Simplex> >,
293 mpl::identity<detail::OrthonormalPolynomialSet<N,Order,RealDim,PolySetType,T,TheTAG,Hypercube> > >::type::type result_type;
294 typedef result_type type;
297 template<u
int16_type TheNewTAG>
300 typedef OrthonormalPolynomialSet<Order,PolySetType,TheNewTAG> type;
303 typedef OrthonormalPolynomialSet<Order,Scalar,TheTAG> component_basis_type;
305 static const uint16_type nOrder = Order;
306 static const uint16_type TAG = TheTAG;