27#ifndef BLOCK_LATTICE_INTERACTION_HH
28#define BLOCK_LATTICE_INTERACTION_HH
38template <
typename T,
unsigned D>
55 bool validRange =
true;
56 for (
unsigned iDim = 0; iDim < D; ++iDim) {
60 int intersectionRange = end[iDim] - start[iDim];
61 validRange = validRange && (intersectionRange >= 0.);
67template <
typename T,
unsigned D>
75 surfaceOutOfGeometry =
false;
76 for (
unsigned i = 0; i < D; ++i) {
77 T
const particleMin = position[i] - circumRadius;
78 T
const particleMax = position[i] + circumRadius;
79 if (particleMin < cellMin[i] && periodic[i]) {
80 surfaceOutOfGeometry =
true;
82 position[i], cellMax[i], cellMin[i]);
84 else if (particleMax > cellMax[i] && periodic[i]) {
85 surfaceOutOfGeometry =
true;
87 position[i], cellMax[i], cellMin[i]);
94template <
typename T,
typename DESCRIPTOR,
typename PARTICLETYPE>
103 constexpr unsigned D = DESCRIPTOR::d;
104 using namespace descriptors;
108 bool surfaceOutOfGeometry;
111 cellMax, originalPosition, circumRadius,
115 if (!surfaceOutOfGeometry) {
120 particle.template setField<GENERAL, POSITION>(ghostPos);
123 particle.template setField<GENERAL, POSITION>(originalPosition);
127 if constexpr(!particles::access::providesParallelization<PARTICLETYPE>()) {
129 bool centerOfMassOutOfDomain =
false;
130 for (
unsigned iDim = 0; iDim < D; ++iDim) {
131 centerOfMassOutOfDomain = centerOfMassOutOfDomain ||
132 originalPosition[iDim] < cellMin[iDim] ||
133 originalPosition[iDim] > cellMax[iDim];
138 if (centerOfMassOutOfDomain) {
139 particle.template setField<GENERAL, POSITION>(ghostPos);
145template <
typename T,
typename DESCRIPTOR,
typename PARTICLETYPE,
146 typename PARTICLECONTACTTYPE,
typename WALLCONTACTTYPE,
typename F>
158 F getSetupPeriodicity)
160 constexpr unsigned D = DESCRIPTOR::d;
161 using namespace descriptors;
164 bool surfaceOutOfGeometry;
168 surfaceOutOfGeometry, ghostPos, cellMin, cellMax, originalPosition,
169 circumRadius, getSetupPeriodicity());
172 if (!surfaceOutOfGeometry) {
174 particleSystem, contactContainer, iP, particle,
179 particle.template setField<GENERAL, POSITION>(ghostPos);
181 particleSystem, contactContainer, iP, particle,
182 solidBoundaries, cellMin, cellMax,
183 getSetupPeriodicity);
185 particle.template setField<GENERAL, POSITION>(originalPosition);
187 particleSystem, contactContainer, iP, particle,
188 solidBoundaries, cellMin, cellMax,
189 getSetupPeriodicity);
194template <
typename T,
typename DESCRIPTOR,
typename F>
200 constexpr unsigned D = DESCRIPTOR::d;
204 T {1} / blockGeometry.
getDeltaR(), start,
205 end, position, circumRadius)) {
219template <
typename T,
typename DESCRIPTOR,
typename PARTICLETYPE>
225 constexpr unsigned D = DESCRIPTOR::d;
227 using namespace descriptors;
229 particle.template getField<SURFACE, SINDICATOR>()->getCircumRadius();
230 auto position = particle.template getField<GENERAL, POSITION>();
235 blockGeometry, blockLattice, padding, position, circumRadius,
239 blockGeometry.
getPhysR(physR, latticeR);
246 auto material = blockGeometry.
get(latticeR);
253 for (
unsigned iDim = 0; iDim != D; ++iDim) {
257 T solidVolumeFraction = 1. - porosity[0];
260 auto cell = blockLattice.
get(latticeR);
261 if constexpr (!DESCRIPTOR::template provides<descriptors::VELOCITY_SOLID>()) {
262 cell.template getFieldPointer<
264 cell.template getFieldPointer<
268 cell.template getFieldPointer<
271 cell.template getFieldPointer<descriptors::POROSITY>() *=
279template <
typename T,
typename DESCRIPTOR,
typename PARTICLETYPE,
280 typename PARTICLECONTACTTYPE,
typename WALLCONTACTTYPE,
typename F>
293 using namespace descriptors;
294 constexpr unsigned D = DESCRIPTOR::d;
295 const T deltaX = blockGeometry.
getDeltaR();
296 auto sIndicator = particle.template getField<SURFACE, SINDICATOR>();
299 detectionDistance, sIndicator->getCircumRadius(), sIndicator->getEpsilon());
300 auto position = particle.template getField<GENERAL, POSITION>();
307 blockGeometry, blockLattice, padding, position, circumRadius,
311 blockGeometry.
getPhysR(physR, latticeR);
317 sIndicator->getEpsilon());
322 auto material = blockGeometry.
getMaterial(latticeR);
323 if (material != 0 && wallMaterials.find(material)==wallMaterials.end()) {
328 for (
unsigned iDim = 0; iDim != D; ++iDim) {
332 T solidVolumeFraction = 1. - porosity[0];
335 auto cell = blockLattice.
get(latticeR);
336 if constexpr (!DESCRIPTOR::template provides<descriptors::VELOCITY_SOLID>()) {
337 cell.template getFieldPointer<descriptors::VELOCITY_NUMERATOR>() +=
338 porosity[0] * velocity;
339 cell.template getFieldPointer<
343 cell.template getFieldPointer<descriptors::VELOCITY_SOLID>() += velocity;
345 cell.template getFieldPointer<descriptors::POROSITY>() *=
353 if (signedDist < detectionDistance) {
355 bool ignoreCell =
false;
356 for (
unsigned i = 0; i < D; ++i) {
357 if (latticeR[i] < 0 ||
358 latticeR[i] > blockGeometry.
getExtent()[i] - 1) {
365 auto contactHelper = blockLattice.
get(latticeR)
366 .template getFieldPointer<
369 std::size_t iP2 = contactHelper[0] - 1;
370 if (contactHelper[0] > 0 && iP2 != iP) {
371 std::size_t localiP2 = iP2;
374 for (std::size_t localiP = 0; localiP < particleSystem.
size();
376 const std::size_t globalParticleID =
377 particleSystem.
get(localiP)
378 .template getField<PARALLELIZATION, ID>();
379 if (globalParticleID == iP2) {
385 auto particle2 = particleSystem.
get(localiP2);
388 bool detectParticleContacts =
true;
389 if constexpr (access::providesComputeContact<PARTICLETYPE>()) {
390 detectParticleContacts =
394 if (detectParticleContacts) {
396 contactContainer, std::array<std::size_t, 2>({iP2, iP}),
397 PhysR<T, D>(physR), particle2, particle, cellMin, cellMax,
398 getSetupPeriodicity, deltaX);
402 contactHelper[0] = iP + 1;
406 bool detectWallContacts =
true;
407 if constexpr (access::providesComputeContact<PARTICLETYPE>()) {
412 if (detectWallContacts) {
415 unsigned solidBoundaryID = 0;
418 const T solidBoundaryDetectionDistance =
420 solidBoundary.getEnlargementForContact();
421 if (solidBoundary.getIndicator()->signedDistance(
422 PhysR<T, D>(physR)) < solidBoundaryDetectionDistance) {
424 contactContainer, iP, solidBoundaryID,
PhysR<T, D>(physR),
425 particle, cellMin, cellMax, getSetupPeriodicity, deltaX);
435template <
typename T,
typename DESCRIPTOR>
440 constexpr unsigned D = DESCRIPTOR::d;
445 using namespace descriptors;
447 blockMin, blockMax, [&](
const LatticeR<D>& latticeR) {
448 auto cell = blockLattice.
get(latticeR);
449 particles::resetAllParticleRelatedFields<DESCRIPTOR, decltype(cell), typename decltype(cell)::value_t>(cell);
454template <
typename T,
typename DESCRIPTOR>
459 constexpr unsigned D = DESCRIPTOR::d;
464 using namespace descriptors;
467 auto cell = blockLattice.
get(latticeR);
468 particles::resetParticleContactRelatedFields<DESCRIPTOR, decltype(cell), typename decltype(cell)::value_t>(cell);
Representation of a block geometry.
Vector< T, D > getPhysR(LatticeR< D > latticeR)
int getMaterial(LatticeR< D > latticeR) const
returns the (iX,iY) entry in the 2D scalar field
T getDeltaR() const
Read only access to the voxel size given in SI units (meter)
Vector< T, D > getOrigin() const
Read only access to the origin position given in SI units (meter)
Vector< int, D > getExtent() const
Returns the extend of the block in lattice units.
std::enable_if_t< sizeof...(L)==D, int > get(L... latticeR) const
Read-only access to a material number.
Platform-abstracted block lattice for external access and inter-block interaction.
Cell< T, DESCRIPTOR > get(CellID iCell)
Get Cell interface for index iCell.
void forSpatialLocations(F f) const
int getPadding() const
Read only access to padding.
Conversion between physical and lattice units, as well as discretization.
constexpr T getLatticeVelocity(T physVelocity) const
conversion from physical to lattice velocity
auto & get()
Expose container.
constexpr std::size_t size()
Size of ParticleSystem.
constexpr bool providesParallelization()
Provides group PARALLELIZATION.
T getRadius(Particle< T, PARTICLETYPE > &particle)
bool isContactComputationEnabled(Particle< T, PARTICLETYPE > &particle)
Check if contact should be regarded (specification for a single particle)
Vector< T, PARTICLETYPE::d > getPosition(Particle< T, PARTICLETYPE > particle)
T movePositionToEnd(const T position, const T max, const T min)
T movePositionToStart(const T position, const T max, const T min)
T evalCircumRadius(T const contactDetectionDistance, T const circumRadius, T const epsilon)
void updateContacts(ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, std::array< std::size_t, 2 > &&ids, const PhysR< T, PARTICLETYPE::d > &pos, Particle< T, PARTICLETYPE > &particle1, Particle< T, PARTICLETYPE > &particle2, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
T evalContactDetectionDistance(Particle< T, PARTICLETYPE > &particle, T const physDeltaX)
constexpr Vector< T, PARTICLETYPE::d > calculateLocalVelocity(Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &input)
const S signedDistanceToParticle(Particle< T, PARTICLETYPE > &particle, const PhysR< S, PARTICLETYPE::d > &input)
bool evalSolidVolumeFraction(T output[], const S input[], Particle< T, PARTICLETYPE > &particle)
void resetBlockParticleField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice)
bool getBlockParticleIntersection(const BlockGeometry< T, D > &blockGeometry, T invDeltaX, LatticeR< D > &start, LatticeR< D > &end, Vector< T, D > position, T circumRadius)
void setBlockParticleField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice, UnitConverter< T, DESCRIPTOR > const &converter, Particle< T, PARTICLETYPE > &particle)
void resetBlockContactField(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice)
void checkSmoothIndicatorOutOfGeometry(bool &outOfGeometry, Vector< T, D > &ghostPos, const PhysR< T, D > &cellMin, const PhysR< T, D > &cellMax, const Vector< T, D > &position, T circumRadius, const Vector< bool, D > &periodic)
void forSpatialLocationsInBlockParticleIntersection(const BlockGeometry< T, DESCRIPTOR::d > &blockGeometry, BlockLattice< T, DESCRIPTOR > &blockLattice, Vector< T, DESCRIPTOR::d > position, T circumRadius, F f)
bool evalSolidVolumeFraction(T output[], T signedDist, T eps) any_platform
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > ceil(const ADf< T, DIM > &a)
ADf< T, DIM > floor(const ADf< T, DIM > &a)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
bool nearZero(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.
std::unordered_set< int > getLatticeMaterials(const std::vector< SolidBoundary< T, D > > &solidBoundaries)
Get material numbers of multiple solid boundaries in std::vector.