10 #include <dune/common/fvector.hh>
29 template <
class F,
int dimD,
unsigned int deriv>
39 typedef Dune::FieldVector<F,size>
Block;
91 template <
class F,
unsigned int deriv>
103 template <
class F,
int dimD>
111 typedef Dune::FieldVector<F,size>
Block;
166 template <
class F,
int dimD,
int dimR,
unsigned int deriv,
171 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
172 struct Derivatives<F,dimD,dimR,deriv,
value>
173 :
public Derivatives<F,dimD,dimR,deriv-1,value>
183 static const unsigned int dimDomain = dimD;
184 static const unsigned int dimRange = dimR;
186 enum { size = Base::size+ThisLFETensor::size*dimR };
187 typedef Dune::FieldVector<F,size>
Block;
194 This &operator=(
const Dune::FieldVector<ThisLFETensor,dimR> &t)
199 template <
unsigned int dorder>
202 tensor<dorder>() = t;
217 void axpy(
const F &a,
const This &y)
224 void assign(
const Derivatives<Fy,dimD,dimR,deriv,value> &y)
230 void assign(
const Derivatives<Fy,dimD,dimR,deriv,derivative> &y)
233 for (
int rr=0; rr<dimR; ++rr)
234 tensor_[rr] = y[rr].
template tensor<deriv>()[0];
237 template <
class Fy,
int dimRy>
238 void assign(
const Derivatives<Fy,dimD,dimRy,deriv,value> &y,
unsigned int r)
240 assign<Fy,dimRy>(y.block(),r);
244 void assign(
unsigned int r,
const Derivatives<Fy,dimD,1,deriv,value> &y)
249 void assign(
unsigned int r,
const Derivatives<Fy,dimD,1,deriv,derivative> &y)
256 return reinterpret_cast<Block&
>(*this);
260 return reinterpret_cast<const Block&
>(*this);
263 template <
unsigned int dorder>
264 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
const
267 const integral_constant<int,dorder> a = {};
270 template <
unsigned int dorder>
271 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
274 return tensor(integral_constant<int,dorder>());
276 template <
unsigned int dorder>
277 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
const
280 const integral_constant<int,dorder> a = {};
281 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
283 template <
unsigned int dorder>
284 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
287 const integral_constant<int,dorder> a = {};
288 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
297 template <
class Fy,
int dimRy>
298 void assign(
const FieldVector<Fy,size*dimRy> &y,
unsigned int r)
300 Base::template assign<Fy,dimRy>(
reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&
>(y),r);
301 tensor_[0] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[Base::size*dimRy+r*ThisLFETensor::size]);
304 void assign(
unsigned int r,
const FieldVector<Fy,size/dimR> &y)
306 Base::assign(r,
reinterpret_cast<const FieldVector<Fy,Base::size/dimR
>&>(y));
307 tensor_[r] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[Base::size/dimR]);
310 template <
class Fy,
unsigned int dy>
311 void assign(
const Derivatives<Fy,dimD,dimR,dy,derivative> &y)
314 for (
int rr=0; rr<dimR; ++rr)
315 tensor_[rr] = y[rr].
template tensor<deriv>()[0];
318 template <
int dorder>
319 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
320 tensor(
const integral_constant<int,dorder> &dorderVar)
const
322 return Base::tensor(dorderVar);
324 const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
325 tensor(
const integral_constant<int,deriv> &dorderVar)
const
329 template <
int dorder>
330 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
331 tensor(
const integral_constant<int,dorder> &dorderVar)
333 return Base::tensor(dorderVar);
335 Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
336 tensor(
const integral_constant<int,deriv> &dorderVar)
340 Dune::FieldVector<ThisLFETensor,dimR>
tensor_;
343 template <
class F,
int dimD,
int dimR>
344 struct Derivatives<F,dimD,dimR,0,
value>
353 static const unsigned int dimDomain = dimD;
354 static const unsigned int dimRange = dimR;
356 enum { size = ThisLFETensor::size*dimR };
357 typedef Dune::FieldVector<F,size>
Block;
362 for (
int r=0; r<dimR; ++r)
363 tensor_[r] = field_cast<F>(f);
366 This &operator=(
const Dune::FieldVector<ThisLFETensor,dimR> &t)
384 void axpy(
const F &a,
const This &y)
389 void assign(
const Derivatives<Fy,dimD,dimR,0,value> &y)
394 void assign(
const Derivatives<Fy,dimD,dimR,0,derivative> &y)
396 for (
int rr=0; rr<dimR; ++rr)
397 tensor_[rr] = y[rr].
template tensor<0>()[0];
399 template <
class Fy,
int dimRy>
400 void assign(
const Derivatives<Fy,dimD,dimRy,0,value> &y,
unsigned int r)
402 assign<Fy,dimRy>(y.block(),r);
405 void assign(
unsigned int r,
const Derivatives<Fy,dimD,1,0,value> &y)
407 tensor_[r].assign(y[0]);
410 void assign(
unsigned int r,
const Derivatives<Fy,dimD,1,0,derivative> &y)
412 tensor_[r].assign(y[0][0]);
417 return reinterpret_cast<Block&
>(*this);
421 return reinterpret_cast<const Block&
>(*this);
430 template <
int dorder>
431 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
const
435 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
439 template <
unsigned int dorder>
440 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
const
443 const integral_constant<int,dorder> a = {};
444 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
446 template <
unsigned int dorder>
447 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
450 const integral_constant<int,dorder> a = {};
451 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR
>&>(tensor(a));
455 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
456 tensor(
const integral_constant<int,0> &dorderVar)
const
460 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
461 tensor(
const integral_constant<int,0> &dorderVar)
465 template <
class Fy,
unsigned int dy>
466 void assign(
const Derivatives<Fy,dimD,dimR,dy,derivative> &y)
468 for (
int rr=0; rr<dimR; ++rr)
469 tensor_[rr] = y[rr].
template tensor<0>()[0];
471 template <
class Fy,
int dimRy>
472 void assign(
const FieldVector<Fy,size*dimRy> &y,
unsigned int r)
474 tensor_[0] =
reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&
>(y[r*ThisLFETensor::size]);
477 void assign(
unsigned int r,
const FieldVector<Fy,size/dimR> &y)
481 Dune::FieldVector<ThisLFETensor,dimR>
tensor_;
485 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
495 static const unsigned int dimDomain = dimD;
496 static const unsigned int dimRange = dimR;
498 enum { size = ScalarDeriv::size*dimR };
499 typedef Dune::FieldVector<F,size>
Block;
520 void axpy(
const FF &a,
const This &y)
522 block().
axpy(field_cast<F>(a),y.
block());
526 void assign(
const Derivatives<Fy,dimD,dimR,deriv,derivative> &y)
532 void assign(
const Derivatives<Fy,dimD,dimR,deriv,value> &y)
534 for (
unsigned int rr=0; rr<dimR; ++rr)
535 deriv_[rr].assign(y,rr);
538 template <
class Fy,DerivativeLayout layouty>
539 void assign(
unsigned int r,
const Derivatives<Fy,dimD,1,deriv,layouty> &y)
541 deriv_[r].assign(r,y);
546 return reinterpret_cast<Block&
>(*this);
550 return reinterpret_cast<const Block&
>(*this);
560 Dune::FieldVector<ScalarDeriv,dimR>
deriv_;
566 template <
class Vec1,
class Vec2,
unsigned int deriv>
569 template <
class Field>
570 static void apply(
unsigned int r,
const Field &a,
571 const Vec1 &x, Vec2 &y)
576 template <
class F1,
int dimD,
int dimR,
582 typedef Derivatives<F1,dimD,dimR,d,value>
Vec1;
583 template <
class Field>
584 static void apply(
unsigned int r,
const Field &a,
585 const Vec1 &x, Vec2 &y)
587 const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
588 for (
int i=0; i<y.size; ++i)
592 template <
class F1,
int dimD,
int dimR,
598 typedef Derivatives<F1,dimD,dimR,d,derivative>
Vec1;
599 template <
class Field>
600 static void apply(
unsigned int r,
const Field &a,
601 const Vec1 &x, Vec2 &y)
603 for (
int rr=0; rr<dimR; ++rr)
605 Vec2,deriv>::
apply(rr,a,x[rr],y);
608 template <
class F1,
int dimD,
614 typedef Derivatives<F1,dimD,1,d,derivative>
Vec1;
615 template <
class Field>
616 static void apply(
unsigned int r,
const Field &a,
617 const Vec1 &x, Vec2 &y)
620 Vec2,deriv>
::apply(r,a,x[0],y);
623 template <
class F1,
int dimD,
629 typedef Derivatives<F1,dimD,1,d,value>
Vec1;
630 template <
class Field>
631 static void apply(
unsigned int r,
const Field &a,
632 const Vec1 &x, Vec2 &y)
635 const unsigned int rr = r*LFETensorType::size;
636 const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
637 for (
int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
645 template <
class Vec1,
class Vec2>
648 static void apply(
unsigned int r,
const Vec1 &vec1,Vec2 &vec2)
656 Derivatives<F2,dimD,dimR,deriv,layout> >
658 typedef Derivatives<F1,dimD,dimR,deriv,layout>
Vec1;
659 typedef Derivatives<F2,dimD,dimR,deriv,layout>
Vec2;
665 template <
int dimD,
int dimR,
unsigned int deriv,
670 typedef Derivatives<F1,dimD,dimR,deriv,value>
Vec1;
671 typedef Derivatives<F2,dimD,dimR,deriv,derivative>
Vec2;
677 template <
int dimD,
int dimR,
unsigned int deriv,
680 Derivatives<F2,dimD,dimR,deriv,
value> >
682 typedef Derivatives<F1,dimD,dimR,deriv,derivative>
Vec1;
683 typedef Derivatives<F2,dimD,dimR,deriv,value>
Vec2;
692 Derivatives<F2,dimD,dimR,deriv,
value> >
694 typedef Derivatives<F1,dimD,1,deriv,layout>
Vec1;
695 typedef Derivatives<F2,dimD,dimR,deriv,value>
Vec2;
706 typedef Derivatives<F1,dimD,1,deriv,layout>
Vec1;
707 typedef Derivatives<F2,dimD,dimR,deriv,derivative>
Vec2;
713 template <
int dimD,
unsigned int deriv,
716 Derivatives<F2,dimD,1,deriv,
value> >
718 typedef Derivatives<F1,dimD,1,deriv,value>
Vec1;
719 typedef Derivatives<F2,dimD,1,deriv,value>
Vec2;
725 template <
int dimD,
unsigned int deriv,
730 typedef Derivatives<F1,dimD,1,deriv,derivative>
Vec1;
731 typedef Derivatives<F2,dimD,1,deriv,derivative>
Vec2;
737 template <
int dimD,
unsigned int deriv,
740 Derivatives<F2,dimD,1,deriv,
value> >
742 typedef Derivatives<F1,dimD,1,deriv,derivative>
Vec1;
743 typedef Derivatives<F2,dimD,1,deriv,value>
Vec2;
749 template <
int dimD,
unsigned int deriv,
754 typedef Derivatives<F1,dimD,1,deriv,value>
Vec1;
755 typedef Derivatives<F2,dimD,1,deriv,derivative>
Vec2;
766 typedef Derivatives<F1,dimD,1,deriv,layout>
Vec1;
773 template <
int dimD,
int dimR,
774 class F1,
unsigned int deriv,
778 typedef Derivatives<F1,dimD,dimR,deriv,value>
Vec1;
779 typedef FieldVector<F2,dimR>
Vec2;
785 template <
int dimD,
int dimR,
786 class F1,
unsigned int deriv,
790 typedef Derivatives<F1,dimD,dimR,deriv,derivative>
Vec1;
791 typedef FieldVector<F2,dimR>
Vec2;
794 for (
int rr=0; rr<dimR; ++rr)
795 field_cast(vec1[rr].
template tensor<0>()[0].block(),vec2[rr]);
799 class F1,
unsigned int deriv,
803 typedef Derivatives<F1,dimD,1,deriv,value>
Vec1;
804 typedef FieldVector<F2,dimR>
Vec2;
807 field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
811 class F1,
unsigned int deriv,
815 typedef Derivatives<F1,dimD,1,deriv,derivative>
Vec1;
816 typedef FieldVector<F2,dimR>
Vec2;
819 field_cast(vec1[0].
template tensor<0>()[0].block(),vec2[r]);
823 class F1,
unsigned int deriv,
827 typedef Derivatives<F1,dimD,1,deriv,value>
Vec1;
828 typedef FieldVector<F2,1>
Vec2;
831 field_cast(vec1.template tensor<0>()[0].block(),vec2);
835 class F1,
unsigned int deriv,
839 typedef Derivatives<F1,dimD,1,deriv,derivative>
Vec1;
840 typedef FieldVector<F2,1>
Vec2;
843 field_cast(vec1[0].
template tensor<0>()[0].block(),vec2);
850 template <
class F,
int dimD,
unsigned int deriv>
851 std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
853 return out << tensor.block();
856 template <
class F,
int dimD,
unsigned int deriv>
857 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
859 out <<
static_cast<const ScalarDerivatives< F,dimD,deriv-1
> &>(d);
860 out <<
" , " << d.tensor() << std::endl;
863 template <
class F,
int dimD>
864 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
866 out << d.tensor() << std::endl;
870 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
871 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,derivative > &d )
875 for (
int r=1; r<dimR; ++r)
877 out <<
" , " << d[r];
879 out <<
" ) " << std::endl;
882 template <
class F,
int dimD,
int dimR,
unsigned int deriv>
883 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,value > &d )
885 out <<
static_cast<const Derivatives< F,dimD,dimR,deriv-1,
value > &>(d);
888 for (
int r=1; r<dimR; ++r)
890 out <<
" , " << d[r];
892 out <<
" ) " << std::endl;
895 template <
class F,
int dimD,
int dimR>
896 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,derivative > &d )
900 for (
int r=1; r<dimR; ++r)
902 out <<
" , " << d[r];
904 out <<
" ) " << std::endl;
907 template <
class F,
int dimD,
int dimR>
908 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,value > &d )
912 for (
int r=1; r<dimR; ++r)
914 out <<
" , " << d[r];
916 out <<
" ) " << std::endl;
919 template <
class F,
int dimD,
int dimR,
unsigned int deriv,DerivativeLayout layout>
920 std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
922 out <<
"Number of basis functions: " << y.size() << std::endl;
923 for (
unsigned int i=0; i<y.size(); ++i)
925 out <<
"Base " << i <<
" : " << std::endl;
932 #endif // DUNE_TENSOR_HH