18 const Size boundaryThreshold)
19 : scheduler(scheduler)
21 , boundaryThreshold(boundaryThreshold) {
22 kernel = Factory::getKernel<3>(settings);
30 gravity->
build(scheduler, storage);
33 gravity->
evalAll(scheduler, dv, stats);
37 for (
Size i = 0; i < r.
size(); ++i) {
42 finder->
build(scheduler, r);
43 Float maxRadius = 0._f;
44 for (
Size i = 0; i < r.
size(); ++i) {
45 maxRadius =
max(maxRadius, r[i][
H]);
55 for (
Size i = 0; i < r.
size(); ++i) {
66 const bool isBoundary = neighs.
size() < boundaryThreshold;
71 const Size j = n.index;
72 const Float hbar = 0.5_f * (r[i][
H] + r[j][
H]);
74 if (i == j || n.distanceSqr >=
sqr(kernel.
radius() * hbar)) {
80 Aii += m[j] * grad /
sqr(rho[i]);
81 const Vector Aij = m[j] * grad /
sqr(rho[j]);
83 for (
int k = 0; k < 3; ++k) {
84 matrix.insert(3 * i + k, j, Aij[k] - lambda *
weight);
87 for (
int k = 0; k < 3; ++k) {
88 matrix.insert(3 * i + k, j, -lambda *
weight);
95 for (
int k = 0; k < 3; ++k) {
96 matrix.insert(3 * i + k, i, Aii[k] + weightSum * lambda);
97 b[3 * i + k] = dv[i][k];
100 for (
int k = 0; k < 3; ++k) {
101 matrix.insert(3 * i + k, i, 1.f + weightSum * lambda);
117 for (
Size i = 0; i < r.
size(); ++i) {
124 const Size j = n.index;
125 const Float hbar = 0.5_f * (r[i][
H] + r[j][
H]);
127 if (i == j || n.distanceSqr >=
sqr(kernel.
radius() * hbar)) {
133 Aii -= m[j] * lapl /
sqr(rho[i]);
134 const Float Aij = m[j] * lapl /
sqr(rho[j]);
135 matrix.insert(i, j, Aij);
137 divDv -= m[j] / rho[j] *
dot(dv[j] - dv[i], grad);
161 if (neighs.size() < boundaryThreshold) {
163 Aii += 100._f * m[i] /
sqr(rho[i]) * kernel.
value(r[i], r[i]) /
sqr(r[i][
H]);
167 matrix.insert(i, i, Aii);
229 template <
bool Symmetrize>
237 p[i] += m[j] / rho[j] * tr3;
238 s[i] += m[j] / rho[j] * ds;
240 p[j] += m[i] / rho[i] * tr3;
241 s[j] += m[i] / rho[i] * ds;
249 derivatives.
require(makeAuto<DisplacementGradient>(settings));
271 return additional + makeTerm<DisplacementTerm>() + makeTerm<ConstSmoothingLength>();
277 : scheduler(scheduler)
278 , equationSolver(scheduler, settings, getEquations(equations)) {
279 kernel = Factory::getKernel<3>(settings);
281 boundaryThreshold = 18;
291 finder->
build(scheduler, r);
295 equationSolver.
integrate(storage, stats);
306 for (
Size i = 0; i < dv.size(); ++i) {
307 for (
Size j = 0; j < 3; ++j) {
308 const Float x = -rho[i] * dv[i][j];
313 b_avg /= dv.size() * 3;
328 for (
Size i = 0; i < r.
size(); ++i) {
331 for (
Size k = 0; k < neighs.size(); ++k) {
332 const Size j = neighs[k].index;
334 const Vector dr = r[i] - r[j];
344 for (
Size a = 0; a < 3; ++a) {
345 for (
Size b = 0; b < 3; ++b) {
346 matrix.
insert(3 * i + a, 3 * i + b, mij(a, b));
347 matrix.
insert(3 * i + a, 3 * j + b, -mij(a, b));
348 matrix.
insert(3 * j + a, 3 * j + b, mji(a, b));
349 matrix.
insert(3 * j + a, 3 * i + b, -mji(a, b));
372 for (
Size i = 0; i < u.
size(); ++i) {
373 if (neighCnts[i] < boundaryThreshold) {
376 for (
Size j = 0; j < 3; ++j) {
377 u[i][j] = a.
value()[3 * i + j];
383 equationSolver.
integrate(storage, stats);
402 equationSolver.
create(storage, material);
@ UNIQUE
Only a single derivative accumulates to this buffer.
INLINE bool isReal(const AntisymmetricTensor &t)
#define SPH_ASSERT(x,...)
Right-hand side term of SPH equations.
Computes quantities to get equilibrium state.
uint32_t Size
Integral type used to index arrays (by default).
double Float
Precision used withing the code. Use Float instead of float or double where precision is important.
Base class for solvers of gravity.
INLINE T laplacian(const T &value, const Vector &grad, const Vector &dr)
SPH approximation of laplacian, computed from a kernel gradient.
INLINE Float weight(const Vector &r1, const Vector &r2)
SPH-specific implementation of particle materials.
constexpr INLINE T max(const T &f1, const T &f2)
constexpr INLINE T sqr(const T &f) noexcept
Return a squared value.
INLINE T sqrt(const T f)
Return a squared root of a value.
INLINE auto abs(const T &f)
#define INLINE
Macros for conditional compilation based on selected compiler.
#define NAMESPACE_SPH_END
const SuccessTag SUCCESS
Global constant for successful outcome.
INLINE Outcome makeFailed(TArgs &&... args)
Constructs failed object with error message.
@ DISPLACEMENT
Displacement vector, a solution of the stress analysis.
@ DEVIATORIC_STRESS
Deviatoric stress tensor, always a traceless tensor.
@ PRESSURE
Pressure, affected by yielding and fragmentation model, always a scalar quantity.
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
@ ENERGY
Specific internal energy, always a scalar quantity.
@ DENSITY
Density, always a scalar quantity.
@ MASS
Paricles masses, always a scalar quantity.
@ NEIGHBOUR_CNT
Number of neighbouring particles (in radius h * kernel.radius)
@ ZERO
Quantity without derivatives, or "zero order" of quantity.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
INLINE SymmetricTensor symmetricOuter(const Vector &v1, const Vector &v2)
SYMMETRIZED outer product of two vectors.
INLINE Float getSqrLength(const Vector &v)
BasicVector< Float > Vector
INLINE float dot(const BasicVector< float > &v1, const BasicVector< float > &v2)
Make sure the vector is trivially constructible and destructible, needed for fast initialization of a...
INLINE Vector getNormalized(const Vector &v)
Storage for accumulating derivatives.
Array< TValue > & getBuffer(const QuantityId id, const OrderEnum order)
Returns the buffer of given quantity and given order.
void insert(const QuantityId id, const OrderEnum order, const BufferSource source)
Creates a new storage with given ID.
Object providing safe access to continuous memory of data.
INLINE TCounter size() const
void clear()
Removes all elements from the array, but does NOT release the memory.
INLINE TCounter size() const noexcept
Container holding derivatives and the storage they accumulate to.
virtual void require(AutoPtr< IDerivative > &&derivative)
Adds derivative if not already present.
Helper template for derivatives that define both the symmetrized and asymmetric variant.
Derivative computing components of stress tensor from known displacement vectors.
DisplacementGradient(const RunSettings &settings)
INLINE void additionalCreate(Accumulated &results)
INLINE void additionalInitialize(const Storage &input, Accumulated &results)
INLINE void eval(const Size i, const Size j, const Vector &grad)
INLINE bool additionalEquals(const DisplacementGradient &UNUSED(other)) const
virtual void initialize(IScheduler &UNUSED(scheduler), Storage &UNUSED(storage), const Float UNUSED(t)) override
virtual void finalize(IScheduler &UNUSED(scheduler), Storage &UNUSED(storage), const Float UNUSED(t)) override
virtual void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings) override
Sets derivatives required by this term.
virtual void create(Storage &storage, IMaterial &UNUSED(material)) const override
Material holding equation of state.
const IEos & getEos() const
Returns the equation of state.
Container holding equation terms.
EquilibriumEnergySolver(IScheduler &scheduler, const RunSettings &settings, AutoPtr< IGravity > &&gravity, const Size boundaryThreshold)
~EquilibriumEnergySolver()
Outcome solve(Storage &storage, Statistics &stats)
EquilibriumStressSolver(IScheduler &scheduler, const RunSettings &settings, const EquationHolder &equations)
Constructs the solver.
void create(Storage &storage, IMaterial &material)
Creates all the necessary quantities in the storage.
~EquilibriumStressSolver()
Outcome solve(Storage &storage, Statistics &stats)
Computed pressure and deviatoric stress are placed into the storage.
Wrapper of type that either contains a value of given type, or an error message.
const Error & error() const
Returns the error message.
Type & value()
Returns the reference to expected value.
void build(IScheduler &scheduler, ArrayView< const Vector > points)
Constructs the struct with an array of vectors.
virtual Size findAll(const Size index, const Float radius, Array< NeighbourRecord > &neighbours) const =0
Finds all neighbours within given radius from the point given by index.
Base class for equations of state.
virtual Float getInternalEnergy(const Float rho, const Float p) const =0
Inverted function; computes specific internal energy u from given density rho and pressure p.
Represents a term or terms appearing in evolutionary equations.
virtual void build(IScheduler &scheduler, const Storage &storage)=0
Builds the accelerating structure.
virtual void evalAll(IScheduler &scheduler, ArrayView< Vector > dv, Statistics &stats) const =0
Evaluates the gravitational acceleration concurrently.
Material settings and functions specific for one material.
INLINE TValue getParam(const BodySettingsId paramIdx) const
Returns a parameter associated with given particle.
Interface that allows unified implementation of sequential and parallelized versions of algorithms.
void build(IScheduler &scheduler, ArrayView< const Vector > points, Flags< FinderFlag > flags=FinderFlag::MAKE_RANK)
Constructs the struct with an array of vectors.
virtual Size findLowerRank(const Size index, const Float radius, Array< NeighbourRecord > &neighbours) const =0
Finds all points within radius that have a lower rank in smoothing length.
Non-owning wrapper of a material and particles with this material.
INLINE IMaterial & material()
Returns the reference to the material of the particles.
INLINE IndexSequence sequence()
Returns iterable index sequence.
Sparse representation of matrix of arbitrary dimension.
void insert(const Size i, const Size j, const Float value)
Adds a values to given element of the matrix.
void resize(const Size rows, const Size cols)
Changes the size of the matrix, removing all previous entries.
Expected< Array< Float > > solve(const Array< Float > &values, const Solver solver, const Float tolerance=0.)
@ BICGSTAB
Stabilized bi-conjugate gradient method, can be used for any square matrix.
@ LSCG
Least-square conjugate gradient, can be used for any matrix.
Object holding various statistics about current run.
Container storing all quantities used within the simulations.
Size getMaterialCnt() const
Return the number of materials in the storage.
Quantity & insert(const QuantityId key, const OrderEnum order, const TValue &defaultValue)
Creates a quantity in the storage, given its key, value type and order.
StaticArray< Array< TValue > &, 3 > getAll(const QuantityId key)
Retrieves quantity buffers from the storage, given its key and value type.
auto getValues(const QuantityId first, const QuantityId second, const TArgs... others)
Retrieves an array of quantities from the key.
Array< TValue > & getD2t(const QuantityId key)
Retrieves a quantity second derivative from the storage, given its key and value type.
MaterialView getMaterial(const Size matIdx) const
Returns an object containing a reference to given material.
void zeroHighestDerivatives()
Sets all highest-level derivatives of quantities to zero.
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
virtual void create(Storage &storage, IMaterial &material) const override
Initializes all quantities needed by the solver in the storage.
virtual void integrate(Storage &storage, Statistics &stats) override
Computes derivatives of all time-dependent quantities.
Symmetric tensor of 2nd order.
INLINE Float trace() const
Return the trace of the tensor.
static INLINE SymmetricTensor identity()
Returns an identity tensor.
INLINE Float value(const Vector &r1, const Vector &r2) const
INLINE Float radius() const
INLINE Vector grad(const Vector &r1, const Vector &r2) const
Symmetric traceless 2nd order tensor.
static INLINE TracelessTensor null()
Returns a tensor with all zeros.
Creating code components based on values from settings.
@ ELASTIC_MODULUS
Elastic modulus lambda (a.k.a Lame's first parameter) of the material.
@ SHEAR_MODULUS
Shear modulus mu (a.k.a Lame's second parameter) of the material.
constexpr Float gravity
Gravitational constant (CODATA 2014)
AutoPtr< ISymmetricFinder > getFinder(const RunSettings &settings)
Overload of std::swap for Sph::Array.
Holds information about a neighbour particles.