20 : domain(
std::move(domain)) {
27 : domain(
std::move(domain)) {
29 params.searchRadius = Factory::getKernel<3>(settings).radius();
48 domain->
addGhosts(r, ghosts, params.searchRadius, params.minimalDist);
52 for (
Ghost& g : ghosts) {
60 for (
Size i = 0; i < ghosts.size(); ++i) {
61 const Size ghostIdx = ghostIdxs[i];
62 r[ghostIdx] = ghosts[i].position;
68 for (
Size i = 0; i < ghosts.size(); ++i) {
69 const Size ghostIdx = ghostIdxs[i];
71 const Vector deltaR = r[ghosts[i].index] - ghosts[i].position;
74 const Float perp =
dot(normal, v[ghosts[i].index]);
76 const Vector v0 = v[ghostIdx] - 2._f * normal * perp;
78 v[ghostIdx] = ghostVelocity(ghosts[i].
position).valueOr(v0);
89 for (
Size i = 0; i < ghosts.size(); ++i) {
90 const Size ghostIdx = ghostIdxs[i];
91 flag[ghostIdx] =
FLAG;
101 ghostVelocity = newGhostVelocity;
106 "Solver changed the number of particles. This is currently not consistent with the implementation of "
110 storage.
remove(ghostIdxs);
113 storage.
setUserData(makeShared<GhostParticlesData>(std::move(ghosts)));
121 : fixedParticles(
std::move(params.material)) {
123 Box box = params.domain->getBoundingBox();
139 for (
Size i = 0; i < dummies.
size();) {
140 if (params.domain->contains(dummies[i])) {
159 solver->create(fixedParticles, mat);
161 mat->create(fixedParticles, context);
176 for (
Size i = 0; i < flag.
size(); ++i) {
177 if (flag[i] ==
Size(-1)) {
195 : domain(
std::move(domain))
220 for (
Size i = 0; i < r.
size(); ++i) {
231 for (
Size i = 0; i < r.
size(); ++i) {
239 iterate<VisitorEnum::HIGHEST_DERIVATIVES>(storage, [
this](
QuantityId,
auto& d2f) {
240 using T =
typename std::decay_t<decltype(d2f)>::Type;
261 for (
Size i = 0; i < r.
size(); ++i) {
266 iterate<VisitorEnum::ALL_BUFFERS>(storage, [&toRemove](
auto& v) { v.
remove(toRemove); });
325 for (
Size i = 0; i < r.
size(); ++i) {
328 int(r[i][
Y] < domain.
lower()[
Y]),
329 int(r[i][
Z] < domain.
lower()[
Z]));
331 int(r[i][
Y] > domain.
upper()[
Y]),
332 int(r[i][
Z] > domain.
upper()[
Z]));
335 r[i] += domain.
size() * (lowerFlags - upperFlags);
339 for (
Size j = 0; j < 3; ++j) {
354 for (
Size i = 0; i < ghostIdxs.
size(); ++i) {
355 r[ghostIdxs[i]] = ghosts[i].position;
356 r[ghostIdxs[i]][
H] = r[ghosts[i].index][
H];
362 for (
Size i = 0; i < ghosts.
size(); ++i) {
363 flag[ghostIdxs[i]] =
Size(-1);
369 storage.
remove(ghostIdxs);
372 storage.
setUserData(makeShared<GhostParticlesData>(std::move(ghosts)));
387 for (
Size i = 0; i < r.
size(); ++i) {
388 if (r[i][
Z] < 0._f) {
389 r[i][
Z] = 0.1_f * r[i][
H];
393 for (
Size j = 0; j < 3; ++j) {
404 for (
Size i = 0; i < ghostIdxs.
size(); ++i) {
405 r[ghostIdxs[i]] = ghosts[i].position;
410 storage.
remove(ghostIdxs);
413 storage.
setUserData(makeShared<GhostParticlesData>(std::move(ghosts)));
426 for (
Size i = 0; i < r.
size(); ++i) {
446 for (
Size i = 0; i < r.size(); ++i) {
448 r[i] =
Vector(domain.
clamp(r[i][0]), 0._f, 0._f, r[i][
H]);
449 v[i] =
Vector(v[i][0], 0._f, 0._f);
454 using Type =
typename std::decay_t<decltype(dv)>::Type;
456 for (
Size i : { 0u, 1u, 2u, 3u, 4u, s - 4, s - 3, s - 2, s - 1 }) {
460 iterate<VisitorEnum::SECOND_ORDER>(
462 using Type =
typename std::decay_t<decltype(dv)>::Type;
464 for (
Size i : { 0u, 1u, 2u, 3u, 4u, s - 4, s - 3, s - 2, s - 1 }) {
INLINE bool isReal(const AntisymmetricTensor &t)
#define SPH_ASSERT(x,...)
Filling spatial domain with SPH particles.
Object defining computational domain.
@ OUTSIDE
Marks all vectors outside of the domain.
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 all particle materials.
Base interface for all solvers.
Functions for iterating over individual quantities in Storage.
Ordinary iterator for custom containers.
constexpr INLINE Float pow< 3 >(const Float v)
#define NAMESPACE_SPH_END
QuantityId
Unique IDs of basic quantities of SPH particles.
@ FLAG
ID of original body, used to implement discontinuities between bodies in SPH.
@ DEVIATORIC_STRESS
Deviatoric stress tensor, always a traceless tensor.
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
@ MASS
Paricles masses, always a scalar quantity.
@ SECOND
Quantity with 1st and 2nd derivative.
@ ZERO
Quantity without derivatives, or "zero order" of quantity.
SequentialScheduler SEQUENTIAL
Global instance of the sequential scheduler.
Interface for executing tasks (potentially) asynchronously.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
INLINE Float getLength(const Vector &v)
Returns the length of the vector. Enabled only for vectors of floating-point precision.
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)
INLINE TCounter size() const
Generic dynamically allocated resizable storage.
INLINE void push(U &&u)
Adds new element to the end of the array, resizing the array if necessary.
void remove(const TCounter idx)
Removes an element with given index from the array.
ArrayView< T, TCounter > view() noexcept
Explicit conversion to arrayview.
void clear()
Removes all elements from the array, but does NOT release the memory.
INLINE TCounter size() const noexcept
INLINE bool empty() const noexcept
Block aligned with coordinate axes, defined by its center and length of each side.
Helper object defining three-dimensional interval (box).
INLINE void extend(const Vector &v)
Enlarges the box to contain the vector.
INLINE const Vector & lower() const
Returns lower bounds of the box.
INLINE const Vector & upper() const
Returns upper bounds of the box.
INLINE Vector center() const
Returns the center of the box.
INLINE Vector size() const
Returns box dimensions.
INLINE Float volume() const
Returns the volume of the box.
FixedParticles(const RunSettings &settings, Params &¶ms)
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
virtual void initialize(Storage &storage) override
Applies the boundary conditions before the derivatives are computed.
Boundary condition that nulls all highest derivates of selected particles.
void freeze(const Size flag)
Adds body ID particles of which shall be frozen by boundary conditions.
FrozenParticles()
Constructs boundary conditions with no particles frozen.
SharedPtr< IDomain > domain
void thaw(const Size flag)
Remove a body from the list of frozen bodies.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
virtual void initialize(Storage &storage) override
Applies the boundary conditions before the derivatives are computed.
GhostParticles(SharedPtr< IDomain > domain, const Float searchRadius, const Float minimalDist)
Creates new boundary conditions.
void setVelocityOverride(Function< Optional< Vector >(const Vector &r)> ghostVelocity)
Specifies a functor that overrides the default velocity assinged to each ghost.
static constexpr Size FLAG
Special flag denoting ghost particles.
virtual void getSubset(ArrayView< const Vector > vs, Array< Size > &output, const SubsetType type) const =0
Returns an array of indices, marking vectors with given property by their index.
virtual void getDistanceToBoundary(ArrayView< const Vector > vs, Array< Float > &distances) const =0
Returns distances of particles lying close to the boundary.
virtual void project(ArrayView< Vector > vs, Optional< ArrayView< Size >> indices=NOTHING) const =0
Projects vectors outside of the domain onto its boundary.
virtual void addGhosts(ArrayView< const Vector > vs, Array< Ghost > &ghosts, const Float eta=2._f, const Float eps=0.05_f) const =0
Duplicates positions located close to the boundary, placing copies ("ghosts") symmetrically to the ot...
virtual bool contains(const Vector &v) const =0
Checks if the given point lies inside the domain.
Object representing a 1D interval of real numbers.
INLINE Float clamp(const Float &value) const
Clamps the given value by the interval.
KillEscapersBoundary(SharedPtr< IDomain > domain)
virtual void initialize(Storage &storage) override
Applies the boundary conditions before the derivatives are computed.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
Non-owning wrapper of a material and particles with this material.
Wrapper of type value of which may or may not be present.
PeriodicBoundary(const Box &domain)
virtual void initialize(Storage &storage) override
Applies the boundary conditions before the derivatives are computed.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
Projection1D(const Interval &domain)
Constructs using range as 1D domain.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
TValue get(const TEnum idx, std::enable_if_t<!std::is_enum< std::decay_t< TValue >>::value, int >=0) const
Returns a value of given type from the settings.
Array with fixed number of allocated elements.
Container storing all quantities used within the simulations.
void merge(Storage &&other)
Merges another storage into this object.
Outcome isValid(const Flags< ValidFlag > flags=ValidFlag::COMPLETE) const
Checks whether the storage is in valid state.
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.
Array< TValue > & getDt(const QuantityId key)
Retrieves a quantity derivative from the storage, given its key and value type.
void setUserData(SharedPtr< IStorageUserData > newData)
Stores new user data into the storage.
Size getParticleCnt() const
Returns the number of particles.
Array< Size > duplicate(ArrayView< const Size > idxs, const Flags< IndicesFlag > flags=EMPTY_FLAGS)
Duplicates some particles in the storage.
@ INDICES_SORTED
Use if the given array is already sorted (optimization)
void remove(ArrayView< const Size > idxs, const Flags< IndicesFlag > flags=EMPTY_FLAGS)
Removes specified particles from the storage.
MaterialView getMaterial(const Size matIdx) const
Returns an object containing a reference to given material.
bool has(const QuantityId key) const
Checks if the storage contains quantity with given key.
Storage clone(const Flags< VisitorEnum > buffers) const
Clones specified buffers of the storage.
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
virtual void initialize(Storage &storage) override
Applies the boundary conditions before the derivatives are computed.
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
Symmetric traceless 2nd order tensor.
WindTunnel(SharedPtr< IDomain > domain, const Float radius)
virtual void finalize(Storage &storage) override
Applies the boundary conditions after the derivatives are computed.
Creating code components based on values from settings.
@ DENSITY
Density at zero pressure.
@ DOMAIN_GHOST_MIN_DIST
Minimal distance between a particle and its ghost, in units of smoothing length.
AutoPtr< ISolver > getSolver(IScheduler &scheduler, const RunSettings &settings)
Vector position(const Float a, const Float e, const Float u)
Computes the position on the elliptic trajectory.
Overload of std::swap for Sph::Array.
Shared data used when creating all bodies in the simulation.