3636#ifdef MP_UNITS_IMPORT_STD
3737import std;
3838#else
39- #include < cmath>
4039#include < concepts>
4140#include < cstddef>
4241#include < type_traits>
42+ #include < utility>
4343#if MP_UNITS_HOSTED
4444#include < ostream>
4545#endif
@@ -48,22 +48,75 @@ import std;
4848
4949namespace mp_units {
5050
51- MP_UNITS_EXPORT template <detail::Scalar T, std::size_t R, std::size_t C>
51+ template <class X >
52+ struct is_cartesian_vector : std::false_type {};
53+ template <class U >
54+ struct is_cartesian_vector <cartesian_vector<U>> : std::true_type {};
55+
56+ MP_UNITS_EXPORT template <detail::Scalar U, std::size_t R, std::size_t C>
5257class cartesian_tensor ;
5358
59+ template <class X >
60+ struct is_cartesian_tensor : std::false_type {};
61+ template <class U , std::size_t R, std::size_t C>
62+ struct is_cartesian_tensor <cartesian_tensor<U, R, C>> : std::true_type {};
63+
5464MP_UNITS_EXPORT template <detail::Scalar T, std::size_t R, std::size_t C>
5565class cartesian_tensor {
5666 static_assert (R >= 1 && R <= 3 && C >= 1 && C <= 3 , " cartesian_tensor supports sizes up to 3x3" );
5767
5868public:
59- T _data_[R * C];
6069 using value_type = T;
6170 static constexpr std::size_t rows_v = R;
6271 static constexpr std::size_t cols_v = C;
6372
73+ // aggregate on purpose
74+ T _data_[R * C];
75+
6476 [[nodiscard]] constexpr T& operator ()(std::size_t r, std::size_t c) { return _data_[r * C + c]; }
6577 [[nodiscard]] constexpr const T& operator ()(std::size_t r, std::size_t c) const { return _data_[r * C + c]; }
6678
79+ [[nodiscard]] constexpr cartesian_tensor operator +() const { return *this ; }
80+ [[nodiscard]] constexpr cartesian_tensor operator -() const
81+ {
82+ cartesian_tensor out{};
83+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = -_data_[i];
84+ return out;
85+ }
86+
87+ template <class U >
88+ requires requires (T& t, const U& u) { t += u; }
89+ constexpr cartesian_tensor& operator +=(const cartesian_tensor<U, R, C>& rhs)
90+ {
91+ for (std::size_t i = 0 ; i < R * C; ++i) _data_[i] += rhs._data_ [i];
92+ return *this ;
93+ }
94+
95+ template <class U >
96+ requires requires (T& t, const U& u) { t -= u; }
97+ constexpr cartesian_tensor& operator -=(const cartesian_tensor<U, R, C>& rhs)
98+ {
99+ for (std::size_t i = 0 ; i < R * C; ++i) _data_[i] -= rhs._data_ [i];
100+ return *this ;
101+ }
102+
103+ template <class S >
104+ requires detail::Scalar<S> && (!is_cartesian_vector<std::remove_cvref_t <S>>::value) &&
105+ (!is_cartesian_tensor<std::remove_cvref_t <S>>::value) && requires (T& t, const S& s) { t *= s; }
106+ constexpr cartesian_tensor& operator *=(const S& scalar)
107+ {
108+ for (std::size_t i = 0 ; i < R * C; ++i) _data_[i] *= scalar;
109+ return *this ;
110+ }
111+
112+ template <class S >
113+ requires detail::Scalar<S> && (!is_cartesian_vector<std::remove_cvref_t <S>>::value) &&
114+ (!is_cartesian_tensor<std::remove_cvref_t <S>>::value) && requires (T& t, const S& s) { t /= s; }
115+ constexpr cartesian_tensor& operator /=(const S& scalar)
116+ {
117+ for (std::size_t i = 0 ; i < R * C; ++i) _data_[i] /= scalar;
118+ return *this ;
119+ }
67120
68121#if MP_UNITS_HOSTED
69122 friend std::ostream& operator <<(std::ostream& os, const cartesian_tensor& A)
@@ -79,130 +132,164 @@ class cartesian_tensor {
79132 return os;
80133 }
81134#endif
135+
136+ template <class U >
137+ requires requires (const T& t, const U& u) { t + u; }
138+ [[nodiscard]] friend constexpr auto operator +(const cartesian_tensor& A, const cartesian_tensor<U, R, C>& B)
139+ {
140+ using CT = std::common_type_t <T, U>;
141+ cartesian_tensor<CT, R, C> out{};
142+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = static_cast <CT>(A._data_ [i]) + static_cast <CT>(B._data_ [i]);
143+ return out;
144+ }
145+
146+ template <class U >
147+ requires requires (const T& t, const U& u) { t - u; }
148+ [[nodiscard]] friend constexpr auto operator -(const cartesian_tensor& A, const cartesian_tensor<U, R, C>& B)
149+ {
150+ using CT = std::common_type_t <T, U>;
151+ cartesian_tensor<CT, R, C> out{};
152+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = static_cast <CT>(A._data_ [i]) - static_cast <CT>(B._data_ [i]);
153+ return out;
154+ }
155+
156+ template <class U >
157+ requires (!treat_as_floating_point<T> && !treat_as_floating_point<U> && requires (const T& t, const U& u) { t % u; })
158+ [[nodiscard]] friend constexpr auto operator %(const cartesian_tensor& A, const cartesian_tensor<U, R, C>& B)
159+ {
160+ using CT = std::common_type_t <T, U>;
161+ cartesian_tensor<CT, R, C> out{};
162+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = static_cast <CT>(A._data_ [i] % B._data_ [i]);
163+ return out;
164+ }
165+
166+ template <class S >
167+ requires detail::Scalar<S> && (!is_cartesian_vector<std::remove_cvref_t <S>>::value) &&
168+ (!is_cartesian_tensor<std::remove_cvref_t <S>>::value) && requires (const T& t, const S& s) { t * s; }
169+ [[nodiscard]] friend constexpr auto operator *(const cartesian_tensor& A, const S& scalar)
170+ {
171+ using CT = std::common_type_t <T, S>;
172+ cartesian_tensor<CT, R, C> out{};
173+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = static_cast <CT>(A._data_ [i]) * static_cast <CT>(scalar);
174+ return out;
175+ }
176+
177+ template <class S >
178+ requires detail::Scalar<S> && (!is_cartesian_vector<std::remove_cvref_t <S>>::value) &&
179+ (!is_cartesian_tensor<std::remove_cvref_t <S>>::value) && requires (const S& s, const T& t) { s * t; }
180+ [[nodiscard]] friend constexpr auto operator *(const S& scalar, const cartesian_tensor& A)
181+ {
182+ return A * scalar;
183+ }
184+
185+ template <class S >
186+ requires detail::Scalar<S> && (!is_cartesian_vector<std::remove_cvref_t <S>>::value) &&
187+ (!is_cartesian_tensor<std::remove_cvref_t <S>>::value) && requires (const T& t, const S& s) { t / s; }
188+ [[nodiscard]] friend constexpr auto operator /(const cartesian_tensor& A, const S& scalar)
189+ {
190+ using CT = std::common_type_t <T, S>;
191+ cartesian_tensor<CT, R, C> out{};
192+ for (std::size_t i = 0 ; i < R * C; ++i) out._data_ [i] = static_cast <CT>(A._data_ [i]) / static_cast <CT>(scalar);
193+ return out;
194+ }
195+
196+ template <class U >
197+ requires std::equality_comparable_with<T, U>
198+ [[nodiscard]] friend constexpr bool operator ==(const cartesian_tensor& A, const cartesian_tensor<U, R, C>& B)
199+ {
200+ for (std::size_t i = 0 ; i < R * C; ++i)
201+ if (!(A._data_ [i] == B._data_ [i])) return false ;
202+ return true ;
203+ }
82204};
83205
206+ template <class T , std::size_t R, std::size_t C>
207+ inline constexpr bool treat_as_floating_point<cartesian_tensor<T, R, C>> = treat_as_floating_point<T>;
208+
209+ template <class S , class T , std::size_t R, std::size_t C>
210+ inline constexpr bool is_value_preserving<S, cartesian_tensor<T, R, C>> = is_value_preserving<S, T>;
211+
84212
85- template <typename T, typename U, std::size_t R, std::size_t K, std::size_t C>
213+ template <class T , class U , std::size_t R, std::size_t K, std::size_t C>
86214[[nodiscard]] constexpr auto matmul (const cartesian_tensor<T, R, K>& A, const cartesian_tensor<U, K, C>& B)
87215{
88216 using CT = std::common_type_t <T, U>;
89- cartesian_tensor<CT, R, C> Rm {};
90- for (std::size_t r = 0 ; r < R; ++r)
217+ cartesian_tensor<CT, R, C> out {};
218+ for (std::size_t r = 0 ; r < R; ++r) {
91219 for (std::size_t c = 0 ; c < C; ++c) {
92220 CT acc{};
93- for (std::size_t k = 0 ; k < K; ++k) acc += static_cast <CT>(A (r, k)) * static_cast <CT>(B (k, c));
94- Rm (r, c) = acc;
221+ for (std::size_t k = 0 ; k < K; ++k) acc = acc + static_cast <CT>(A (r, k)) * static_cast <CT>(B (k, c));
222+ out (r, c) = acc;
95223 }
96- return Rm;
224+ }
225+ return out;
97226}
98227
99- template <typename T, typename U>
100- [[nodiscard]] constexpr auto matvec (const cartesian_tensor<T, 3 , 3 >& tensor , const cartesian_vector<U>& vector )
228+ template <class T , class U >
229+ [[nodiscard]] constexpr auto matvec (const cartesian_tensor<T, 3 , 3 >& M , const cartesian_vector<U>& x )
101230{
102231 using CT = std::common_type_t <T, U>;
103232 cartesian_vector<CT> y{};
104- for (std::size_t r = 0 ; r < 3 ; ++r) {
105- CT acc{};
106- for (std::size_t c = 0 ; c < 3 ; ++c) acc += static_cast <CT>(tensor (r, c)) * static_cast <CT>(vector[c]);
107- y[r] = acc;
108- }
233+ y[0 ] = static_cast <CT>(M (0 , 0 )) * static_cast <CT>(x[0 ]) + static_cast <CT>(M (0 , 1 )) * static_cast <CT>(x[1 ]) +
234+ static_cast <CT>(M (0 , 2 )) * static_cast <CT>(x[2 ]);
235+ y[1 ] = static_cast <CT>(M (1 , 0 )) * static_cast <CT>(x[0 ]) + static_cast <CT>(M (1 , 1 )) * static_cast <CT>(x[1 ]) +
236+ static_cast <CT>(M (1 , 2 )) * static_cast <CT>(x[2 ]);
237+ y[2 ] = static_cast <CT>(M (2 , 0 )) * static_cast <CT>(x[0 ]) + static_cast <CT>(M (2 , 1 )) * static_cast <CT>(x[1 ]) +
238+ static_cast <CT>(M (2 , 2 )) * static_cast <CT>(x[2 ]);
109239 return y;
110240}
111241
112- template <typename T, typename U, std::size_t R, std::size_t C>
242+ template <class T , class U , std::size_t R, std::size_t C>
113243[[nodiscard]] constexpr auto double_contraction (const cartesian_tensor<T, R, C>& A, const cartesian_tensor<U, R, C>& B)
114244{
115245 using CT = std::common_type_t <T, U>;
116246 CT acc{};
117- for (std::size_t i = 0 ; i < R * C; ++i) acc += static_cast <CT>(A._data_ [i]) * static_cast <CT>(B._data_ [i]);
247+ for (std::size_t i = 0 ; i < R * C; ++i) acc = acc + static_cast <CT>(A._data_ [i]) * static_cast <CT>(B._data_ [i]);
118248 return acc;
119249}
120250
121-
122- template <typename T, typename U>
123- [[nodiscard]] constexpr auto outer_numeric (const cartesian_vector<T>& lhs, const cartesian_vector<U>& rhs)
251+ template <class T , class U >
252+ [[nodiscard]] constexpr auto outer_numeric (const cartesian_vector<T>& a, const cartesian_vector<U>& b)
124253{
125254 using CT = std::common_type_t <T, U>;
126- cartesian_tensor<CT, 3 , 3 > Rm {};
255+ cartesian_tensor<CT, 3 , 3 > out {};
127256 for (std::size_t i = 0 ; i < 3 ; ++i)
128- for (std::size_t j = 0 ; j < 3 ; ++j) Rm (i, j) = static_cast <CT>(lhs[i]) * static_cast <CT>(rhs[j]);
129- return Rm;
130- }
131-
132- template <typename T, typename U, std::size_t R, std::size_t C>
133- requires requires (const T& t, const U& u) { t + u; }
134- [[nodiscard]] constexpr auto operator +(const cartesian_tensor<T, R, C>& A, const cartesian_tensor<U, R, C>& B)
135- {
136- using CT = std::common_type_t <T, U>;
137- cartesian_tensor<CT, R, C> Rm{};
138- for (std::size_t i = 0 ; i < R * C; ++i) Rm._data_ [i] = static_cast <CT>(A._data_ [i]) + static_cast <CT>(B._data_ [i]);
139- return Rm;
140- }
141-
142- template <typename T, typename U, std::size_t R, std::size_t C>
143- requires requires (const T& t, const U& u) { t - u; }
144- [[nodiscard]] constexpr auto operator -(const cartesian_tensor<T, R, C>& A, const cartesian_tensor<U, R, C>& B)
145- {
146- using CT = std::common_type_t <T, U>;
147- cartesian_tensor<CT, R, C> Rm{};
148- for (std::size_t i = 0 ; i < R * C; ++i) Rm._data_ [i] = static_cast <CT>(A._data_ [i]) - static_cast <CT>(B._data_ [i]);
149- return Rm;
150- }
151-
152- template <typename T, typename U, std::size_t R, std::size_t C>
153- requires (!treat_as_floating_point<T> && !treat_as_floating_point<U> && requires (const T& t, const U& u) { t % u; })
154- [[nodiscard]] constexpr auto operator %(const cartesian_tensor<T, R, C>& A, const cartesian_tensor<U, R, C>& B)
155- {
156- using CT = std::common_type_t <T, U>;
157- cartesian_tensor<CT, R, C> Rm{};
158- for (std::size_t i = 0 ; i < R * C; ++i) Rm._data_ [i] = static_cast <CT>(A._data_ [i] % B._data_ [i]);
159- return Rm;
160- }
161-
162- template <typename T, typename S, std::size_t R, std::size_t C>
163- requires requires (const T& t, const S& s) { t * s; }
164- [[nodiscard]] constexpr auto operator *(const cartesian_tensor<T, R, C>& tensor, const S& scalar)
165- {
166- using CT = std::common_type_t <T, S>;
167- cartesian_tensor<CT, R, C> Rm{};
168- for (std::size_t i = 0 ; i < R * C; ++i) Rm._data_ [i] = static_cast <CT>(tensor._data_ [i]) * static_cast <CT>(scalar);
169- return Rm;
170- }
171-
172- template <typename S, typename U, std::size_t R, std::size_t C>
173- requires requires (const S& s, const U& u) { s * u; }
174- [[nodiscard]] constexpr auto operator *(const S& scalar, const cartesian_tensor<U, R, C>& tensor)
175- {
176- return tensor * scalar;
177- }
178-
179- template <typename T, typename S, std::size_t R, std::size_t C>
180- requires requires (const T& t, const S& s) { t / s; }
181- [[nodiscard]] constexpr auto operator /(const cartesian_tensor<T, R, C>& tensor, const S& scalar)
182- {
183- using CT = std::common_type_t <T, S>;
184- cartesian_tensor<CT, R, C> Rm{};
185- for (std::size_t i = 0 ; i < R * C; ++i) Rm._data_ [i] = static_cast <CT>(tensor._data_ [i]) / static_cast <CT>(scalar);
186- return Rm;
257+ for (std::size_t j = 0 ; j < 3 ; ++j) out (i, j) = static_cast <CT>(a[i]) * static_cast <CT>(b[j]);
258+ return out;
187259}
188260
189261} // namespace mp_units
190262
263+ namespace std {
264+ template <class T , class U , std::size_t R, std::size_t C>
265+ struct common_type <mp_units::cartesian_tensor<T, R, C>, mp_units::cartesian_tensor<U, R, C>> {
266+ using type = mp_units::cartesian_tensor<common_type_t <T, U>, R, C>;
267+ };
268+ template <class T , class U , std::size_t R, std::size_t C>
269+ struct common_type <mp_units::cartesian_tensor<T, R, C>, U> {
270+ using type = mp_units::cartesian_tensor<common_type_t <T, U>, R, C>;
271+ };
272+ template <class T , class U , std::size_t R, std::size_t C>
273+ struct common_type <U, mp_units::cartesian_tensor<T, R, C>> {
274+ using type = mp_units::cartesian_tensor<common_type_t <U, T>, R, C>;
275+ };
276+ } // namespace std
277+
191278#if MP_UNITS_HOSTED
192279template <typename T, std::size_t R, std::size_t C, typename Char>
193280struct MP_UNITS_STD_FMT ::formatter<mp_units::cartesian_tensor<T, R, C>, Char> :
194281 formatter<std::basic_string_view<Char>, Char> {
195- template <typename Ctx >
196- auto format (const mp_units::cartesian_tensor<T, R, C>& A, Ctx & ctx) const
282+ template <typename FormatContext >
283+ auto format (const mp_units::cartesian_tensor<T, R, C>& A, FormatContext & ctx) const
197284 {
198285 auto out = ctx.out ();
199286 for (std::size_t r = 0 ; r < R; ++r) {
200- out = format_to (out, " {}" , (r == 0 ? " [[" : " [" ));
287+ out = MP_UNITS_STD_FMT:: format_to (out, " {}" , (r == 0 ? " [[" : " [" ));
201288 for (std::size_t c = 0 ; c < C; ++c) {
202- out = format_to (out, " {}" , A (r, c));
203- if (c + 1 != C) out = format_to (out, " {}" , " , " );
289+ out = MP_UNITS_STD_FMT:: format_to (out, " {}" , A (r, c));
290+ if (c + 1 != C) out = MP_UNITS_STD_FMT:: format_to (out, " {}" , " , " );
204291 }
205- out = format_to (out, " {}" , (r + 1 == R ? " ]]" : " ]\n " ));
292+ out = MP_UNITS_STD_FMT:: format_to (out, " {}" , (r + 1 == R ? " ]]" : " ]\n " ));
206293 }
207294 return out;
208295 }
0 commit comments