Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
grav_moments.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
18#include "shambase/aliases_int.hpp"
19#include "shambase/type_traits.hpp"
20#include "shammath/symtensor_collections.hpp"
21
22namespace shamphys {
23
25 template<class T, u32 low_order, u32 high_order>
26 inline shammath::SymTensorCollection<T, low_order, high_order> get_M_mat(
27 const shammath::SymTensorCollection<T, low_order, high_order> &D,
28 const shammath::SymTensorCollection<T, low_order, high_order> &Q) {
29
30 using namespace shammath;
31
32 if constexpr (low_order == 0 && high_order == 5) {
33
34 const T &TD0 = D.t0;
35 const SymTensor3d_1<T> &TD1 = D.t1;
36 const SymTensor3d_2<T> &TD2 = D.t2;
37 const SymTensor3d_3<T> &TD3 = D.t3;
38 const SymTensor3d_4<T> &TD4 = D.t4;
39 const SymTensor3d_5<T> &TD5 = D.t5;
40
41 const T &TQ0 = Q.t0;
42 const SymTensor3d_1<T> &TQ1 = Q.t1;
43 const SymTensor3d_2<T> &TQ2 = Q.t2;
44 const SymTensor3d_3<T> &TQ3 = Q.t3;
45 const SymTensor3d_4<T> &TQ4 = Q.t4;
46 const SymTensor3d_5<T> &TQ5 = Q.t5;
47
48 auto M_0 = (TD0 * TQ0) + (TD1 * TQ1) + ((TD2 * TQ2)) * T(1 / 2.)
49 + ((TD3 * TQ3)) * T(1 / 6.) + ((TD4 * TQ4)) * T(1 / 24.)
50 + ((TD5 * TQ5)) * T(1 / 120.);
51 auto M_1 = T(-1.) * (TD1 * TQ0) - (TD2 * TQ1) - ((TD3 * TQ2)) * T(1 / 2.)
52 - ((TD4 * TQ3)) * T(1 / 6.) - ((TD5 * TQ4)) * T(1 / 24.);
53 auto M_2
54 = (TD2 * TQ0) + (TD3 * TQ1) + ((TD4 * TQ2)) * T(1 / 2.) + ((TD5 * TQ3)) * T(1 / 6.);
55 auto M_3 = T(-1.) * (TD3 * TQ0) - (TD4 * TQ1) - ((TD5 * TQ2)) * T(1 / 2.);
56 auto M_4 = (TD4 * TQ0) + (TD5 * TQ1);
57 auto M_5 = T(-1.) * (TD5 * TQ0);
58
59 return SymTensorCollection<T, 0, 5>{M_0, M_1, M_2, M_3, M_4, M_5};
60 } else if constexpr (low_order == 0 && high_order == 4) {
61
62 const T &TD0 = D.t0;
63 const SymTensor3d_1<T> &TD1 = D.t1;
64 const SymTensor3d_2<T> &TD2 = D.t2;
65 const SymTensor3d_3<T> &TD3 = D.t3;
66 const SymTensor3d_4<T> &TD4 = D.t4;
67
68 const T &TQ0 = Q.t0;
69 const SymTensor3d_1<T> &TQ1 = Q.t1;
70 const SymTensor3d_2<T> &TQ2 = Q.t2;
71 const SymTensor3d_3<T> &TQ3 = Q.t3;
72 const SymTensor3d_4<T> &TQ4 = Q.t4;
73
74 auto M_0 = (TD0 * TQ0) + (TD1 * TQ1) + ((TD2 * TQ2)) * T(1 / 2.)
75 + ((TD3 * TQ3)) * T(1 / 6.) + ((TD4 * TQ4)) * T(1 / 24.);
76 auto M_1 = T(-1.) * (TD1 * TQ0) - (TD2 * TQ1) - ((TD3 * TQ2)) * T(1 / 2.)
77 - ((TD4 * TQ3)) * T(1 / 6.);
78 auto M_2 = (TD2 * TQ0) + (TD3 * TQ1) + ((TD4 * TQ2)) * T(1 / 2.);
79 auto M_3 = T(-1.) * (TD3 * TQ0) - (TD4 * TQ1);
80 auto M_4 = (TD4 * TQ0);
81
82 return SymTensorCollection<T, 0, 4>{M_0, M_1, M_2, M_3, M_4};
83 } else if constexpr (low_order == 0 && high_order == 3) {
84
85 const T &TD0 = D.t0;
86 const SymTensor3d_1<T> &TD1 = D.t1;
87 const SymTensor3d_2<T> &TD2 = D.t2;
88 const SymTensor3d_3<T> &TD3 = D.t3;
89
90 const T &TQ0 = Q.t0;
91 const SymTensor3d_1<T> &TQ1 = Q.t1;
92 const SymTensor3d_2<T> &TQ2 = Q.t2;
93 const SymTensor3d_3<T> &TQ3 = Q.t3;
94
95 auto M_0
96 = (TD0 * TQ0) + (TD1 * TQ1) + ((TD2 * TQ2)) * T(1 / 2.) + ((TD3 * TQ3)) * T(1 / 6.);
97 auto M_1 = T(-1.) * (TD1 * TQ0) - (TD2 * TQ1) - ((TD3 * TQ2)) * T(1 / 2.);
98 auto M_2 = (TD2 * TQ0) + (TD3 * TQ1);
99 auto M_3 = T(-1.) * (TD3 * TQ0);
100
101 return SymTensorCollection<T, 0, 3>{M_0, M_1, M_2, M_3};
102 } else if constexpr (low_order == 0 && high_order == 2) {
103
104 const T &TD0 = D.t0;
105 const SymTensor3d_1<T> &TD1 = D.t1;
106 const SymTensor3d_2<T> &TD2 = D.t2;
107
108 const T &TQ0 = Q.t0;
109 const SymTensor3d_1<T> &TQ1 = Q.t1;
110 const SymTensor3d_2<T> &TQ2 = Q.t2;
111
112 auto M_0 = (TD0 * TQ0) + (TD1 * TQ1) + ((TD2 * TQ2)) * T(1 / 2.);
113 auto M_1 = T(-1.) * (TD1 * TQ0) - (TD2 * TQ1);
114 auto M_2 = TD2 * TQ0;
115
116 return SymTensorCollection<T, 0, 2>{M_0, M_1, M_2};
117 } else if constexpr (low_order == 0 && high_order == 1) {
118
119 const T &TD0 = D.t0;
120 const SymTensor3d_1<T> &TD1 = D.t1;
121
122 const T &TQ0 = Q.t0;
123 const SymTensor3d_1<T> &TQ1 = Q.t1;
124
125 auto M_0 = (TD0 * TQ0) + (TD1 * TQ1);
126 auto M_1 = T(-1.) * (TD1 * TQ0);
127
128 return SymTensorCollection<T, 0, 1>{M_0, M_1};
129 } else if constexpr (low_order == 0 && high_order == 0) {
130
131 const T &TD0 = D.t0;
132
133 const T &TQ0 = Q.t0;
134
135 auto M_0 = TD0 * TQ0;
136
137 return SymTensorCollection<T, 0, 0>{M_0};
138 } else {
139 static_assert(shambase::always_false_v<T>, "This combination of orders is not valid");
140 }
141 }
142
145 template<class T, u32 high_order>
146 inline shammath::SymTensorCollection<T, 1, high_order + 1> get_dM_mat(
147 const shammath::SymTensorCollection<T, 1, high_order + 1> &D,
148 const shammath::SymTensorCollection<T, 0, high_order> &Q) {
149
150 using namespace shammath;
151
152 if constexpr (high_order == 4) {
153 // T & TD0 = D.t0;
154 const SymTensor3d_1<T> &TD1 = D.t1;
155 const SymTensor3d_2<T> &TD2 = D.t2;
156 const SymTensor3d_3<T> &TD3 = D.t3;
157 const SymTensor3d_4<T> &TD4 = D.t4;
158 const SymTensor3d_5<T> &TD5 = D.t5;
159
160 const T &TQ0 = Q.t0;
161 const SymTensor3d_1<T> &TQ1 = Q.t1;
162 const SymTensor3d_2<T> &TQ2 = Q.t2;
163 const SymTensor3d_3<T> &TQ3 = Q.t3;
164 const SymTensor3d_4<T> &TQ4 = Q.t4;
165 // SymTensor3d_5<T> & TQ5 = Q.t5;
166
167 constexpr T _1i2 = T(1. / 2.);
168 constexpr T _1i6 = T(1. / 6.);
169 constexpr T _1i24 = T(1. / 24.);
170
171 auto M_1 = TD1 * TQ0 + TD2 * TQ1 + (TD3 * TQ2) * _1i2 + (TD4 * TQ3) * _1i6
172 + (TD5 * TQ4) * _1i24;
173 auto M_2 = (T(-1.) * TD2 * TQ0) - TD3 * TQ1 - (TD4 * TQ2) * _1i2 - (TD5 * TQ3) * _1i6;
174 auto M_3 = TD3 * TQ0 + TD4 * TQ1 + (TD5 * TQ2) * _1i2;
175 auto M_4 = (T(-1.) * (TD4 * TQ0)) - TD5 * TQ1;
176 auto M_5 = TD5 * TQ0;
177
178 return SymTensorCollection<T, 1, 5>{M_1, M_2, M_3, M_4, M_5};
179 } else if constexpr (high_order == 3) {
180
181 // T & TD0 = D.t0;
182 const SymTensor3d_1<T> &TD1 = D.t1;
183 const SymTensor3d_2<T> &TD2 = D.t2;
184 const SymTensor3d_3<T> &TD3 = D.t3;
185 const SymTensor3d_4<T> &TD4 = D.t4;
186
187 const T &TQ0 = Q.t0;
188 const SymTensor3d_1<T> &TQ1 = Q.t1;
189 const SymTensor3d_2<T> &TQ2 = Q.t2;
190 const SymTensor3d_3<T> &TQ3 = Q.t3;
191
192 auto M_1 = TD1 * TQ0 + TD2 * TQ1 + (TD3 * TQ2) * T(1. / 2.) + (TD4 * TQ3) * T(1. / 6.);
193 auto M_2 = (T(-1.) * TD2 * TQ0) - TD3 * TQ1 - (TD4 * TQ2) * T(1. / 2.);
194 auto M_3 = TD3 * TQ0 + TD4 * TQ1;
195 auto M_4 = T(-1.) * (TD4 * TQ0);
196
197 return SymTensorCollection<T, 1, 4>{M_1, M_2, M_3, M_4};
198 } else if constexpr (high_order == 2) {
199
200 // T & TD0 = D.t0;
201 const SymTensor3d_1<T> &TD1 = D.t1;
202 const SymTensor3d_2<T> &TD2 = D.t2;
203 const SymTensor3d_3<T> &TD3 = D.t3;
204
205 const T &TQ0 = Q.t0;
206 const SymTensor3d_1<T> &TQ1 = Q.t1;
207 const SymTensor3d_2<T> &TQ2 = Q.t2;
208
209 auto M_1 = TD1 * TQ0 + TD2 * TQ1 + (TD3 * TQ2) * (1. / 2.);
210 auto M_2 = (T(-1.) * TD2 * TQ0) - TD3 * TQ1;
211 auto M_3 = TD3 * TQ0;
212
213 return SymTensorCollection<T, 1, 3>{M_1, M_2, M_3};
214 } else if constexpr (high_order == 1) {
215
216 // T & TD0 = D.t0;
217 const SymTensor3d_1<T> &TD1 = D.t1;
218 const SymTensor3d_2<T> &TD2 = D.t2;
219
220 const T &TQ0 = Q.t0;
221 const SymTensor3d_1<T> &TQ1 = Q.t1;
222
223 auto M_1 = TD1 * TQ0 + TD2 * TQ1;
224 auto M_2 = T(-1.) * TD2 * TQ0;
225
226 return SymTensorCollection<T, 1, 2>{M_1, M_2};
227 } else if constexpr (high_order == 0) {
228 // T & TD0 = D.t0;
229 const SymTensor3d_1<T> &TD1 = D.t1;
230
231 const T &TQ0 = Q.t0;
232
233 auto M_1 = TD1 * TQ0;
234
235 return SymTensorCollection<T, 1, 1>{M_1};
236 } else {
237 static_assert(shambase::always_false_v<T>, "This combination of orders is not valid");
238 }
239 }
240} // namespace shamphys
aliases_int.hpp
shamphys::get_dM_mat
shammath::SymTensorCollection< T, 1, high_order+1 > get_dM_mat(const shammath::SymTensorCollection< T, 1, high_order+1 > &D, const shammath::SymTensorCollection< T, 0, high_order > &Q)
Definition grav_moments.hpp:146
shamphys::get_M_mat
shammath::SymTensorCollection< T, low_order, high_order > get_M_mat(const shammath::SymTensorCollection< T, low_order, high_order > &D, const shammath::SymTensorCollection< T, low_order, high_order > &Q)
Contraction of the green function derivatives (D_n) and the multipole moments (Q_n)
Definition grav_moments.hpp:26
shambase::throw_with_loc
void throw_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Throw an exception and append the source location to it.
Definition exception.hpp:132
shammath
namespace for math utility
Definition AABB.hpp:26
type_traits.hpp
Traits for C++ types.
shammath::SymTensor3d_1
Definition symtensors.hpp:23
shammath::SymTensor3d_2
Definition symtensors.hpp:91
shammath::SymTensor3d_3
Definition symtensors.hpp:192
shammath::SymTensor3d_4
Definition symtensors.hpp:354
shammath::SymTensor3d_5
Definition symtensors.hpp:589
symtensor_collections.hpp
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8