10 #include <dune/common/fvector.hh>
11 #include <dune/common/misc.hh>
12 #include <dune/common/typetraits.hh>
31 template <
class F,
int dimD,
unsigned int deriv>
41 typedef Dune::FieldVector<F,size>
Block;
93 template <
class F,
unsigned int deriv>
105 template <
class F,
int dimD>
113 typedef Dune::FieldVector<F,size>
Block;
168 template <
class F,
int dimD,
int dimR,
unsigned int deriv,
173 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
185 static const unsigned int dimDomain = dimD;
186 static const unsigned int dimRange = dimR;
188 enum { size = Base::size+ThisLFETensor::size*dimR };
189 typedef Dune::FieldVector<F,size>
Block;
196 This &operator=(
const Dune::FieldVector<ThisLFETensor,dimR> &t)
201 template <
unsigned int dorder>
204 tensor<dorder>() = t;
219 void axpy(
const F &a,
const This &y)
235 for (
int rr=0;rr<dimR;++rr)
236 tensor_[rr] = y[rr].
template tensor<deriv>()[0];
239 template <
class Fy,
int dimRy>
242 assign<Fy,dimRy>(y.block(),r);
258 return reinterpret_cast<Block&
>(*this);
262 return reinterpret_cast<const Block&
>(*this);
265 template <
unsigned int dorder>
266 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
const
269 const integral_constant<int,dorder> a = {};
272 template <
unsigned int dorder>
273 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
276 return tensor(integral_constant<int,dorder>());
278 template <
unsigned int dorder>
279 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
const
282 const integral_constant<int,dorder> a = {};
283 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
285 template <
unsigned int dorder>
286 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
289 const integral_constant<int,dorder> a = {};
290 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
299 template <
class Fy,
int dimRy>
300 void assign(
const FieldVector<Fy,size*dimRy> &y,
unsigned int r)
302 Base::template assign<Fy,dimRy>(
reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&
>(y),r);
303 tensor_[0] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[Base::size*dimRy+r*ThisLFETensor::size]);
306 void assign(
unsigned int r,
const FieldVector<Fy,size/dimR> &y)
308 Base::assign(r,
reinterpret_cast<const FieldVector<Fy,Base::size/dimR
>&>(y));
309 tensor_[r] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[Base::size/dimR]);
312 template <
class Fy,
unsigned int dy>
316 for (
int rr=0;rr<dimR;++rr)
317 tensor_[rr] = y[rr].
template tensor<deriv>()[0];
320 template <
int dorder>
321 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
322 tensor(
const integral_constant<int,dorder> &dorderVar)
const
324 return Base::tensor(dorderVar);
326 const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
327 tensor(
const integral_constant<int,deriv> &dorderVar)
const
331 template <
int dorder>
332 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
333 tensor(
const integral_constant<int,dorder> &dorderVar)
335 return Base::tensor(dorderVar);
337 Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
338 tensor(
const integral_constant<int,deriv> &dorderVar)
342 Dune::FieldVector<ThisLFETensor,dimR>
tensor_;
345 template <
class F,
int dimD,
int dimR>
355 static const unsigned int dimDomain = dimD;
356 static const unsigned int dimRange = dimR;
358 enum { size = ThisLFETensor::size*dimR };
359 typedef Dune::FieldVector<F,size>
Block;
364 for (
int r=0;r<dimR;++r)
365 tensor_[r] = field_cast<F>(f);
368 This &operator=(
const Dune::FieldVector<ThisLFETensor,dimR> &t)
386 void axpy(
const F &a,
const This &y)
398 for (
int rr=0;rr<dimR;++rr)
399 tensor_[rr] = y[rr].
template tensor<0>()[0];
401 template <
class Fy,
int dimRy>
404 assign<Fy,dimRy>(y.block(),r);
409 tensor_[r].assign(y[0]);
414 tensor_[r].assign(y[0][0]);
419 return reinterpret_cast<Block&
>(*this);
423 return reinterpret_cast<const Block&
>(*this);
432 template <
int dorder>
433 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
const
437 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
441 template <
unsigned int dorder>
442 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
const
445 const integral_constant<int,dorder> a = {};
446 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
448 template <
unsigned int dorder>
449 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
452 const integral_constant<int,dorder> a = {};
453 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
457 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
458 tensor(
const integral_constant<int,0> &dorderVar)
const
462 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
463 tensor(
const integral_constant<int,0> &dorderVar)
467 template <
class Fy,
unsigned int dy>
470 for (
int rr=0;rr<dimR;++rr)
471 tensor_[rr] = y[rr].
template tensor<0>()[0];
473 template <
class Fy,
int dimRy>
474 void assign(
const FieldVector<Fy,size*dimRy> &y,
unsigned int r)
476 tensor_[0] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[r*ThisLFETensor::size]);
479 void assign(
unsigned int r,
const FieldVector<Fy,size/dimR> &y)
483 Dune::FieldVector<ThisLFETensor,dimR>
tensor_;
487 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
497 static const unsigned int dimDomain = dimD;
498 static const unsigned int dimRange = dimR;
500 enum { size = ScalarDeriv::size*dimR };
501 typedef Dune::FieldVector<F,size>
Block;
522 void axpy(
const FF &a,
const This &y)
524 block().
axpy(field_cast<F>(a),y.
block());
536 for (
unsigned int rr=0;rr<dimR;++rr)
537 deriv_[rr].assign(y,rr);
540 template <
class Fy,DerivativeLayout layouty>
543 deriv_[r].assign(r,y);
548 return reinterpret_cast<Block&
>(*this);
552 return reinterpret_cast<const Block&
>(*this);
562 Dune::FieldVector<ScalarDeriv,dimR>
deriv_;
568 template <
class Vec1,
class Vec2,
unsigned int deriv>
571 template <
class Field>
572 static void apply(
unsigned int r,
const Field &a,
573 const Vec1 &x, Vec2 &y)
578 template <
class F1,
int dimD,
int dimR,
585 template <
class Field>
586 static void apply(
unsigned int r,
const Field &a,
587 const Vec1 &x, Vec2 &y)
589 const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
590 for (
int i=0;i<y.size;++i)
594 template <
class F1,
int dimD,
int dimR,
601 template <
class Field>
602 static void apply(
unsigned int r,
const Field &a,
603 const Vec1 &x, Vec2 &y)
605 for (
int rr=0;rr<dimR;++rr)
607 Vec2,deriv>::
apply(rr,a,x[rr],y);
610 template <
class F1,
int dimD,
617 template <
class Field>
618 static void apply(
unsigned int r,
const Field &a,
619 const Vec1 &x, Vec2 &y)
622 Vec2,deriv>
::apply(r,a,x[0],y);
625 template <
class F1,
int dimD,
632 template <
class Field>
633 static void apply(
unsigned int r,
const Field &a,
634 const Vec1 &x, Vec2 &y)
637 const unsigned int rr = r*LFETensorType::size;
638 const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
639 for (
int i=0;i<FieldVector<F1,LFETensorType::size>::dimension;++i)
647 template <
class Vec1,
class Vec2>
650 static void apply(
unsigned int r,
const Vec1 &vec1,Vec2 &vec2)
667 template <
int dimD,
int dimR,
unsigned int deriv,
679 template <
int dimD,
int dimR,
unsigned int deriv,
715 template <
int dimD,
unsigned int deriv,
727 template <
int dimD,
unsigned int deriv,
739 template <
int dimD,
unsigned int deriv,
751 template <
int dimD,
unsigned int deriv,
775 template <
int dimD,
int dimR,
776 class F1,
unsigned int deriv,
781 typedef FieldVector<F2,dimR>
Vec2;
787 template <
int dimD,
int dimR,
788 class F1,
unsigned int deriv,
793 typedef FieldVector<F2,dimR>
Vec2;
796 for (
int rr=0;rr<dimR;++rr)
797 field_cast(vec1[rr].
template tensor<0>()[0].block(),vec2[rr]);
801 class F1,
unsigned int deriv,
806 typedef FieldVector<F2,dimR>
Vec2;
809 field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
813 class F1,
unsigned int deriv,
818 typedef FieldVector<F2,dimR>
Vec2;
821 field_cast(vec1[0].
template tensor<0>()[0].block(),vec2[r]);
825 class F1,
unsigned int deriv,
830 typedef FieldVector<F2,1>
Vec2;
833 field_cast(vec1.template tensor<0>()[0].block(),vec2);
837 class F1,
unsigned int deriv,
842 typedef FieldVector<F2,1>
Vec2;
845 field_cast(vec1[0].
template tensor<0>()[0].block(),vec2);
852 template <
class F,
int dimD,
unsigned int deriv>
853 std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
855 return out << tensor.block();
858 template <
class F,
int dimD,
unsigned int deriv>
859 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
861 out <<
static_cast<const ScalarDerivatives< F,dimD,deriv-1
> &>(d);
862 out <<
" , " << d.tensor() << std::endl;
865 template <
class F,
int dimD>
866 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
868 out << d.tensor() << std::endl;
872 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
873 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,derivative > &d )
877 for (
int r=1;r<dimR;++r)
879 out <<
" , " << d[r];
881 out <<
" ) " << std::endl;
884 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
885 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,value > &d )
890 for (
int r=1;r<dimR;++r)
892 out <<
" , " << d[r];
894 out <<
" ) " << std::endl;
897 template <
class F,
int dimD,
int dimR>
898 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,derivative > &d )
902 for (
int r=1;r<dimR;++r)
904 out <<
" , " << d[r];
906 out <<
" ) " << std::endl;
909 template <
class F,
int dimD,
int dimR>
910 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,value > &d )
914 for (
int r=1;r<dimR;++r)
916 out <<
" , " << d[r];
918 out <<
" ) " << std::endl;
921 template <
class F,
int dimD,
int dimR,
unsigned int deriv,DerivativeLayout layout>
922 std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
924 out <<
"Number of basis functions: " << y.size() << std::endl;
925 for (
unsigned int i=0;i<y.size();++i)
927 out <<
"Base " << i <<
" : " << std::endl;
934 #endif // DUNE_TENSOR_HH