1 #include "physics/Integrals.h"
2 #include "post/Analysis.h"
3 #include "quantities/Storage.h"
4 #include "system/Factory.h"
5
6 NAMESPACE_SPH_BEGIN
7
8
evaluate(const Storage & storage) const9 Float TotalMass::evaluate(const Storage& storage) const {
10 Float total(0._f);
11 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
12 SPH_ASSERT(!m.empty());
13 for (Size i = 0; i < m.size(); ++i) {
14 total += m[i];
15 }
16 SPH_ASSERT(isReal(total));
17 return total;
18 }
19
TotalMomentum(const Float omega)20 TotalMomentum::TotalMomentum(const Float omega)
21 : omega(0._f, 0._f, omega) {}
22
23
evaluate(const Storage & storage) const24 Vector TotalMomentum::evaluate(const Storage& storage) const {
25 BasicVector<double> total(0.); // compute in double precision to avoid round-off error during accumulation
26 ArrayView<const Vector> r, v, dv;
27 tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
28 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
29 SPH_ASSERT(!v.empty());
30 for (Size i = 0; i < v.size(); ++i) {
31 total += vectorCast<double>(m[i] * (v[i] + cross(omega, r[i])));
32 }
33 SPH_ASSERT(isReal(total));
34 return vectorCast<Float>(total);
35 }
36
TotalAngularMomentum(const Float omega)37 TotalAngularMomentum::TotalAngularMomentum(const Float omega)
38 : frameFrequency(0._f, 0._f, omega) {}
39
evaluate(const Storage & storage) const40 Vector TotalAngularMomentum::evaluate(const Storage& storage) const {
41 BasicVector<double> total(0.);
42 ArrayView<const Vector> r, v, dv;
43 tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
44 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
45
46 // angular momentum with respect to origin
47 SPH_ASSERT(!v.empty());
48 for (Size i = 0; i < v.size(); ++i) {
49 total += vectorCast<double>(m[i] * cross(r[i], (v[i] + cross(frameFrequency, r[i]))));
50 }
51
52 // local angular momentum
53 /// \todo consolidate the angular momentum here - always in local frame? introduce physical quantity?
54 /*if (storage.has(QuantityId::ANGULAR_VELOCITY)) {
55 ArrayView<const Vector> omega = storage.getValue<Vector>(QuantityId::ANGULAR_VELOCITY);
56 ArrayView<const SymmetricTensor> I = storage.getValue<SymmetricTensor>(QuantityId::MOMENT_OF_INERTIA);
57 for (Size i = 0; i < v.size(); ++i) {
58 total += I[i] * omega[i];
59 }
60 }*/
61
62 SPH_ASSERT(isReal(total));
63 return vectorCast<Float>(total);
64 }
65
TotalEnergy(const Float omega)66 TotalEnergy::TotalEnergy(const Float omega)
67 : omega(0._f, 0._f, omega) {}
68
evaluate(const Storage & storage) const69 Float TotalEnergy::evaluate(const Storage& storage) const {
70 double total = 0.;
71 // add kinetic energy
72 ArrayView<const Vector> r, v, dv;
73 tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
74 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
75 for (Size i = 0; i < v.size(); ++i) {
76 SPH_ASSERT(!v.empty());
77 total += 0.5 * m[i] * getSqrLength(v[i]);
78 }
79
80 if (storage.has(QuantityId::ENERGY)) {
81 // add internal energy
82 ArrayView<const Float> u = storage.getValue<Float>(QuantityId::ENERGY);
83 for (Size i = 0; i < v.size(); ++i) {
84 SPH_ASSERT(!v.empty());
85 total += m[i] * u[i];
86 }
87 }
88
89 SPH_ASSERT(isReal(total));
90 return Float(total);
91 }
92
TotalKineticEnergy(const Float omega)93 TotalKineticEnergy::TotalKineticEnergy(const Float omega)
94 : omega(0._f, 0._f, omega) {}
95
evaluate(const Storage & storage) const96 Float TotalKineticEnergy::evaluate(const Storage& storage) const {
97 double total = 0.;
98 ArrayView<const Vector> r, v, dv;
99 tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
100 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
101 for (Size i = 0; i < v.size(); ++i) {
102 SPH_ASSERT(!v.empty());
103 total += 0.5 * m[i] * getSqrLength(v[i]);
104 }
105 SPH_ASSERT(isReal(total));
106 return Float(total);
107 }
108
evaluate(const Storage & storage) const109 Float TotalInternalEnergy::evaluate(const Storage& storage) const {
110 double total = 0.;
111 if (!storage.has(QuantityId::ENERGY)) {
112 return Float(total);
113 }
114 ArrayView<const Float> u = storage.getValue<Float>(QuantityId::ENERGY);
115 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
116 SPH_ASSERT(!m.empty());
117 for (Size i = 0; i < m.size(); ++i) {
118 total += double(m[i] * u[i]);
119 }
120 SPH_ASSERT(isReal(total));
121 return Float(total);
122 }
123
CenterOfMass(const Optional<Size> bodyId)124 CenterOfMass::CenterOfMass(const Optional<Size> bodyId)
125 : bodyId(bodyId) {}
126
evaluate(const Storage & storage) const127 Vector CenterOfMass::evaluate(const Storage& storage) const {
128 Vector com(0._f);
129 Float totalMass = 0._f;
130 ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
131 ArrayView<const Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
132 auto accumulate = [&](const Size i) {
133 totalMass += m[i];
134 com += m[i] * r[i];
135 };
136
137 if (bodyId) {
138 ArrayView<const Size> ids = storage.getValue<Size>(QuantityId::FLAG);
139 for (Size i = 0; i < r.size(); ++i) {
140 if (ids[i] == bodyId.value()) {
141 accumulate(i);
142 }
143 }
144 } else {
145 for (Size i = 0; i < r.size(); ++i) {
146 accumulate(i);
147 }
148 }
149 return com / totalMass;
150 }
151
QuantityMeans(const QuantityId id,const Optional<Size> bodyId)152 QuantityMeans::QuantityMeans(const QuantityId id, const Optional<Size> bodyId)
153 : quantity(id)
154 , bodyId(bodyId) {}
155
QuantityMeans(AutoPtr<IUserQuantity> && func,const Optional<Size> bodyId)156 QuantityMeans::QuantityMeans(AutoPtr<IUserQuantity>&& func, const Optional<Size> bodyId)
157 : quantity(std::move(func))
158 , bodyId(bodyId) {}
159
evaluate(const Storage & storage) const160 MinMaxMean QuantityMeans::evaluate(const Storage& storage) const {
161 MinMaxMean means;
162 auto accumulate = [&](const auto& getter) {
163 const Size size = storage.getParticleCnt();
164 if (bodyId) {
165 ArrayView<const Size> ids = storage.getValue<Size>(QuantityId::FLAG);
166 for (Size i = 0; i < size; ++i) {
167 if (ids[i] == bodyId.value()) {
168 means.accumulate(getter(i));
169 }
170 }
171 } else {
172 for (Size i = 0; i < size; ++i) {
173 means.accumulate(getter(i));
174 }
175 }
176 };
177 if (quantity.getTypeIdx() == 0) {
178 const QuantityId id = quantity.get<QuantityId>();
179 SPH_ASSERT(storage.has(id));
180 ArrayView<const Float> q = storage.getValue<Float>(id);
181 accumulate([&](const Size i) { return q[i]; });
182 } else {
183 SPH_ASSERT(quantity.getTypeIdx() == 1);
184 const AutoPtr<IUserQuantity>& func = quantity;
185 func->initialize(storage);
186 accumulate([&](const Size i) { return func->evaluate(i); });
187 }
188 return means;
189 }
190
QuantityValue(const QuantityId id,const Size particleIdx)191 QuantityValue::QuantityValue(const QuantityId id, const Size particleIdx)
192 : id(id)
193 , idx(particleIdx) {}
194
evaluate(const Storage & storage) const195 Float QuantityValue::evaluate(const Storage& storage) const {
196 ArrayView<const Float> values = storage.getValue<Float>(id);
197 return values[idx];
198 }
199
200 NAMESPACE_SPH_END
201