diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 34b8f23..d273c8f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,7 +5,6 @@ set(HEADERS Picasso.hpp Picasso_AdaptiveMesh.hpp Picasso_APIC.hpp - Picasso_BatchedLinearAlgebra.hpp Picasso_BilinearMeshMapping.hpp Picasso_CurvilinearMesh.hpp Picasso_FacetGeometry.hpp diff --git a/src/Picasso_APIC.hpp b/src/Picasso_APIC.hpp index 52dc41a..4877ad0 100644 --- a/src/Picasso_APIC.hpp +++ b/src/Picasso_APIC.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_APIC_HPP #define PICASSO_APIC_HPP -#include #include #include @@ -91,7 +90,7 @@ KOKKOS_INLINE_FUNCTION void p2g( constexpr int ncomp = ParticleField::extent_1; // Affine Matrix - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Matrix B_p; for ( int d = 0; d < ncomp; ++d ) { // B_p is sliced from c_p and is transposed. @@ -105,7 +104,7 @@ KOKKOS_INLINE_FUNCTION void p2g( value_type D_p_inv = inertialScaling( sd ); // Project momentum. - Vec3 distance; + Cabana::Vec3 distance; value_type wm_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -187,7 +186,7 @@ KOKKOS_INLINE_FUNCTION void p2g( constexpr int ncomp = ParticleVelocity::extent_1; // Affine Matrix - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Matrix B_p; for ( int d = 0; d < ncomp; ++d ) { // B_p is sliced from c_p and is transposed. @@ -204,7 +203,7 @@ KOKKOS_INLINE_FUNCTION void p2g( value_type D_p_inv = inertialScaling( sd ); // Project momentum. - Vec3 distance; + Cabana::Vec3 distance; value_type wm_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -273,7 +272,7 @@ KOKKOS_INLINE_FUNCTION void p2g( constexpr int ncomp = ParticleField::extent_1; // Affine Matrix - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Matrix B_p; for ( int d = 0; d < ncomp; ++d ) { // B_p is sliced from c_p and is transposed. @@ -282,7 +281,7 @@ KOKKOS_INLINE_FUNCTION void p2g( B_p( d, 2 ) = c_p( 3, d ); } // Project momentum. - Vec3 gm_ip; + Cabana::Vec3 gm_ip; value_type wm_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -360,7 +359,7 @@ KOKKOS_INLINE_FUNCTION void p2g( "APIC requires 3 space dimensions" ); // Affine Matrix - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Matrix B_p; for ( int d = 0; d < 3; ++d ) { // B_p is sliced from c_p and is transposed. @@ -372,7 +371,7 @@ KOKKOS_INLINE_FUNCTION void p2g( const int dim = SplineDataType::entity_type::dim; // Project momentum. - Vec3 gm_ip; + Cabana::Vec3 gm_ip; value_type wm_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -426,15 +425,15 @@ g2p( const GridField& u_i, ParticleField& c_p, const SplineDataType& sd, constexpr int ncomp = ParticleField::extent_1; // Affine Matrix - LinearAlgebra::Vector u_p; - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Vector u_p; + Cabana::LinearAlgebra::Matrix B_p; // Reset the particle values. u_p = 0.0; B_p = 0.0; // Update particle. - Vec3 distance; + Cabana::Vec3 distance; value_type w_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -496,8 +495,8 @@ g2p( const GridMomentum& u_i, ParticleVelocity& c_p, const SplineDataType& sd, "APIC requires 3 space dimensions" ); // Affine Matrix - LinearAlgebra::Vector u_p; - LinearAlgebra::Matrix B_p; + Cabana::LinearAlgebra::Vector u_p; + Cabana::LinearAlgebra::Matrix B_p; // Get the field dimension we are working on. const int dim = SplineDataType::entity_type::dim; @@ -507,7 +506,7 @@ g2p( const GridMomentum& u_i, ParticleVelocity& c_p, const SplineDataType& sd, B_p = 0.0; // Update particle. - Vec3 distance; + Cabana::Vec3 distance; value_type w_ip; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) diff --git a/src/Picasso_BatchedLinearAlgebra.hpp b/src/Picasso_BatchedLinearAlgebra.hpp index f791dce..37122bf 100644 --- a/src/Picasso_BatchedLinearAlgebra.hpp +++ b/src/Picasso_BatchedLinearAlgebra.hpp @@ -116,18 +116,6 @@ namespace LinearAlgebra //---------------------------------------------------------------------------// // Forward declarations. //---------------------------------------------------------------------------// -template -struct MatrixExpression; -template -struct Matrix; -template -struct MatrixView; -template -struct VectorExpression; -template -struct Vector; -template -struct VectorView; template struct QuaternionExpression; template @@ -150,58 +138,6 @@ struct Tensor4Expression; //---------------------------------------------------------------------------// // Type traits. //---------------------------------------------------------------------------// -// Matrix -template -struct is_matrix_impl : public std::false_type -{ -}; - -template -struct is_matrix_impl> : public std::true_type -{ -}; - -template -struct is_matrix_impl> : public std::true_type -{ -}; - -template -struct is_matrix_impl> : public std::true_type -{ -}; - -template -struct is_matrix : public is_matrix_impl::type>::type -{ -}; - -// Vector -template -struct is_vector_impl : public std::false_type -{ -}; - -template -struct is_vector_impl> : public std::true_type -{ -}; - -template -struct is_vector_impl> : public std::true_type -{ -}; - -template -struct is_vector_impl> : public std::true_type -{ -}; - -template -struct is_vector : public is_vector_impl::type>::type -{ -}; - // Quaternion template struct is_quaternion_impl : public std::false_type @@ -304,22 +240,6 @@ createTensor3Expression( const Func& f ) return Tensor3Expression( f ); } -// Matrix -template -KOKKOS_INLINE_FUNCTION MatrixExpression -createMatrixExpression( const Func& f ) -{ - return MatrixExpression( f ); -} - -// Vector. -template -KOKKOS_INLINE_FUNCTION VectorExpression -createVectorExpression( const Func& f ) -{ - return VectorExpression( f ); -} - // Quaternion template KOKKOS_INLINE_FUNCTION QuaternionExpression @@ -382,7 +302,7 @@ struct Tensor4Expression auto vector( const ALL_INDEX_t, const int n, const int p, const int q ) const { - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression( [=]( const int i ) { return ( *this )( i, n, p, q ); } ); } // Get a row as a vector view. @@ -390,7 +310,7 @@ struct Tensor4Expression auto vector( const int m, const ALL_INDEX_t, const int p, const int q ) const { - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression( [=]( const int i ) { return ( *this )( m, i, p, q ); } ); } // Get a row as a vector view. @@ -398,7 +318,7 @@ struct Tensor4Expression auto vector( const int m, const int n, const ALL_INDEX_t, const int q ) const { - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression( [=]( const int i ) { return ( *this )( m, n, i, q ); } ); } @@ -407,7 +327,7 @@ struct Tensor4Expression auto vector( const int m, const int n, const int p, const ALL_INDEX_t ) const { - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression( [=]( const int i ) { return ( *this )( m, n, p, i ); } ); } @@ -416,7 +336,7 @@ struct Tensor4Expression auto matrix( const ALL_INDEX_t, const ALL_INDEX_t, const int p, const int q ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( i, j, p, q ); } ); } @@ -425,7 +345,7 @@ struct Tensor4Expression auto matrix( const ALL_INDEX_t, const int n, const ALL_INDEX_t, const int q ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( i, n, j, q ); } ); } @@ -434,7 +354,7 @@ struct Tensor4Expression auto matrix( const int m, const ALL_INDEX_t, const ALL_INDEX_t, const int q ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( m, i, j, q ); } ); } @@ -443,7 +363,7 @@ struct Tensor4Expression auto matrix( const ALL_INDEX_t, const int n, const int p, const ALL_INDEX_t ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( i, n, p, j ); } ); } @@ -452,7 +372,7 @@ struct Tensor4Expression auto matrix( const int m, const ALL_INDEX_t, const int p, const ALL_INDEX_t ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( m, i, p, j ); } ); } @@ -461,7 +381,7 @@ struct Tensor4Expression auto matrix( const int m, const int n, const ALL_INDEX_t, const ALL_INDEX_t ) const { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression( [=]( const int i, const int j ) { return ( *this )( m, n, i, j ); } ); } @@ -473,1474 +393,163 @@ struct Tensor4Expression { return createTensor3Expression( [=]( const int i, const int j, const int k ) - { return ( *this )( i, j, k, b ); } ); - } - // Get a tensor3 as a Tensor3 view. - KOKKOS_INLINE_FUNCTION - auto tensor3( const ALL_INDEX_t, const ALL_INDEX_t, const int b, - const ALL_INDEX_t ) - { - return createTensor3Expression( - [=]( const int i, const int j, const int k ) - { return ( *this )( i, j, b, k ); } ); - } - // Get a tensor3 as a Tensor3 view. - KOKKOS_INLINE_FUNCTION - auto tensor3( const ALL_INDEX_t, const int b, const ALL_INDEX_t, - const ALL_INDEX_t ) - { - return createTensor3Expression( - [=]( const int i, const int j, const int k ) - { return ( *this )( i, b, j, k ); } ); - } - // Get a tensor3 as a Tensor3 view. - KOKKOS_INLINE_FUNCTION - auto tensor3( const int b, const ALL_INDEX_t, const ALL_INDEX_t, - const ALL_INDEX_t ) - { - return createTensor3Expression( - [=]( const int i, const int j, const int k ) - { return ( *this )( b, i, j, k ); } ); - } -}; - -// Tensor3 expression container. -template -struct Tensor3Expression -{ - static constexpr int extent_0 = M; - static constexpr int extent_1 = N; - static constexpr int extent_2 = P; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - - using eval_type = Tensor3; - using copy_type = Tensor3; - - Func _f; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Tensor3Expression() = default; - - // Create an expression from a callable object. - KOKKOS_INLINE_FUNCTION - Tensor3Expression( const Func& f ) - : _f( f ) - { - } - - // Extent. - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 2 ? extent_2 - : ( d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ) ); - } - - // Evaluate the expression at an index. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i, const int j, const int k ) const - { - return _f( i, j, k ); - } - - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - auto vector( const ALL_INDEX_t, const int n, const int p ) const - { - return createVectorExpression( [=]( const int i ) - { return ( *this )( i, n, p ); } ); - } - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - auto vector( const int m, const ALL_INDEX_t, const int p ) const - { - return createVectorExpression( [=]( const int i ) - { return ( *this )( m, i, p ); } ); - } - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - auto vector( const int m, const int n, const ALL_INDEX_t ) const - { - return createVectorExpression( [=]( const int i ) - { return ( *this )( m, n, i ); } ); - } - - // Get a matrix as a matrix view. - KOKKOS_INLINE_FUNCTION - auto matrix( const ALL_INDEX_t, const ALL_INDEX_t, const int p ) const - { - return createMatrixExpression( - [=]( const int i, const int j ) { return ( *this )( i, j, p ); } ); - } - // Get a matrix as a matrix view. - KOKKOS_INLINE_FUNCTION - auto matrix( const ALL_INDEX_t, const int n, const ALL_INDEX_t ) const - { - return createMatrixExpression( - [=]( const int i, const int j ) { return ( *this )( i, n, j ); } ); - } - // Get a matrix as a matrix view. - KOKKOS_INLINE_FUNCTION - auto matrix( const int m, const ALL_INDEX_t, const ALL_INDEX_t ) const - { - return createMatrixExpression( - [=]( const int i, const int j ) { return ( *this )( m, i, j ); } ); - } -}; - -// Matrix expression container. -template -struct MatrixExpression -{ - static constexpr int extent_0 = M; - static constexpr int extent_1 = N; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - - using eval_type = Matrix; - using copy_type = Matrix; - - Func _f; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - MatrixExpression() = default; - - // Create an expression from a callable object. - KOKKOS_INLINE_FUNCTION - MatrixExpression( const Func& f ) - : _f( f ) - { - } - - // Extent. - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ); - } - - // Evaluate the expression at an index. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i, const int j ) const - { - return _f( i, j ); - } - - // Get a row as a vector expression. - KOKKOS_INLINE_FUNCTION - auto row( const int n ) const - { - return createVectorExpression( [=]( const int i ) - { return ( *this )( n, i ); } ); - } - - // Get a column as a vector expression. - KOKKOS_INLINE_FUNCTION - auto column( const int n ) const - { - return createVectorExpression( [=]( const int i ) - { return ( *this )( i, n ); } ); - } -}; - -//---------------------------------------------------------------------------// -// Vector expression container. -template -struct VectorExpression -{ - static constexpr int extent_0 = N; - static constexpr int extent_1 = 1; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - - using eval_type = Vector; - using copy_type = Vector; - - Func _f; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - VectorExpression() = default; - - // Create an expression from a callable object. - KOKKOS_INLINE_FUNCTION - VectorExpression( const Func& f ) - : _f( f ) - { - } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int ) const { return N; } - - // Evaluate the expression at an index. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i ) const { return _f( i ); } - - // Evaluate the expression at an index. 2D version for vectors treated as - // matrices. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i, int ) const { return _f( i ); } -}; - -//---------------------------------------------------------------------------// -// Quaternion expression container. -template -struct QuaternionExpression -{ - static constexpr int extent_0 = 4; - static constexpr int extent_1 = 1; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - - using eval_type = Quaternion; - using copy_type = Quaternion; - - Func _f; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - QuaternionExpression() = default; - - // Create an expression from a callable object. - KOKKOS_INLINE_FUNCTION - QuaternionExpression( const Func& f ) - : _f( f ) - { - } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int ) const { return extent_0; } - - // Evaluate the expression at an index. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i ) const { return _f( i ); } - - // Evaluate the expression at an index. 2D version for vectors treated as - // matrices. - KOKKOS_INLINE_FUNCTION - value_type operator()( const int i, int ) const { return _f( i ); } -}; - -//---------------------------------------------------------------------------// -// Matrix -//---------------------------------------------------------------------------// -// Dense matrix. -template -struct Matrix -{ - T _d[M][N]; - - static constexpr int extent_0 = M; - static constexpr int extent_1 = N; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = MatrixView; - using copy_type = Matrix; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Matrix() = default; - - // Initializer list constructor. - KOKKOS_INLINE_FUNCTION - Matrix( const std::initializer_list> data ) - { - int i = 0; - int j = 0; - for ( const auto& row : data ) - { - j = 0; - for ( const auto& value : row ) - { - _d[i][j] = value; - ++j; - } - ++i; - } - } - - // Deep copy constructor. Triggers expression evaluation. - template < - class Expression, - typename std::enable_if::value, int>::type = 0> - KOKKOS_INLINE_FUNCTION Matrix( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = e( i, j ); - } - - // Scalar constructor. - KOKKOS_INLINE_FUNCTION - Matrix( const T value ) - { - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = value; - } - - // Assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = e( i, j ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) += e( i, j ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) -= e( i, j ); - return *this; - } - - // Initializer list assignment operator. - KOKKOS_INLINE_FUNCTION - Matrix& - operator=( const std::initializer_list> data ) - { - int i = 0; - int j = 0; - for ( const auto& row : data ) - { - j = 0; - for ( const auto& value : row ) - { - _d[i][j] = value; - ++j; - } - ++i; - } - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - Matrix& operator=( const T value ) - { - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = value; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return N; } - - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride( const int d ) const { return ( 0 == d ) ? N : 1; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ); - } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i, const int j ) const - { - return _d[i][j]; - } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int i, const int j ) { return _d[i][j]; } - - // Get the raw data. - KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( &_d[0][0] ); } - - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView row( const int n ) const - { - return VectorView( const_cast( &_d[n][0] ), 1 ); - } - - // Get a column as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView column( const int n ) const - { - return VectorView( const_cast( &_d[0][n] ), N ); - } -}; - -//---------------------------------------------------------------------------// -// Scalar overload. -template -struct Matrix -{ - T _d; - - static constexpr int extent_0 = 1; - static constexpr int extent_1 = 1; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = MatrixView; - using copy_type = Matrix; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Matrix() = default; - - // Scalar constructor. - KOKKOS_INLINE_FUNCTION - Matrix( const T value ) - : _d( value ) - { - } - - // Deep copy constructor. Triggers expression evaluation. - template < - class Expression, - typename std::enable_if::value, int>::type = 0> - KOKKOS_INLINE_FUNCTION Matrix( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d = e( 0, 0 ); - } - - // Assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d = e( 0, 0 ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d += e( 0, 0 ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d -= e( 0, 0 ); - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - Matrix& operator=( const T value ) - { - _d = value; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride( const int ) const { return 1; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int ) const { return 1; } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int, const int ) const { return _d; } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int, const int ) { return _d; } - - // Get the raw data. - KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( &_d ); } - - // Scalar conversion operator. - KOKKOS_INLINE_FUNCTION - operator value_type() const { return _d; } -}; - -template -struct Matrix -{ - T _d[3][3]; - - static constexpr int extent_0 = 3; - static constexpr int extent_1 = 3; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = MatrixView; - using copy_type = Matrix; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Matrix() = default; - - // Initializer list constructor. - KOKKOS_INLINE_FUNCTION - Matrix( const std::initializer_list> data ) - { - int i = 0; - int j = 0; - for ( const auto& row : data ) - { - j = 0; - for ( const auto& value : row ) - { - _d[i][j] = value; - ++j; - } - ++i; - } - } - - // Quaternion to matrix explicit constructor - explicit KOKKOS_INLINE_FUNCTION Matrix( const Quaternion& q ) - { - value_type n = q( 0 ) * q( 0 ) + q( 1 ) * q( 1 ) + q( 2 ) * q( 2 ) + - q( 3 ) * q( 3 ); - value_type s = n == 0 ? 0 : 2.0 / n; - - _d[0][0] = 1.0 - s * ( q( 2 ) * q( 2 ) + q( 3 ) * q( 3 ) ); - _d[0][1] = s * ( q( 1 ) * q( 2 ) - q( 0 ) * q( 3 ) ); - _d[0][2] = s * ( q( 1 ) * q( 3 ) + q( 0 ) * q( 2 ) ); - _d[1][0] = s * ( q( 1 ) * q( 2 ) + q( 0 ) * q( 3 ) ); - _d[1][1] = 1.0 - s * ( q( 1 ) * q( 1 ) + q( 3 ) * q( 3 ) ); - _d[1][2] = s * ( q( 2 ) * q( 3 ) - q( 0 ) * q( 1 ) ); - _d[2][0] = s * ( q( 1 ) * q( 3 ) - q( 0 ) * q( 2 ) ); - _d[2][1] = s * ( q( 2 ) * q( 3 ) + q( 0 ) * q( 1 ) ); - _d[2][2] = 1.0 - s * ( q( 1 ) * q( 1 ) + q( 2 ) * q( 2 ) ); - } - - // Deep copy constructor. Triggers expression evaluation. - template < - class Expression, - typename std::enable_if::value, int>::type = 0> - KOKKOS_INLINE_FUNCTION Matrix( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) = e( i, j ); - } - - // Scalar constructor. - KOKKOS_INLINE_FUNCTION - Matrix( const T value ) - { - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) = value; - } - - // Assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) = e( i, j ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) += e( i, j ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Matrix&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) -= e( i, j ); - return *this; - } - - // Initializer list assignment operator. - KOKKOS_INLINE_FUNCTION - Matrix& - operator=( const std::initializer_list> data ) - { - int i = 0; - int j = 0; - for ( const auto& row : data ) - { - j = 0; - for ( const auto& value : row ) - { - _d[i][j] = value; - ++j; - } - ++i; - } - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - Matrix& operator=( const T value ) - { - for ( int i = 0; i < extent_0; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < extent_1; ++j ) - ( *this )( i, j ) = value; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return extent_1; } - - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride( const int d ) const { return ( 0 == d ) ? extent_1 : 1; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ); - } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i, const int j ) const - { - return _d[i][j]; - } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int i, const int j ) { return _d[i][j]; } - - // Get the raw data. - KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( &_d[0][0] ); } - - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView row( const int n ) const - { - return VectorView( const_cast( &_d[n][0] ), 1 ); - } - - // Get a column as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView column( const int n ) const - { - return VectorView( const_cast( &_d[0][n] ), extent_1 ); - } -}; - -//---------------------------------------------------------------------------// -// View for wrapping matrix data. -// -// NOTE: Data in this view may be non-contiguous. -template -struct MatrixView -{ - T* _d; - int _stride[2]; - - static constexpr int extent_0 = M; - static constexpr int extent_1 = N; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = MatrixView; - using copy_type = Matrix; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - MatrixView() = default; - - // Matrix constructor. - KOKKOS_INLINE_FUNCTION - MatrixView( const Matrix& m ) - : _d( m.data() ) - { - _stride[0] = m.stride_0(); - _stride[1] = m.stride_1(); - } - - // Pointer constructor. - KOKKOS_INLINE_FUNCTION - MatrixView( T* data, const int stride_0, const int stride_1 ) - : _d( data ) - { - _stride[0] = stride_0; - _stride[1] = stride_1; - } - - // Assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, MatrixView&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = e( i, j ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, MatrixView&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) += e( i, j ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, MatrixView&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) -= e( i, j ); - return *this; - } - - // Initializer list assignment operator. - KOKKOS_INLINE_FUNCTION - MatrixView& - operator=( const std::initializer_list> data ) - { - int i = 0; - int j = 0; - for ( const auto& row : data ) - { - j = 0; - for ( const auto& value : row ) - { - ( *this )( i, j ) = value; - ++j; - } - ++i; - } - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - MatrixView& operator=( const T value ) - { - for ( int i = 0; i < M; ++i ) -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int j = 0; j < N; ++j ) - ( *this )( i, j ) = value; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return _stride[0]; } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_1() const { return _stride[1]; } - - KOKKOS_INLINE_FUNCTION - int stride( const int d ) const { return _stride[d]; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ); - } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i, const int j ) const - { - return _d[_stride[0] * i + _stride[1] * j]; - } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int i, const int j ) - { - return _d[_stride[0] * i + _stride[1] * j]; - } - - // Get a row as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView row( const int n ) const - { - return VectorView( const_cast( &_d[_stride[0] * n] ), - _stride[1] ); - } - - // Get a column as a vector view. - KOKKOS_INLINE_FUNCTION - VectorView column( const int n ) const - { - return VectorView( const_cast( &_d[_stride[1] * n] ), - _stride[0] ); - } - - // Get the raw data. - KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( _d ); } -}; - -//---------------------------------------------------------------------------// -// Vector -//---------------------------------------------------------------------------// -// Dense vector. -template -struct Vector -{ - T _d[N]; - - static constexpr int extent_0 = N; - static constexpr int extent_1 = 1; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = VectorView; - using copy_type = Vector; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Vector() = default; - - // Initializer list constructor. - KOKKOS_INLINE_FUNCTION - Vector( const std::initializer_list data ) - { - int i = 0; - for ( const auto& value : data ) - { - _d[i] = value; - ++i; - } - } - - // Deep copy constructor. Triggers expression evaluation. - template < - class Expression, - typename std::enable_if::value, int>::type = 0> - KOKKOS_INLINE_FUNCTION Vector( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = e( i ); - } - - // Scalar constructor. - KOKKOS_INLINE_FUNCTION - Vector( const T value ) - { -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = value; - } - - // Deep copy assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = e( i ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) += e( i ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) -= e( i ); - return *this; - } - - // Initializer list assignment operator. - KOKKOS_INLINE_FUNCTION - Vector& operator=( const std::initializer_list data ) - { - int i = 0; - for ( const auto& value : data ) - { - _d[i] = value; - ++i; - } - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - Vector& operator=( const T value ) - { -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = value; - return *this; - } - - // Normalization - KOKKOS_INLINE_FUNCTION - Vector& normalize() - { -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - auto norm2 = Kokkos::sqrt( ~( *this ) * ( *this ) ); - for ( int i = 0; i < N; ++i ) - ( *this )( i ) /= norm2; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return 1; } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 0; } - - KOKKOS_INLINE_FUNCTION - int stride( const int d ) const { return ( 0 == d ) ? 1 : 0; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? 1 : 0 ); - } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i ) const { return _d[i]; } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int i ) { return _d[i]; } - - // Access an individual element. 2D version for vectors treated as matrices. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i, const int ) const { return _d[i]; } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int i, const int ) { return _d[i]; } - - // Get the raw data. - KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( &_d[0] ); } -}; - -//---------------------------------------------------------------------------// -// Scalar overload. -template -struct Vector -{ - T _d; - - static constexpr int extent_0 = 1; - static constexpr int extent_1 = 1; - - using value_type = T; - using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - - using eval_type = VectorView; - using copy_type = Vector; - - // Default constructor. - KOKKOS_DEFAULTED_FUNCTION - Vector() = default; - - // Scalar constructor. - KOKKOS_INLINE_FUNCTION - Vector( const T value ) - : _d( value ) - { - } - - // Deep copy constructor. Triggers expression evaluation. - template < - class Expression, - typename std::enable_if::value, int>::type = 0> - KOKKOS_INLINE_FUNCTION Vector( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d = e( 0 ); - } - - // Deep copy assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d = e( 0 ); - return *this; - } - - // Addition assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator+=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d += e( 0 ); - return *this; - } - - // Subtraction assignment operator. Triggers expression evaluation. - template - KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, Vector&>::type - operator-=( const Expression& e ) - { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - - _d -= e( 0 ); - return *this; - } - - // Scalar value assignment. - KOKKOS_INLINE_FUNCTION - Vector& operator=( const T value ) - { - _d = value; - return *this; - } - - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_0() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 1; } - - KOKKOS_INLINE_FUNCTION - int stride( const int ) const { return 1; } - - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int ) const { return 1; } - - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int ) const { return _d; } - - KOKKOS_INLINE_FUNCTION - reference operator()( const int ) { return _d; } - - // Access an individual element. 2D version for vectors treated as matrices. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int, const int ) const { return _d; } - + { return ( *this )( i, j, k, b ); } ); + } + // Get a tensor3 as a Tensor3 view. KOKKOS_INLINE_FUNCTION - reference operator()( const int, const int ) { return _d; } - - // Get the raw data. + auto tensor3( const ALL_INDEX_t, const ALL_INDEX_t, const int b, + const ALL_INDEX_t ) + { + return createTensor3Expression( + [=]( const int i, const int j, const int k ) + { return ( *this )( i, j, b, k ); } ); + } + // Get a tensor3 as a Tensor3 view. KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( &_d ); } - - // Scalar conversion operator. + auto tensor3( const ALL_INDEX_t, const int b, const ALL_INDEX_t, + const ALL_INDEX_t ) + { + return createTensor3Expression( + [=]( const int i, const int j, const int k ) + { return ( *this )( i, b, j, k ); } ); + } + // Get a tensor3 as a Tensor3 view. KOKKOS_INLINE_FUNCTION - operator value_type() const { return _d; } + auto tensor3( const int b, const ALL_INDEX_t, const ALL_INDEX_t, + const ALL_INDEX_t ) + { + return createTensor3Expression( + [=]( const int i, const int j, const int k ) + { return ( *this )( b, i, j, k ); } ); + } }; -//---------------------------------------------------------------------------// -// View for wrapping vector data. -// -// NOTE: Data in this view may be non-contiguous. -template -struct VectorView +// Tensor3 expression container. +template +struct Tensor3Expression { - T* _d; - int _stride; - - static constexpr int extent_0 = N; - static constexpr int extent_1 = 1; + static constexpr int extent_0 = M; + static constexpr int extent_1 = N; + static constexpr int extent_2 = P; using value_type = T; using non_const_value_type = typename std::remove_cv::type; - using pointer = T*; - using reference = T&; - using const_reference = const T&; - using eval_type = VectorView; - using copy_type = Vector; + using eval_type = Tensor3; + using copy_type = Tensor3; + + Func _f; // Default constructor. KOKKOS_DEFAULTED_FUNCTION - VectorView() = default; + Tensor3Expression() = default; - // Vector constructor. + // Create an expression from a callable object. KOKKOS_INLINE_FUNCTION - VectorView( const Vector& v ) - : _d( v.data() ) - , _stride( v.stride_0() ) + Tensor3Expression( const Func& f ) + : _f( f ) { } - // Pointer constructor. + // Extent. KOKKOS_INLINE_FUNCTION - VectorView( T* data, const int stride ) - : _d( data ) - , _stride( stride ) + constexpr int extent( const int d ) const { + return d == 2 ? extent_2 + : ( d == 0 ? extent_0 : ( d == 1 ? extent_1 : 0 ) ); } - // Assignment operator. Triggers expression evaluation. - template + // Evaluate the expression at an index. KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, VectorView&>::type - operator=( const Expression& e ) + value_type operator()( const int i, const int j, const int k ) const { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = e( i ); - return *this; + return _f( i, j, k ); } - // Addition assignment operator. Triggers expression evaluation. - template + // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, VectorView&>::type - operator+=( const Expression& e ) + auto vector( const ALL_INDEX_t, const int n, const int p ) const { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) += e( i ); - return *this; + return Cabana::LinearAlgebra::createVectorExpression( + [=]( const int i ) { return ( *this )( i, n, p ); } ); } - - // Subtraction assignment operator. Triggers expression evaluation. - template + // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - typename std::enable_if::value, VectorView&>::type - operator-=( const Expression& e ) + auto vector( const int m, const ALL_INDEX_t, const int p ) const { - static_assert( Expression::extent_0 == extent_0, "Extents must match" ); - static_assert( Expression::extent_1 == extent_1, "Extents must match" ); - -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) -= e( i ); - return *this; + return Cabana::LinearAlgebra::createVectorExpression( + [=]( const int i ) { return ( *this )( m, i, p ); } ); } - - // Initializer list assignment operator. + // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView& operator=( const std::initializer_list data ) + auto vector( const int m, const int n, const ALL_INDEX_t ) const { - int i = 0; - for ( const auto& value : data ) - { - ( *this )( i ) = value; - ++i; - } - return *this; + return Cabana::LinearAlgebra::createVectorExpression( + [=]( const int i ) { return ( *this )( m, n, i ); } ); } - // Scalar value assignment. + // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - VectorView& operator=( const T value ) + auto matrix( const ALL_INDEX_t, const ALL_INDEX_t, const int p ) const { -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - for ( int i = 0; i < N; ++i ) - ( *this )( i ) = value; - return *this; + return Cabana::LinearAlgebra::createMatrixExpression( + [=]( const int i, const int j ) { return ( *this )( i, j, p ); } ); } - - // Normalization + // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - VectorView& normalize() + auto matrix( const ALL_INDEX_t, const int n, const ALL_INDEX_t ) const { -#if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) -#pragma unroll -#endif - auto norm2 = Kokkos::sqrt( ~( *this ) * ( *this ) ); - for ( int i = 0; i < N; ++i ) - ( *this )( i ) /= norm2; - return *this; + return Cabana::LinearAlgebra::createMatrixExpression( + [=]( const int i, const int j ) { return ( *this )( i, n, j ); } ); } - - // Strides. + // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - int stride_0() const { return _stride; } + auto matrix( const int m, const ALL_INDEX_t, const ALL_INDEX_t ) const + { + return Cabana::LinearAlgebra::createMatrixExpression( + [=]( const int i, const int j ) { return ( *this )( m, i, j ); } ); + } +}; - // Strides. - KOKKOS_INLINE_FUNCTION - int stride_1() const { return 0; } +//---------------------------------------------------------------------------// +// Quaternion expression container. +template +struct QuaternionExpression +{ + static constexpr int extent_0 = 4; + static constexpr int extent_1 = 1; - KOKKOS_INLINE_FUNCTION - int stride( const int d ) const { return ( 0 == d ) ? _stride : 0; } + using value_type = T; + using non_const_value_type = typename std::remove_cv::type; - // Extent - KOKKOS_INLINE_FUNCTION - constexpr int extent( const int d ) const - { - return d == 0 ? extent_0 : ( d == 1 ? 1 : 0 ); - } + using eval_type = Quaternion; + using copy_type = Quaternion; - // Access an individual element. - KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i ) const { return _d[_stride * i]; } + Func _f; - KOKKOS_INLINE_FUNCTION - reference operator()( const int i ) { return _d[_stride * i]; } + // Default constructor. + KOKKOS_DEFAULTED_FUNCTION + QuaternionExpression() = default; - // Access an individual element. 2D version for vectors treated as matrices. + // Create an expression from a callable object. KOKKOS_INLINE_FUNCTION - const_reference operator()( const int i, const int ) const + QuaternionExpression( const Func& f ) + : _f( f ) { - return _d[_stride * i]; } + // Extent KOKKOS_INLINE_FUNCTION - reference operator()( const int i, const int ) { return _d[_stride * i]; } + constexpr int extent( const int ) const { return extent_0; } - // Get the raw data. + // Evaluate the expression at an index. KOKKOS_INLINE_FUNCTION - pointer data() const { return const_cast( _d ); } + value_type operator()( const int i ) const { return _f( i ); } + + // Evaluate the expression at an index. 2D version for vectors treated as + // matrices. + KOKKOS_INLINE_FUNCTION + value_type operator()( const int i, int ) const { return _f( i ); } }; //---------------------------------------------------------------------------// @@ -2184,43 +793,52 @@ struct Tensor3 // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const ALL_INDEX_t, const int n, const int p ) const + Cabana::LinearAlgebra::VectorView + vector( const ALL_INDEX_t, const int n, const int p ) const { - return VectorView( const_cast( &_d[0][n][p] ), N * P ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[0][n][p] ), N * P ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const ALL_INDEX_t, const int p ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const ALL_INDEX_t, const int p ) const { - return VectorView( const_cast( &_d[m][0][p] ), P ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[m][0][p] ), P ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const ALL_INDEX_t ) const + Cabana::LinearAlgebra::VectorView vector( const int m, const int n, + const ALL_INDEX_t ) const { - return VectorView( const_cast( &_d[m][n][0] ), 1 ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[m][n][0] ), 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const ALL_INDEX_t, - const int p ) const + Cabana::LinearAlgebra::MatrixView + matrix( const ALL_INDEX_t, const ALL_INDEX_t, const int p ) const { - return MatrixView( const_cast( &_d[0][0][p] ), N * P, P ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[0][0][p] ), N * P, P ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView + matrix( const ALL_INDEX_t, const int n, const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[0][n][0] ), N * P, 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[0][n][0] ), N * P, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView + matrix( const int m, const ALL_INDEX_t, const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[m][0][0] ), P, 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[m][0][0] ), P, 1 ); } }; @@ -2423,49 +1041,52 @@ struct Tensor3View // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const ALL_INDEX_t, const int n, const int p ) const + Cabana::LinearAlgebra::VectorView + vector( const ALL_INDEX_t, const int n, const int p ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[1] * n + _stride[2] * p] ), N * P ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const ALL_INDEX_t, const int p ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const ALL_INDEX_t, const int p ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[0] * m + _stride[2] * p] ), P ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const ALL_INDEX_t ) const + Cabana::LinearAlgebra::VectorView vector( const int m, const int n, + const ALL_INDEX_t ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[0] * m + _stride[1] * n] ), 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const ALL_INDEX_t, - const int p ) const + Cabana::LinearAlgebra::MatrixView + matrix( const ALL_INDEX_t, const ALL_INDEX_t, const int p ) const { - return MatrixView( const_cast( &_d[_stride[2] * p] ), - N * P, P ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[_stride[2] * p] ), N * P, P ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView + matrix( const ALL_INDEX_t, const int n, const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[_stride[1] * n] ), - N * P, 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[_stride[1] * n] ), N * P, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView + matrix( const int m, const ALL_INDEX_t, const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[_stride[0] * m] ), P, - 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[_stride[0] * m] ), P, 1 ); } // Get the raw data. @@ -2756,77 +1377,92 @@ struct Tensor4 // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const ALL_INDEX_t, const int n, const int p, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const ALL_INDEX_t, const int n, const int p, const int q ) const { - return VectorView( const_cast( &_d[0][n][p][q] ), N * P * Q ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[0][n][p][q] ), N * P * Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const ALL_INDEX_t, const int p, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const ALL_INDEX_t, const int p, const int q ) const { - return VectorView( const_cast( &_d[m][0][p][q] ), P * Q ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[m][0][p][q] ), P * Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const ALL_INDEX_t, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const int n, const ALL_INDEX_t, const int q ) const { - return VectorView( const_cast( &_d[m][n][0][q] ), Q ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[m][n][0][q] ), Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const int n, const int p, const ALL_INDEX_t ) const { - return VectorView( const_cast( &_d[m][n][p][0] ), 1 ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[m][n][p][0] ), 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const ALL_INDEX_t, - const int p, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const ALL_INDEX_t, + const int p, + const int q ) const { - return MatrixView( const_cast( &_d[0][0][p][q] ), - N * P * Q, P * Q ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[0][0][p][q] ), N * P * Q, P * Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, - const ALL_INDEX_t, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const int n, + const ALL_INDEX_t, + const int q ) const { - return MatrixView( const_cast( &_d[0][n][0][q] ), - N * P * Q, Q ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[0][n][0][q] ), N * P * Q, Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, - const ALL_INDEX_t, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, + const ALL_INDEX_t, + const ALL_INDEX_t, + const int q ) const { - return MatrixView( const_cast( &_d[m][0][0][q] ), P * Q, - Q ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[m][0][0][q] ), P * Q, Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const int n, const int p, + const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[0][n][p][0] ), - N * P * Q, 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[0][n][p][0] ), N * P * Q, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, + const ALL_INDEX_t, + const int p, + const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[m][0][p][0] ), P * Q, - 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[m][0][p][0] ), P * Q, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const int n, const ALL_INDEX_t, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, const int n, + const ALL_INDEX_t, + const ALL_INDEX_t ) const { - return MatrixView( const_cast( &_d[m][n][0][0] ), Q, 1 ); + return Cabana::LinearAlgebra::MatrixView( + const_cast( &_d[m][n][0][0] ), Q, 1 ); } // Get a tensor3 as a Tensor3 view. @@ -3096,40 +1732,40 @@ struct Tensor4View // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const ALL_INDEX_t, const int n, const int p, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const ALL_INDEX_t, const int n, const int p, const int q ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[1] * n + _stride[2] * p + _stride[3] * q] ), N * P * Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const ALL_INDEX_t, const int p, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const ALL_INDEX_t, const int p, const int q ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[0] * m + _stride[2] * p + _stride[3] * q] ), P * Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const ALL_INDEX_t, - const int q ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const int n, const ALL_INDEX_t, const int q ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[0] * m + _stride[1] * n + _stride[3] * q] ), Q ); } // Get a row as a vector view. KOKKOS_INLINE_FUNCTION - VectorView vector( const int m, const int n, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::VectorView + vector( const int m, const int n, const int p, const ALL_INDEX_t ) const { - return VectorView( + return Cabana::LinearAlgebra::VectorView( const_cast( &_d[_stride[0] * m + _stride[1] * n + _stride[2] * p] ), 1 ); @@ -3137,53 +1773,63 @@ struct Tensor4View // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const ALL_INDEX_t, - const int p, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const ALL_INDEX_t, + const int p, + const int q ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[2] * p + _stride[3] * q] ), N * P * Q, P * Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, - const ALL_INDEX_t, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const int n, + const ALL_INDEX_t, + const int q ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[1] * n + _stride[3] * q] ), N * P * Q, Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, - const ALL_INDEX_t, const int q ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, + const ALL_INDEX_t, + const ALL_INDEX_t, + const int q ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[0] * m + _stride[3] * q] ), P * Q, Q ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const ALL_INDEX_t, const int n, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const ALL_INDEX_t, + const int n, const int p, + const ALL_INDEX_t ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[1] * n + _stride[2] * p] ), N * P * Q, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const ALL_INDEX_t, const int p, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, + const ALL_INDEX_t, + const int p, + const ALL_INDEX_t ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[0] * m + _stride[2] * p] ), P * Q, 1 ); } // Get a matrix as a matrix view. KOKKOS_INLINE_FUNCTION - MatrixView matrix( const int m, const int n, const ALL_INDEX_t, - const ALL_INDEX_t ) const + Cabana::LinearAlgebra::MatrixView matrix( const int m, const int n, + const ALL_INDEX_t, + const ALL_INDEX_t ) const { - return MatrixView( + return Cabana::LinearAlgebra::MatrixView( const_cast( &_d[_stride[0] * m + _stride[1] * n] ), Q, 1 ); } @@ -3291,7 +1937,8 @@ struct Quaternion // Scalar + vector constructor. KOKKOS_INLINE_FUNCTION - Quaternion( const T value, const VectorView vec ) + Quaternion( const T value, + const Cabana::LinearAlgebra::VectorView vec ) { ( *this )( 0 ) = value; #if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) @@ -3425,16 +2072,18 @@ struct Quaternion // Get the vector part KOKKOS_INLINE_FUNCTION - VectorView vector() const + Cabana::LinearAlgebra::VectorView vector() const { - return VectorView( const_cast( &_d[1] ), 1 ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[1] ), 1 ); } // Get the vector part KOKKOS_INLINE_FUNCTION - VectorView vector() + Cabana::LinearAlgebra::VectorView vector() { - return VectorView( const_cast( &_d[1] ), 1 ); + return Cabana::LinearAlgebra::VectorView( + const_cast( &_d[1] ), 1 ); } }; @@ -3659,10 +2308,11 @@ KOKKOS_INLINE_FUNCTION void deepCopy( Tensor3A& a, const ExpressionB& b ) //---------------------------------------------------------------------------// // Matrix-matrix deep copy. //---------------------------------------------------------------------------// -template < - class MatrixA, class ExpressionB, - typename std::enable_if_t< - is_matrix::value && is_matrix::value, int> = 0> +template ::value && + Cabana::LinearAlgebra::is_matrix::value, + int> = 0> KOKKOS_INLINE_FUNCTION void deepCopy( MatrixA& a, const ExpressionB& b ) { static_assert( std::is_same::value && is_vector::value, int> = 0> +template ::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION void deepCopy( VectorX& x, const ExpressionY& y ) { static_assert( std::is_same::value, int> = 0> -KOKKOS_INLINE_FUNCTION auto operator~( const Expression& e ) -{ - return createMatrixExpression( - [=]( const int i, const int j ) { return e( j, i ); } ); -} - -//---------------------------------------------------------------------------// -// Vector operator. -template ::value, int> = 0> -KOKKOS_INLINE_FUNCTION auto operator~( const Expression& e ) -{ - return createMatrixExpression( - [=]( const int, const int j ) { return e( j ); } ); -} - //---------------------------------------------------------------------------// // Quaternion conjugate template ::value && - is_matrix::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator+( const ExpressionA& a, const ExpressionB& b ) { @@ -3775,18 +2405,20 @@ KOKKOS_INLINE_FUNCTION auto operator+( const ExpressionA& a, "extent_0 must match" ); static_assert( ExpressionA::extent_1 == ExpressionB::extent_1, "extent_1 must match" ); - return createMatrixExpression( - [=]( const int i, const int j ) { return a( i, j ) + b( i, j ); } ); + return Cabana::LinearAlgebra::createMatrixExpression< + typename ExpressionA::value_type, ExpressionA::extent_0, + ExpressionA::extent_1>( [=]( const int i, const int j ) + { return a( i, j ) + b( i, j ); } ); } //---------------------------------------------------------------------------// // Matrix-matrix subtraction. //---------------------------------------------------------------------------// template ::value && - is_matrix::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator-( const ExpressionA& a, const ExpressionB& b ) { @@ -3797,9 +2429,10 @@ KOKKOS_INLINE_FUNCTION auto operator-( const ExpressionA& a, "extent_0 must match" ); static_assert( ExpressionA::extent_1 == ExpressionB::extent_1, "extent_1 must_match" ); - return createMatrixExpression( - [=]( const int i, const int j ) { return a( i, j ) - b( i, j ); } ); + return Cabana::LinearAlgebra::createMatrixExpression< + typename ExpressionA::value_type, ExpressionA::extent_0, + ExpressionA::extent_1>( [=]( const int i, const int j ) + { return a( i, j ) - b( i, j ); } ); } //---------------------------------------------------------------------------// @@ -3807,9 +2440,10 @@ KOKKOS_INLINE_FUNCTION auto operator-( const ExpressionA& a, //---------------------------------------------------------------------------// template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && is_matrix::value, - Matrix> + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value, + Cabana::LinearAlgebra::Matrix> operator*( const ExpressionA& a, const ExpressionB& b ) { static_assert( std::is_same + Cabana::LinearAlgebra::Matrix c = static_cast( 0 ); for ( int i = 0; i < ExpressionA::extent_0; ++i ) @@ -3840,8 +2474,10 @@ operator*( const ExpressionA& a, const ExpressionB& b ) //---------------------------------------------------------------------------// template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && is_vector::value, - Vector> + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_vector::value, + Cabana::LinearAlgebra::Vector> operator*( const ExpressionA& a, const ExpressionX& x ) { static_assert( std::is_same y = - static_cast( 0 ); + Cabana::LinearAlgebra::Vector + y = static_cast( 0 ); for ( int i = 0; i < ExpressionA::extent_0; ++i ) #if defined( KOKKOS_ENABLE_PRAGMA_UNROLL ) @@ -3870,9 +2507,10 @@ operator*( const ExpressionA& a, const ExpressionX& x ) //---------------------------------------------------------------------------// template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && is_vector::value, - Matrix> + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_vector::value, + Cabana::LinearAlgebra::Matrix> operator*( const ExpressionX& x, const ExpressionA& a ) { static_assert( std::is_same + Cabana::LinearAlgebra::Matrix y; for ( int i = 0; i < ExpressionX::extent_0; ++i ) @@ -3900,9 +2538,10 @@ operator*( const ExpressionX& x, const ExpressionA& a ) // Vector-vector addition. //---------------------------------------------------------------------------// template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator+( const ExpressionX& x, const ExpressionY& y ) { @@ -3911,8 +2550,8 @@ KOKKOS_INLINE_FUNCTION auto operator+( const ExpressionX& x, "value_type must match" ); static_assert( ExpressionX::extent_0 == ExpressionY::extent_0, "extent must match" ); - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression< + typename ExpressionX::value_type, ExpressionX::extent_0>( [=]( const int i ) { return x( i ) + y( i ); } ); } @@ -3920,9 +2559,10 @@ KOKKOS_INLINE_FUNCTION auto operator+( const ExpressionX& x, // Vector-vector subtraction. //---------------------------------------------------------------------------// template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator-( const ExpressionX& x, const ExpressionY& y ) { @@ -3931,8 +2571,8 @@ KOKKOS_INLINE_FUNCTION auto operator-( const ExpressionX& x, "value_type must match" ); static_assert( ExpressionX::extent_0 == ExpressionY::extent_0, "extent must match" ); - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression< + typename ExpressionX::value_type, ExpressionX::extent_0>( [=]( const int i ) { return x( i ) - y( i ); } ); } @@ -4011,9 +2651,10 @@ KOKKOS_INLINE_FUNCTION auto operator&( const ExpressionX& x, // Quaternion-matrix conjugation //---------------------------------------------------------------------------// template ::value && - is_quaternion::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + is_quaternion::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator&( const ExpressionX& X, const ExpressionY& q ) { @@ -4099,11 +2740,11 @@ KOKKOS_INLINE_FUNCTION auto operator|( const ExpressionX& x, //---------------------------------------------------------------------------// // Cross product template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - is_vector::value, - Vector> - operator%( const ExpressionX& x, const ExpressionY& y ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value && + Cabana::LinearAlgebra::is_vector::value, + Cabana::LinearAlgebra::Vector> +operator%( const ExpressionX& x, const ExpressionY& y ) { static_assert( std::is_same::value, @@ -4114,7 +2755,7 @@ KOKKOS_INLINE_FUNCTION "cross product is for 3-vectors" ); typename ExpressionX::eval_type x_eval = x; typename ExpressionY::eval_type y_eval = y; - return Vector{ + return Cabana::LinearAlgebra::Vector{ x_eval( 1 ) * y_eval( 2 ) - x_eval( 2 ) * y_eval( 1 ), x_eval( 2 ) * y_eval( 0 ) - x_eval( 0 ) * y_eval( 2 ), x_eval( 0 ) * y_eval( 1 ) - x_eval( 1 ) * y_eval( 0 ) }; @@ -4123,9 +2764,10 @@ KOKKOS_INLINE_FUNCTION //---------------------------------------------------------------------------// // Element-wise multiplication. template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator&( const ExpressionX& x, const ExpressionY& y ) { @@ -4134,17 +2776,18 @@ KOKKOS_INLINE_FUNCTION auto operator&( const ExpressionX& x, "value_type must match" ); static_assert( ExpressionX::extent_0 == ExpressionY::extent_0, "extent_0 must match" ); - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression< + typename ExpressionX::value_type, ExpressionX::extent_0>( [=]( const int i ) { return x( i ) * y( i ); } ); } //---------------------------------------------------------------------------// // Element-wise division. template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto operator|( const ExpressionX& x, const ExpressionY& y ) { @@ -4153,8 +2796,8 @@ KOKKOS_INLINE_FUNCTION auto operator|( const ExpressionX& x, "value_type must match" ); static_assert( ExpressionX::extent_0 == ExpressionY::extent_0, "extent must match" ); - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression< + typename ExpressionX::value_type, ExpressionX::extent_0>( [=]( const int i ) { return x( i ) / y( i ); } ); } @@ -4205,17 +2848,20 @@ operator*( const ExpressionA& a, const typename ExpressionA::value_type& s ) // Matrix. template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator*( const typename ExpressionA::value_type& s, const ExpressionA& a ) { - return createMatrixExpression( - [=]( const int i, const int j ) { return s * a( i, j ); } ); + return Cabana::LinearAlgebra::createMatrixExpression< + typename ExpressionA::value_type, ExpressionA::extent_0, + ExpressionA::extent_1>( [=]( const int i, const int j ) + { return s * a( i, j ); } ); } template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator*( const ExpressionA& a, const typename ExpressionA::value_type& s ) { @@ -4225,17 +2871,19 @@ operator*( const ExpressionA& a, const typename ExpressionA::value_type& s ) //---------------------------------------------------------------------------// // Vector. template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator*( const typename ExpressionX::value_type& s, const ExpressionX& x ) { - return createVectorExpression( + return Cabana::LinearAlgebra::createVectorExpression< + typename ExpressionX::value_type, ExpressionX::extent_0>( [=]( const int i ) { return s * x( i ); } ); } template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator*( const ExpressionX& x, const typename ExpressionX::value_type& s ) { @@ -4286,7 +2934,8 @@ operator/( const ExpressionA& a, const typename ExpressionA::value_type& s ) // Matrix. template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator/( const ExpressionA& a, const typename ExpressionA::value_type& s ) { @@ -4297,7 +2946,8 @@ operator/( const ExpressionA& a, const typename ExpressionA::value_type& s ) //---------------------------------------------------------------------------// // Vector. template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION auto operator/( const ExpressionX& x, const typename ExpressionX::value_type& s ) { @@ -4321,12 +2971,11 @@ operator/( const ExpressionX& x, const typename ExpressionX::value_type& s ) //---------------------------------------------------------------------------// // 2x2 specialization template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - Expression::extent_0 == 2 && - Expression::extent_1 == 2, - typename Expression::value_type> - operator!( const Expression& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + Expression::extent_0 == 2 && Expression::extent_1 == 2, + typename Expression::value_type> +operator!( const Expression& a ) { return a( 0, 0 ) * a( 1, 1 ) - a( 0, 1 ) * a( 1, 0 ); } @@ -4334,12 +2983,11 @@ KOKKOS_INLINE_FUNCTION //---------------------------------------------------------------------------// // 3x3 specialization template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - Expression::extent_0 == 3 && - Expression::extent_1 == 3, - typename Expression::value_type> - operator!( const Expression& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + Expression::extent_0 == 3 && Expression::extent_1 == 3, + typename Expression::value_type> +operator!( const Expression& a ) { return a( 0, 0 ) * a( 1, 1 ) * a( 2, 2 ) + a( 0, 1 ) * a( 1, 2 ) * a( 2, 0 ) + @@ -4353,10 +3001,10 @@ KOKKOS_INLINE_FUNCTION // LU decomposition. //---------------------------------------------------------------------------// template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value, - typename ExpressionA::copy_type> - LU( const ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, + typename ExpressionA::copy_type> +LU( const ExpressionA& a ) { using value_type = typename ExpressionA::value_type; constexpr int m = ExpressionA::extent_0; @@ -4390,10 +3038,10 @@ KOKKOS_INLINE_FUNCTION // Matrix trace //---------------------------------------------------------------------------// template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value, - typename ExpressionA::value_type> - trace( const ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, + typename ExpressionA::value_type> +trace( const ExpressionA& a ) { static_assert( ExpressionA::extent_1 == ExpressionA::extent_0, "matrix must be square" ); @@ -4411,9 +3059,9 @@ KOKKOS_INLINE_FUNCTION // Identity //---------------------------------------------------------------------------// template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value, void> - identity( ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, void> +identity( ExpressionA& a ) { static_assert( ExpressionA::extent_1 == ExpressionA::extent_0, "matrix must be square" ); @@ -4447,11 +3095,13 @@ KOKKOS_INLINE_FUNCTION // Diagonal matrix. //---------------------------------------------------------------------------// template ::value, int> = 0> + typename std::enable_if_t< + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION auto diagonal( const ExpressionX& x ) { - return createMatrixExpression( + return Cabana::LinearAlgebra::createMatrixExpression< + typename ExpressionX::value_type, ExpressionX::extent_0, + ExpressionX::extent_0>( [=]( const int i, const int j ) { return ( i == j ) @@ -4466,9 +3116,10 @@ KOKKOS_INLINE_FUNCTION auto diagonal( const ExpressionX& x ) //---------------------------------------------------------------------------// template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + is_tensor3::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, const ExpressionV& v, std::integral_constant ) @@ -4478,8 +3129,8 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, typename ExpressionT::eval_type t_eval = t; typename ExpressionV::eval_type v_eval = v; - Matrix + Cabana::LinearAlgebra::Matrix res = static_cast( 0 ); for ( int i = 0; i < ExpressionT::extent_1; ++i ) @@ -4494,9 +3145,10 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, } template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + is_tensor3::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, const ExpressionV& v, std::integral_constant ) @@ -4506,8 +3158,8 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, typename ExpressionT::eval_type t_eval = t; typename ExpressionV::eval_type v_eval = v; - Matrix + Cabana::LinearAlgebra::Matrix res = static_cast( 0 ); for ( int i = 0; i < ExpressionT::extent_0; ++i ) @@ -4522,9 +3174,10 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, } template ::value && - is_vector::value, - int> = 0> + typename std::enable_if_t< + is_tensor3::value && + Cabana::LinearAlgebra::is_vector::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, const ExpressionV& v, std::integral_constant ) @@ -4534,8 +3187,8 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, typename ExpressionT::eval_type t_eval = t; typename ExpressionV::eval_type v_eval = v; - Matrix + Cabana::LinearAlgebra::Matrix res = static_cast( 0 ); for ( int i = 0; i < ExpressionT::extent_0; ++i ) @@ -4550,9 +3203,10 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, } template ::value && - is_matrix::value, - int> = 0> + typename std::enable_if_t< + is_tensor4::value && + Cabana::LinearAlgebra::is_matrix::value, + int> = 0> KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, const ExpressionM& m ) { @@ -4563,8 +3217,8 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, typename ExpressionT::eval_type t_eval = t; typename ExpressionM::eval_type m_eval = m; - Matrix + Cabana::LinearAlgebra::Matrix res = static_cast( 0 ); for ( int i = 0; i < ExpressionT::extent_0; ++i ) @@ -4584,19 +3238,17 @@ KOKKOS_INLINE_FUNCTION auto contract( const ExpressionT& t, //---------------------------------------------------------------------------// // 2x2 specialization with determinant given. template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - ExpressionA::extent_0 == 2 && - ExpressionA::extent_1 == 2, - typename ExpressionA::copy_type> - inverse( const ExpressionA& a, - const typename ExpressionA::value_type& a_det ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + ExpressionA::extent_0 == 2 && ExpressionA::extent_1 == 2, + typename ExpressionA::copy_type> +inverse( const ExpressionA& a, const typename ExpressionA::value_type& a_det ) { typename ExpressionA::eval_type a_eval = a; auto a_det_inv = static_cast( 1 ) / a_det; - Matrix a_inv; + Cabana::LinearAlgebra::Matrix a_inv; a_inv( 0, 0 ) = a_eval( 1, 1 ) * a_det_inv; a_inv( 0, 1 ) = -a_eval( 0, 1 ) * a_det_inv; @@ -4609,12 +3261,11 @@ KOKKOS_INLINE_FUNCTION //---------------------------------------------------------------------------// // 2x2 specialization. template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - ExpressionA::extent_0 == 2 && - ExpressionA::extent_1 == 2, - typename ExpressionA::copy_type> - inverse( const ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + ExpressionA::extent_0 == 2 && ExpressionA::extent_1 == 2, + typename ExpressionA::copy_type> +inverse( const ExpressionA& a ) { typename ExpressionA::eval_type a_eval = a; return inverse( a_eval, !a_eval ); @@ -4623,19 +3274,17 @@ KOKKOS_INLINE_FUNCTION //---------------------------------------------------------------------------// // 3x3 specialization with determinant given. template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - ExpressionA::extent_0 == 3 && - ExpressionA::extent_1 == 3, - typename ExpressionA::copy_type> - inverse( const ExpressionA& a, - const typename ExpressionA::value_type& a_det ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3, + typename ExpressionA::copy_type> +inverse( const ExpressionA& a, const typename ExpressionA::value_type& a_det ) { typename ExpressionA::eval_type a_eval = a; auto a_det_inv = static_cast( 1 ) / a_det; - Matrix a_inv; + Cabana::LinearAlgebra::Matrix a_inv; a_inv( 0, 0 ) = ( a_eval( 1, 1 ) * a_eval( 2, 2 ) - a_eval( 1, 2 ) * a_eval( 2, 1 ) ) * @@ -4673,12 +3322,11 @@ KOKKOS_INLINE_FUNCTION //---------------------------------------------------------------------------// // 3x3 specialization. template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - ExpressionA::extent_0 == 3 && - ExpressionA::extent_1 == 3, - typename ExpressionA::copy_type> - inverse( const ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3, + typename ExpressionA::copy_type> +inverse( const ExpressionA& a ) { typename ExpressionA::eval_type a_eval = a; return inverse( a_eval, !a_eval ); @@ -4688,14 +3336,14 @@ KOKKOS_INLINE_FUNCTION // General case. template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value && !( ExpressionA::extent_0 == 2 && ExpressionA::extent_1 == 2 ) && !( ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3 ), typename ExpressionA::copy_type> inverse( const ExpressionA& a ) { - Matrix + Cabana::LinearAlgebra::Matrix ident; identity( ident ); return a ^ ident; @@ -4707,15 +3355,15 @@ inverse( const ExpressionA& a ) // https://arxiv.org/abs/1606.08395 //---------------------------------------------------------------------------// template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value, - typename ExpressionA::copy_type> - exponential( const ExpressionA& a ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value, + typename ExpressionA::copy_type> +exponential( const ExpressionA& a ) { static_assert( ExpressionA::extent_1 == ExpressionA::extent_0, "matrix must be square" ); - Matrix + Cabana::LinearAlgebra::Matrix ident; identity( ident ); @@ -4758,8 +3406,9 @@ KOKKOS_INLINE_FUNCTION // 2x2 specialization. Single and multiple RHS supported. template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && - ( is_matrix::value || is_vector::value ) && + Cabana::LinearAlgebra::is_matrix::value && + ( Cabana::LinearAlgebra::is_matrix::value || + Cabana::LinearAlgebra::is_vector::value ) && ExpressionA::extent_0 == 2 && ExpressionA::extent_1 == 2, typename ExpressionB::copy_type> operator^( const ExpressionA& a, const ExpressionB& b ) @@ -4771,8 +3420,9 @@ operator^( const ExpressionA& a, const ExpressionB& b ) // 3x3 specialization. Single and multiple RHS supported. template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && - ( is_matrix::value || is_vector::value ) && + Cabana::LinearAlgebra::is_matrix::value && + ( Cabana::LinearAlgebra::is_matrix::value || + Cabana::LinearAlgebra::is_vector::value ) && ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3, typename ExpressionB::copy_type> operator^( const ExpressionA& a, const ExpressionB& b ) @@ -4784,8 +3434,9 @@ operator^( const ExpressionA& a, const ExpressionB& b ) // General case. Single and multiple RHS supported. template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && - ( is_matrix::value || is_vector::value ) && + Cabana::LinearAlgebra::is_matrix::value && + ( Cabana::LinearAlgebra::is_matrix::value || + Cabana::LinearAlgebra::is_vector::value ) && !( ExpressionA::extent_0 == 2 && ExpressionA::extent_1 == 2 ) && !( ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3 ), typename ExpressionB::copy_type> @@ -4869,10 +3520,12 @@ KOKKOS_INLINE_FUNCTION void condNegSwap( bool b, ValueType& lv, ValueType& rv ) template <> KOKKOS_INLINE_FUNCTION void -condNegSwap>( bool b, VectorView& lv, - VectorView& rv ) +condNegSwap>( + bool b, Cabana::LinearAlgebra::VectorView& lv, + Cabana::LinearAlgebra::VectorView& rv ) { - using vector_type = typename VectorView::copy_type; + using vector_type = + typename Cabana::LinearAlgebra::VectorView::copy_type; vector_type tmpv; if ( b ) @@ -4962,12 +3615,11 @@ Quaternion givensQR( double a1, double a2, double tol, } template -KOKKOS_INLINE_FUNCTION - typename std::enable_if_t::value && - ExpressionA::extent_0 == 3 && - ExpressionA::extent_1 == 3, - Quaternion> - cyclicJacobi( const ExpressionA& A, const int max_iter ) +KOKKOS_INLINE_FUNCTION typename std::enable_if_t< + Cabana::LinearAlgebra::is_matrix::value && + ExpressionA::extent_0 == 3 && ExpressionA::extent_1 == 3, + Quaternion> +cyclicJacobi( const ExpressionA& A, const int max_iter ) { // Form the symmetric positive-definite matrix from A auto S = ~A * A; @@ -4989,7 +3641,7 @@ KOKKOS_INLINE_FUNCTION Quaternion q_total = givensQuaternion( s_pp, s_pq, s_qq, pq ); // Convert to rotation matrix for conjugation - Matrix Q_1{ q_total }; + Cabana::LinearAlgebra::Matrix Q_1{ q_total }; S = ~Q_1 * S * Q_1; @@ -5007,7 +3659,7 @@ KOKKOS_INLINE_FUNCTION auto q = givensQuaternion( s_pp, s_pq, s_qq, pq ); // Convert to rotation matrix for conjugation - Matrix Q{ q }; + Cabana::LinearAlgebra::Matrix Q{ q }; S = ~Q * S * Q; // Update to (k+1) q_total = q_total & q; @@ -5063,8 +3715,10 @@ KOKKOS_INLINE_FUNCTION void sortSingularValues( MatrixType& B, MatrixType& V ) template KOKKOS_INLINE_FUNCTION typename std::enable_if_t< - is_matrix::value && is_matrix::value && - is_matrix::value && is_matrix::value, + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value && + Cabana::LinearAlgebra::is_matrix::value, void> svd( const ExpressionA& A, EigenU& U, Diagonal& D, EigenV& V ) { @@ -5078,9 +3732,9 @@ svd( const ExpressionA& A, EigenU& U, Diagonal& D, EigenV& V ) double q_mag = Kokkos::sqrt( q( 0 ) * q( 0 ) + q( 1 ) * q( 1 ) + q( 2 ) * q( 2 ) + q( 3 ) * q( 3 ) ); - Matrix V_rot{ q / q_mag }; + Cabana::LinearAlgebra::Matrix V_rot{ q / q_mag }; - Matrix Q = 0.0; + Cabana::LinearAlgebra::Matrix Q = 0.0; // Perform a QR factorization of the matrix B = AV auto B = A * V_rot; @@ -5088,21 +3742,22 @@ svd( const ExpressionA& A, EigenU& U, Diagonal& D, EigenV& V ) sortSingularValues( B, V_rot ); auto q21 = givensQR( B( 0, 0 ), B( 1, 0 ), tol, { 1, 0 } ); - auto Q1 = static_cast>( q21 ); + auto Q1 = static_cast>( q21 ); auto B1 = ~Q1 * B; auto q31 = givensQR( B1( 0, 0 ), B1( 2, 0 ), tol, { 2, 0 } ); - auto Q2 = static_cast>( q31 ); + auto Q2 = static_cast>( q31 ); auto B2 = ~Q2 * B1; auto q32 = givensQR( B2( 1, 1 ), B2( 2, 1 ), tol, { 2, 1 } ); - auto Q3 = static_cast>( q32 ); + auto Q3 = static_cast>( q32 ); auto B3 = ~Q3 * B2; - Q = static_cast>( q21 & q31 & q32 ); + Q = static_cast>( q21 & q31 & + q32 ); U = Q; D = B3; @@ -5114,9 +3769,11 @@ svd( const ExpressionA& A, EigenU& U, Diagonal& D, EigenV& V ) //---------------------------------------------------------------------------// // template // KOKKOS_INLINE_FUNCTION -// typename std::enable_if_t::value && -// is_vector::value && -// is_matrix::value, +// typename +// std::enable_if_t::value && +// Cabana::LinearAlgebra::Cabana::LinearAlgebra::is_vector::value +// && +// Cabana::LinearAlgebra::is_matrix::value, // typename ExpressionA::copy_type> // eigendecomposition( const ExpressionA& a, // Eigenvalues& e_real, @@ -5149,22 +3806,22 @@ svd( const ExpressionA& A, EigenU& U, Diagonal& D, EigenV& V ) //---------------------------------------------------------------------------// template -using Mat2 = LinearAlgebra::Matrix; +using Mat2 = Cabana::LinearAlgebra::Matrix; template -using MatView2 = LinearAlgebra::MatrixView; +using MatView2 = Cabana::LinearAlgebra::MatrixView; template -using Vec2 = LinearAlgebra::Vector; +using Vec2 = Cabana::LinearAlgebra::Vector; template -using VecView2 = LinearAlgebra::VectorView; +using VecView2 = Cabana::LinearAlgebra::VectorView; template -using Mat3 = LinearAlgebra::Matrix; +using Mat3 = Cabana::LinearAlgebra::Matrix; template -using MatView3 = LinearAlgebra::MatrixView; +using MatView3 = Cabana::LinearAlgebra::MatrixView; template -using Vec3 = LinearAlgebra::Vector; +using Vec3 = Cabana::LinearAlgebra::Vector; template -using VecView3 = LinearAlgebra::VectorView; +using VecView3 = Cabana::LinearAlgebra::VectorView; //---------------------------------------------------------------------------// diff --git a/src/Picasso_BilinearMeshMapping.hpp b/src/Picasso_BilinearMeshMapping.hpp index e6bbef0..db5ff9e 100644 --- a/src/Picasso_BilinearMeshMapping.hpp +++ b/src/Picasso_BilinearMeshMapping.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_BILINEARMESHMAPPING_HPP #define PICASSO_BILINEARMESHMAPPING_HPP -#include #include #include #include @@ -137,11 +136,11 @@ class CurvilinearMeshMapping> template static KOKKOS_INLINE_FUNCTION void transformationMetrics( const mesh_mapping& mapping, const ReferenceCoords& reference_coords, - LinearAlgebra::Matrix& - jacobian, + Cabana::LinearAlgebra::Matrix& jacobian, typename ReferenceCoords::value_type& jacobian_det, - LinearAlgebra::Matrix& - jacobian_inv ) + Cabana::LinearAlgebra::Matrix& jacobian_inv ) { using value_type = typename ReferenceCoords::value_type; @@ -186,7 +185,7 @@ class CurvilinearMeshMapping> jacobian_det = !jacobian; - jacobian_inv = LinearAlgebra::inverse( jacobian, jacobian_det ); + jacobian_inv = Cabana::LinearAlgebra::inverse( jacobian, jacobian_det ); } // Reverse mapping. Given coordinates in the physical frame compute the @@ -262,11 +261,11 @@ struct CurvilinearMeshMapping> template static KOKKOS_INLINE_FUNCTION void transformationMetrics( const mesh_mapping& mapping, const ReferenceCoords& reference_coords, - LinearAlgebra::Matrix& - jacobian, + Cabana::LinearAlgebra::Matrix& jacobian, typename ReferenceCoords::value_type& jacobian_det, - LinearAlgebra::Matrix& - jacobian_inv ) + Cabana::LinearAlgebra::Matrix& jacobian_inv ) { using value_type = typename ReferenceCoords::value_type; @@ -299,7 +298,7 @@ struct CurvilinearMeshMapping> jacobian_det = !jacobian; - jacobian_inv = LinearAlgebra::inverse( jacobian, jacobian_det ); + jacobian_inv = Cabana::LinearAlgebra::inverse( jacobian, jacobian_det ); } // Reverse mapping. Given coordinates in the physical frame compute the diff --git a/src/Picasso_CurvilinearMesh.hpp b/src/Picasso_CurvilinearMesh.hpp index ed335d4..863f3e7 100644 --- a/src/Picasso_CurvilinearMesh.hpp +++ b/src/Picasso_CurvilinearMesh.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_CURVILINEARMESH_HPP #define PICASSO_CURVILINEARMESH_HPP -#include #include #include @@ -76,11 +75,12 @@ struct CurvilinearMeshMapping template static KOKKOS_INLINE_FUNCTION void transformationMetrics( const Mapping& mapping, const ReferenceCoords& reference_coords, - LinearAlgebra::Matrix& jacobian, + Cabana::LinearAlgebra::Matrix& jacobian, typename ReferenceCoords::value_type& jacobian_det, - LinearAlgebra::Matrix& jacobian_inv ); + Cabana::LinearAlgebra::Matrix& + jacobian_inv ); // Reverse mapping. Given coordinates in the physical frame compute the // coordinates in the reference frame. The data in reference_coords @@ -123,12 +123,12 @@ struct DefaultCurvilinearMeshMapping int max_iter = 15; // Iteration data. - LinearAlgebra::Vector x_ref_old; - LinearAlgebra::Vector x_phys_new; - LinearAlgebra::Matrix + Cabana::LinearAlgebra::Vector x_ref_old; + Cabana::LinearAlgebra::Vector x_phys_new; + Cabana::LinearAlgebra::Matrix jacobian; value_type jacobian_det; - LinearAlgebra::Matrix + Cabana::LinearAlgebra::Matrix jacobian_inv; // Newton iterations. diff --git a/src/Picasso_FacetGeometry.hpp b/src/Picasso_FacetGeometry.hpp index 52a1319..c9e2542 100644 --- a/src/Picasso_FacetGeometry.hpp +++ b/src/Picasso_FacetGeometry.hpp @@ -12,11 +12,11 @@ #ifndef PICASSO_FACETGEOMETRY_HPP #define PICASSO_FACETGEOMETRY_HPP -#include - #include #include +#include + #include #include @@ -511,8 +511,8 @@ KOKKOS_FUNCTION bool pointFacetProjection( const float x[3], const float r[3], { // Build the system of equations to solve for intersection. Fire the ray // in the Y direction - this choice is arbitrary. - Mat3 A; - Vec3 b; + Cabana::Mat3 A; + Cabana::Vec3 b; for ( int i = 0; i < 3; ++i ) { A( i, 0 ) = facets( f, 1, i ) - facets( f, 0, i ); @@ -528,7 +528,7 @@ KOKKOS_FUNCTION bool pointFacetProjection( const float x[3], const float r[3], return false; // Solve the system. - VecView3 y_view( y, 1 ); + Cabana::VecView3 y_view( y, 1 ); y_view = A ^ b; // Check the solution for inclusion in the triangle. diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 11bc6c0..8be4779 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -205,7 +205,7 @@ struct Scalar : ScalarBase template using field_type = Scalar; template - using gradient_type = LinearAlgebra::Vector; + using gradient_type = Cabana::LinearAlgebra::Vector; }; template @@ -232,11 +232,11 @@ struct Vector : VectorBase static constexpr int size = D; static constexpr int dim0 = D; using data_type = value_type[D]; - using linear_algebra_type = LinearAlgebra::VectorView; + using linear_algebra_type = Cabana::LinearAlgebra::VectorView; template using field_type = Vector; template - using gradient_type = LinearAlgebra::Matrix; + using gradient_type = Cabana::LinearAlgebra::Matrix; }; template @@ -247,11 +247,11 @@ struct Quaternion : VectorBase static constexpr int size = 4; static constexpr int dim0 = 4; using data_type = value_type[4]; - using linear_algebra_type = LinearAlgebra::QuaternionView; + using linear_algebra_type = Picasso::LinearAlgebra::QuaternionView; template using field_type = Quaternion; template - using matrix_type = Mat3; + using matrix_type = Cabana::Mat3; }; template @@ -279,7 +279,7 @@ struct Matrix : MatrixBase static constexpr int dim0 = D0; static constexpr int dim1 = D1; using data_type = value_type[D0][D1]; - using linear_algebra_type = LinearAlgebra::MatrixView; + using linear_algebra_type = Cabana::LinearAlgebra::MatrixView; template using field_type = Matrix; }; @@ -310,7 +310,8 @@ struct Tensor3 : Tensor3Base static constexpr int dim1 = D1; static constexpr int dim2 = D2; using data_type = value_type[D0][D1][D2]; - using linear_algebra_type = LinearAlgebra::Tensor3View; + using linear_algebra_type = + Picasso::LinearAlgebra::Tensor3View; template using field_type = Tensor3; }; @@ -342,7 +343,8 @@ struct Tensor4 : Tensor4Base static constexpr int dim2 = D2; static constexpr int dim3 = D3; using data_type = value_type[D0][D1][D2][D3]; - using linear_algebra_type = LinearAlgebra::Tensor4View; + using linear_algebra_type = + Picasso::LinearAlgebra::Tensor4View; template using field_type = Tensor4; }; diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index 7a04382..9032687 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_PARTICLEINIT_HPP #define PICASSO_PARTICLEINIT_HPP -#include #include #include @@ -499,32 +498,32 @@ void initializeParticlesSurface( InitRandom, const ExecutionSpace&, auto rand = pool.get_state( f ); // Particle coordinate. - Vec3 px; + Cabana::Vec3 px; // Particle. particle_type particle; - Vec3 a = { surface.facets( f, 0, 0 ), - surface.facets( f, 0, 1 ), - surface.facets( f, 0, 2 ) }; + Cabana::Vec3 a = { surface.facets( f, 0, 0 ), + surface.facets( f, 0, 1 ), + surface.facets( f, 0, 2 ) }; - Vec3 b = { surface.facets( f, 1, 0 ), - surface.facets( f, 1, 1 ), - surface.facets( f, 1, 2 ) }; + Cabana::Vec3 b = { surface.facets( f, 1, 0 ), + surface.facets( f, 1, 1 ), + surface.facets( f, 1, 2 ) }; - Vec3 c = { surface.facets( f, 2, 0 ), - surface.facets( f, 2, 1 ), - surface.facets( f, 2, 2 ) }; + Cabana::Vec3 c = { surface.facets( f, 2, 0 ), + surface.facets( f, 2, 1 ), + surface.facets( f, 2, 2 ) }; - Vec3 ab = b - a; - Vec3 ac = c - a; + Cabana::Vec3 ab = b - a; + Cabana::Vec3 ac = c - a; auto cross = ab % ac; double c2 = ~cross * cross; double sqc2 = sqrt( c2 ); double facet_area = 0.5 * sqc2; - Vec3 pan = + Cabana::Vec3 pan = cross / sqc2; // Create unit normal vector for the facet area // Particle surface area. diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index fa41a0c..0e301a8 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_PARTICLEINTERPOLATION_HPP #define PICASSO_PARTICLEINTERPOLATION_HPP -#include #include #include @@ -112,8 +111,8 @@ namespace G2P //---------------------------------------------------------------------------// // G2P scalar value. Requires SplineValue when constructing the spline data. template ::value, - int> = 0> + typename std::enable_if_t< + !Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, const ViewType& view, Scalar& result ) { @@ -124,7 +123,7 @@ KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, // G2P vector value. Requires SplineValue when constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, const ViewType& view, ResultVector& result ) { @@ -144,7 +143,7 @@ KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, // constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void gradient( const SplineDataType& sd, const ViewType& view, ResultVector& result ) { @@ -164,7 +163,7 @@ gradient( const SplineDataType& sd, const ViewType& view, ResultVector& result ) // constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_matrix::value, int> = 0> KOKKOS_INLINE_FUNCTION void gradient( const SplineDataType& sd, const ViewType& view, ResultMatrix& result ) { @@ -202,8 +201,8 @@ namespace P2G //---------------------------------------------------------------------------// // P2G scalar value. Requires SplineValue when constructing the spline data. template ::value, - int> = 0> + typename std::enable_if_t< + !Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, const Scalar& value, const ScatterViewType& view ) @@ -215,7 +214,7 @@ KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, // P2G vector value. Requires SplineValue when constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void value( const SplineDataType& sd, const ValueVector& value, const ScatterViewType& view ) @@ -243,7 +242,7 @@ KOKKOS_INLINE_FUNCTION void gradient( const SplineDataType& sd, // constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_vector::value, int> = 0> KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd, const ValueVector& value, const ScatterViewType& view ) @@ -260,7 +259,7 @@ KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd, // constructing the spline data. template ::value, int> = 0> + Cabana::LinearAlgebra::is_matrix::value, int> = 0> KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd, const ValueMatrix& value, const ScatterViewType& view ) diff --git a/src/Picasso_ParticleList.hpp b/src/Picasso_ParticleList.hpp index 9eca0f6..60a656d 100644 --- a/src/Picasso_ParticleList.hpp +++ b/src/Picasso_ParticleList.hpp @@ -88,20 +88,22 @@ get( Cabana::ParticleView& particle, FieldTag, // Get a view of a particle member as a vector. (Works for both Particle // and ParticleView) template -KOKKOS_FORCEINLINE_FUNCTION typename std::enable_if< - LinearAlgebra::is_vector::value, - typename FieldTag::linear_algebra_type>::type -get( ParticleType& particle, FieldTag tag ) +KOKKOS_FORCEINLINE_FUNCTION + typename std::enable_if::value, + typename FieldTag::linear_algebra_type>::type + get( ParticleType& particle, FieldTag tag ) { return typename FieldTag::linear_algebra_type( &( Cabana::get( particle, tag, 0 ) ), ParticleType::vector_length ); } template -KOKKOS_FORCEINLINE_FUNCTION typename std::enable_if< - LinearAlgebra::is_vector::value, - const typename FieldTag::linear_algebra_type>::type -get( const ParticleType& particle, FieldTag tag ) +KOKKOS_FORCEINLINE_FUNCTION + typename std::enable_if::value, + const typename FieldTag::linear_algebra_type>::type + get( const ParticleType& particle, FieldTag tag ) { return typename FieldTag::linear_algebra_type( const_cast( @@ -113,10 +115,11 @@ get( const ParticleType& particle, FieldTag tag ) // Get a view of a particle member as a matrix. (Works for both Particle // and ParticleView) template -KOKKOS_FORCEINLINE_FUNCTION typename std::enable_if< - LinearAlgebra::is_matrix::value, - typename FieldTag::linear_algebra_type>::type -get( ParticleType& particle, FieldTag tag ) +KOKKOS_FORCEINLINE_FUNCTION + typename std::enable_if::value, + typename FieldTag::linear_algebra_type>::type + get( ParticleType& particle, FieldTag tag ) { return typename FieldTag::linear_algebra_type( &( Cabana::get( particle, tag, 0, 0 ) ), @@ -125,10 +128,11 @@ get( ParticleType& particle, FieldTag tag ) } template -KOKKOS_FORCEINLINE_FUNCTION typename std::enable_if< - LinearAlgebra::is_matrix::value, - const typename FieldTag::linear_algebra_type>::type -get( const ParticleType& particle, FieldTag tag ) +KOKKOS_FORCEINLINE_FUNCTION + typename std::enable_if::value, + const typename FieldTag::linear_algebra_type>::type + get( const ParticleType& particle, FieldTag tag ) { return typename FieldTag::linear_algebra_type( const_cast( diff --git a/src/Picasso_PolyPIC.hpp b/src/Picasso_PolyPIC.hpp index b2641e9..b6848bd 100644 --- a/src/Picasso_PolyPIC.hpp +++ b/src/Picasso_PolyPIC.hpp @@ -14,7 +14,6 @@ #include -#include #include #include @@ -78,17 +77,17 @@ KOKKOS_INLINE_FUNCTION void p2g( int ncomp = ParticleField::extent_1; // Affine material motion operator. - Mat3 am_p = { + Cabana::Mat3 am_p = { { 1.0 + dt * u_p( 1, 0 ), dt * u_p( 1, 1 ), dt * u_p( 1, 2 ) }, { dt * u_p( 2, 0 ), 1.0 + dt * u_p( 2, 1 ), dt * u_p( 2, 2 ) }, { dt * u_p( 3, 0 ), dt * u_p( 3, 1 ), 1.0 + dt * u_p( 3, 2 ) } }; // Invert the affine operator. - auto am_inv_p = LinearAlgebra::inverse( am_p ); + auto am_inv_p = Cabana::LinearAlgebra::inverse( am_p ); // Project mass and mass-weighted field. - LinearAlgebra::Vector basis; - Vec3 distance; + Cabana::LinearAlgebra::Vector basis; + Cabana::Vec3 distance; value_type wm; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -181,17 +180,17 @@ KOKKOS_INLINE_FUNCTION void p2g( const int dim = ( 3 == ncomp ) ? SplineDataType::entity_type::dim : 0; // Affine material motion operator. - Mat3 am_p = { + Cabana::Mat3 am_p = { { 1.0 + dt * u_p( 1, 0 ), dt * u_p( 1, 1 ), dt * u_p( 1, 2 ) }, { dt * u_p( 2, 0 ), 1.0 + dt * u_p( 2, 1 ), dt * u_p( 2, 2 ) }, { dt * u_p( 3, 0 ), dt * u_p( 3, 1 ), 1.0 + dt * u_p( 3, 2 ) } }; // Invert the affine operator. - auto am_inv_p = LinearAlgebra::inverse( am_p ); + auto am_inv_p = Cabana::LinearAlgebra::inverse( am_p ); // Project mass and mass-weighted field. - LinearAlgebra::Vector basis; - Vec3 distance; + Cabana::LinearAlgebra::Vector basis; + Cabana::Vec3 distance; value_type wm; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) @@ -260,7 +259,7 @@ KOKKOS_INLINE_FUNCTION void g2p( c_p = 0.0; // Update particle. - LinearAlgebra::Vector coeff; + Cabana::LinearAlgebra::Vector coeff; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) for ( int k = 0; k < SplineDataType::num_knot; ++k ) @@ -328,7 +327,7 @@ KOKKOS_INLINE_FUNCTION void g2p( c_p.column( dim ) = 0.0; // Update particle. - LinearAlgebra::Vector coeff; + Cabana::LinearAlgebra::Vector coeff; for ( int i = 0; i < SplineDataType::num_knot; ++i ) for ( int j = 0; j < SplineDataType::num_knot; ++j ) for ( int k = 0; k < SplineDataType::num_knot; ++k ) diff --git a/src/Picasso_UniformCartesianMeshMapping.hpp b/src/Picasso_UniformCartesianMeshMapping.hpp index 68f8122..82c97ac 100644 --- a/src/Picasso_UniformCartesianMeshMapping.hpp +++ b/src/Picasso_UniformCartesianMeshMapping.hpp @@ -12,7 +12,6 @@ #ifndef PICASSO_UNIFORMCARTESIANMESHMAPPING_HPP #define PICASSO_UNIFORMCARTESIANMESHMAPPING_HPP -#include #include #include #include @@ -91,11 +90,11 @@ class CurvilinearMeshMapping< template static KOKKOS_INLINE_FUNCTION void transformationMetrics( const mesh_mapping& mapping, const ReferenceCoords&, - LinearAlgebra::Matrix& jacobian, + Cabana::LinearAlgebra::Matrix& jacobian, typename ReferenceCoords::value_type& jacobian_det, - LinearAlgebra::Matrix& jacobian_inv ) + Cabana::LinearAlgebra::Matrix& jacobian_inv ) { for ( std::size_t i = 0; i < num_space_dim; ++i ) for ( std::size_t j = 0; j < num_space_dim; ++j ) diff --git a/unit_test/CMakeLists.txt b/unit_test/CMakeLists.txt index 869a95b..8dd3af7 100644 --- a/unit_test/CMakeLists.txt +++ b/unit_test/CMakeLists.txt @@ -142,7 +142,6 @@ macro(Picasso_add_tests) endmacro() Picasso_add_tests(NAMES - BatchedLinearAlgebra PolyPIC APIC ) diff --git a/unit_test/tstAPIC.hpp b/unit_test/tstAPIC.hpp index 3419d19..4d10905 100644 --- a/unit_test/tstAPIC.hpp +++ b/unit_test/tstAPIC.hpp @@ -609,7 +609,7 @@ void collocatedTest() createSpline( Location(), InterpolationOrder(), local_mesh, x, SplineValue(), SplineDistance() ); - LinearAlgebra::Matrix aff = 0.0; + Cabana::LinearAlgebra::Matrix aff = 0.0; APIC::g2p( gv_wrapper, aff, sd ); @@ -651,7 +651,7 @@ void collocatedTest() local_mesh, x, SplineValue(), SplineGradient(), SplineDistance(), SplineCellSize() ); - LinearAlgebra::Matrix aff = 0.0; + Cabana::LinearAlgebra::Matrix aff = 0.0; for ( int d = 0; d < 3; ++d ) aff( 0, d ) = pb( 0, d ); @@ -777,7 +777,7 @@ void staggeredTest() InterpolationOrder(), local_mesh, x, SplineValue(), SplineDistance() ); - LinearAlgebra::Matrix aff = 0.0; + Cabana::LinearAlgebra::Matrix aff = 0.0; APIC::g2p( gs_wrapper, aff, sd ); @@ -815,7 +815,7 @@ void staggeredTest() SplineValue(), SplineGradient(), SplineDistance(), SplineCellSize() ); - LinearAlgebra::Matrix aff = 0.0; + Cabana::LinearAlgebra::Matrix aff = 0.0; for ( int d = 0; d < 3; ++d ) aff( 0, d ) = pb( 0, d ); diff --git a/unit_test/tstBilinearMeshMapping.hpp b/unit_test/tstBilinearMeshMapping.hpp index 6f25791..d975a56 100644 --- a/unit_test/tstBilinearMeshMapping.hpp +++ b/unit_test/tstBilinearMeshMapping.hpp @@ -9,7 +9,6 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ -#include #include #include #include @@ -116,18 +115,18 @@ void mappingTest3d() "check_mapping", TEST_EXECSPACE{}, ghosted_cells, KOKKOS_LAMBDA( const int i, const int j, const int k ) { // Map to physical frame. - LinearAlgebra::Vector cell_coords = 0.0; + Cabana::LinearAlgebra::Vector cell_coords = 0.0; int ijk[3] = { i, j, k }; local_mesh.coordinates( Cabana::Grid::Cell{}, ijk, cell_coords.data() ); - LinearAlgebra::VectorView phys_coords( + Cabana::LinearAlgebra::VectorView phys_coords( &cell_forward_map( i, j, k, 0 ), cell_forward_map.stride( 3 ) ); CurvilinearMeshMapping::mapToPhysicalFrame( mapping, cell_coords, phys_coords ); // Map to reference frame. Use the bad guess to test search. - LinearAlgebra::VectorView ref_coords( + Cabana::LinearAlgebra::VectorView ref_coords( &cell_reverse_map( i, j, k, 0 ), cell_reverse_map.stride( 3 ) ); ref_coords = { guess[Dim::I], guess[Dim::J], guess[Dim::K] }; cell_map_success( i, j, k ) = @@ -135,7 +134,7 @@ void mappingTest3d() mapping, phys_coords, ref_coords ); // Default map to reference frame. Use a good guess in the cell. - LinearAlgebra::VectorView default_ref_coords( + Cabana::LinearAlgebra::VectorView default_ref_coords( &default_cell_reverse_map( i, j, k, 0 ), default_cell_reverse_map.stride( 3 ) ); default_ref_coords = { cell_coords( Dim::I ) - 0.1, @@ -285,17 +284,17 @@ void mappingTest2d() "check_mapping", TEST_EXECSPACE{}, ghosted_cells, KOKKOS_LAMBDA( const int i, const int j ) { // Map to physical frame. - LinearAlgebra::Vector cell_coords = 0.0; + Cabana::LinearAlgebra::Vector cell_coords = 0.0; int ijk[2] = { i, j }; local_mesh.coordinates( Cabana::Grid::Cell{}, ijk, cell_coords.data() ); - LinearAlgebra::VectorView phys_coords( + Cabana::LinearAlgebra::VectorView phys_coords( &cell_forward_map( i, j, 0 ), cell_forward_map.stride( 2 ) ); CurvilinearMeshMapping::mapToPhysicalFrame( mapping, cell_coords, phys_coords ); // Map to reference frame. Use a bad guess to test the search. - LinearAlgebra::VectorView ref_coords( + Cabana::LinearAlgebra::VectorView ref_coords( &cell_reverse_map( i, j, 0 ), cell_reverse_map.stride( 2 ) ); ref_coords = { guess[Dim::I], guess[Dim::J] }; cell_map_success( i, j ) = @@ -304,7 +303,7 @@ void mappingTest2d() // Default map to reference frame. Use a good guess that is in the // cell. - LinearAlgebra::VectorView default_ref_coords( + Cabana::LinearAlgebra::VectorView default_ref_coords( &default_cell_reverse_map( i, j, 0 ), default_cell_reverse_map.stride( 2 ) ); default_ref_coords = { cell_coords( Dim::I ) - 0.1, diff --git a/unit_test/tstGridOperator3d.hpp b/unit_test/tstGridOperator3d.hpp index 77e6426..9878dbc 100644 --- a/unit_test/tstGridOperator3d.hpp +++ b/unit_test/tstGridOperator3d.hpp @@ -119,13 +119,14 @@ struct ParticleFunc auto bar_out = scatter_deps.get( FieldLocation::Cell(), BarOut() ); // Get particle data. - auto foop = get( particle, FooP() ); - auto& barp = Picasso::get( particle, BarP() ); + auto foop = Cabana::get( particle, FooP() ); + auto& barp = Cabana::get( particle, BarP() ); // Zero-order cell interpolant. auto spline = createSpline( FieldLocation::Cell(), InterpolationOrder<0>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue() ); + Cabana::get( particle, Field::LogicalPosition<3>() ), + SplineValue() ); // Interpolate to the particles. G2P::value( spline, foo_in, foop ); @@ -232,7 +233,7 @@ struct GridTensor3Func boo( i, j, k ) = levi_civita; const int index[3] = { i, j, k }; - auto boo_t_foo_in = LinearAlgebra::contract( + auto boo_t_foo_in = Picasso::LinearAlgebra::contract( boo( index ), foo_in( index ), SpaceDim<2>{} ); auto boo_t_foo_in_t_fez_in = boo_t_foo_in * fez_in( index ); @@ -247,12 +248,12 @@ struct GridTensor3Func { { 9, 10, 11 }, { 12, 13, 14 }, { 15, 16, 17 } }, { { 18, 19, 20 }, { 21, 22, 23 }, { 24, 25, 26 } } }; - mat_i_out( index ) = - LinearAlgebra::contract( tensor, foo_in( index ), SpaceDim<0>() ); - mat_j_out( index ) = - LinearAlgebra::contract( tensor, foo_in( index ), SpaceDim<1>() ); - mat_k_out( index ) = - LinearAlgebra::contract( tensor, foo_in( index ), SpaceDim<2>() ); + mat_i_out( index ) = Picasso::LinearAlgebra::contract( + tensor, foo_in( index ), SpaceDim<0>() ); + mat_j_out( index ) = Picasso::LinearAlgebra::contract( + tensor, foo_in( index ), SpaceDim<1>() ); + mat_k_out( index ) = Picasso::LinearAlgebra::contract( + tensor, foo_in( index ), SpaceDim<2>() ); } }; @@ -304,7 +305,7 @@ struct GridTensor4Func { 0.0, lam, 0.0 }, { 0.0, 0.0, lam + 2 * mu } } } }; - Picasso::Mat3 strain = { + Cabana::Mat3 strain = { { 0.5, 1.0, 0 }, { 1.0, 0, 0 }, { 0, 0, 0 } }; const int index[3] = { i, j, k }; @@ -313,7 +314,8 @@ struct GridTensor4Func // "generalized" Hook's law, which is just a double contraction of a // fourth-order stiffness tensor with a strain tensor. baz_out is the // stress tensor in this example - baz_out( index ) = LinearAlgebra::contract( cam( index ), strain ); + baz_out( index ) = + Picasso::LinearAlgebra::contract( cam( index ), strain ); } }; @@ -353,7 +355,7 @@ void gatherScatterTest() const double, particle_type& p ) { for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Field::LogicalPosition<3>(), d ) = x[d]; + Cabana::get( p, Field::LogicalPosition<3>(), d ) = x[d]; return true; }; diff --git a/unit_test/tstMarchingCubes.hpp b/unit_test/tstMarchingCubes.hpp index 22c3044..ba7c770 100644 --- a/unit_test/tstMarchingCubes.hpp +++ b/unit_test/tstMarchingCubes.hpp @@ -102,13 +102,13 @@ void runTest( const Phi0& phi_0, const std::string& stl_filename ) mc_data->facets ); for ( std::size_t f = 0; f < facets.extent( 0 ); ++f ) { - Vec3 ab = { facets( f, 1, 0 ) - facets( f, 0, 0 ), - facets( f, 1, 1 ) - facets( f, 0, 1 ), - facets( f, 1, 2 ) - facets( f, 0, 2 ) }; + Cabana::Vec3 ab = { facets( f, 1, 0 ) - facets( f, 0, 0 ), + facets( f, 1, 1 ) - facets( f, 0, 1 ), + facets( f, 1, 2 ) - facets( f, 0, 2 ) }; - Vec3 ac = { facets( f, 2, 0 ) - facets( f, 0, 0 ), - facets( f, 2, 1 ) - facets( f, 0, 1 ), - facets( f, 2, 2 ) - facets( f, 0, 2 ) }; + Cabana::Vec3 ac = { facets( f, 2, 0 ) - facets( f, 0, 0 ), + facets( f, 2, 1 ) - facets( f, 0, 1 ), + facets( f, 2, 2 ) - facets( f, 0, 2 ) }; auto cross = ab % ac; diff --git a/unit_test/tstParticleInterpolation.hpp b/unit_test/tstParticleInterpolation.hpp index b790792..afa2668 100644 --- a/unit_test/tstParticleInterpolation.hpp +++ b/unit_test/tstParticleInterpolation.hpp @@ -73,12 +73,13 @@ struct ScalarValueP2G scatter_deps.get( FieldLocation::Node(), NodeScalar() ); // Get particle data. - auto particle_scalar = Picasso::get( particle, ParticleScalar() ); + auto particle_scalar = Cabana::get( particle, ParticleScalar() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue() ); + Cabana::get( particle, Field::LogicalPosition<3>() ), + SplineValue() ); // Interpolate to grid. P2G::value( spline, particle_scalar, node_scalar ); @@ -101,12 +102,13 @@ struct VectorValueP2G scatter_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto particle_vector = get( particle, ParticleVector() ); + auto particle_vector = Cabana::get( particle, ParticleVector() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue() ); + Cabana::get( particle, Field::LogicalPosition<3>() ), + SplineValue() ); // Interpolate to grid. P2G::value( spline, particle_vector, node_vector ); @@ -129,12 +131,12 @@ struct ScalarGradientP2G scatter_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto particle_scalar = Picasso::get( particle, ParticleScalar() ); + auto particle_scalar = Cabana::get( particle, ParticleScalar() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -158,12 +160,12 @@ struct VectorDivergenceP2G scatter_deps.get( FieldLocation::Node(), NodeScalar() ); // Get particle data. - auto particle_vector = get( particle, ParticleVector() ); + auto particle_vector = Cabana::get( particle, ParticleVector() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -187,12 +189,12 @@ struct TensorDivergenceP2G scatter_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto particle_tensor = get( particle, ParticleTensor() ); + auto particle_tensor = Cabana::get( particle, ParticleTensor() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -217,12 +219,13 @@ struct ScalarValueG2P gather_deps.get( FieldLocation::Node(), NodeScalar() ); // Get particle data. - auto& particle_scalar = Picasso::get( particle, ParticleScalar() ); + auto& particle_scalar = Cabana::get( particle, ParticleScalar() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue() ); + Cabana::get( particle, Field::LogicalPosition<3>() ), + SplineValue() ); // Interpolate to grid. G2P::value( spline, node_scalar, particle_scalar ); @@ -246,12 +249,13 @@ struct VectorValueG2P gather_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto particle_vector = get( particle, ParticleVector() ); + auto particle_vector = Cabana::get( particle, ParticleVector() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue() ); + Cabana::get( particle, Field::LogicalPosition<3>() ), + SplineValue() ); // Interpolate to grid. G2P::value( spline, node_vector, particle_vector ); @@ -275,12 +279,12 @@ struct ScalarGradientG2P gather_deps.get( FieldLocation::Node(), NodeScalar() ); // Get particle data. - auto particle_vector = get( particle, ParticleVector() ); + auto particle_vector = Cabana::get( particle, ParticleVector() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -305,12 +309,12 @@ struct VectorGradientG2P gather_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto particle_tensor = get( particle, ParticleTensor() ); + auto particle_tensor = Cabana::get( particle, ParticleTensor() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -335,12 +339,12 @@ struct VectorDivergenceG2P gather_deps.get( FieldLocation::Node(), NodeVector() ); // Get particle data. - auto& particle_scalar = Picasso::get( particle, ParticleScalar() ); + auto& particle_scalar = Cabana::get( particle, ParticleScalar() ); // Node Interpolant auto spline = createSpline( FieldLocation::Node(), InterpolationOrder<1>(), local_mesh, - get( particle, Field::LogicalPosition<3>() ), SplineValue(), + Cabana::get( particle, Field::LogicalPosition<3>() ), SplineValue(), SplineGradient() ); // Interpolate to grid. @@ -389,7 +393,7 @@ void interpolationTest() const double, particle_type& p ) { for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Field::LogicalPosition<3>(), d ) = x[d]; + Cabana::get( p, Field::LogicalPosition<3>(), d ) = x[d]; return true; }; diff --git a/unit_test/tstParticleList.hpp b/unit_test/tstParticleList.hpp index ec554a1..58e6afc 100644 --- a/unit_test/tstParticleList.hpp +++ b/unit_test/tstParticleList.hpp @@ -104,15 +104,15 @@ void linearAlgebraTest() KOKKOS_LAMBDA( const int p ) { auto particle = particles.getParticleView( p ); - auto px_v = get( particle, Field::LogicalPosition<3>() ); + auto px_v = Cabana::get( particle, Field::LogicalPosition<3>() ); for ( int d = 0; d < 3; ++d ) px_v( d ) += p + d; - Picasso::get( particle, Foo() ) += p; + Cabana::get( particle, Foo() ) += p; - Picasso::get( particle, Field::Color() ) += p; + Cabana::get( particle, Field::Color() ) += p; - auto pf_m = Picasso::get( particle, Bar() ); + auto pf_m = Cabana::get( particle, Bar() ); for ( int i = 0; i < 3; ++i ) for ( int j = 0; j < 3; ++j ) pf_m( i, j ) += p + i + j; diff --git a/unit_test/tstPolyPIC.hpp b/unit_test/tstPolyPIC.hpp index 0b72697..11f5895 100644 --- a/unit_test/tstPolyPIC.hpp +++ b/unit_test/tstPolyPIC.hpp @@ -9,7 +9,6 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ -#include #include #include #include @@ -266,12 +265,12 @@ void collocatedTest() Kokkos::parallel_for( "g2p", Kokkos::RangePolicy( 0, 1 ), KOKKOS_LAMBDA( const int ) { - Vec3 x = { px, py, pz }; + Cabana::Vec3 x = { px, py, pz }; auto sd = createSpline( Location(), InterpolationOrder(), local_mesh, x, SplineValue(), SplineGradient() ); - LinearAlgebra::Matrix modes; + Cabana::LinearAlgebra::Matrix modes; PolyPIC::g2p( gv_wrapper, modes, sd ); for ( int r = 0; r < num_mode; ++r ) @@ -303,12 +302,12 @@ void collocatedTest() Kokkos::parallel_for( "p2g", Kokkos::RangePolicy( 0, 1 ), KOKKOS_LAMBDA( const int ) { - Vec3 x = { px, py, pz }; + Cabana::Vec3 x = { px, py, pz }; auto sd = createSpline( Location(), InterpolationOrder(), local_mesh, x, SplineValue(), SplineGradient(), SplineDistance() ); - LinearAlgebra::Matrix modes; + Cabana::LinearAlgebra::Matrix modes; for ( int r = 0; r < num_mode; ++r ) for ( int d = 0; d < 3; ++d ) modes( r, d ) = pc( r, d ); @@ -431,12 +430,12 @@ void staggeredTest() Kokkos::parallel_for( "g2p", Kokkos::RangePolicy( 0, 1 ), KOKKOS_LAMBDA( const int ) { - Vec3 x = { px, py, pz }; + Cabana::Vec3 x = { px, py, pz }; auto sd = createSpline( Location(), InterpolationOrder(), local_mesh, x, SplineValue(), SplineGradient() ); - LinearAlgebra::Matrix modes = 0.0; + Cabana::LinearAlgebra::Matrix modes = 0.0; PolyPIC::g2p( gs_wrapper, modes, sd ); for ( int r = 0; r < num_mode; ++r ) @@ -472,12 +471,12 @@ void staggeredTest() Kokkos::parallel_for( "p2g", Kokkos::RangePolicy( 0, 1 ), KOKKOS_LAMBDA( const int ) { - Vec3 x = { px, py, pz }; + Cabana::Vec3 x = { px, py, pz }; auto sd = createSpline( Location(), InterpolationOrder(), local_mesh, x, SplineValue(), SplineGradient(), SplineDistance() ); - LinearAlgebra::Matrix modes; + Cabana::LinearAlgebra::Matrix modes; for ( int r = 0; r < num_mode; ++r ) for ( int d = 0; d < 3; ++d ) modes( r, d ) = pc( r, d ); diff --git a/unit_test/tstUniformCartesianMeshMapping.hpp b/unit_test/tstUniformCartesianMeshMapping.hpp index 55de3a9..ceeffb6 100644 --- a/unit_test/tstUniformCartesianMeshMapping.hpp +++ b/unit_test/tstUniformCartesianMeshMapping.hpp @@ -9,7 +9,6 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ -#include #include #include #include @@ -109,17 +108,17 @@ void mappingTest3d() "check_mapping", TEST_EXECSPACE{}, ghosted_cells, KOKKOS_LAMBDA( const int i, const int j, const int k ) { // Map to physical frame. - LinearAlgebra::Vector cell_coords = 0.0; + Cabana::LinearAlgebra::Vector cell_coords = 0.0; int ijk[3] = { i, j, k }; local_mesh.coordinates( Cabana::Grid::Cell{}, ijk, cell_coords.data() ); - LinearAlgebra::VectorView phys_coords( + Cabana::LinearAlgebra::VectorView phys_coords( &cell_forward_map( i, j, k, 0 ), cell_forward_map.stride( 3 ) ); CurvilinearMeshMapping::mapToPhysicalFrame( mapping, cell_coords, phys_coords ); // Map to reference frame. - LinearAlgebra::VectorView ref_coords( + Cabana::LinearAlgebra::VectorView ref_coords( &cell_reverse_map( i, j, k, 0 ), cell_reverse_map.stride( 3 ) ); ref_coords = { cell_coords( Dim::I ) - 0.1, cell_coords( Dim::J ) - 0.1, @@ -129,7 +128,7 @@ void mappingTest3d() mapping, phys_coords, ref_coords ); // Default map to reference frame. - LinearAlgebra::VectorView default_ref_coords( + Cabana::LinearAlgebra::VectorView default_ref_coords( &default_cell_reverse_map( i, j, k, 0 ), default_cell_reverse_map.stride( 3 ) ); default_ref_coords = { cell_coords( Dim::I ) - 0.1, @@ -272,17 +271,17 @@ void mappingTest2d() "check_mapping", TEST_EXECSPACE{}, ghosted_cells, KOKKOS_LAMBDA( const int i, const int j ) { // Map to physical frame. - LinearAlgebra::Vector cell_coords = 0.0; + Cabana::LinearAlgebra::Vector cell_coords = 0.0; int ijk[2] = { i, j }; local_mesh.coordinates( Cabana::Grid::Cell{}, ijk, cell_coords.data() ); - LinearAlgebra::VectorView phys_coords( + Cabana::LinearAlgebra::VectorView phys_coords( &cell_forward_map( i, j, 0 ), cell_forward_map.stride( 2 ) ); CurvilinearMeshMapping::mapToPhysicalFrame( mapping, cell_coords, phys_coords ); // Map to reference frame. - LinearAlgebra::VectorView ref_coords( + Cabana::LinearAlgebra::VectorView ref_coords( &cell_reverse_map( i, j, 0 ), cell_reverse_map.stride( 2 ) ); ref_coords = { cell_coords( Dim::I ) - 0.1, cell_coords( Dim::J ) - 0.1 }; @@ -291,7 +290,7 @@ void mappingTest2d() mapping, phys_coords, ref_coords ); // Default map to reference frame. - LinearAlgebra::VectorView default_ref_coords( + Cabana::LinearAlgebra::VectorView default_ref_coords( &default_cell_reverse_map( i, j, 0 ), default_cell_reverse_map.stride( 2 ) ); default_ref_coords = { cell_coords( Dim::I ) - 0.1,