1#ifndef DUNE_FEM_SHAPEFUNCTIONSET_VECTORIAL_HH
2#define DUNE_FEM_SHAPEFUNCTIONSET_VECTORIAL_HH
9#include <dune/fem/common/fmatrixcol.hh>
10#include <dune/fem/common/utility.hh>
11#include <dune/fem/space/basisfunctionset/functor.hh>
12#include <dune/fem/space/common/functionspace.hh>
24 template<
class Scalar,
class Vectorial >
25 struct MakeVectorialTraits;
27 template<
class K,
int dimR >
28 struct MakeVectorialTraits< FieldVector< K, 1 >, FieldVector< K, dimR > >
30 typedef FieldVector< K, 1 > ScalarType;
31 typedef FieldVector< K, dimR > VectorialType;
33 typedef typename FieldTraits< VectorialType >::field_type field_type;
34 typedef typename VectorialType::size_type ComponentType;
35 typedef typename VectorialType::size_type size_type;
37 static const size_type factor = dimR;
39 static ComponentType begin () {
return ComponentType( 0 ); }
40 static ComponentType end () {
return ComponentType( factor ); }
42 static VectorialType zeroVectorial () {
return VectorialType( K( field_type( 0 ) ) ); }
44 static const K &access (
const ScalarType &x ) {
return x[ 0 ]; }
45 static K &access ( ScalarType &x ) {
return x[ 0 ]; }
47 static const K &access (
const VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
48 static K &access ( VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
50 static size_type index (
const ComponentType &i ) {
return i; }
53 template<
class K,
int dimR >
54 struct MakeVectorialTraits< ExplicitFieldVector< K, 1 >, ExplicitFieldVector< K, dimR > >
56 typedef FieldVector< K, 1 > ScalarType;
57 typedef ExplicitFieldVector< K, dimR > VectorialType;
59 typedef typename FieldTraits< VectorialType >::field_type field_type;
60 typedef typename VectorialType::size_type ComponentType;
61 typedef typename VectorialType::size_type size_type;
63 static const size_type factor = dimR;
65 static ComponentType begin () {
return ComponentType( 0 ); }
66 static ComponentType end () {
return ComponentType( factor ); }
68 static VectorialType zeroVectorial () {
return VectorialType( K( field_type( 0 ) ) ); }
70 static const K &access (
const ScalarType &x ) {
return x[ 0 ]; }
71 static K &access ( ScalarType &x ) {
return x[ 0 ]; }
73 static const K &access (
const VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
74 static K &access ( VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
76 static size_type index (
const ComponentType &i ) {
return i; }
79 template<
class K,
int dimR,
int dimD >
80 struct MakeVectorialTraits< FieldMatrix< K, 1, dimD >, FieldMatrix< K, dimR, dimD > >
82 typedef FieldMatrix< K, 1, dimD > ScalarType;
83 typedef FieldMatrix< K, dimR, dimD > VectorialType;
85 typedef typename FieldTraits< VectorialType >::field_type field_type;
86 typedef typename VectorialType::size_type ComponentType;
87 typedef typename VectorialType::size_type size_type;
89 static const size_type factor = dimR;
91 static ComponentType begin () {
return ComponentType( 0 ); }
92 static ComponentType end () {
return ComponentType( factor ); }
94 static VectorialType zeroVectorial () {
return VectorialType( K( field_type( 0 ) ) ); }
96 static const FieldVector< K, dimD > &access (
const ScalarType &x ) {
return x[ 0 ]; }
97 static FieldVector< K, dimD > &access ( ScalarType &x ) {
return x[ 0 ]; }
99 static const FieldVector< K, dimD > &access (
const VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
100 static FieldVector< K, dimD > &access ( VectorialType &x,
const ComponentType &i ) {
return x[ i ]; }
102 static size_type index (
const ComponentType &i ) {
return i; }
109 template<
class Scalar,
class Vectorial >
110 class BasicMakeVectorialExpression
112 typedef BasicMakeVectorialExpression< Scalar, Vectorial > ThisType;
114 typedef MakeVectorialTraits< Scalar, Vectorial > Traits;
117 typedef typename Traits::ScalarType ScalarType;
118 typedef typename Traits::VectorialType VectorialType;
120 typedef typename Traits::field_type field_type;
121 typedef typename Traits::ComponentType ComponentType;
122 typedef typename Traits::size_type size_type;
124 BasicMakeVectorialExpression (
const ComponentType &component,
const ScalarType &scalar )
125 : component_( component ),
129 operator VectorialType ()
const
131 VectorialType vectorial = Traits::zeroVectorial();
132 Traits::access( vectorial, component() ) = Traits::access( scalar() );
136 const ThisType &operator*= (
const field_type &s )
142 const ThisType &operator/= (
const field_type &s )
148 const ComponentType &component ()
const {
return component_; }
150 const ScalarType &scalar ()
const {
return scalar_; }
151 ScalarType &scalar () {
return scalar_; }
154 ComponentType component_;
163 template<
class Scalar,
class Vectorial >
164 class MakeVectorialExpression
165 :
public BasicMakeVectorialExpression< Scalar, Vectorial >
167 typedef MakeVectorialExpression< Scalar, Vectorial > ThisType;
168 typedef BasicMakeVectorialExpression< Scalar, Vectorial > BaseType;
171 typedef typename BaseType::ComponentType ComponentType;
172 typedef typename BaseType::ScalarType ScalarType;
174 MakeVectorialExpression (
const ComponentType &component,
const ScalarType &scalar )
175 : BaseType( component, scalar )
179 template<
class K,
int dimR >
180 class MakeVectorialExpression< FieldVector< K, 1 >, FieldVector< K, dimR > >
181 :
public BasicMakeVectorialExpression< FieldVector< K, 1 >, FieldVector< K, dimR > >
183 typedef MakeVectorialExpression< FieldVector< K, 1 >, FieldVector< K, dimR > > ThisType;
184 typedef BasicMakeVectorialExpression< FieldVector< K, 1 >, FieldVector< K, dimR > > BaseType;
187 typedef typename BaseType::ScalarType ScalarType;
188 typedef typename BaseType::VectorialType VectorialType;
190 typedef typename BaseType::field_type field_type;
191 typedef typename BaseType::ComponentType ComponentType;
192 typedef typename BaseType::size_type size_type;
194 using BaseType::component;
195 using BaseType::scalar;
197 MakeVectorialExpression (
const ComponentType &component,
const ScalarType &scalar )
198 : BaseType( component, scalar )
201 field_type operator* (
const ThisType &other )
const
203 return (component() == other.component() ? scalar() * other.scalar() : field_type( 0 ));
206 field_type operator* (
const VectorialType &other )
const
208 return (scalar()[ 0 ] * other[ component() ]);
211 field_type one_norm ()
const {
return scalar().one_norm(); }
212 field_type two_norm ()
const {
return scalar().two_norm(); }
213 field_type two_norm2 ()
const {
return scalar().two_norm2(); }
214 field_type infinity_norm ()
const {
return scalar().infinity_norm(); }
216 size_type
size ()
const {
return dimR; }
218 friend field_type operator* (
const VectorialType &a, ThisType &b ) {
return b*a; }
221 template<
class K,
int dimR,
int dimD >
222 class MakeVectorialExpression< FieldMatrix< K, 1, dimD >, FieldMatrix< K, dimR, dimD > >
223 :
public BasicMakeVectorialExpression< FieldMatrix< K, 1, dimD >, FieldMatrix< K, dimR, dimD > >
225 typedef MakeVectorialExpression< FieldMatrix< K, 1, dimD >, FieldMatrix< K, dimR, dimD > > ThisType;
226 typedef BasicMakeVectorialExpression< FieldMatrix< K, 1, dimD >, FieldMatrix< K, dimR, dimD > > BaseType;
229 typedef typename BaseType::ScalarType ScalarType;
230 typedef typename BaseType::VectorialType VectorialType;
232 typedef typename BaseType::field_type field_type;
233 typedef typename BaseType::ComponentType ComponentType;
234 typedef typename BaseType::size_type size_type;
236 using BaseType::component;
237 using BaseType::scalar;
239 MakeVectorialExpression (
const ComponentType &component,
const ScalarType &scalar )
240 : BaseType( component, scalar )
243 template<
class X,
class Y >
244 void mv (
const X &x, Y &y )
const
246 for( size_type i = 0; i < rows(); ++i )
247 y[ i ] = field_type( 0 );
248 for( size_type j= 0; j < cols(); ++j )
249 y[ component() ] += scalar()[ component() ][ j ] * x[ j ];
252 template<
class X,
class Y >
253 void mtv (
const X &x, Y &y )
const
255 for( size_type i = 0; i < rows(); ++i )
256 y[ i ] = scalar()[ i ][ component() ] * x[ component() ];
259 template<
class X,
class Y >
260 void umv (
const X &x, Y &y )
const
262 for( size_type j= 0; j < cols(); ++j )
263 y[ component() ] += scalar()[ component() ][ j ] * x[ j ];
266 template<
class X,
class Y >
267 void umtv (
const X &x, Y &y )
const
269 for( size_type i = 0; i < rows(); ++i )
270 y[ i ] += scalar()[ i ][ component() ] * x[ component() ];
273 template<
class X,
class Y >
274 void mmv (
const X &x, Y &y )
const
276 for( size_type j= 0; j < cols(); ++j )
277 y[ component() ] -= scalar()[ component() ][ j ] * x[ j ];
280 template<
class X,
class Y >
281 void mmtv (
const X &x, Y &y )
const
283 for( size_type i = 0; i < rows(); ++i )
284 y[ i ] -= scalar()[ i ][ component() ] * x[ component() ];
287 field_type frobenius_norm ()
const {
return scalar().frobenius_norm(); }
288 field_type frobenius_norm2 ()
const {
return scalar().frobenius_norm2(); }
289 field_type infinity_norm ()
const {
return scalar().infinity_norm(); }
291 field_type determinant ()
const {
return (dimR == 1 ? scalar().determinant() : field_type( 0 )); }
293 size_type N ()
const {
return rows(); }
294 size_type M ()
const {
return cols(); }
296 size_type rows ()
const {
return dimR; }
297 size_type cols ()
const {
return scalar().cols(); }
305 template<
class Scalar,
class Vectorial >
307 operator== (
const MakeVectorialExpression< Scalar, Vectorial > &a,
308 const MakeVectorialExpression< Scalar, Vectorial > &b )
310 return ((a.component() == b.component()) && (a.scalar() == b.scalar()));
313 template<
class Scalar,
class Vectorial >
315 operator!= (
const MakeVectorialExpression< Scalar, Vectorial > &a,
316 const MakeVectorialExpression< Scalar, Vectorial > &b )
318 return ((a.component() != b.component()) || (a.scalar() != b.scalar()));
321 template<
class Scalar,
class Vectorial >
324 const MakeVectorialExpression< Scalar, Vectorial > &b )
326 return (a ==
static_cast< Vectorial
>( b ));
329 template<
class Scalar,
class Vectorial >
332 const MakeVectorialExpression< Scalar, Vectorial > &b )
334 return (a !=
static_cast< Vectorial
>( b ));
337 template<
class Scalar,
class Vectorial >
339 operator== (
const MakeVectorialExpression< Scalar, Vectorial > &a,
342 return (
static_cast< Vectorial
>( a ) == b);
345 template<
class Scalar,
class Vectorial >
347 operator!= (
const MakeVectorialExpression< Scalar, Vectorial > &a,
350 return (
static_cast< Vectorial
>( a ) != b);
353 template<
class Scalar,
class Vectorial >
355 axpy (
const typename MakeVectorialTraits< Scalar, Vectorial >::field_type &a,
356 const MakeVectorialExpression< Scalar, Vectorial > &x,
357 typename MakeVectorialTraits< Scalar, Vectorial >::VectorialType &y )
359 typedef MakeVectorialTraits< Scalar, Vectorial > Traits;
360 axpy( a, Traits::access( x.scalar() ), Traits::access( y, x.component() ) );
363 template<
class GeometryJacobianInverseTransposed,
class K,
int ROWS >
364 void jacobianTransformation (
const GeometryJacobianInverseTransposed &gjit,
365 const MakeVectorialExpression< FieldMatrix< K, 1, GeometryJacobianInverseTransposed::cols >, FieldMatrix< K, ROWS, GeometryJacobianInverseTransposed::cols > > &a,
366 FieldMatrix< K, ROWS, GeometryJacobianInverseTransposed::rows > &b )
368 typedef MakeVectorialTraits< FieldMatrix< K, 1, GeometryJacobianInverseTransposed::cols >, FieldMatrix< K, ROWS, GeometryJacobianInverseTransposed::cols > > Traits;
369 typedef MakeVectorialTraits< FieldMatrix< K, 1, GeometryJacobianInverseTransposed::rows >, FieldMatrix< K, ROWS, GeometryJacobianInverseTransposed::rows > > RgTraits;
370 b = RgTraits::zeroVectorial();
371 gjit.mv( Traits::access( a.scalar() ), b[ a.component() ] );
374 template<
class GeometryJacobianInverseTransposed,
class K,
int SIZE >
375 void hessianTransformation (
const GeometryJacobianInverseTransposed &gjit,
376 const MakeVectorialExpression< ExplicitFieldVector< FieldMatrix< K, GeometryJacobianInverseTransposed::cols, GeometryJacobianInverseTransposed::cols >, 1 >, ExplicitFieldVector< FieldMatrix< K, GeometryJacobianInverseTransposed::cols, GeometryJacobianInverseTransposed::cols >, SIZE > > &a,
377 ExplicitFieldVector< FieldMatrix< K, GeometryJacobianInverseTransposed::rows, GeometryJacobianInverseTransposed::rows >, SIZE > &b )
379 const int dimLocal = GeometryJacobianInverseTransposed::cols;
380 const int dimGlobal = GeometryJacobianInverseTransposed::rows;
381 typedef MakeVectorialTraits< FieldVector< FieldMatrix< K, dimLocal, dimLocal >, 1 >, FieldVector< FieldMatrix< K, dimLocal, dimLocal >, SIZE > > Traits;
382 typedef MakeVectorialTraits< FieldVector< FieldMatrix< K, dimGlobal, dimGlobal >, 1 >, FieldVector< FieldMatrix< K, dimGlobal, dimGlobal >, SIZE > > RgTraits;
384 b = RgTraits::zeroVectorial();
387 FieldMatrix< K, dimLocal, dimGlobal > c;
388 for(
int i = 0; i < dimLocal; ++i )
389 gjit.mv( Traits::access( a.scalar() )[ i ], c[ i ] );
392 for(
int i = 0; i < dimGlobal; ++i )
394 FieldMatrixColumn< const FieldMatrix< K, dimLocal, dimGlobal > > ci( c, i );
395 FieldMatrixColumn< FieldMatrix< K, dimGlobal, dimGlobal > > bi( RgTraits::access( b, a.component() ), i );
401 template<
class Scalar,
class Vectorial >
402 inline typename MakeVectorialTraits< Scalar, Vectorial >::field_type
403 scalarProduct (
const MakeVectorialExpression< Scalar, Vectorial > &a,
406 return scalarProduct( a.scalar()[ 0 ], b[ a.component() ] );
409 template<
class Scalar,
class Vectorial >
410 inline typename MakeVectorialTraits< Scalar, Vectorial >::field_type
411 scalarProduct (
const Vectorial &a,
412 const MakeVectorialExpression< Scalar, Vectorial > &b )
414 return scalarProduct( a[ b.component() ], b.scalar()[ 0 ] );
417 template<
class Scalar,
class Vectorial >
418 inline typename MakeVectorialTraits< Scalar, Vectorial >::field_type
419 scalarProduct (
const MakeVectorialExpression< Scalar, Vectorial > &a,
420 const MakeVectorialExpression< Scalar, Vectorial > &b )
422 typedef typename MakeVectorialTraits< Scalar, Vectorial >::field_type field_type;
423 return (a.component() == b.component() ? scalarProduct( a.scalar(), b.scalar() ) : field_type( 0 ));
431 template<
class ScalarFunctionSpace,
class RangeVector >
434 template<
class DomainField,
class RangeField,
int dimD,
int dimR >
435 struct ToNewRange< FunctionSpace< DomainField, RangeField, dimD, 1 >, FieldVector< RangeField, dimR > >
437 typedef FunctionSpace< DomainField, RangeField, dimD, dimR > Type;
445 template<
class ScalarShapeFunctionSet,
class RangeVector >
446 class VectorialShapeFunctionSet
448 typedef VectorialShapeFunctionSet< ScalarShapeFunctionSet, RangeVector > ThisType;
451 typedef ScalarShapeFunctionSet ScalarShapeFunctionSetType;
455 static constexpr bool codegenShapeFunctionSet = detail::IsCodegenShapeFunctionSet< ScalarShapeFunctionSetType >::value;
456 static const int pointSetId = detail::SelectPointSetId< ScalarShapeFunctionSetType >::value;
459 typedef typename ScalarShapeFunctionSetType::FunctionSpaceType ScalarFunctionSpaceType;
461 static const std::size_t dimRangeFactor = MakeVectorialTraits< typename ScalarFunctionSpaceType::RangeType, RangeVector >::factor;
463 template<
class Functor,
class Vectorial >
464 struct VectorialFunctor;
467 typedef typename ToNewRange< ScalarFunctionSpaceType, RangeVector >::Type
FunctionSpaceType;
473 template<
class ... Args >
474 VectorialShapeFunctionSet ( Args&& ... args )
475 : scalarShapeFunctionSet_(
std::forward< Args > ( args ) ... )
478 explicit VectorialShapeFunctionSet (
const ScalarShapeFunctionSetType &scalarShapeFunctionSet )
479 : scalarShapeFunctionSet_( scalarShapeFunctionSet )
482 const ScalarShapeFunctionSetType &scalarShapeFunctionSet ()
const {
return scalarShapeFunctionSet_; }
484 int order ()
const {
return scalarShapeFunctionSet().order(); }
487 std::size_t
size ()
const {
return dimRangeFactor * scalarShapeFunctionSet().size(); }
489 template<
class Po
int,
class Functor >
490 void evaluateEach (
const Point &x, Functor functor )
const;
492 template<
class Po
int,
class Functor >
493 void jacobianEach (
const Point &x, Functor functor )
const;
495 template<
class Po
int,
class Functor >
496 void hessianEach (
const Point &x, Functor functor )
const;
499 ScalarShapeFunctionSet scalarShapeFunctionSet_;
507 template<
class ScalarShapeFunctionSet,
class RangeVector >
508 template<
class Functor,
class Vectorial >
509 struct VectorialShapeFunctionSet< ScalarShapeFunctionSet, RangeVector >::VectorialFunctor
511 explicit VectorialFunctor (
const Functor &functor )
512 : functor_( functor )
515 template<
class Scalar >
516 void operator() (
const std::size_t i,
const Scalar &value )
518 typedef MakeVectorialTraits< Scalar, Vectorial > Traits;
519 typedef MakeVectorialExpression< Scalar, Vectorial > Expression;
520 const typename Traits::ComponentType end = Traits::end();
521 for(
typename Traits::ComponentType k = Traits::begin(); k != end; ++k )
522 functor_( i*Traits::factor + Traits::index( k ), Expression( k, value ) );
534 template<
class ScalarShapeFunctionSet,
class RangeVector >
535 template<
class Po
int,
class Functor >
536 inline void VectorialShapeFunctionSet< ScalarShapeFunctionSet, RangeVector >
537 ::evaluateEach (
const Point &x, Functor functor )
const
540 scalarShapeFunctionSet().evaluateEach( x, VectorialFunctor< Functor, VectorialType >( functor ) );
544 template<
class ScalarShapeFunctionSet,
class RangeVector >
545 template<
class Po
int,
class Functor >
546 inline void VectorialShapeFunctionSet< ScalarShapeFunctionSet, RangeVector >
547 ::jacobianEach (
const Point &x, Functor functor )
const
550 scalarShapeFunctionSet().jacobianEach( x, VectorialFunctor< Functor, VectorialType >( functor ) );
554 template<
class ScalarShapeFunctionSet,
class RangeVector >
555 template<
class Po
int,
class Functor >
556 inline void VectorialShapeFunctionSet< ScalarShapeFunctionSet, RangeVector >
557 ::hessianEach (
const Point &x, Functor functor )
const
560 scalarShapeFunctionSet().hessianEach( x, VectorialFunctor< Functor, VectorialType >( functor ) );
Definition: explicitfieldvector.hh:75
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
FunctionSpaceTraits::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
A vector valued function space.
Definition: functionspace.hh:60
EnableIfInterOperable< T1, T2, bool >::type operator==(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for equality.
Definition: iteratorfacades.hh:238
EnableIfInterOperable< T1, T2, bool >::type operator!=(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for inequality.
Definition: iteratorfacades.hh:260
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75