Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines classes for Tensor
6 :
7 : #pragma once
8 :
9 : #include <concepts>
10 : #include <cstddef>
11 : #include <pup.h>
12 : #include <pup_stl.h>
13 :
14 : #include "DataStructures/Tensor/Expressions/AddSubtract.hpp"
15 : #include "DataStructures/Tensor/Expressions/Contract.hpp"
16 : #include "DataStructures/Tensor/Expressions/DataTypeSupport.hpp"
17 : #include "DataStructures/Tensor/Expressions/Divide.hpp"
18 : #include "DataStructures/Tensor/Expressions/Evaluate.hpp"
19 : #include "DataStructures/Tensor/Expressions/IndexPropertyCheck.hpp"
20 : #include "DataStructures/Tensor/Expressions/LhsTensorSymmAndIndices.hpp"
21 : #include "DataStructures/Tensor/Expressions/Negate.hpp"
22 : #include "DataStructures/Tensor/Expressions/NumberAsExpression.hpp"
23 : #include "DataStructures/Tensor/Expressions/Product.hpp"
24 : #include "DataStructures/Tensor/Expressions/SpatialSpacetimeIndex.hpp"
25 : #include "DataStructures/Tensor/Expressions/SquareRoot.hpp"
26 : #include "DataStructures/Tensor/Expressions/TensorAsExpression.hpp"
27 : #include "DataStructures/Tensor/Expressions/TensorExpression.hpp"
28 : #include "DataStructures/Tensor/Expressions/TensorIndex.hpp"
29 : #include "DataStructures/Tensor/Expressions/TensorIndexTransformation.hpp"
30 : #include "DataStructures/Tensor/Expressions/TimeIndex.hpp"
31 : #include "DataStructures/Tensor/IndexType.hpp"
32 : #include "DataStructures/Tensor/Structure.hpp"
33 : #include "DataStructures/Tensor/TypeAliases.hpp"
34 : #include "Utilities/Autodiff/Autodiff.hpp"
35 : #include "Utilities/ErrorHandling/Error.hpp"
36 : #include "Utilities/ForceInline.hpp"
37 : #include "Utilities/Gsl.hpp"
38 : #include "Utilities/MakeArray.hpp"
39 : #include "Utilities/MakeWithValue.hpp"
40 : #include "Utilities/PrettyType.hpp"
41 : #include "Utilities/SetNumberOfGridPoints.hpp"
42 : #include "Utilities/Simd/Simd.hpp"
43 : #include "Utilities/StdHelpers.hpp"
44 : #include "Utilities/TMPL.hpp"
45 : #include "Utilities/TypeTraits.hpp"
46 : #include "Utilities/TypeTraits/IsStreamable.hpp"
47 : // Not all the TensorExpression includes are necessary, but we include them so
48 : // that Tensor.hpp provides a uniform interface to Tensors and TensorExpressions
49 :
50 : /// \cond
51 : template <typename X, typename Symm = Symmetry<>,
52 : typename IndexList = index_list<>>
53 : class Tensor;
54 : /// \endcond
55 :
56 : namespace Tensor_detail {
57 : template <typename T, typename = void>
58 : struct is_tensor : std::false_type {};
59 : template <typename X, typename Symm, typename IndexList>
60 : struct is_tensor<Tensor<X, Symm, IndexList>> : std::true_type {};
61 : } // namespace Tensor_detail
62 :
63 : /*!
64 : * \ingroup TensorGroup
65 : * \brief Represents an object with multiple components
66 : *
67 : * \details
68 : * Tensor is a container that represents indexable geometric objects. Each index
69 : * has its own dimension, valence, and frame and must be either spatial or
70 : * spacetime. Note that the dimension passed to `SpatialIndex` and
71 : * `SpacetimeIndex` is always the spatial dimension of the index. Tensors with
72 : * symmetric indices are stored only once and must be of the same
73 : * type. A list of available type aliases can be found in the ::tnsr namespace
74 : * where the adopted conventions are:
75 : *
76 : * 1. Upper case for contravariant or upper indices, lower case for covariant or
77 : * lower indices.
78 : *
79 : * 2. `a, b, c, d` are used for spacetime indices while `i, j, k, l` are used
80 : * for spatial indices.
81 : *
82 : * 3. `::Scalar` is not inside the `::tnsr` namespace but is used to represent
83 : * a scalar with no indices.
84 : *
85 : * \example
86 : * \snippet Test_Tensor.cpp scalar
87 : * \snippet Test_Tensor.cpp spatial_vector
88 : * \snippet Test_Tensor.cpp spacetime_vector
89 : * \snippet Test_Tensor.cpp rank_3_122
90 : *
91 : * \tparam X the type held
92 : * \tparam Symm the ::Symmetry of the indices
93 : * \tparam IndexList a typelist of \ref SpacetimeIndex "TensorIndexType"'s
94 : */
95 : template <typename X, typename Symm, template <typename...> class IndexList,
96 : typename... Indices>
97 1 : class Tensor<X, Symm, IndexList<Indices...>> {
98 : static_assert(sizeof...(Indices) < 5,
99 : "If you are sure you need rank 5 or higher Tensor's please "
100 : "file an issue on GitHub or discuss with a core developer of "
101 : "SpECTRE.");
102 : public:
103 : /// The type of the sequence that holds the data
104 1 : using storage_type =
105 : std::array<X, Tensor_detail::Structure<Symm, Indices...>::size()>;
106 : /// The type that is stored by the Tensor
107 1 : using type = X;
108 : /// Typelist of the symmetry of the Tensor
109 : ///
110 : /// \details
111 : /// For a rank-3 tensor symmetric in the last two indices,
112 : /// \f$T_{a(bc)}\f$, the ::Symmetry is `<2, 1, 1>`. For a non-symmetric rank-2
113 : /// tensor the ::Symmetry is `<2, 1>`.
114 1 : using symmetry = Symm;
115 : /// Typelist of the \ref SpacetimeIndex "TensorIndexType"'s that the
116 : /// Tensor has
117 1 : using index_list = tmpl::list<Indices...>;
118 : /// The number of indices that the Tensor has
119 1 : static constexpr size_t num_tensor_indices = sizeof...(Indices);
120 : /// The Tensor_detail::Structure for the particular tensor index structure
121 : ///
122 : /// Each tensor index structure, e.g. \f$T_{ab}\f$, \f$T_a{}^b\f$ or
123 : /// \f$T^{ab}\f$ has its own Tensor_detail::TensorStructure that holds
124 : /// information about how the data is stored, what the multiplicity of the
125 : /// stored indices are, the number of (independent) components, etc.
126 1 : using structure = Tensor_detail::Structure<Symm, Indices...>;
127 : /// The type of the TensorExpression that would represent this Tensor in a
128 : /// tensor expression.
129 : template <typename ArgsList>
130 1 : using TE = tenex::TensorAsExpression<Tensor<X, Symm, IndexList<Indices...>>,
131 : ArgsList>;
132 :
133 0 : Tensor() = default;
134 0 : ~Tensor() = default;
135 0 : Tensor(const Tensor&) = default;
136 0 : Tensor(Tensor&&) = default;
137 0 : Tensor& operator=(const Tensor&) = default;
138 0 : Tensor& operator=(Tensor&&) = default;
139 :
140 : /// Initialize a vector or scalar from an array
141 : ///
142 : /// \example
143 : /// \snippet Test_Tensor.cpp init_vector
144 : /// \param data the values of the individual components of the Vector
145 1 : explicit constexpr Tensor(storage_type data)
146 : requires(sizeof...(Indices) <= 1);
147 :
148 : /// Constructor that passes "args" to constructor of X and initializes each
149 : /// component to be the same
150 : template <typename Arg0, typename... Args>
151 1 : explicit constexpr Tensor(Arg0&& arg0, Args&&... args)
152 : // NOLINTNEXTLINE(readability-simplify-boolean-expr)
153 : requires(not(sizeof...(Args) == 0 and
154 : std::same_as<std::decay_t<Arg0>, Tensor>) and
155 : std::constructible_from<X, Arg0, Args...>);
156 :
157 0 : using value_type = typename storage_type::value_type;
158 0 : using reference = typename storage_type::reference;
159 0 : using const_reference = typename storage_type::const_reference;
160 0 : using iterator = typename storage_type::iterator;
161 0 : using const_iterator = typename storage_type::const_iterator;
162 0 : using pointer = typename storage_type::pointer;
163 0 : using const_pointer = typename storage_type::const_pointer;
164 0 : using reverse_iterator = typename storage_type::reverse_iterator;
165 0 : using const_reverse_iterator = typename storage_type::const_reverse_iterator;
166 :
167 0 : constexpr iterator begin() { return data_.begin(); }
168 0 : constexpr const_iterator begin() const { return data_.begin(); }
169 0 : constexpr const_iterator cbegin() const { return data_.begin(); }
170 :
171 0 : constexpr iterator end() { return data_.end(); }
172 0 : constexpr const_iterator end() const { return data_.end(); }
173 0 : constexpr const_iterator cend() const { return data_.end(); }
174 :
175 0 : constexpr reverse_iterator rbegin() { return data_.rbegin(); }
176 0 : constexpr const_reverse_iterator rbegin() const { return data_.rbegin(); }
177 0 : constexpr const_reverse_iterator crbegin() const { return data_.rbegin(); }
178 :
179 0 : constexpr reverse_iterator rend() { return data_.rend(); }
180 0 : constexpr const_reverse_iterator rend() const { return data_.rend(); }
181 0 : constexpr const_reverse_iterator crend() const { return data_.rend(); }
182 :
183 : /// @{
184 : /// Get data entry using an array representing a tensor index
185 : ///
186 : /// \details
187 : /// Let \f$T_{abc}\f$ be a Tensor.
188 : /// Then `get({{0, 2, 1}})` returns the \f$T_{0 2 1}\f$ component.
189 : /// \param tensor_index the index at which to get the data
190 : template <typename T>
191 1 : SPECTRE_ALWAYS_INLINE constexpr reference get(
192 : const std::array<T, sizeof...(Indices)>& tensor_index) {
193 : return gsl::at(data_, structure::get_storage_index(tensor_index));
194 : }
195 : template <typename T>
196 1 : SPECTRE_ALWAYS_INLINE constexpr const_reference get(
197 : const std::array<T, sizeof...(Indices)>& tensor_index) const {
198 : return gsl::at(data_, structure::get_storage_index(tensor_index));
199 : }
200 : /// @}
201 : /// @{
202 : /// Get data entry using a list of integers representing a tensor index
203 : ///
204 : /// \details
205 : /// Let \f$T_{abc}\f$ be a Tensor.
206 : /// Then `get(0, 2, 1)` returns the \f$T_{0 2 1}\f$ component.
207 : /// \param n the index at which to get the data
208 : template <typename... N>
209 1 : constexpr reference get(N... n) {
210 : static_assert(
211 : sizeof...(Indices) == sizeof...(N),
212 : "the number of tensor indices specified must match the rank of "
213 : "the tensor");
214 : return gsl::at(data_, structure::get_storage_index(n...));
215 : }
216 : template <typename... N>
217 1 : constexpr const_reference get(N... n) const {
218 : static_assert(
219 : sizeof...(Indices) == sizeof...(N),
220 : "the number of tensor indices specified must match the rank of "
221 : "the tensor");
222 : return gsl::at(data_, structure::get_storage_index(n...));
223 : }
224 : /// @}
225 :
226 : /// @{
227 : /// Retrieve the index `N...` by computing the storage index at compile time
228 : // clang-tidy: redundant declaration (bug in clang-tidy)
229 : template <int... N, typename... Args>
230 : friend SPECTRE_ALWAYS_INLINE constexpr typename Tensor<Args...>::reference
231 1 : get(Tensor<Args...>& t); // NOLINT
232 : // clang-tidy: redundant declaration (bug in clang-tidy)
233 : template <int... N, typename... Args>
234 : friend SPECTRE_ALWAYS_INLINE constexpr
235 : typename Tensor<Args...>::const_reference
236 1 : get(const Tensor<Args...>& t); // NOLINT
237 : /// @}
238 :
239 : /// @{
240 : /// Retrieve a TensorExpression object with the index structure passed in
241 : template <typename... TensorIndices>
242 1 : SPECTRE_ALWAYS_INLINE constexpr auto operator()(
243 : TensorIndices... /*meta*/) const {
244 : static_assert((... and tt::is_tensor_index<TensorIndices>::value),
245 : "The tensor expression must be created using TensorIndex "
246 : "objects to represent generic indices, e.g. ti::a, ti::b, "
247 : "etc.");
248 : static_assert(
249 : tenex::tensorindex_list_is_valid<tmpl::list<TensorIndices...>>::value,
250 : "Cannot create a tensor expression with a repeated generic index. "
251 : "(Note that the concrete time indices, ti::T and ti::t, can be "
252 : "repeated.) If you intend to contract, ensure that the indices to "
253 : "contract have opposite valences.");
254 : static_assert(
255 : std::is_same_v<tmpl::integral_list<UpLo, TensorIndices::valence...>,
256 : tmpl::integral_list<UpLo, Indices::ul...>>,
257 : "The valences of the generic indices in the expression do "
258 : "not match the valences of the indices in the Tensor.");
259 : static_assert((... and (not(TensorIndices::is_spacetime and
260 : (Indices::index_type == IndexType::Spatial)))),
261 : "Cannot use a spacetime index for a spatial index. e.g. "
262 : "Cannot do R(ti::a), where R's index is spatial, because "
263 : "ti::a denotes a generic spacetime index.");
264 : return tenex::contract(TE<tmpl::list<TensorIndices...>>{*this});
265 : }
266 : /// @}
267 :
268 : /// @{
269 : /// Return i'th component of storage vector
270 1 : constexpr reference operator[](const size_t storage_index) {
271 : return gsl::at(data_, storage_index);
272 : }
273 1 : constexpr const_reference operator[](const size_t storage_index) const {
274 : return gsl::at(data_, storage_index);
275 : }
276 : /// @}
277 :
278 : /// Return the number of independent components of the Tensor
279 : ///
280 : /// \details
281 : /// Returns the number of independent components of the Tensor taking into
282 : /// account symmetries. For example, let \f$T_{ab}\f$ be a n-dimensional
283 : /// rank-2 symmetric tensor, then the number of independent components is
284 : /// \f$n(n+1)/2\f$.
285 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t size() {
286 : return structure::size();
287 : }
288 :
289 : /// Returns the rank of the Tensor
290 : ///
291 : /// \details
292 : /// The rank of a tensor is the number of indices it has. For example, the
293 : /// tensor \f$v^a\f$ is rank-1, the tensor \f$\phi\f$ is rank-0, and the
294 : /// tensor \f$T_{abc}\f$ is rank-3.
295 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t rank() {
296 : return sizeof...(Indices);
297 : }
298 :
299 : /// @{
300 : /// Given an iterator or storage index, get the canonical tensor index.
301 : /// For scalars this is defined to be std::array<int, 1>{{0}}
302 : SPECTRE_ALWAYS_INLINE constexpr std::array<size_t, sizeof...(Indices)>
303 1 : get_tensor_index(const const_iterator& iter) const {
304 : return structure::get_canonical_tensor_index(
305 : static_cast<size_t>(iter - begin()));
306 : }
307 : SPECTRE_ALWAYS_INLINE static constexpr std::array<size_t, sizeof...(Indices)>
308 1 : get_tensor_index(const size_t storage_index) {
309 : return structure::get_canonical_tensor_index(storage_index);
310 : }
311 : /// @}
312 :
313 : /// @{
314 : /// Get the storage index of the tensor index. Should only be used when
315 : /// optimizing code in which computing the storage index is a bottleneck.
316 : template <typename... N>
317 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t get_storage_index(
318 : const N... args) {
319 : return structure::get_storage_index(args...);
320 : }
321 : template <typename I>
322 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t get_storage_index(
323 : const std::array<I, sizeof...(Indices)>& tensor_index) {
324 : return structure::get_storage_index(tensor_index);
325 : }
326 : /// @}
327 :
328 : /// @{
329 : /// Given an iterator or storage index, get the multiplicity of an index
330 : ///
331 : /// \see TensorMetafunctions::compute_multiplicity
332 1 : SPECTRE_ALWAYS_INLINE constexpr size_t multiplicity(
333 : const const_iterator& iter) const {
334 : return structure::multiplicity(static_cast<size_t>(iter - begin()));
335 : }
336 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t multiplicity(
337 : const size_t storage_index) {
338 : return structure::multiplicity(storage_index);
339 : }
340 : /// @}
341 :
342 : /// @{
343 : /// Get dimensionality of i'th tensor index
344 : ///
345 : /// \snippet Test_Tensor.cpp index_dim
346 : /// \see ::index_dim
347 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t index_dim(const size_t i) {
348 : static_assert(sizeof...(Indices),
349 : "A scalar does not have any indices from which you can "
350 : "retrieve the dimensionality.");
351 : return gsl::at(structure::dims(), i);
352 : }
353 : /// @}
354 :
355 : /// @{
356 : /// Return an array corresponding to the ::Symmetry of the Tensor
357 : SPECTRE_ALWAYS_INLINE static constexpr std::array<int, sizeof...(Indices)>
358 1 : symmetries() {
359 : return structure::symmetries();
360 : }
361 : /// @}
362 :
363 : /// @{
364 : /// Return array of the ::IndexType's (spatial or spacetime)
365 : SPECTRE_ALWAYS_INLINE static constexpr std::array<IndexType,
366 : sizeof...(Indices)>
367 1 : index_types() {
368 : return structure::index_types();
369 : }
370 : /// @}
371 :
372 : /// @{
373 : /// Return array of dimensionality of each index
374 : ///
375 : /// \snippet Test_Tensor.cpp index_dim
376 : /// \see index_dim ::index_dim
377 : SPECTRE_ALWAYS_INLINE static constexpr std::array<size_t, sizeof...(Indices)>
378 1 : index_dims() {
379 : return structure::dims();
380 : }
381 : /// @}
382 :
383 : /// @{
384 : /// Return array of the valence of each index (::UpLo)
385 : SPECTRE_ALWAYS_INLINE static constexpr std::array<UpLo, sizeof...(Indices)>
386 1 : index_valences() {
387 : return structure::index_valences();
388 : }
389 : /// @}
390 :
391 : /// @{
392 : /// Returns std::tuple of the ::Frame of each index
393 1 : SPECTRE_ALWAYS_INLINE static constexpr auto index_frames() {
394 : return Tensor_detail::Structure<Symm, Indices...>::index_frames();
395 : }
396 : /// @}
397 :
398 : /// @{
399 : /// \brief Given a tensor index, get the canonical label associated with the
400 : /// canonical \ref SpacetimeIndex "TensorIndexType"
401 : ///
402 : /// \param tensor_index The index of the tensor component to label
403 : /// \param axis_labels The labels for the indices. Defaults to "t", "x", "y"
404 : /// and "z" for spacetime indices and "x", "y" and "z" for spatial indices.
405 : /// Note that a tensor can have indices of different types, so we specify
406 : /// labels for each index individually.
407 : template <typename T = int>
408 1 : static std::string component_name(
409 : const std::array<T, rank()>& tensor_index = std::array<T, rank()>{},
410 : const std::array<std::string, rank()>& axis_labels =
411 : make_array<rank()>(std::string(""))) {
412 : return structure::component_name(tensor_index, axis_labels);
413 : }
414 : /// @}
415 :
416 : /// @{
417 : /// \brief Suffix to append to the tensor name that indicates the component
418 : ///
419 : /// The suffix is empty for scalars, otherwise it is an underscore followed by
420 : /// the `Tensor::component_name` of either the `tensor_index` or the canonical
421 : /// tensor index obtained from the `storage_index`. Use `axis_labels` to
422 : /// overwrite the default labels for each component (see
423 : /// `Tensor::component_name`).
424 : ///
425 : /// An example use case for the suffix is to label tensor components in
426 : /// data files.
427 : ///
428 : /// \see Tensor::component_name
429 : template <typename IndexType = int>
430 1 : static std::string component_suffix(
431 : const std::array<IndexType, rank()>& tensor_index =
432 : std::array<IndexType, rank()>{},
433 : const std::array<std::string, rank()>& axis_labels =
434 : make_array<rank()>(std::string(""))) {
435 : return rank() == 0 ? "" : "_" + component_name(tensor_index, axis_labels);
436 : }
437 :
438 1 : static std::string component_suffix(
439 : const size_t storage_index,
440 : const std::array<std::string, rank()>& axis_labels) {
441 : return component_suffix(get_tensor_index(storage_index), axis_labels);
442 : }
443 :
444 1 : static std::string component_suffix(size_t storage_index);
445 : /// @}
446 :
447 : /// Copy tensor data into an `std::vector<X>` along with the
448 : /// component names into a `std::vector<std::string>`
449 : /// \requires `std::is_same<X, DataVector>::%value` is true
450 1 : std::pair<std::vector<std::string>, std::vector<X>> get_vector_of_data()
451 : const;
452 :
453 : /// \cond HIDDEN_SYMBOLS
454 : /// Serialization function used by Charm++
455 : // NOLINTNEXTLINE(google-runtime-references)
456 : void pup(PUP::er& p) { p | data_; }
457 : /// \endcond
458 :
459 : private:
460 : // clang-tidy: redundant declaration
461 : /// \cond
462 : template <int I, class... Ts>
463 : friend SPECTRE_ALWAYS_INLINE constexpr size_t index_dim( // NOLINT
464 : const Tensor<Ts...>& /*t*/);
465 : /// \endcond
466 :
467 0 : storage_type data_;
468 : };
469 :
470 : // ================================================================
471 : // Template Definitions - Variadic templates must be in header
472 : // ================================================================
473 :
474 : template <typename X, typename Symm, template <typename...> class IndexList,
475 : typename... Indices>
476 : // Implementation note: we explicitly prevent inlining and mark the function
477 : // as used so that GDB's pretty printing facilities have access to this
478 : // function.
479 0 : __attribute__((noinline)) __attribute__((used)) std::string
480 : Tensor<X, Symm, IndexList<Indices...>>::component_suffix(
481 : const size_t storage_index) {
482 : return component_suffix(get_tensor_index(storage_index),
483 : make_array<rank()>(std::string("")));
484 : }
485 :
486 : template <typename X, typename Symm, template <typename...> class IndexList,
487 : typename... Indices>
488 : constexpr Tensor<X, Symm, IndexList<Indices...>>::Tensor(storage_type data)
489 : requires(sizeof...(Indices) <= 1)
490 : : data_(std::move(data)) {}
491 :
492 : template <typename X, typename Symm, template <typename...> class IndexList,
493 : typename... Indices>
494 : template <typename Arg0, typename... Args>
495 : constexpr Tensor<X, Symm, IndexList<Indices...>>::Tensor(Arg0&& arg0,
496 : Args&&... args)
497 : // NOLINTNEXTLINE(readability-simplify-boolean-expr)
498 : requires(not(sizeof...(Args) == 0 and
499 : std::same_as<std::decay_t<Arg0>, Tensor>) and
500 : std::constructible_from<X, Arg0, Args...>)
501 : : data_(make_array<size(), X>(std::forward<Arg0>(arg0),
502 : std::forward<Args>(args)...)) {}
503 :
504 : template <typename X, typename Symm, template <typename...> class IndexList,
505 : typename... Indices>
506 : std::pair<std::vector<std::string>, std::vector<X>>
507 : Tensor<X, Symm, IndexList<Indices...>>::get_vector_of_data() const {
508 : std::vector<value_type> serialized_tensor(size());
509 : std::vector<std::string> component_names(size());
510 : for (size_t i = 0; i < data_.size(); ++i) {
511 : component_names[i] = component_name(get_tensor_index(i));
512 : serialized_tensor[i] = gsl::at(data_, i);
513 : }
514 : return std::make_pair(component_names, serialized_tensor);
515 : }
516 :
517 : template <int... N, typename... Args>
518 0 : SPECTRE_ALWAYS_INLINE constexpr typename Tensor<Args...>::reference get(
519 : Tensor<Args...>& t) {
520 : static_assert(Tensor<Args...>::rank() == sizeof...(N),
521 : "the number of tensor indices specified must match the rank "
522 : "of the tensor");
523 : return gsl::at(
524 : t.data_, Tensor<Args...>::structure::template get_storage_index<N...>());
525 : }
526 :
527 : template <int... N, typename... Args>
528 0 : SPECTRE_ALWAYS_INLINE constexpr typename Tensor<Args...>::const_reference get(
529 : const Tensor<Args...>& t) {
530 : static_assert(Tensor<Args...>::rank() == sizeof...(N),
531 : "the number of tensor indices specified must match the rank "
532 : "of the tensor");
533 : return gsl::at(
534 : t.data_, Tensor<Args...>::structure::template get_storage_index<N...>());
535 : }
536 :
537 : template <typename X, typename Symm, template <typename...> class IndexList,
538 : typename... Indices>
539 0 : bool operator==(const Tensor<X, Symm, IndexList<Indices...>>& lhs,
540 : const Tensor<X, Symm, IndexList<Indices...>>& rhs) {
541 : return std::equal(lhs.begin(), lhs.end(), rhs.begin());
542 : }
543 : template <typename X, typename Symm, template <typename...> class IndexList,
544 : typename... Indices>
545 0 : bool operator!=(const Tensor<X, Symm, IndexList<Indices...>>& lhs,
546 : const Tensor<X, Symm, IndexList<Indices...>>& rhs) {
547 : return not(lhs == rhs);
548 : }
549 :
550 : /// \ingroup TensorGroup
551 : /// Get dimensionality of i'th tensor index
552 : ///
553 : /// \snippet Test_Tensor.cpp index_dim
554 : template <int I, class... Ts>
555 1 : SPECTRE_ALWAYS_INLINE constexpr size_t index_dim(const Tensor<Ts...>& /*t*/) {
556 : return Tensor<Ts...>::structure::template dim<I>();
557 : }
558 :
559 : // We place the stream operators in the header file so they do not need to be
560 : // explicitly instantiated.
561 : template <typename X, typename Symm, template <typename...> class IndexList,
562 : typename... Indices>
563 0 : std::ostream& operator<<(std::ostream& os,
564 : const Tensor<X, Symm, IndexList<Indices...>>& x) {
565 : static_assert(tt::is_streamable_v<decltype(os), X>,
566 : "operator<< is not defined for the type you are trying to "
567 : "stream in Tensor");
568 : for (size_t i = 0; i < x.size() - 1; ++i) {
569 : os << "T" << x.get_tensor_index(i) << "=" << x[i] << "\n";
570 : }
571 : size_t i = x.size() - 1;
572 : os << "T" << x.get_tensor_index(i) << "=" << x[i];
573 : return os;
574 : }
575 :
576 : namespace MakeWithValueImpls {
577 : template <typename T, typename... Structure>
578 0 : struct NumberOfPoints<Tensor<T, Structure...>> {
579 : static SPECTRE_ALWAYS_INLINE size_t
580 0 : apply(const Tensor<T, Structure...>& input) {
581 : return number_of_points(*input.begin());
582 : }
583 : };
584 :
585 : template <typename T, typename... Structure>
586 0 : struct MakeWithSize<Tensor<T, Structure...>> {
587 : template <typename U>
588 0 : static SPECTRE_ALWAYS_INLINE Tensor<T, Structure...> apply(const size_t size,
589 : const U value) {
590 : return Tensor<T, Structure...>(make_with_value<T>(size, value));
591 : }
592 : };
593 :
594 : template <typename... Structure, typename T>
595 0 : struct MakeWithValueImpl<Tensor<double, Structure...>, T> {
596 0 : static SPECTRE_ALWAYS_INLINE Tensor<double, Structure...> apply(
597 : const T& /*input*/, const double value) {
598 : return Tensor<double, Structure...>(value);
599 : }
600 : };
601 :
602 : template <typename... Structure, typename T>
603 : struct MakeWithValueImpl<Tensor<std::complex<double>, Structure...>, T> {
604 : static SPECTRE_ALWAYS_INLINE Tensor<std::complex<double>, Structure...> apply(
605 : const T& /*input*/, const std::complex<double> value) {
606 : return Tensor<std::complex<double>, Structure...>(value);
607 : }
608 : };
609 :
610 : #ifdef SPECTRE_AUTODIFF
611 : template <typename... Structure, typename T>
612 0 : struct MakeWithValueImpl<
613 : Tensor<autodiff::HigherOrderDual<2, double>, Structure...>, T> {
614 : static SPECTRE_ALWAYS_INLINE
615 : Tensor<autodiff::HigherOrderDual<2, double>, Structure...>
616 0 : apply(const T& /*input*/, const double value) {
617 : return Tensor<autodiff::HigherOrderDual<2, double>, Structure...>(value);
618 : }
619 : };
620 :
621 : template <typename... Structure, typename T>
622 0 : struct MakeWithValueImpl<
623 : Tensor<autodiff::HigherOrderDual<2, simd::batch<double>>, Structure...>,
624 : T> {
625 : static SPECTRE_ALWAYS_INLINE
626 : Tensor<autodiff::HigherOrderDual<2, simd::batch<double>>, Structure...>
627 0 : apply(const T& /*input*/, const double value) {
628 : return Tensor<autodiff::HigherOrderDual<2, simd::batch<double>>,
629 : Structure...>(simd::batch<double>::broadcast(value));
630 : }
631 : };
632 : #endif // SPECTRE_AUTODIFF
633 : } // namespace MakeWithValueImpls
634 :
635 : template <typename T, typename... Structure>
636 0 : struct SetNumberOfGridPointsImpls::SetNumberOfGridPointsImpl<
637 : Tensor<T, Structure...>> {
638 0 : static constexpr bool is_trivial = SetNumberOfGridPointsImpl<T>::is_trivial;
639 0 : static SPECTRE_ALWAYS_INLINE void apply(
640 : const gsl::not_null<Tensor<T, Structure...>*> result, const size_t size) {
641 : for (auto& component : *result) {
642 : set_number_of_grid_points(make_not_null(&component), size);
643 : }
644 : }
645 : };
|
