olb::particles::contact Namespace Reference

Classes
struct	ContactContainer
struct	ContactProperties
	Object that stores properties which are necessary for the computation of contact forces N = number of different materials The material here is an identifier of a solid material with certain (mechanical) properties This material is something completely different from the lattice's material number, which is used to assign boundary conditions and dynamics. More...
struct	ContactProperty
struct	MaterialProperties
	Class storing properties that are necessary for the computation of the contact orce N = number of different materials [0]: modulus of elasticity [1]: Poisson's ratio. More...
struct	particle_particle
struct	particle_particle< T, PARTICLETYPE, ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX >, WALLCONTACTTYPE, BBCORRECTIONMETHOD, CONVEX, useSDF >
struct	particle_wall
struct	particle_wall< T, PARTICLETYPE, PARTICLECONTACTTYPE, WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX >, BBCORRECTIONMETHOD, CONVEX, useSDF >
struct	ParticleContact
struct	ParticleContactArbitraryFromOverlapVolume
	An object holding data for a contact which is described analog to Nassauer and Kuna (2013) More...
struct	WallContact
struct	WallContactArbitraryFromOverlapVolume

Functions
template<typename T >
constexpr T	evalEffectiveYoungModulus (T E1, T E2, T nu1, T nu2)
template<typename T >
constexpr T	evalDampingFactor (const T coefficientOfRestitution, const T initialRelativeVelocityMagnitude)
	Calculates the damping factor according to Carvalho & Martins (2019) (10.1016/j.mechmachtheory.2019.03.028)
template<typename T , typename PARTICLETYPE >
T	evalContactDetectionDistance (Particle< T, PARTICLETYPE > &particle, T const physDeltaX)
template<typename T >
T	evalCircumRadius (T const contactDetectionDistance, T const circumRadius, T const epsilon)
template<typename T , typename PARTICLETYPE >
T	evalCircumRadius (Particle< T, PARTICLETYPE > &particle, T const physDeltaX, T const circumRadius, T const epsilon)
template<typename T , typename PARTICLETYPE >
T	evalCircumRadius (Particle< T, PARTICLETYPE > &particle, T const physDeltaX)
template<typename CONTACTTYPE >
void	cleanContacts (std::vector< CONTACTTYPE > &contacts)
template<typename CONTACTTYPE >
void	resetResponsibleRank (CONTACTTYPE &contact)
template<typename CONTACTTYPE >
bool	isResponsibleRankSet (CONTACTTYPE &contact)
template<typename T , unsigned D>
void	updateMinMax (PhysR< T, D > &min, PhysR< T, D > &max, const PhysR< T, D > &pos)
template<typename T , unsigned D>
unsigned short	domainDimension (PhysR< T, D > &min, PhysR< T, D > &max)
template<typename T , unsigned N, typename F >
constexpr ContactProperties< T, N >	createContactProperties (F f)
template<typename T , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >
void	communicateContacts (ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >
void	accountForPeriodicParticleBoundary (SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, const SuperGeometry< T, PARTICLETYPE::d > &sGeometry, F getSetupPeriodicity)
	update positions in contact to account for periodic boundaries
template<typename T , typename PARTICLETYPE , bool CONVEX>
std::unordered_set< int >	evalDestRanksForDetectionCommunication (ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, const Vector< bool, PARTICLETYPE::d > &periodicity)
	evaluate responsible ranks for particle-particle contact
template<typename T , typename PARTICLETYPE , bool CONVEX>
std::unordered_set< int >	evalDestRanksForDetectionCommunication (WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, const Vector< bool, PARTICLETYPE::d > &periodicity)
	evaluate responsible rank for solid boundary contact
template<typename T , typename PARTICLETYPE , bool CONVEX>
std::unordered_set< int >	evalDestRanksForPostContactTreatmentCommunication (ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, const Vector< bool, PARTICLETYPE::d > &periodicity)
	evaluate ranks that touch both particles of a particle-particle contact
template<typename T , typename PARTICLETYPE , bool CONVEX>
std::unordered_set< int >	evalDestRanksForPostContactTreatmentCommunication (WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, const Vector< bool, PARTICLETYPE::d > &periodicity)
	evaluate ranks that touch the particle of a particle-wall contact
template<typename T , typename PARTICLETYPE , typename CONTACTTYPE , bool CONVEX, bool SANITYCHECK = false>
void	receiveContact (ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &newcontact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, std::vector< CONTACTTYPE > &contacts, int rankOrig)
template<typename T , typename PARTICLETYPE , typename CONTACTTYPE , bool CONVEX, bool SANITYCHECK = false>
void	receiveContact (WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &newcontact, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, std::vector< CONTACTTYPE > &contacts, int rankOrig)
template<typename T , typename PARTICLETYPE , typename CONTACTTYPE >
void	communicateParallelContacts (std::vector< CONTACTTYPE > &contacts, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, MPI_Comm contactDetectionComm, const Vector< bool, PARTICLETYPE::d > &periodicity)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >
void	communicateContactsParallel (ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, MPI_Comm particleContactDetectionComm, MPI_Comm wallContactDetectionComm, const Vector< bool, PARTICLETYPE::d > &periodicity)
template<typename T , typename PARTICLETYPE , typename CONTACTTYPE >
void	communicatePostContactTreatmentContacts (std::vector< CONTACTTYPE > &contacts, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, MPI_Comm postContactTreatmentComm, const Vector< bool, PARTICLETYPE::d > &periodicity)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >
void	communicatePostContactTreatmentContacts (ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, const T deltaX, MPI_Comm particlePostContactTreatmentComm, MPI_Comm wallPostContactTreatmentComm, const Vector< bool, PARTICLETYPE::d > &periodicity)
	It is necessary to communicate contacts so that no information is lost.
template<typename T , unsigned D>
int	getContactTreatmentResultsSerialSize ()
template<typename T , unsigned D>
void	extendContactTreatmentResultsDataMap (const std::size_t &globalParticleID, Vector< T, D > &force, Vector< T, utilities::dimensions::convert< D >::rotation > &torque, const int destRank, std::multimap< int, std::unique_ptr< std::uint8_t[]> > &dataMap)
template<typename T , typename PARTICLETYPE >
void	receiveContactTreatmentResults (SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, MPI_Comm contactTreatmentComm, singleton::MpiNonBlockingHelper &mpiNbHelper)
template<typename T , typename PARTICLETYPE >
void	communicateContactTreatmentResults (SuperParticleSystem< T, PARTICLETYPE > &particleSystem, std::multimap< int, std::unique_ptr< std::uint8_t[]> > &dataMap, MPI_Comm contactTreatmentComm)
template<typename CONTACTTYPE >
auto	getContactIterator (std::vector< CONTACTTYPE > &contacts, const std::function< bool(const CONTACTTYPE &)> condition)
std::array< std::size_t, 2 >	sortParticleIDs (const std::array< std::size_t, 2 > &ids)
template<typename T , typename PARTICLETYPE , typename F >
void	unifyPositions (Particle< T, PARTICLETYPE > &particle1, Particle< T, PARTICLETYPE > &particle2, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, std::array< PhysR< T, PARTICLETYPE::d >, 2 > &positions, T deltaX)
template<typename T , typename PARTICLETYPE , typename F >
PhysR< T, PARTICLETYPE::d >	unifyPosition (Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
template<typename T , typename PARTICLETYPE , typename F >
PhysR< T, PARTICLETYPE::d >	evalContactPosition (Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &particlePos, const PhysR< T, PARTICLETYPE::d > &contactPos, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
template<typename PARTICLECONTACTTYPE , bool IS_INPUT_SORTED = false>
bool	particleContactConsistsOfIDs (PARTICLECONTACTTYPE &particleContact, const std::array< size_t, 2 > &ids)
template<typename T , typename PARTICLETYPE , bool CONVEX, typename F >
void	updateContact (ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, Particle< T, PARTICLETYPE > &particle1, Particle< T, PARTICLETYPE > &particle2, const PhysR< T, PARTICLETYPE::d > &contactPos, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
template<typename T , typename PARTICLETYPE , bool CONVEX, typename F >
void	updateContact (WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &contact, Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &contactPos, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >
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)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >
void	updateContacts (ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, size_t particleID, unsigned wallID, const PhysR< T, PARTICLETYPE::d > &pos, Particle< T, PARTICLETYPE > &particle, const PhysR< T, PARTICLETYPE::d > &cellMin, const PhysR< T, PARTICLETYPE::d > &cellMax, F getSetupPeriodicity, T deltaX)
template<typename T , unsigned D, typename CONTACTTYPE >
constexpr T	evalCurrentDampingFactor (CONTACTTYPE &contact, const T coefficientOfRestitution, const Vector< T, D > &initialRelativeNormalVelocity)
template<typename T , typename PARTICLETYPE >
constexpr Vector< T, PARTICLETYPE::d >	evalRelativeVelocity (Particle< T, PARTICLETYPE > &particleA, Particle< T, PARTICLETYPE > &particleB, const Vector< T, PARTICLETYPE::d > &position)
	Returns the relative velocity of two particles at a defined position.
template<typename T , typename PARTICLETYPE >
constexpr Vector< T, PARTICLETYPE::d >	evalRelativeVelocity (Particle< T, PARTICLETYPE > &particle, const Vector< T, PARTICLETYPE::d > &position)
template<typename T , unsigned D>
constexpr Vector< T, D >	evalRelativeNormalVelocity (const Vector< T, D > &contactNormal, const Vector< T, D > &relVel)
template<typename T , typename F >
constexpr void	forEachPosition (Vector< T, 3 > startPos, Vector< T, 3 > endPos, F f)
template<typename T , typename F >
constexpr void	forEachPosition (Vector< T, 2 > startPos, Vector< T, 2 > endPos, F f)
template<typename T , typename F >
constexpr void	forEachPositionWithBreak (Vector< T, 3 > startPos, Vector< T, 3 > endPos, F f)
template<typename T , typename F >
constexpr void	forEachPositionWithBreak (Vector< T, 2 > startPos, Vector< T, 2 > endPos, F f)
template<typename T , unsigned D>
constexpr Vector< T, D >	evalContactDeltaX (const Vector< T, D > &min, const Vector< T, D > &max, const T physDeltaX, const unsigned contactBoxResolutionPerDirection)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , bool CONVEX>
void	prepareForceCalculation (ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer)
template<typename T , unsigned D>
Vector< T, D >	calcElasticAndViscousForce (T effectiveE, T k, T overlapVolume, T indentation, T dampingFactor, const Vector< T, D > &relVelNormal, const Vector< T, D > &contactNormal)
	Elastic and viscous force from overlap volume according to Nassauer & Kuna (2013)
template<typename T , unsigned D>
T	evalApproxSurface (const Vector< T, D > &normalizedNormal, const PhysR< T, D > &contactPhysDeltaX)
	Takes a normalized normal and the voxel sizes and approximates the corresponding surface.
template<typename T , unsigned D>
Vector< T, D >	getTangentialForceFromOverlapVolume (const Vector< T, D > &contactNormal, const T coefficientStaticFriction, const T coefficientKineticFriction, const T staticKineticTransitionVelocity, const T k, const Vector< T, D > &relVel, const Vector< T, D > &normalForce)
template<typename T , unsigned D>
Vector< T, D >	getNormalForceFromOverlapVolume (const Vector< T, D > &contactNormal, const T indentation, const T overlapVolume, const T effectiveE, const T dampingFactor, const T k, const Vector< T, D > &relVel)
template<typename T , typename PARTICLETYPE >
void	applyForceOnParticle (std::multimap< int, std::unique_ptr< std::uint8_t[]> > &dataMap, Particle< T, PARTICLETYPE > &particle, XParticleSystem< T, PARTICLETYPE > &particleSystem, const Vector< T, PARTICLETYPE::d > &contactPoint, const Vector< T, PARTICLETYPE::d > &normalForce, const Vector< T, PARTICLETYPE::d > &tangentialForce)
template<typename T , typename PARTICLETYPE , typename CONTACTPROPERTIES , typename F = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>
void	applyForceFromOverlapVolume (std::multimap< int, std::unique_ptr< std::uint8_t[]> > &dataMap, Particle< T, PARTICLETYPE > &particleA, Particle< T, PARTICLETYPE > &particleB, XParticleSystem< T, PARTICLETYPE > &particleSystem, const PhysR< T, PARTICLETYPE::d > &contactPoint, const Vector< T, PARTICLETYPE::d > &contactNormal, const T indentation, const T overlapVolume, const Vector< T, PARTICLETYPE::d > &relVel, const T dampingFactor, const CONTACTPROPERTIES &contactProperties, T k, const std::size_t &particleAID, const std::size_t &particleBID, F processContactForce=defaults::processContactForce< T, PARTICLETYPE::d >)
template<typename T , typename PARTICLETYPE , typename CONTACTPROPERTIES , typename F = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>
void	applyForceFromOverlapVolume (std::multimap< int, std::unique_ptr< std::uint8_t[]> > &dataMap, Particle< T, PARTICLETYPE > &particle, XParticleSystem< T, PARTICLETYPE > &particleSystem, SolidBoundary< T, PARTICLETYPE::d > &solidBoundary, const PhysR< T, PARTICLETYPE::d > &contactPoint, const Vector< T, PARTICLETYPE::d > &contactNormal, const T indentation, const T overlapVolume, const Vector< T, PARTICLETYPE::d > &relVel, const T dampingFactor, const CONTACTPROPERTIES &contactProperties, T k, const std::size_t &particleID, const std::size_t &wallID, F processContactForce=defaults::processContactForce< T, PARTICLETYPE::d >)
template<typename T , unsigned D>
void	evalStartAndEndPos (Vector< long, D > &startPos, Vector< long, D > &endPos, const PhysR< T, D > &max, const PhysR< T, D > &min, const PhysR< T, D > &contactPhysDeltaX)
template<typename T , unsigned D>
T	evalForceViaVolumeCoefficient (const Vector< T, D > &min, const Vector< T, D > &max, T overlapVolume)
template<typename T >
constexpr T	evalForceViaVolumeCoefficient (const T contactVolume, const T contactSurface, const T indentationDepth)
template<typename T , unsigned D, typename F1 , typename F2 , typename F3 , typename F4 , typename F5 >
void	processCell (const PhysR< T, D > &pos, const PhysR< T, D > &contactPhysDeltaX, PhysR< T, D > &center, Vector< T, D > &contactNormal, unsigned &volumeCellCount, T &surfaceCellCount, F1 isInside, F2 isOnSurface, F3 calcNormalA, F4 calcNormalB, F5 updateMinMax)
template<typename T , unsigned D, typename F1 , typename F2 , typename F3 , typename F4 >
void	processContactViaVolume (const Vector< T, D > &min, const Vector< T, D > &max, T physDeltaX, unsigned contactBoxResolutionPerDirection, F1 resetContactBox, F2 processCellWrapped, F3 calculateIndentation, F4 applyForce)
template<typename T , unsigned D, typename F1 >
bool	evalStartingPoint (const Vector< long, D > &startPos, const Vector< long, D > &endPos, const Vector< T, D > &contactPhysDeltaX, Vector< long, D > &currPos, F1 isInsideContact)
template<typename T , unsigned D, typename F1 , typename F2 >
void	evalBoundingBoxRecursive (const Vector< long, D > &startPos, const Vector< long, D > &endPos, const Vector< T, D > &contactPhysDeltaX, const Vector< long, D > &currPos, const Vector< short, D > &dirIn, F1 isInsideContact, F2 updateMinMax)
template<typename T , unsigned D, typename F1 , typename F2 , typename F3 >
void	correctBoundingBox (const PhysR< T, D > &min, const PhysR< T, D > &max, const T physDeltaX, const unsigned contactBoxResolutionPerDirection, F1 isInsideContact, F2 resetContactBox, F3 updateMinMax)
template<typename T , unsigned D, typename F1 , typename F2 , typename F3 >
void	correctBoundingBoxNewContact (PhysR< T, D > min, PhysR< T, D > max, T physDeltaX, unsigned contactBoxResolutionPerDirection, F1 signedDistanceToIntersection, F2 resetContactBox, F3 updateMinMax)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename CONTACTPROPERTIES , unsigned BBCORRECTIONMETHOD = 0, typename F1 = decltype(defaults::periodicity<PARTICLETYPE::d>), typename F2 = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>
void	processContacts (SuperParticleSystem< T, PARTICLETYPE > &particleSystem, std::vector< SolidBoundary< T, PARTICLETYPE::d > > &solidBoundaries, ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, const CONTACTPROPERTIES &contactProperties, const SuperGeometry< T, PARTICLETYPE::d > &sGeometry, MPI_Comm contactTreatmentComm, const unsigned contactBoxResolutionPerDirection=8, const T k=static_cast< T >(4./(3 *util::sqrt(M_PI))), F1 getSetupPeriodicity=defaults::periodicity< PARTICLETYPE::d >, F2 processContactForce=defaults::processContactForce< T, PARTICLETYPE::d >)
template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename CONTACTPROPERTIES , unsigned BBCORRECTIONMETHOD = 0, typename F1 = decltype(defaults::periodicity<PARTICLETYPE::d>), typename F2 = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>
void	processContacts (ParticleSystem< T, PARTICLETYPE > &particleSystem, std::vector< SolidBoundary< T, PARTICLETYPE::d > > &solidBoundaries, ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &contactContainer, const CONTACTPROPERTIES &contactProperties, const SuperGeometry< T, PARTICLETYPE::d > &sGeometry, const unsigned contactBoxResolutionPerDirection=8, const T k=static_cast< T >(4./(3 *util::sqrt(M_PI))), F1 getSetupPeriodicity=defaults::periodicity< PARTICLETYPE::d >, F2 processContactForce=defaults::processContactForce< T, PARTICLETYPE::d >)

Function Documentation

accountForPeriodicParticleBoundary()

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >

void olb::particles::contact::accountForPeriodicParticleBoundary	(	SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		const SuperGeometry< T, PARTICLETYPE::d > &	sGeometry,
		F	getSetupPeriodicity )

update positions in contact to account for periodic boundaries

Definition at line 176 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
 
  // if a periodic setup is used, we have to update min and max in rare occasions
  // i.e. when the first/only particle "jumps" from the end to the beginning
  // because then the min and max from the step before isn't valid anymore
  const T physDeltaX =
      sGeometry.getCuboidGeometry().getMotherCuboid().getDeltaR();
  const PhysR<T, PARTICLETYPE::d> cellMin =
      particles::communication::getCuboidMin<T, PARTICLETYPE::d>(
          sGeometry.getCuboidGeometry());
  const PhysR<T, PARTICLETYPE::d> cellMax =
      particles::communication::getCuboidMax<T, PARTICLETYPE::d>(
          sGeometry.getCuboidGeometry(), cellMin);
 
  communication::forParticlesInSuperParticleSystem(
      sParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        const std::size_t currentGlobalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
 
        for (WALLCONTACTTYPE& wallContact : contactContainer.wallContacts) {
          if (!wallContact.isEmpty() && !wallContact.isNew() &&
              wallContact.getParticleID() == currentGlobalParticleID &&
              wallContact.getResponsibleRank() == singleton::mpi().getRank()) {
            const PhysR<T, PARTICLETYPE::d> min = wallContact.getMin();
            const PhysR<T, PARTICLETYPE::d> max = wallContact.getMax();
            wallContact.resetMinMax();
            const PhysR<T, PARTICLETYPE::d> unifiedPos = contact::unifyPosition(
                particle, cellMin, cellMax, getSetupPeriodicity, physDeltaX);
            wallContact.setParticlePosition(unifiedPos);
            const PhysR<T, PARTICLETYPE::d> newMin =
                contact::evalContactPosition(
                    particle, wallContact.getParticlePosition(), min, cellMin,
                    cellMax, getSetupPeriodicity, physDeltaX);
            const PhysR<T, PARTICLETYPE::d> newMax =
                contact::evalContactPosition(
                    particle, wallContact.getParticlePosition(), max, cellMin,
                    cellMax, getSetupPeriodicity, physDeltaX);
            wallContact.updateMinMax(newMin);
            wallContact.updateMinMax(newMax);
          }
        }
        for (PARTICLECONTACTTYPE& particleContact :
             contactContainer.particleContacts) {
          if (!particleContact.isEmpty() && !particleContact.isNew() &&
              particleContact.getIDs()[0] == currentGlobalParticleID &&
              particleContact.getResponsibleRank() ==
                  singleton::mpi().getRank()) {
            const PhysR<T, PARTICLETYPE::d> min = particleContact.getMin();
            const PhysR<T, PARTICLETYPE::d> max = particleContact.getMax();
            particleContact.resetMinMax();
            std::array<PhysR<T, PARTICLETYPE::d>, 2> unifiedPos;
            auto particleB = sParticleSystem.get(particleContact.getIDs()[1]);
            contact::unifyPositions(particle, particleB, cellMin, cellMax,
                                    getSetupPeriodicity, unifiedPos,
                                    physDeltaX);
            particleContact.setParticlePositions(unifiedPos);
            const PhysR<T, PARTICLETYPE::d> newMin =
                contact::evalContactPosition(
                    particle, particleContact.getParticlePositions()[0], min,
                    cellMin, cellMax, getSetupPeriodicity, physDeltaX);
            const PhysR<T, PARTICLETYPE::d> newMax =
                contact::evalContactPosition(
                    particle, particleContact.getParticlePositions()[0], max,
                    cellMin, cellMax, getSetupPeriodicity, physDeltaX);
            particleContact.updateMinMax(newMin);
            particleContact.updateMinMax(newMax);
          }
        }
      });
}

References evalContactPosition(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::get(), olb::SuperStructure< T, D >::getCuboidGeometry(), olb::singleton::MpiManager::getRank(), olb::singleton::mpi(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::particleContacts, unifyPosition(), unifyPositions(), and olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::wallContacts.

Here is the call graph for this function:

applyForceFromOverlapVolume() [1/2]

template<typename T , typename PARTICLETYPE , typename CONTACTPROPERTIES , typename F = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>

void olb::particles::contact::applyForceFromOverlapVolume	(	std::multimap< int, std::unique_ptr< std::uint8_t[]> > &	dataMap,
		Particle< T, PARTICLETYPE > &	particle,
		XParticleSystem< T, PARTICLETYPE > &	particleSystem,
		SolidBoundary< T, PARTICLETYPE::d > &	solidBoundary,
		const PhysR< T, PARTICLETYPE::d > &	contactPoint,
		const Vector< T, PARTICLETYPE::d > &	contactNormal,
		const T	indentation,
		const T	overlapVolume,
		const Vector< T, PARTICLETYPE::d > &	relVel,
		const T	dampingFactor,
		const CONTACTPROPERTIES &	contactProperties,
		T	k,
		const std::size_t &	particleID,
		const std::size_t &	wallID,
		F	processContactForce = defaults::processContactForce<T, PARTICLETYPE::d> )

Definition at line 448 of file particleContactForceFunctions.h.

{
  using namespace descriptors;
  constexpr unsigned D = PARTICLETYPE::d;
 
  const unsigned particleMaterial =
      particle.template getField<MECHPROPERTIES, MATERIAL>();
  const unsigned wallMaterial = solidBoundary.getContactMaterial();
 
  // Calculation of contact force
  const Vector<T, D> normalForce = getNormalForceFromOverlapVolume(
      contactNormal, indentation, overlapVolume,
      contactProperties.getEffectiveYoungsModulus(particleMaterial,
                                                  wallMaterial),
      dampingFactor, k, relVel);
  const Vector<T, D> tangentialForce = getTangentialForceFromOverlapVolume(
      contactNormal,
      contactProperties.getStaticFrictionCoefficient(particleMaterial,
                                                     wallMaterial),
      contactProperties.getKineticFrictionCoefficient(particleMaterial,
                                                      wallMaterial),
      contactProperties.getStaticKineticTransitionVelocity(particleMaterial,
                                                           wallMaterial),
      k, relVel, normalForce);
 
  processContactForce(contactPoint, contactNormal, normalForce, tangentialForce,
                      std::array<std::size_t, 2>({particleID, wallID}),
                      overlapVolume, indentation, true);
 
  // Applying contact forces on the particle
  applyForceOnParticle(dataMap, particle, particleSystem, contactPoint,
                       normalForce, tangentialForce);
}

References applyForceOnParticle(), olb::SolidBoundary< T, D >::getContactMaterial(), getNormalForceFromOverlapVolume(), and getTangentialForceFromOverlapVolume().

Here is the call graph for this function:

applyForceFromOverlapVolume() [2/2]

template<typename T , typename PARTICLETYPE , typename CONTACTPROPERTIES , typename F = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>

void olb::particles::contact::applyForceFromOverlapVolume	(	std::multimap< int, std::unique_ptr< std::uint8_t[]> > &	dataMap,
		Particle< T, PARTICLETYPE > &	particleA,
		Particle< T, PARTICLETYPE > &	particleB,
		XParticleSystem< T, PARTICLETYPE > &	particleSystem,
		const PhysR< T, PARTICLETYPE::d > &	contactPoint,
		const Vector< T, PARTICLETYPE::d > &	contactNormal,
		const T	indentation,
		const T	overlapVolume,
		const Vector< T, PARTICLETYPE::d > &	relVel,
		const T	dampingFactor,
		const CONTACTPROPERTIES &	contactProperties,
		T	k,
		const std::size_t &	particleAID,
		const std::size_t &	particleBID,
		F	processContactForce = defaults::processContactForce<T, PARTICLETYPE::d> )

Definition at line 403 of file particleContactForceFunctions.h.

{
  using namespace descriptors;
  constexpr unsigned D = PARTICLETYPE::d;
 
  const unsigned materialA =
      particleA.template getField<MECHPROPERTIES, MATERIAL>();
  const unsigned materialB =
      particleB.template getField<MECHPROPERTIES, MATERIAL>();
 
  // Calculation of contact force
  const Vector<T, D> normalForce = getNormalForceFromOverlapVolume(
      contactNormal, indentation, overlapVolume,
      contactProperties.getEffectiveYoungsModulus(materialA, materialB),
      dampingFactor, k, relVel);
  const Vector<T, D> tangentialForce = getTangentialForceFromOverlapVolume(
      contactNormal,
      contactProperties.getStaticFrictionCoefficient(materialA, materialB),
      contactProperties.getKineticFrictionCoefficient(materialA, materialB),
      contactProperties.getStaticFrictionCoefficient(materialA, materialB), k,
      relVel, normalForce);
 
  processContactForce(contactPoint, contactNormal, normalForce, tangentialForce,
                      std::array<std::size_t, 2>({particleAID, particleBID}),
                      overlapVolume, indentation, false);
 
  // Applying forces on the particle
  applyForceOnParticle(dataMap, particleA, particleSystem, contactPoint,
                       normalForce, tangentialForce);
  applyForceOnParticle(dataMap, particleB, particleSystem, contactPoint,
                       -1 * normalForce, -1 * tangentialForce);
}

References applyForceOnParticle(), getNormalForceFromOverlapVolume(), and getTangentialForceFromOverlapVolume().

Here is the call graph for this function:

Here is the caller graph for this function:

applyForceOnParticle()

template<typename T , typename PARTICLETYPE >

void olb::particles::contact::applyForceOnParticle	(	std::multimap< int, std::unique_ptr< std::uint8_t[]> > &	dataMap,
		Particle< T, PARTICLETYPE > &	particle,
		XParticleSystem< T, PARTICLETYPE > &	particleSystem,
		const Vector< T, PARTICLETYPE::d > &	contactPoint,
		const Vector< T, PARTICLETYPE::d > &	normalForce,
		const Vector< T, PARTICLETYPE::d > &	tangentialForce )

Definition at line 310 of file particleContactForceFunctions.h.

{
  using namespace descriptors;
  constexpr unsigned D            = PARTICLETYPE::d;
  constexpr unsigned Drot         = utilities::dimensions::convert<D>::rotation;
  Vector<T, D>       contactForce = normalForce + tangentialForce;
  const PhysR<T, D>  position = particle.template getField<GENERAL, POSITION>();
  const PhysR<T, D>  contactPointDistance = contactPoint - position;
  Vector<T, Drot>    torqueFromContact(
      particles::dynamics::torque_from_force<D, T>::calculate(
          contactForce, contactPointDistance));
 
  if (access::isContactComputationEnabled(particle)) {
    bool applyNow = true;
 
    if constexpr (access::providesParallelization<PARTICLETYPE>()) {
      auto& loadBalancer = particleSystem.getSuperStructure().getLoadBalancer();
      const std::size_t globalParticleIC =
          particle.template getField<PARALLELIZATION, IC>();
      int responsibleRank = loadBalancer.rank(globalParticleIC);
      if (responsibleRank != singleton::mpi().getRank()) {
        applyNow = false;
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
        extendContactTreatmentResultsDataMap(globalParticleID, contactForce,
                                             torqueFromContact, responsibleRank,
                                             dataMap);
      }
    }
 
    if (applyNow) {
      if constexpr (access::providesParallelization<PARTICLETYPE>()) {
        // Apply force only to the main particle
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
        particles::communication::forParticlesInSuperParticleSystem<
            T, PARTICLETYPE, conditions::valid_particles>(
            particleSystem,
            [&](Particle<T, PARTICLETYPE>&       particle,
                ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
              const int globalParticleIC =
                  particle.template getField<PARALLELIZATION, IC>();
              const std::size_t currGlobalParticleID =
                  particle.template getField<PARALLELIZATION, ID>();
              if (currGlobalParticleID == globalParticleID &&
                  globiC == globalParticleIC) {
                // Applying contact force on the particle
                const Vector<T, D> force =
                    particle.template getField<FORCING, FORCE>();
                particle.template setField<FORCING, FORCE>(force +
                                                           contactForce);
 
                // Torque calculation and application
                const Vector<T, Drot> torque(
                    particle.template getField<FORCING, TORQUE>());
                particle.template setField<FORCING, TORQUE>(
                    utilities::dimensions::convert<D>::serialize_rotation(
                        torque + torqueFromContact));
#ifdef VERBOSE_CONTACT_COMMUNICATION
                std::cout << "Rank " << singleton::mpi().getRank()
                          << " procsses data of particle " << globalParticleID
                          << " (contactForce = " << contactForce + force
                          << "; torque = " << torqueFromContact + torque << ")"
                          << std::endl;
#endif
              }
            });
      }
      else {
        // Applying contact force on the particle
        const Vector<T, D> force = particle.template getField<FORCING, FORCE>();
        particle.template setField<FORCING, FORCE>(force + contactForce);
 
        // Torque calculation and application
        const Vector<T, Drot> torque(
            particle.template getField<FORCING, TORQUE>());
        particle.template setField<FORCING, TORQUE>(
            utilities::dimensions::convert<D>::serialize_rotation(
                torque + torqueFromContact));
      }
    }
  }
}

References extendContactTreatmentResultsDataMap(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::singleton::MpiManager::getRank(), olb::particles::access::isContactComputationEnabled(), and olb::singleton::mpi().

Here is the call graph for this function:

Here is the caller graph for this function:

calcElasticAndViscousForce()

template<typename T , unsigned D>

Vector< T, D > olb::particles::contact::calcElasticAndViscousForce	(	T	effectiveE,
		T	k,
		T	overlapVolume,
		T	indentation,
		T	dampingFactor,
		const Vector< T, D > &	relVelNormal,
		const Vector< T, D > &	contactNormal )

Elastic and viscous force from overlap volume according to Nassauer & Kuna (2013)

Definition at line 208 of file particleContactForceFunctions.h.

{
  // Sum is necessary to get the sign of the damping since velocity and normal have the exact same/opposite direction
  T sum = relVelNormal[0] * contactNormal[0];
  for (unsigned iD = 1; iD < D; ++iD) {
    sum += relVelNormal[iD] * contactNormal[iD];
  }
 
  const T forceMagnitude =
      effectiveE * k * util::sqrt(overlapVolume * indentation) *
      (1 + dampingFactor * util::sign(-sum) * norm(relVelNormal));
  return util::max(forceMagnitude, T {0}) * contactNormal;
}

References olb::util::max(), olb::norm(), olb::util::sign(), and olb::util::sqrt().

Here is the call graph for this function:

Here is the caller graph for this function:

cleanContacts()

template<typename CONTACTTYPE >

void olb::particles::contact::cleanContacts

(

std::vector< CONTACTTYPE > &

contacts

)

Definition at line 107 of file contactFunctions.h.

{
  contacts.erase(std::remove_if(contacts.begin(), contacts.end(),
                                [](const CONTACTTYPE& contact) {
                                  return contact.isEmpty();
                                }),
                 contacts.end());
}

Here is the caller graph for this function:

communicateContacts()

template<typename T , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >

void olb::particles::contact::communicateContacts

(

ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &

contactContainer

)

Definition at line 38 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  OstreamManager clout(std::cout, "contactCommunication");
 
  if (singleton::mpi().getSize() > 1) {
    int sumParticleContact = 0;
    int sumWallContact     = 0;
    int countParticleContactsPerRank[singleton::mpi().getSize()];
    int countWallContactsPerRank[singleton::mpi().getSize()];
    int localCountParticleContacts = 0;
    int localCountWallContacts     = 0;
 
    //clout << "communicate number of contacts ..." << std::endl;
    if (singleton::mpi().isMainProcessor()) {
      countParticleContactsPerRank[singleton::mpi().bossId()] =
          contactContainer.particleContacts.size();
      sumParticleContact +=
          countParticleContactsPerRank[singleton::mpi().bossId()];
      countWallContactsPerRank[singleton::mpi().bossId()] =
          contactContainer.wallContacts.size();
      sumWallContact += countWallContactsPerRank[singleton::mpi().bossId()];
      for (int rank = 0; rank < singleton::mpi().getSize(); ++rank) {
        if (rank != singleton::mpi().bossId()) {
          singleton::mpi().receive(&countParticleContactsPerRank[rank], 1,
                                   rank);
          sumParticleContact += countParticleContactsPerRank[rank];
          singleton::mpi().receive(&countWallContactsPerRank[rank], 1, rank);
          sumWallContact += countWallContactsPerRank[rank];
        }
      }
    }
    else {
      localCountParticleContacts =
          contactContainer.particleContacts
              .size(); // size of std::vector (called contacts in container)
      singleton::mpi().send(&localCountParticleContacts, 1,
                            singleton::mpi().bossId()); // 0 is rank of master
      localCountWallContacts = contactContainer.wallContacts.size();
      singleton::mpi().send(&localCountWallContacts, 1,
                            singleton::mpi().bossId());
    }
    //clout << "communicate number of contacts ... OK" << std::endl;
 
    //clout << "broadcasting number of contacts ..." << std::endl;
    singleton::mpi().bCast(&sumParticleContact, 1, singleton::mpi().bossId());
    singleton::mpi().bCast(&sumWallContact, 1, singleton::mpi().bossId());
    //std::cout << sumParticleContact << std::endl;
    //singleton::mpi().barrier();
    //clout << "broadcasting number of contacts ... OK" << std::endl;
 
    if (sumParticleContact > 0 || sumWallContact > 0) {
      ContactContainer<T, PARTICLECONTACTTYPE, WALLCONTACTTYPE> container(
          sumParticleContact, sumWallContact);
      int particleContactObjSize = sizeof(PARTICLECONTACTTYPE);
      int wallContactObjSize     = sizeof(WALLCONTACTTYPE);
 
      //clout << "Bundling contact objects ..." << std::endl;
      if (singleton::mpi().isMainProcessor()) {
        PARTICLECONTACTTYPE* currParticle = container.particleContacts.data();
        WALLCONTACTTYPE*     currWall     = container.wallContacts.data();
        for (int rank = 0; rank < singleton::mpi().getSize(); ++rank) {
          std::uint8_t* bytesParticle = (std::uint8_t*)currParticle;
          std::uint8_t* bytesWall     = (std::uint8_t*)currWall;
          if (rank != singleton::mpi().bossId()) {
            singleton::mpi().receive(bytesParticle,
                                     countParticleContactsPerRank[rank] *
                                         particleContactObjSize,
                                     rank);
            singleton::mpi().receive(
                bytesWall, countWallContactsPerRank[rank] * wallContactObjSize,
                rank);
          }
          else {
            std::memcpy(
                bytesParticle,
                (std::uint8_t*)contactContainer.particleContacts.data(),
                countParticleContactsPerRank[singleton::mpi().bossId()] *
                    particleContactObjSize);
            std::memcpy(bytesWall,
                        (std::uint8_t*)contactContainer.wallContacts.data(),
                        countWallContactsPerRank[singleton::mpi().bossId()] *
                            wallContactObjSize);
          }
          currParticle += countParticleContactsPerRank[rank];
          currWall += countWallContactsPerRank[rank];
        }
      }
      else {
        PARTICLECONTACTTYPE* currParticle =
            contactContainer.particleContacts.data();
        WALLCONTACTTYPE* currWall      = contactContainer.wallContacts.data();
        std::uint8_t*    bytesParticle = (std::uint8_t*)currParticle;
        std::uint8_t*    bytesWall     = (std::uint8_t*)currWall;
        singleton::mpi().send(
            bytesParticle, localCountParticleContacts * particleContactObjSize,
            singleton::mpi().bossId());
        singleton::mpi().send(bytesWall,
                              localCountWallContacts * wallContactObjSize,
                              singleton::mpi().bossId());
      }
      //clout << "Bundling contact objects ... OK" << std::endl;
 
      if (singleton::mpi().isMainProcessor()) {
        container.combineContacts();
      }
 
      //clout << "Broadcasting contact objects ..." << std::endl;
      std::uint8_t* containerParticlesBytes =
          (std::uint8_t*)container.particleContacts.data();
      singleton::mpi().bCast(containerParticlesBytes,
                             sumParticleContact * particleContactObjSize,
                             singleton::mpi().bossId());
      std::uint8_t* containerWallBytes =
          (std::uint8_t*)container.wallContacts.data();
      singleton::mpi().bCast(containerWallBytes,
                             sumWallContact * wallContactObjSize,
                             singleton::mpi().bossId());
      //singleton::mpi().barrier();
      //clout << "Broadcasting contact objects ... OK" << std::endl;
      /*
            clout << "#######" << std::endl;
            clout << "id: " << container.wallContacts[3].id << std::endl;
            clout << "#######" << std::endl;
      */
      //std::cout << "size before: " << container.particleContacts.size() << std::endl;
      container.cleanContacts();
      //std::cout << "size after: " << container.particleContacts.size() << std::endl;
      contactContainer = container;
    }
  }
#endif
}

References olb::singleton::MpiManager::bCast(), olb::singleton::MpiManager::bossId(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::cleanContacts(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::combineContacts(), olb::singleton::MpiManager::getSize(), olb::singleton::MpiManager::isMainProcessor(), olb::singleton::mpi(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::particleContacts, olb::singleton::MpiManager::receive(), olb::singleton::MpiManager::send(), and olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::wallContacts.

Here is the call graph for this function:

Here is the caller graph for this function:

communicateContactsParallel()

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >

void olb::particles::contact::communicateContactsParallel	(	ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		MPI_Comm	particleContactDetectionComm,
		MPI_Comm	wallContactDetectionComm,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

Definition at line 743 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  OstreamManager clout(std::cout, "contactCommunication");
 
  communicateParallelContacts(contactContainer.wallContacts, sParticleSystem,
                              deltaX, wallContactDetectionComm, periodicity);
 
  communicateParallelContacts(contactContainer.particleContacts,
                              sParticleSystem, deltaX,
                              particleContactDetectionComm, periodicity);
#endif
}

References communicateParallelContacts(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::particleContacts, and olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::wallContacts.

Here is the call graph for this function:

communicateContactTreatmentResults()

template<typename T , typename PARTICLETYPE >

void olb::particles::contact::communicateContactTreatmentResults	(	SuperParticleSystem< T, PARTICLETYPE > &	particleSystem,
		std::multimap< int, std::unique_ptr< std::uint8_t[]> > &	dataMap,
		MPI_Comm	contactTreatmentComm )

Definition at line 1052 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  if constexpr (access::providesParallelization<PARTICLETYPE>()) {
 
    std::size_t serialSize =
        getContactTreatmentResultsSerialSize<T, PARTICLETYPE::d>();
 
    //Create non blocking mpi helper
    singleton::MpiNonBlockingHelper mpiNbHelper;
 
    std::map<int, std::vector<std::uint8_t>> rankDataMapSorted;
    communication::fillSendBuffer(dataMap, rankDataMapSorted, serialSize);
 
    // Send mapped data
    communication::sendMappedData(
        rankDataMapSorted, particleSystem.getNeighbourRanks(), serialSize,
        contactTreatmentComm, mpiNbHelper);
 
    receiveContactTreatmentResults(particleSystem, contactTreatmentComm,
                                   mpiNbHelper);
  }
#endif
}

References olb::particles::communication::fillSendBuffer(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getNeighbourRanks(), receiveContactTreatmentResults(), and olb::particles::communication::sendMappedData().

Here is the call graph for this function:

communicateParallelContacts()

template<typename T , typename PARTICLETYPE , typename CONTACTTYPE >

void olb::particles::contact::communicateParallelContacts	(	std::vector< CONTACTTYPE > &	contacts,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		MPI_Comm	contactDetectionComm,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

Definition at line 667 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  OstreamManager clout(std::cout, "contactCommunication");
 
  // Retrieve rank of direct neighbours
  auto& listNeighbourRanks         = sParticleSystem.getNeighbourRanks();
  constexpr std::size_t serialSize = CONTACTTYPE::getSerialSize();
 
  std::multimap<int, std::unique_ptr<std::uint8_t[]>> dataMap;
 
  for (CONTACTTYPE& contact : contacts) {
    // Reset responsible rank because we evaluate the new one in the next step
    resetResponsibleRank(contact);
 
    if (!contact.isEmpty()) {
      std::unordered_set<int> destRanks {evalDestRanksForDetectionCommunication(
          contact, sParticleSystem, deltaX, periodicity)};
 
      // Only new contacts need communication because they originate from the on-lattice
      // contact detection. Known contacts are not updated during that step.
      if (particles::contact::isResponsibleRankSet(contact)
          // not new contacts need to be communicated as the responsible processors for each particle must know the responsible rank for contact treatment, which is evaluated above
          && (contact.isNew() ||
              contact.getResponsibleRank() == singleton::mpi().getRank())) {
        for (const int destRank : destRanks) {
          if (destRank != singleton::mpi().getRank()) {
            std::unique_ptr<std::uint8_t[]> buffer(
                new std::uint8_t[serialSize] {});
            contact.serialize(buffer.get());
            dataMap.insert(std::make_pair(destRank, std::move(buffer)));
          }
        }
      }
 
      // Reset contact on other ranks to avoid unnecessary communication
      //if (destRanks.find(singleton::mpi().getRank()) == destRanks.end()) {
      //  contact.resetMinMax();
      //}
    }
  }
 
  //Create non blocking mpi helper
  singleton::MpiNonBlockingHelper mpiNbHelper;
 
  std::map<int, std::vector<std::uint8_t>> rankDataMapSorted;
  communication::fillSendBuffer(dataMap, rankDataMapSorted, serialSize);
 
  //Send mapped data
  communication::sendMappedData(rankDataMapSorted, listNeighbourRanks,
                                serialSize, contactDetectionComm, mpiNbHelper);
 
  //Receive data and iterate buffer
  communication::receiveAndExecuteForData(
      listNeighbourRanks, serialSize, contactDetectionComm, mpiNbHelper,
      [&](int rankOrig, std::uint8_t* buffer) {
        CONTACTTYPE newcontact;
        newcontact.deserialize(buffer);
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
        std::cout << "From Rank " << rankOrig << ": ";
#endif
 
        receiveContact(newcontact, sParticleSystem, contacts, rankOrig);
      });
#endif
}

References evalDestRanksForDetectionCommunication(), olb::particles::communication::fillSendBuffer(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getNeighbourRanks(), olb::singleton::MpiManager::getRank(), isResponsibleRankSet(), olb::singleton::mpi(), olb::particles::communication::receiveAndExecuteForData(), receiveContact(), resetResponsibleRank(), and olb::particles::communication::sendMappedData().

Here is the call graph for this function:

Here is the caller graph for this function:

communicatePostContactTreatmentContacts() [1/2]

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE >

void olb::particles::contact::communicatePostContactTreatmentContacts	(	ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		MPI_Comm	particlePostContactTreatmentComm,
		MPI_Comm	wallPostContactTreatmentComm,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

It is necessary to communicate contacts so that no information is lost.

Definition at line 832 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  OstreamManager clout(std::cout, "postContactTreatmentCommunication");
 
  communicatePostContactTreatmentContacts(
      contactContainer.wallContacts, sParticleSystem, deltaX,
      wallPostContactTreatmentComm, periodicity);
 
  communicatePostContactTreatmentContacts(
      contactContainer.particleContacts, sParticleSystem, deltaX,
      particlePostContactTreatmentComm, periodicity);
#endif
}

References communicatePostContactTreatmentContacts(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::particleContacts, and olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::wallContacts.

Here is the call graph for this function:

communicatePostContactTreatmentContacts() [2/2]

template<typename T , typename PARTICLETYPE , typename CONTACTTYPE >

void olb::particles::contact::communicatePostContactTreatmentContacts	(	std::vector< CONTACTTYPE > &	contacts,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		MPI_Comm	postContactTreatmentComm,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

Definition at line 764 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  OstreamManager clout(std::cout, "postContactTreatmentCommunication");
 
  // Retrieve rank of direct neighbours
  auto& listNeighbourRanks         = sParticleSystem.getNeighbourRanks();
  constexpr std::size_t serialSize = CONTACTTYPE::getSerialSize();
 
  std::multimap<int, std::unique_ptr<std::uint8_t[]>> dataMap;
 
  for (CONTACTTYPE& contact : contacts) {
    if (contact.getResponsibleRank() == singleton::mpi().getRank()) {
      if (!contact.isEmpty()) {
        std::unordered_set<int> destRanks {
            evalDestRanksForPostContactTreatmentCommunication(
                contact, sParticleSystem, deltaX, periodicity)};
        resetResponsibleRank(contact);
        for (const int destRank : destRanks) {
          if (destRank != singleton::mpi().getRank()) {
            std::unique_ptr<std::uint8_t[]> buffer(
                new std::uint8_t[serialSize] {});
            contact.serialize(buffer.get());
            dataMap.insert(std::make_pair(destRank, std::move(buffer)));
          }
        }
      }
    }
    // Otherwise reset min and max as they are not accurate anymore
    else {
      contact.resetMinMax();
    }
  }
 
  //Create non blocking mpi helper
  singleton::MpiNonBlockingHelper mpiNbHelper;
 
  std::map<int, std::vector<std::uint8_t>> rankDataMapSorted;
  communication::fillSendBuffer(dataMap, rankDataMapSorted, serialSize);
 
  //Send mapped data
  communication::sendMappedData(rankDataMapSorted, listNeighbourRanks,
                                serialSize, postContactTreatmentComm,
                                mpiNbHelper);
 
  // Remove all empty contacts (= all contacts that were empty and all contacts where the process wasn't responsible) to ensure that 100% of the data send by the responsible rank is used
  cleanContacts(contacts);
 
  //Receive data and iterate buffer
  communication::receiveAndExecuteForData(
      listNeighbourRanks, serialSize, postContactTreatmentComm, mpiNbHelper,
      [&](int rankOrig, std::uint8_t* buffer) {
        CONTACTTYPE newcontact;
        newcontact.deserialize(buffer);
        receiveContact(newcontact, sParticleSystem, contacts, rankOrig);
      });
#endif
}

References cleanContacts(), evalDestRanksForPostContactTreatmentCommunication(), olb::particles::communication::fillSendBuffer(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getNeighbourRanks(), olb::singleton::MpiManager::getRank(), olb::singleton::mpi(), olb::particles::communication::receiveAndExecuteForData(), receiveContact(), resetResponsibleRank(), and olb::particles::communication::sendMappedData().

Here is the call graph for this function:

Here is the caller graph for this function:

correctBoundingBox()

template<typename T , unsigned D, typename F1 , typename F2 , typename F3 >

void olb::particles::contact::correctBoundingBox	(	const PhysR< T, D > &	min,
		const PhysR< T, D > &	max,
		const T	physDeltaX,
		const unsigned	contactBoxResolutionPerDirection,
		F1	isInsideContact,
		F2	resetContactBox,
		F3	updateMinMax )

Definition at line 731 of file particleContactForceFunctions.h.

{
  const Vector<T, D> contactPhysDeltaX =
      evalContactDeltaX(min, max, physDeltaX, contactBoxResolutionPerDirection);
  Vector<long, D> startPos;
  Vector<long, D> endPos;
  evalStartAndEndPos(startPos, endPos, max, min, contactPhysDeltaX);
  resetContactBox();
 
  Vector<long, D> origin;
  const bool      originFound = evalStartingPoint(
      startPos, endPos, contactPhysDeltaX, origin, isInsideContact);
 
  if (originFound) {
    PhysR<T, D> physOrigin;
    for (unsigned iD = 0; iD < D; ++iD) {
      physOrigin[iD] = contactPhysDeltaX[iD] * origin[iD];
    }
    updateMinMax(physOrigin);
 
    for (unsigned i = 0; i < utilities::dimensions::convert<D>::neighborsCount;
         ++i) {
      const Vector<short, D> direction =
          utilities::dimensions::convert<D>::neighborDirections[i];
      evalBoundingBoxRecursive(startPos, endPos, contactPhysDeltaX,
                               origin + direction, direction, isInsideContact,
                               updateMinMax);
    }
  }
}

References evalBoundingBoxRecursive(), evalContactDeltaX(), evalStartAndEndPos(), evalStartingPoint(), and updateMinMax().

Here is the call graph for this function:

Here is the caller graph for this function:

correctBoundingBoxNewContact()

template<typename T , unsigned D, typename F1 , typename F2 , typename F3 >

void olb::particles::contact::correctBoundingBoxNewContact	(	PhysR< T, D >	min,
		PhysR< T, D >	max,
		T	physDeltaX,
		unsigned	contactBoxResolutionPerDirection,
		F1	signedDistanceToIntersection,
		F2	resetContactBox,
		F3	updateMinMax )

Definition at line 766 of file particleContactForceFunctions.h.

{
  Vector<T, D> contactPhysDeltaX =
      evalContactDeltaX(min, max, physDeltaX, contactBoxResolutionPerDirection);
 
  for (unsigned iD = 0; iD < D; ++iD) {
    if (util::nearZero(min[iD] - max[iD])) {
      min[iD] -= contactPhysDeltaX[iD];
      /* T {0.5} * contactBoxResolutionPerDirection * contactPhysDeltaX[iD]; */
      max[iD] += contactPhysDeltaX[iD];
      /* T {0.5} * contactBoxResolutionPerDirection * contactPhysDeltaX[iD]; */
    }
  }
 
  Vector<long, D> startPos;
  Vector<long, D> endPos;
  evalStartAndEndPos(startPos, endPos, max, min, contactPhysDeltaX);
  resetContactBox();
 
  forEachPosition(startPos, endPos, [&](const Vector<long, D>& pos) {
    bool isOnContactSurface = false;
    for (unsigned iD = 0; iD < D; ++iD) {
      if (pos[iD] == startPos[iD] || pos[iD] == endPos[iD]) {
        isOnContactSurface = true;
        break;
      }
    }
 
    if (isOnContactSurface) {
      // Determine physical position
      PhysR<T, D> physPos;
      for (unsigned iD = 0; iD < D; ++iD) {
        physPos[iD] = pos[iD] * contactPhysDeltaX[iD];
      }
 
      for (unsigned i = 0;
           i < utilities::dimensions::convert<D>::neighborsCount; ++i) {
        // do not regard direction if the neighbor is also on the surface
        bool regardDirection = true;
        //const Vector<long,D> neighborPos = pos + direction;
        const Vector<long, D> neighborPos =
            pos + Vector<long, D>(
                      utilities::dimensions::convert<D>::neighborDirections[i]);
        for (unsigned iD = 0; iD < D; ++iD) {
          if (neighborPos[iD] <= startPos[iD] ||
              neighborPos[iD] >= endPos[iD]) {
            regardDirection = false;
            break;
          }
        }
 
        if (regardDirection) {
          const PhysR<T, D> normalizedDirection =
              PhysR<T, D>(utilities::dimensions::convert<
                          D>::template normalizedNeighborDirections<T>[i]);
          T distance;
 
          // TODO: the sdf from intersections is not exact - replace with another option
          // Only change min & max if boundary was found
          /* if (util::distance( */
          /*         distance, physPos, direction, physDeltaX * T {1e-2}, */
          /*         [&](const Vector<T, D>& pos) { */
          /*           return signedDistanceToIntersection(pos); */
          /*         }, */
          /*         [&](const Vector<T, D>& pos) { */
          /*           return signedDistanceToIntersection(pos) < physDeltaX; */
          /*         })) { */
          if (util::distance<T, D, false>(
                  distance, physPos, normalizedDirection, physDeltaX * T {1e-2},
                  [&](const Vector<T, D>& pos) {
                    return signedDistanceToIntersection(pos);
                  },
                  physDeltaX)) {
            // Determine position of boundary in given direction
            const PhysR<T, D> posOnBoundary =
                physPos + distance * normalizedDirection;
 
            // Attempt to update min and max with found position on contact boundary
            updateMinMax(posOnBoundary);
          }
        }
      }
    }
  });
}

References evalContactDeltaX(), evalStartAndEndPos(), forEachPosition(), olb::util::nearZero(), and updateMinMax().

Here is the call graph for this function:

Here is the caller graph for this function:

createContactProperties()

template<typename T , unsigned N, typename F >

constexpr ContactProperties< T, N > olb::particles::contact::createContactProperties ( F f )

constexpr

Definition at line 98 of file contactProperties.h.

{
  ContactProperties<T, N> contactProperties {};
  f(contactProperties);
  return contactProperties;
}

domainDimension()

template<typename T , unsigned D>

unsigned short olb::particles::contact::domainDimension	(	PhysR< T, D > &	min,
		PhysR< T, D > &	max )

Definition at line 44 of file contactHelpers.h.

{
  unsigned short dimension=0;
  for (unsigned iD = 0; iD < D; ++iD) {
    if(max[iD]>min[iD]){
      dimension+=1;
    }
  }
  return dimension;
}

evalApproxSurface()

template<typename T , unsigned D>

T olb::particles::contact::evalApproxSurface	(	const Vector< T, D > &	normalizedNormal,
		const PhysR< T, D > &	contactPhysDeltaX )

Takes a normalized normal and the voxel sizes and approximates the corresponding surface.

Definition at line 227 of file particleContactForceFunctions.h.

{
  T approxSurface;
  if constexpr (D == 3) {
    approxSurface =
        contactPhysDeltaX[1] * contactPhysDeltaX[2] *
        /* util::abs(normalizedNormal[0]); // * normalizedNormal[0]; */
        normalizedNormal[0] * normalizedNormal[0];
    approxSurface +=
        contactPhysDeltaX[0] * contactPhysDeltaX[2] *
        /* util::abs(normalizedNormal[1]); // * normalizedNormal[1]; */
        normalizedNormal[1] * normalizedNormal[1];
    approxSurface +=
        contactPhysDeltaX[0] * contactPhysDeltaX[1] *
        /* util::abs(normalizedNormal[2]); // * normalizedNormal[2]; */
        normalizedNormal[2] * normalizedNormal[2];
  }
  else {
    approxSurface =
        contactPhysDeltaX[0] *
        /* util::abs(normalizedNormal[1]); // * normalizedNormal[1]; */
        normalizedNormal[1] * normalizedNormal[1];
    approxSurface +=
        contactPhysDeltaX[1] *
        /* util::abs(normalizedNormal[0]); // * normalizedNormal[0]; */
        normalizedNormal[0] * normalizedNormal[0];
  }
  return approxSurface;
  /* return approxSurface / norm(contactPhysDeltaX); */
}

Here is the caller graph for this function:

evalBoundingBoxRecursive()

template<typename T , unsigned D, typename F1 , typename F2 >

void olb::particles::contact::evalBoundingBoxRecursive	(	const Vector< long, D > &	startPos,
		const Vector< long, D > &	endPos,
		const Vector< T, D > &	contactPhysDeltaX,
		const Vector< long, D > &	currPos,
		const Vector< short, D > &	dirIn,
		F1	isInsideContact,
		F2	updateMinMax )

Definition at line 687 of file particleContactForceFunctions.h.

{
  PhysR<T, D> currPhysPos;
  for (unsigned iD = 0; iD < D; ++iD) {
    currPhysPos[iD] = contactPhysDeltaX[iD] * currPos[iD];
  }
  bool isInside   = isInsideContact(currPhysPos);
  bool isInsideBB = true;
  for (unsigned iD = 0; iD < D; ++iD) {
    if (currPos[iD] < startPos[iD] || currPos[iD] > endPos[iD]) {
      isInsideBB = false;
      break;
    }
  }
 
  // if point is inside, then we check the neighboring points too, otherwise they cannot be inside for convex shapes
  if (isInside || isInsideBB) {
    if (isInside) {
      updateMinMax(currPhysPos);
    }
 
    forEachPosition(Vector<unsigned short, D>(0),
                    Vector<unsigned short, D>(abs(dirIn)),
                    [&](const Vector<unsigned short, D>& dirOut) {
                      if (norm(dirOut) >= 1) {
                        Vector<short, D> signedDirOut;
                        for (unsigned iD = 0; iD < D; ++iD) {
                          signedDirOut[iD] = dirOut[iD] * dirIn[iD];
                        }
                        const Vector<long, D> nextPos = signedDirOut + currPos;
                        evalBoundingBoxRecursive(
                            startPos, endPos, contactPhysDeltaX, nextPos,
                            signedDirOut, isInsideContact, updateMinMax);
                      }
                    });
  }
}

References olb::abs(), evalBoundingBoxRecursive(), forEachPosition(), olb::norm(), and updateMinMax().

Here is the call graph for this function:

Here is the caller graph for this function:

evalCircumRadius() [1/3]

template<typename T , typename PARTICLETYPE >

T olb::particles::contact::evalCircumRadius	(	Particle< T, PARTICLETYPE > &	particle,
		T const	physDeltaX )

Definition at line 98 of file contactFunctions.h.

{
  using namespace descriptors;
  auto sIndicator = particle.template getField<SURFACE, SINDICATOR>();
  return evalCircumRadius(particle, physDeltaX, sIndicator->getCircumRadius(),
                          sIndicator->getEpsilon());
}

References evalCircumRadius().

Here is the call graph for this function:

evalCircumRadius() [2/3]

template<typename T , typename PARTICLETYPE >

T olb::particles::contact::evalCircumRadius	(	Particle< T, PARTICLETYPE > &	particle,
		T const	physDeltaX,
		T const	circumRadius,
		T const	epsilon )

Definition at line 89 of file contactFunctions.h.

{
  const T contactDetectionDistance =
      evalContactDetectionDistance(particle, physDeltaX);
  return evalCircumRadius(contactDetectionDistance, circumRadius, epsilon);
}

References evalCircumRadius(), and evalContactDetectionDistance().

Here is the call graph for this function:

evalCircumRadius() [3/3]

template<typename T >

T olb::particles::contact::evalCircumRadius	(	T const	contactDetectionDistance,
		T const	circumRadius,
		T const	epsilon )

Definition at line 82 of file contactFunctions.h.

{
  return circumRadius - epsilon + util::max(epsilon, contactDetectionDistance);
}

References olb::util::max().

Here is the call graph for this function:

Here is the caller graph for this function:

evalContactDeltaX()

template<typename T , unsigned D>

constexpr Vector< T, D > olb::particles::contact::evalContactDeltaX ( const Vector< T, D > & min, const Vector< T, D > & max, const T physDeltaX, const unsigned contactBoxResolutionPerDirection )

constexpr

Definition at line 148 of file particleContactForceFunctions.h.

{
  // Adaptiv resolution of found contact
  const Vector<T, D>  boxSize           = max - min;
  Vector<T, D>        contactPhysDeltaX = boxSize;
  std::array<bool, D> fallbackUsed {[]() {
    static_assert(D == 2 || D == 3, "Only D=2 and D=3 are supported");
    if constexpr (D == 3) {
      return std::array<bool, D> {false, false, false};
    }
    else {
      return std::array<bool, D> {false, false};
    }
  }()};
 
  for (unsigned iD = 0; iD < D; ++iD) {
    contactPhysDeltaX[iD] /= contactBoxResolutionPerDirection;
    if (util::nearZero(contactPhysDeltaX[iD])) {
      contactPhysDeltaX[iD] = physDeltaX / contactBoxResolutionPerDirection;
      fallbackUsed[iD]      = true;
    }
  }
 
  // If fallback was used, check if in another dimension the deltaX is known.
  // If yes, then use the smallest known deltaX. Otherwise, keep the fallback.
  // This is possible because any dimenison for which we cannot find a proper deltaX has probably even smaller lenghts than the others.
  for (unsigned iD = 0; iD < D; ++iD) {
    if (fallbackUsed[iD]) {
      for (unsigned iD2 = 0; iD2 < D; ++iD2) {
        if (iD != iD2) {
          contactPhysDeltaX[iD] =
              util::min(contactPhysDeltaX[iD], contactPhysDeltaX[iD2]);
        }
      }
    }
  }
 
  return contactPhysDeltaX;
}

References olb::util::min(), and olb::util::nearZero().

Here is the call graph for this function:

Here is the caller graph for this function:

evalContactDetectionDistance()

template<typename T , typename PARTICLETYPE >

T olb::particles::contact::evalContactDetectionDistance	(	Particle< T, PARTICLETYPE > &	particle,
		T const	physDeltaX )

Definition at line 59 of file contactFunctions.h.

{
  constexpr unsigned D = PARTICLETYPE::d;
 
  if constexpr (particles::access::providesContactMaterial<PARTICLETYPE>()) {
    constexpr T factor = []() {
      static_assert(D == 2 || D == 3, "Only D=2 and D=3 are supported");
      // TODO: Use with c++20
      //return T{0.5} * (D == 3 ? std::numbers::sqrt3_v<T> : std::numbers::sqrt2_v<T>);
      return T {0.5} * (D == 3 ? 1.7320508075688772935 : 1.4142135623730950488);
    }();
    return factor * physDeltaX +
           particles::access::getEnlargementForContact(particle);
  }
  else {
    return T {0};
  }
 
  __builtin_unreachable();
}

References olb::particles::access::getEnlargementForContact().

Here is the call graph for this function:

Here is the caller graph for this function:

evalContactPosition()

template<typename T , typename PARTICLETYPE , typename F >

PhysR< T, PARTICLETYPE::d > olb::particles::contact::evalContactPosition	(	Particle< T, PARTICLETYPE > &	particle,
		const PhysR< T, PARTICLETYPE::d > &	particlePos,
		const PhysR< T, PARTICLETYPE::d > &	contactPos,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		T	deltaX )

Definition at line 162 of file particleContactDetectionFunctions.h.

{
  using namespace descriptors;
 
  if constexpr (!(isPeriodic(getSetupPeriodicity()))) {
    return contactPos;
  }
  else {
    auto    sIndicator   = particle.template getField<SURFACE, SINDICATOR>();
    const T circumRadius = sIndicator->getCircumRadius() -
                           sIndicator->getEpsilon() +
                           util::max(sIndicator->getEpsilon(), 0.5 * deltaX);
 
    PhysR<T, PARTICLETYPE::d> ghostPos = contactPos;
    for (unsigned i = 0; i < PARTICLETYPE::d; ++i) {
      T const particleMin = particlePos[i] - circumRadius;
      T const particleMax = particlePos[i] + circumRadius;
      if (particleMin > contactPos[i] && getSetupPeriodicity()[i]) {
        ghostPos[i] = communication::movePositionToEnd(
            contactPos[i], cellMax[i], cellMin[i]);
      }
      else if (particleMax < contactPos[i] && getSetupPeriodicity()[i]) {
        ghostPos[i] = communication::movePositionToStart(
            contactPos[i], cellMax[i], cellMin[i]);
      }
    }
    return ghostPos;
  }
}

References olb::particles::isPeriodic(), olb::util::max(), olb::particles::communication::movePositionToEnd(), and olb::particles::communication::movePositionToStart().

Here is the call graph for this function:

Here is the caller graph for this function:

evalCurrentDampingFactor()

template<typename T , unsigned D, typename CONTACTTYPE >

constexpr T olb::particles::contact::evalCurrentDampingFactor ( CONTACTTYPE & contact, const T coefficientOfRestitution, const Vector< T, D > & initialRelativeNormalVelocity )

constexpr

Definition at line 52 of file particleContactForceFunctions.h.

{
  if (contact.getDampingFactor() < 0) {
    contact.setDampingFactorFromInitialVelocity(
        coefficientOfRestitution, norm(initialRelativeNormalVelocity));
  }
  return contact.getDampingFactor();
}

References olb::norm().

Here is the call graph for this function:

Here is the caller graph for this function:

evalDampingFactor()

template<typename T >

constexpr T olb::particles::contact::evalDampingFactor ( const T coefficientOfRestitution, const T initialRelativeVelocityMagnitude )

constexpr

Calculates the damping factor according to Carvalho & Martins (2019) (10.1016/j.mechmachtheory.2019.03.028)

Definition at line 46 of file contactFunctions.h.

{
  // This case should never happen, but could for example due to initial conditions
  if (initialRelativeVelocityMagnitude <= 0) {
    return 0;
  }
  return 1.5 * (1 - coefficientOfRestitution) *
         (11 - coefficientOfRestitution) / (1 + 9 * coefficientOfRestitution) /
         initialRelativeVelocityMagnitude;
}

Here is the caller graph for this function:

evalDestRanksForDetectionCommunication() [1/2]

template<typename T , typename PARTICLETYPE , bool CONVEX>

std::unordered_set< int > olb::particles::contact::evalDestRanksForDetectionCommunication	(	ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

evaluate responsible ranks for particle-particle contact

Definition at line 257 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
  auto& superStructure = sParticleSystem.getSuperStructure();
  auto& loadBalancer   = superStructure.getLoadBalancer();
  /* auto& cuboidGeometry = superStructure.getCuboidGeometry(); */
 
  // Set instead of unordered_set so that it's sorted
  // Therefore, it is easier to find the the rank with smallest ID
  std::set<int>                intersectingRanks {};
  std::array<int, 2>           responsibleRank = {-1, -1};
  std::array<std::set<int>, 2> touchedRanks {};
  std::unordered_set<int>      destRanks {};
 
  // Function that is called on failure
  const auto errorTreatment = [&contact, &destRanks]() {
    contact.resetMinMax();
    resetResponsibleRank(contact);
    destRanks.clear();
  };
 
  // Iterate over all particles to populate responsible and touched ranks of each particle in the contact
  communication::forParticlesInSuperParticleSystem<T, PARTICLETYPE,
                                                   conditions::valid_particles>(
      sParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
 
        for (unsigned short index = 0; index < 2; ++index) {
          if (globalParticleID == contact.getIDs()[index]) {
            const int globiCcenter =
                particle.template getField<PARALLELIZATION, IC>();
            responsibleRank[index] = loadBalancer.rank(globiCcenter);
            touchedRanks[index].insert(responsibleRank[index]);
 
            const PhysR<T, PARTICLETYPE::d> position =
                access::getPosition(particle);
            const T circumRadius = evalCircumRadius(particle, deltaX);
            communication::getSurfaceTouchingICs(
                sParticleSystem, position, circumRadius, periodicity,
                globiCcenter, [&](const int iC) {
                  touchedRanks[index].insert(loadBalancer.rank(iC));
                });
          }
        }
      });
 
  // If at least one particle has no touched cuboids then cancel function call
  // (this should be caused by an invalid particle that somehow still has an existing contact)
  if (touchedRanks[0].empty() || touchedRanks[1].empty() ||
      responsibleRank[0] < 0 || responsibleRank[1] < 0) {
    errorTreatment();
  }
  // If both particles have the same responsible rank then let that rank take care of the contact
  else if (responsibleRank[0] == responsibleRank[1]) {
    destRanks.insert(responsibleRank[0]);
    contact.setResponsibleRank(responsibleRank[0]);
  }
  // If not, then determine treatment rank from intersection of ranks that touch the particles
  else {
    std::set_intersection(
        touchedRanks[0].begin(), touchedRanks[0].end(), touchedRanks[1].begin(),
        touchedRanks[1].end(),
        std::inserter(intersectingRanks, intersectingRanks.begin()));
 
    /*
     * It's possible that no intersection is found, if a rank loses responsibility
     * for the contact, but still holds the old information. At the same time,
     * the contact of both particles breaks, and they are also on different ranks,
     * without any processor that knows both. This is very unlikely, but it can
     * happen.
     * However, the following simple if-condition shouldn't be expensive as it is
     * almost always true.
     */
    if (!intersectingRanks.empty()) {
      destRanks.insert(responsibleRank[0]);
      destRanks.insert(responsibleRank[1]);
 
      // Sort responsible ranks so that the following is consistent
      if (responsibleRank[0] > responsibleRank[1]) {
        std::swap(responsibleRank[0], responsibleRank[1]);
      }
 
      // First use one of the responsible ranks if they know both particles
      if (intersectingRanks.find(responsibleRank[0]) !=
          intersectingRanks.end()) {
        contact.setResponsibleRank(responsibleRank[0]);
      }
      else if (intersectingRanks.find(responsibleRank[1]) !=
               intersectingRanks.end()) {
        contact.setResponsibleRank(responsibleRank[1]);
      }
      // Otherwise use smallest rank of intersection
      else {
        const int responsibleRankContact = *intersectingRanks.begin();
        contact.setResponsibleRank(responsibleRankContact);
        destRanks.insert(responsibleRankContact);
      }
    }
    else {
      errorTreatment();
    }
  }
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
  std::cout << "Contact of ids " << contact.getIDs()[0] << " & "
            << contact.getIDs()[1] << " goes to " << destRanks << " and "
            << contact.getResponsibleRank()
            << " is responsible for the calculation (Rank "
            << singleton::mpi().getRank() << ")"
            << " the intersection was " << intersectingRanks
            << " from touched ranks " << touchedRanks[0] << " and "
            << touchedRanks[1] << std::endl;
#endif
 
  return destRanks;
}

References evalCircumRadius(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getIDs(), olb::SuperStructure< T, D >::getLoadBalancer(), olb::particles::access::getPosition(), olb::singleton::MpiManager::getRank(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getResponsibleRank(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getSuperStructure(), olb::particles::communication::getSurfaceTouchingICs(), olb::singleton::mpi(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::resetMinMax(), resetResponsibleRank(), and olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::setResponsibleRank().

Here is the call graph for this function:

Here is the caller graph for this function:

evalDestRanksForDetectionCommunication() [2/2]

template<typename T , typename PARTICLETYPE , bool CONVEX>

std::unordered_set< int > olb::particles::contact::evalDestRanksForDetectionCommunication	(	WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

evaluate responsible rank for solid boundary contact

Definition at line 383 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
  auto& superStructure = sParticleSystem.getSuperStructure();
  auto& loadBalancer   = superStructure.getLoadBalancer();
 
  bool                    isSuccessful = false;
  std::unordered_set<int> destRanks {};
  int                     destRank;
  // Find responsible rank
  communication::forParticlesInSuperParticleSystem<T, PARTICLETYPE,
                                                   conditions::valid_particles>(
      sParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
        if (globalParticleID == contact.getParticleID()) {
          std::size_t globiCcenter =
              particle.template getField<PARALLELIZATION, IC>();
          destRank     = loadBalancer.rank(globiCcenter);
          isSuccessful = true;
        }
      });
 
  if (isSuccessful) {
    contact.setResponsibleRank(destRank);
    if (destRank != singleton::mpi().getRank()) {
      destRanks.insert(destRank);
    }
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
    std::cout << "Contact of id " << contact.getParticleID()
              << " & solid boundary " << contact.getWallID() << " goes to "
              << destRanks << std::endl;
#endif
  }
 
  return destRanks;
}

References olb::particles::communication::forParticlesInSuperParticleSystem(), olb::SuperStructure< T, D >::getLoadBalancer(), olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::getParticleID(), olb::singleton::MpiManager::getRank(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getSuperStructure(), olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::getWallID(), olb::singleton::mpi(), and olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::setResponsibleRank().

Here is the call graph for this function:

evalDestRanksForPostContactTreatmentCommunication() [1/2]

template<typename T , typename PARTICLETYPE , bool CONVEX>

std::unordered_set< int > olb::particles::contact::evalDestRanksForPostContactTreatmentCommunication	(	ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

evaluate ranks that touch both particles of a particle-particle contact

Definition at line 429 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
  auto& superStructure = sParticleSystem.getSuperStructure();
  auto& loadBalancer   = superStructure.getLoadBalancer();
  auto& cuboidGeometry = superStructure.getCuboidGeometry();
 
  std::unordered_set<int>      destRanks;
  std::array<std::set<int>, 2> touchedRanks;
 
  // Iterate over all particles to populate touched ranks of each particle in the contact
  communication::forParticlesInSuperParticleSystem<T, PARTICLETYPE,
                                                   conditions::valid_particles>(
      sParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
 
        for (unsigned short index = 0; index < 2; ++index) {
          if (globalParticleID == contact.getIDs()[index]) {
            const PhysR<T, PARTICLETYPE::d> position =
                access::getPosition(particle);
            // the particle moved, therefore, the iC may have changed
            int        globiCcenter;
            const bool cuboidFound = communication::getCuboid(
                cuboidGeometry, periodicity, position, globiCcenter);
            if (cuboidFound) {
              touchedRanks[index].insert(loadBalancer.rank(globiCcenter));
 
              const T circumRadius = evalCircumRadius(particle, deltaX);
              communication::getSurfaceTouchingICs(
                  sParticleSystem, position, circumRadius, periodicity,
                  globiCcenter, [&](const int iC) {
                    touchedRanks[index].insert(loadBalancer.rank(iC));
                  });
            }
          }
        }
      });
 
  std::set_intersection(touchedRanks[0].begin(), touchedRanks[0].end(),
                        touchedRanks[1].begin(), touchedRanks[1].end(),
                        std::inserter(destRanks, destRanks.begin()));
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
  std::cout << "Contact of ids " << contact.getIDs()[0] << " & "
            << contact.getIDs()[1] << " goes to " << destRanks << " and "
            << contact.getResponsibleRank()
            << " is responsible for the calculation (Rank "
            << singleton::mpi().getRank() << ")" << '\n'
            << "Touched ranks of 1. particle: " << touchedRanks[0]
            << " and touched ranks of 2. particle: " << touchedRanks[1]
            << std::endl;
#endif
 
  return destRanks;
}

References evalCircumRadius(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::particles::communication::getCuboid(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getIDs(), olb::SuperStructure< T, D >::getLoadBalancer(), olb::particles::access::getPosition(), olb::singleton::MpiManager::getRank(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getResponsibleRank(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getSuperStructure(), olb::particles::communication::getSurfaceTouchingICs(), and olb::singleton::mpi().

Here is the call graph for this function:

Here is the caller graph for this function:

evalDestRanksForPostContactTreatmentCommunication() [2/2]

template<typename T , typename PARTICLETYPE , bool CONVEX>

std::unordered_set< int > olb::particles::contact::evalDestRanksForPostContactTreatmentCommunication	(	WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		const T	deltaX,
		const Vector< bool, PARTICLETYPE::d > &	periodicity )

evaluate ranks that touch the particle of a particle-wall contact

Definition at line 494 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
  auto& superStructure = sParticleSystem.getSuperStructure();
  auto& loadBalancer   = superStructure.getLoadBalancer();
  auto& cuboidGeometry = superStructure.getCuboidGeometry();
 
  std::unordered_set<int> destRanks;
 
  // Find responsible rank
  communication::forParticlesInSuperParticleSystem<T, PARTICLETYPE,
                                                   conditions::valid_particles>(
      sParticleSystem,
      [&](Particle<T, PARTICLETYPE>&       particle,
          ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
        const std::size_t globalParticleID =
            particle.template getField<PARALLELIZATION, ID>();
        if (globalParticleID == contact.getParticleID()) {
          const PhysR<T, PARTICLETYPE::d> position =
              access::getPosition(particle);
          // the particle moved, therefore, the iC may have changed
          int        globiCcenter;
          const bool cuboidFound = communication::getCuboid(
              cuboidGeometry, periodicity, position, globiCcenter);
          if (cuboidFound) {
            destRanks.insert(loadBalancer.rank(globiCcenter));
 
            const T circumRadius {evalCircumRadius(particle, deltaX)};
            communication::getSurfaceTouchingICs(
                sParticleSystem, position, circumRadius, periodicity,
                globiCcenter, [&](const int iC) {
                  destRanks.insert(loadBalancer.rank(iC));
                });
          }
        }
      });
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
  std::cout << "Contact of id " << contact.getParticleID()
            << " & solid boundary " << contact.getWallID() << " goes to "
            << destRanks << " and rank " << contact.getResponsibleRank()
            << " is responsible for the calculation" << std::endl;
#endif
 
  return destRanks;
}

References evalCircumRadius(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::particles::communication::getCuboid(), olb::SuperStructure< T, D >::getLoadBalancer(), olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::getParticleID(), olb::particles::access::getPosition(), olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::getResponsibleRank(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getSuperStructure(), olb::particles::communication::getSurfaceTouchingICs(), and olb::particles::contact::WallContactArbitraryFromOverlapVolume< T, D, CONVEX >::getWallID().

Here is the call graph for this function:

evalEffectiveYoungModulus()

template<typename T >

constexpr T olb::particles::contact::evalEffectiveYoungModulus ( T E1, T E2, T nu1, T nu2 )

constexpr

Definition at line 38 of file contactFunctions.h.

{
  const T denominator = (1 - nu1 * nu1) / E1 + (1 - nu2 * nu2) / E2;
  return T {1} / denominator;
}

evalForceViaVolumeCoefficient() [1/2]

template<typename T >

constexpr T olb::particles::contact::evalForceViaVolumeCoefficient ( const T contactVolume, const T contactSurface, const T indentationDepth )

constexpr

Definition at line 522 of file particleContactForceFunctions.h.

{
  constexpr T a = 0.6118229;
  constexpr T b = 1.900093;
  constexpr T c = 0.651237;
 
  const T x = indentationDepth * contactSurface / contactVolume;
  const T k = a + b * util::exp(-c * x);
  return k;
}

References olb::util::exp().

Here is the call graph for this function:

evalForceViaVolumeCoefficient() [2/2]

template<typename T , unsigned D>

T olb::particles::contact::evalForceViaVolumeCoefficient	(	const Vector< T, D > &	min,
		const Vector< T, D > &	max,
		T	overlapVolume )

Definition at line 504 of file particleContactForceFunctions.h.

{
  const T a = util::log(2. / 3.) / util::log(2.);
  const T m = 2. / util::sqrt(M_PI) * util::pow(4. / M_PI, -a);
 
  Vector<T, D> length            = max - min;
  T            boundingBoxVolume = 1;
  for (unsigned iD = 0; iD < D; ++iD) {
    boundingBoxVolume *= length[iD];
  }
 
  const T k =
      m * util::pow(util::max(boundingBoxVolume / overlapVolume, T {1}), a);
  return k;
}

References olb::util::log(), M_PI, olb::util::max(), olb::util::pow(), and olb::util::sqrt().

Here is the call graph for this function:

evalRelativeNormalVelocity()

template<typename T , unsigned D>

constexpr Vector< T, D > olb::particles::contact::evalRelativeNormalVelocity ( const Vector< T, D > & contactNormal, const Vector< T, D > & relVel )

constexpr

Definition at line 83 of file particleContactForceFunctions.h.

{
  return Vector<T, D>((relVel * contactNormal) * contactNormal);
}

Here is the caller graph for this function:

evalRelativeVelocity() [1/2]

template<typename T , typename PARTICLETYPE >

constexpr Vector< T, PARTICLETYPE::d > olb::particles::contact::evalRelativeVelocity ( Particle< T, PARTICLETYPE > & particle, const Vector< T, PARTICLETYPE::d > & position )

constexpr

Definition at line 75 of file particleContactForceFunctions.h.

{
  return particles::dynamics::calculateLocalVelocity(particle, position);
}

References olb::particles::dynamics::calculateLocalVelocity().

Here is the call graph for this function:

evalRelativeVelocity() [2/2]

template<typename T , typename PARTICLETYPE >

constexpr Vector< T, PARTICLETYPE::d > olb::particles::contact::evalRelativeVelocity ( Particle< T, PARTICLETYPE > & particleA, Particle< T, PARTICLETYPE > & particleB, const Vector< T, PARTICLETYPE::d > & position )

constexpr

Returns the relative velocity of two particles at a defined position.

Definition at line 65 of file particleContactForceFunctions.h.

{
  return particles::dynamics::calculateLocalVelocity(particleA, position) -
         particles::dynamics::calculateLocalVelocity(particleB, position);
}

References olb::particles::dynamics::calculateLocalVelocity().

Here is the call graph for this function:

Here is the caller graph for this function:

evalStartAndEndPos()

template<typename T , unsigned D>

void olb::particles::contact::evalStartAndEndPos	(	Vector< long, D > &	startPos,
		Vector< long, D > &	endPos,
		const PhysR< T, D > &	max,
		const PhysR< T, D > &	min,
		const PhysR< T, D > &	contactPhysDeltaX )

Definition at line 493 of file particleContactForceFunctions.h.

{
  for (unsigned iD = 0; iD < D; ++iD) {
    startPos[iD] = util::floor(min[iD] / contactPhysDeltaX[iD]);
    endPos[iD]   = util::ceil(max[iD] / contactPhysDeltaX[iD]);
  }
}

References olb::util::ceil(), and olb::util::floor().

Here is the call graph for this function:

Here is the caller graph for this function:

evalStartingPoint()

template<typename T , unsigned D, typename F1 >

bool olb::particles::contact::evalStartingPoint	(	const Vector< long, D > &	startPos,
		const Vector< long, D > &	endPos,
		const Vector< T, D > &	contactPhysDeltaX,
		Vector< long, D > &	currPos,
		F1	isInsideContact )

Definition at line 663 of file particleContactForceFunctions.h.

{
  bool pointFound = false;
 
  forEachPositionWithBreak(startPos, endPos,
                           [&](const Vector<long, D>& pos, bool& breakLoops) {
                             PhysR<T, D> physPos;
                             for (unsigned iD = 0; iD < D; ++iD) {
                               physPos[iD] = contactPhysDeltaX[iD] * pos[iD];
                             }
                             if (isInsideContact(physPos)) {
                               currPos    = pos;
                               pointFound = true;
                               breakLoops = true;
                             }
                           });
 
  return pointFound;
}

References forEachPositionWithBreak().

Here is the call graph for this function:

Here is the caller graph for this function:

extendContactTreatmentResultsDataMap()

template<typename T , unsigned D>

void olb::particles::contact::extendContactTreatmentResultsDataMap	(	const std::size_t &	globalParticleID,
		Vector< T, D > &	force,
		Vector< T, utilities::dimensions::convert< D >::rotation > &	torque,
		const int	destRank,
		std::multimap< int, std::unique_ptr< std::uint8_t[]> > &	dataMap )

Definition at line 897 of file particleContactCommunicationFunctions.h.

{
#ifdef PARALLEL_MODE_MPI
  // Create communicatables
  std::array<std::size_t, 1> commGlobalParticleID {globalParticleID};
  ConcreteCommunicatable     communicatableID =
      ConcreteCommunicatable(commGlobalParticleID);
  ConcreteCommunicatable communicatableForce  = ConcreteCommunicatable(force);
  ConcreteCommunicatable communicatableTorque = ConcreteCommunicatable(torque);
 
  // Setting dimension dependent indices
  const std::vector<unsigned> indicesID {0};
  std::vector<unsigned>       indicesForce;
  if constexpr (D == 3) {
    indicesForce = std::vector<unsigned> {0, 1, 2};
  }
  else {
    indicesForce = std::vector<unsigned> {0, 1};
  }
  std::vector<unsigned> indicesTorque;
  if constexpr (D == 3) {
    indicesTorque = indicesForce;
  }
  else {
    indicesTorque = std::vector<unsigned> {0};
  }
 
  // Retrieve serial size
  const std::size_t serialSize = communicatableID.size(indicesID) +
                                 communicatableForce.size(indicesForce) +
                                 communicatableTorque.size(indicesTorque);
 
  std::unique_ptr<std::uint8_t[]> buffer(new std::uint8_t[serialSize] {});
  std::uint8_t*                   bufferRaw = buffer.get();
  std::size_t serialIdx = communicatableID.serialize(indicesID, bufferRaw);
  serialIdx +=
      communicatableForce.serialize(indicesForce, &bufferRaw[serialIdx]);
  serialIdx +=
      communicatableTorque.serialize(indicesTorque, &bufferRaw[serialIdx]);
 
  dataMap.insert(std::make_pair(destRank, std::move(buffer)));
#ifdef VERBOSE_CONTACT_COMMUNICATION
  std::cout << "Rank " << singleton::mpi().getRank()
            << " sends contact treatment results of particle "
            << globalParticleID << " to " << destRank << std::endl;
#endif
#endif
}

References olb::singleton::MpiManager::getRank(), olb::singleton::mpi(), olb::ConcreteCommunicatable< COMMUNICATEE >::serialize(), and olb::ConcreteCommunicatable< COMMUNICATEE >::size().

Here is the call graph for this function:

Here is the caller graph for this function:

forEachPosition() [1/2]

template<typename T , typename F >

constexpr void olb::particles::contact::forEachPosition ( Vector< T, 2 > startPos, Vector< T, 2 > endPos, F f )

constexpr

Definition at line 103 of file particleContactForceFunctions.h.

{
  Vector<T, 2> pos;
  for (pos[0] = startPos[0]; pos[0] <= endPos[0]; ++pos[0]) {
    for (pos[1] = startPos[1]; pos[1] <= endPos[1]; ++pos[1]) {
      f(pos);
    }
  }
}

forEachPosition() [2/2]

template<typename T , typename F >

constexpr void olb::particles::contact::forEachPosition ( Vector< T, 3 > startPos, Vector< T, 3 > endPos, F f )

constexpr

Definition at line 90 of file particleContactForceFunctions.h.

{
  Vector<T, 3> pos;
  for (pos[0] = startPos[0]; pos[0] <= endPos[0]; ++pos[0]) {
    for (pos[1] = startPos[1]; pos[1] <= endPos[1]; ++pos[1]) {
      for (pos[2] = startPos[2]; pos[2] <= endPos[2]; ++pos[2]) {
        f(pos);
      }
    }
  }
}

Here is the caller graph for this function:

forEachPositionWithBreak() [1/2]

template<typename T , typename F >

constexpr void olb::particles::contact::forEachPositionWithBreak ( Vector< T, 2 > startPos, Vector< T, 2 > endPos, F f )

constexpr

Definition at line 132 of file particleContactForceFunctions.h.

{
  bool         breakLoops = false;
  Vector<T, 2> pos;
  for (pos[0] = startPos[0]; pos[0] <= endPos[0]; ++pos[0]) {
    for (pos[1] = startPos[1]; pos[1] <= endPos[1]; ++pos[1]) {
      f(pos, breakLoops);
      if (breakLoops) {
        return;
      }
    }
  }
}

forEachPositionWithBreak() [2/2]

template<typename T , typename F >

constexpr void olb::particles::contact::forEachPositionWithBreak ( Vector< T, 3 > startPos, Vector< T, 3 > endPos, F f )

constexpr

Definition at line 114 of file particleContactForceFunctions.h.

{
  bool         breakLoops = false;
  Vector<T, 3> pos;
  for (pos[0] = startPos[0]; pos[0] <= endPos[0]; ++pos[0]) {
    for (pos[1] = startPos[1]; pos[1] <= endPos[1]; ++pos[1]) {
      for (pos[2] = startPos[2]; pos[2] <= endPos[2]; ++pos[2]) {
        f(pos, breakLoops);
        if (breakLoops) {
          return;
        }
      }
    }
  }
}

Here is the caller graph for this function:

getContactIterator()

template<typename CONTACTTYPE >

auto olb::particles::contact::getContactIterator	(	std::vector< CONTACTTYPE > &	contacts,
		const std::function< bool(const CONTACTTYPE &)>	condition )

Definition at line 38 of file particleContactDetectionFunctions.h.

{
  auto contactIt = std::find_if(contacts.begin(), contacts.end(),
                                [&condition](const auto& contact) -> bool {
                                  return condition(contact);
                                });
  return contactIt;
}

Here is the caller graph for this function:

getContactTreatmentResultsSerialSize()

template<typename T , unsigned D>

int olb::particles::contact::getContactTreatmentResultsSerialSize

(

)

Definition at line 857 of file particleContactCommunicationFunctions.h.

{
  constexpr unsigned Drot = utilities::dimensions::convert<D>::rotation;
 
  // Create communicatables
  std::array<std::size_t, 1> commGlobalParticleID {1};
  ConcreteCommunicatable     communicatableID =
      ConcreteCommunicatable(commGlobalParticleID);
  Vector<T, D>           forceDummy(1);
  ConcreteCommunicatable communicatableForce =
      ConcreteCommunicatable(forceDummy);
  Vector<T, Drot>        torqueDummy(1);
  ConcreteCommunicatable communicatableTorque =
      ConcreteCommunicatable(torqueDummy);
 
  // Setting dimension dependent indices
  const std::vector<unsigned> indicesID {0};
  std::vector<unsigned>       indicesForce;
  if constexpr (D == 3) {
    indicesForce = std::vector<unsigned> {0, 1, 2};
  }
  else {
    indicesForce = std::vector<unsigned> {0, 1};
  }
  std::vector<unsigned> indicesTorque;
  if constexpr (D == 3) {
    indicesTorque = indicesForce;
  }
  else {
    indicesTorque = std::vector<unsigned> {0};
  }
 
  // Retrieve serial size
  return communicatableID.size(indicesID) +
         communicatableForce.size(indicesForce) +
         communicatableTorque.size(indicesTorque);
}

References olb::ConcreteCommunicatable< COMMUNICATEE >::size().

Here is the call graph for this function:

getNormalForceFromOverlapVolume()

template<typename T , unsigned D>

Vector< T, D > olb::particles::contact::getNormalForceFromOverlapVolume	(	const Vector< T, D > &	contactNormal,
		const T	indentation,
		const T	overlapVolume,
		const T	effectiveE,
		const T	dampingFactor,
		const T	k,
		const Vector< T, D > &	relVel )

Definition at line 298 of file particleContactForceFunctions.h.

{
  const Vector<T, D> relVelNormal =
      evalRelativeNormalVelocity(contactNormal, relVel);
  return calcElasticAndViscousForce(effectiveE, k, overlapVolume, indentation,
                                    dampingFactor, relVelNormal, contactNormal);
}

References calcElasticAndViscousForce(), and evalRelativeNormalVelocity().

Here is the call graph for this function:

Here is the caller graph for this function:

getTangentialForceFromOverlapVolume()

template<typename T , unsigned D>

Vector< T, D > olb::particles::contact::getTangentialForceFromOverlapVolume	(	const Vector< T, D > &	contactNormal,
		const T	coefficientStaticFriction,
		const T	coefficientKineticFriction,
		const T	staticKineticTransitionVelocity,
		const T	k,
		const Vector< T, D > &	relVel,
		const Vector< T, D > &	normalForce )

Definition at line 260 of file particleContactForceFunctions.h.

{
  // Only regard if static friction is > 0, because we obtain nan otherwise
  if (coefficientStaticFriction > T {0.}) {
    const Vector<T, D> relVelNormal =
        evalRelativeNormalVelocity(contactNormal, relVel);
    const Vector<T, D> relVelTangential(relVel - relVelNormal);
    T const            magnitudeTangentialVelocity = norm(relVelTangential);
    if (magnitudeTangentialVelocity != T {0}) {
      T const velocityQuotient =
          magnitudeTangentialVelocity / staticKineticTransitionVelocity;
 
      T magnitudeTangentialForce =
          2 * coefficientStaticFriction *
          (1 - T {0.09} * util::pow(coefficientKineticFriction /
                                        coefficientStaticFriction,
                                    4));
      magnitudeTangentialForce -= coefficientKineticFriction;
      magnitudeTangentialForce *= velocityQuotient * velocityQuotient /
                                  (util::pow(velocityQuotient, 4) + 1);
      magnitudeTangentialForce += coefficientKineticFriction;
      magnitudeTangentialForce -= coefficientKineticFriction /
                                  (velocityQuotient * velocityQuotient + 1);
      magnitudeTangentialForce *= norm(normalForce);
 
      return -magnitudeTangentialForce * relVelTangential /
             magnitudeTangentialVelocity;
    }
  }
 
  return Vector<T, D>(T {0.});
}

References evalRelativeNormalVelocity(), olb::norm(), and olb::util::pow().

Here is the call graph for this function:

Here is the caller graph for this function:

isResponsibleRankSet()

template<typename CONTACTTYPE >

bool olb::particles::contact::isResponsibleRankSet

(

CONTACTTYPE &

contact

)

Definition at line 123 of file contactFunctions.h.

{
  return contact.getResponsibleRank() < std::numeric_limits<int>::max();
}

Here is the caller graph for this function:

particleContactConsistsOfIDs()

template<typename PARTICLECONTACTTYPE , bool IS_INPUT_SORTED = false>

bool olb::particles::contact::particleContactConsistsOfIDs	(	PARTICLECONTACTTYPE &	particleContact,
		const std::array< size_t, 2 > &	ids )

Definition at line 198 of file particleContactDetectionFunctions.h.

{
  if constexpr (!IS_INPUT_SORTED) {
    return (particleContact.getIDs()[0] == ids[0] &&
            particleContact.getIDs()[1] == ids[1]) ||
           (particleContact.getIDs()[0] == ids[1] &&
            particleContact.getIDs()[1] == ids[0]);
  }
  else {
    return particleContact.getIDs()[0] == ids[0] &&
           particleContact.getIDs()[1] == ids[1];
  }
  __builtin_unreachable();
}

Here is the caller graph for this function:

prepareForceCalculation()

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , bool CONVEX>

void olb::particles::contact::prepareForceCalculation

(

ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &

contactContainer

)

Definition at line 192 of file particleContactForceFunctions.h.

{
 
  communicateContacts(contactContainer);
 
  if constexpr (!CONVEX) {
    OstreamManager clout(std::cout, "forceCalcPrep");
    // TODO: Split contacts if the points are not connected
    // How to remember properties then? Store old min and max (maybe with a little extra) and then check for to which contact a position belongs?
    clout << "WARNING: Concave particles aren't supported." << std::endl;
  }
}

References communicateContacts().

Here is the call graph for this function:

processCell()

template<typename T , unsigned D, typename F1 , typename F2 , typename F3 , typename F4 , typename F5 >

void olb::particles::contact::processCell	(	const PhysR< T, D > &	pos,
		const PhysR< T, D > &	contactPhysDeltaX,
		PhysR< T, D > &	center,
		Vector< T, D > &	contactNormal,
		unsigned &	volumeCellCount,
		T &	surfaceCellCount,
		F1	isInside,
		F2	isOnSurface,
		F3	calcNormalA,
		F4	calcNormalB,
		F5	updateMinMax )

Definition at line 537 of file particleContactForceFunctions.h.

{
  // outside for normal calculation
  if (!isInside(pos)) {
    bool onSurfaceA = false, onSurfaceB = false;
    // calc normals
    Vector<T, D> normalA, normalB;
    // TODO: Improve mesh size
    T const meshSize = 1e-8 * util::average(contactPhysDeltaX);
    normalA          = calcNormalA(pos, meshSize);
    normalB          = calcNormalB(pos, meshSize);
    if (!util::nearZero(norm(normalA))) {
      normalA = normalize(normalA);
    }
    if (!util::nearZero(norm(normalB))) {
      normalB = normalize(normalB);
    }
 
    if (isOnSurface(pos, contactPhysDeltaX, onSurfaceA, onSurfaceB, normalA,
                    normalB)) {
      T approxSurface;
      if (!util::nearZero(norm(normalA)) && onSurfaceA) {
        approxSurface = evalApproxSurface(normalA, contactPhysDeltaX);
        surfaceCellCount += approxSurface;
        contactNormal += normalA * approxSurface;
      }
      if (!util::nearZero(norm(normalB)) && onSurfaceB) {
        approxSurface = evalApproxSurface(normalB, contactPhysDeltaX);
        surfaceCellCount += approxSurface;
        contactNormal -= normalB * approxSurface;
      }
    }
  }
  // inside volume calculation
  else {
    ++volumeCellCount;
 
    for (unsigned iD = 0; iD < D; ++iD) {
      center[iD] += pos[iD];
    }
    // Updating min and max for next sub-timestep
    updateMinMax(pos);
  }
}

References olb::util::average(), evalApproxSurface(), olb::util::nearZero(), olb::norm(), olb::normalize(), and updateMinMax().

Here is the call graph for this function:

Here is the caller graph for this function:

processContacts() [1/2]

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename CONTACTPROPERTIES , unsigned BBCORRECTIONMETHOD = 0, typename F1 = decltype(defaults::periodicity<PARTICLETYPE::d>), typename F2 = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>

void olb::particles::contact::processContacts	(	ParticleSystem< T, PARTICLETYPE > &	particleSystem,
		std::vector< SolidBoundary< T, PARTICLETYPE::d > > &	solidBoundaries,
		ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		const CONTACTPROPERTIES &	contactProperties,
		const SuperGeometry< T, PARTICLETYPE::d > &	sGeometry,
		const unsigned	contactBoxResolutionPerDirection = 8,
		const T	k = static_cast<T>(4. / (3 * util::sqrt(M_PI))),
		F1	getSetupPeriodicity = defaults::periodicity<PARTICLETYPE::d>,
		F2	processContactForce = defaults::processContactForce<T, PARTICLETYPE::d> )

Definition at line 1657 of file particleContactForceFunctions.h.

{
  std::multimap<int, std::unique_ptr<std::uint8_t[]>> dataMap;
 
  particle_particle<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE,
                    BBCORRECTIONMETHOD>::apply(particleSystem, contactContainer,
                                               contactProperties, sGeometry,
                                               dataMap,
                                               contactBoxResolutionPerDirection,
                                               k, getSetupPeriodicity,
                                               processContactForce);
 
  particle_wall<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE,
                BBCORRECTIONMETHOD>::apply(particleSystem, solidBoundaries,
                                           contactContainer, contactProperties,
                                           sGeometry, dataMap,
                                           contactBoxResolutionPerDirection, k,
                                           getSetupPeriodicity,
                                           processContactForce);
}

processContacts() [2/2]

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename CONTACTPROPERTIES , unsigned BBCORRECTIONMETHOD = 0, typename F1 = decltype(defaults::periodicity<PARTICLETYPE::d>), typename F2 = decltype(defaults::processContactForce<T, PARTICLETYPE::d>)>

void olb::particles::contact::processContacts	(	SuperParticleSystem< T, PARTICLETYPE > &	particleSystem,
		std::vector< SolidBoundary< T, PARTICLETYPE::d > > &	solidBoundaries,
		ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		const CONTACTPROPERTIES &	contactProperties,
		const SuperGeometry< T, PARTICLETYPE::d > &	sGeometry,
		MPI_Comm	contactTreatmentComm,
		const unsigned	contactBoxResolutionPerDirection = 8,
		const T	k = static_cast<T>(4. / (3 * util::sqrt(M_PI))),
		F1	getSetupPeriodicity = defaults::periodicity<PARTICLETYPE::d>,
		F2	processContactForce = defaults::processContactForce<T, PARTICLETYPE::d> )

Definition at line 1609 of file particleContactForceFunctions.h.

{
  std::multimap<int, std::unique_ptr<std::uint8_t[]>> dataMap;
 
  particle_particle<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE,
                    BBCORRECTIONMETHOD>::apply(particleSystem, contactContainer,
                                               contactProperties, sGeometry,
                                               dataMap,
                                               contactBoxResolutionPerDirection,
                                               k, getSetupPeriodicity,
                                               processContactForce);
 
  particle_wall<T, PARTICLETYPE, PARTICLECONTACTTYPE, WALLCONTACTTYPE,
                BBCORRECTIONMETHOD>::apply(particleSystem, solidBoundaries,
                                           contactContainer, contactProperties,
                                           sGeometry, dataMap,
                                           contactBoxResolutionPerDirection, k,
                                           getSetupPeriodicity,
                                           processContactForce);
  if constexpr (access::providesParallelization<PARTICLETYPE>()) {
    communicateContactTreatmentResults<T, PARTICLETYPE>(particleSystem, dataMap
#ifdef PARALLEL_MODE_MPI
                                                        ,
                                                        contactTreatmentComm
#endif
    );
  }
}

processContactViaVolume()

template<typename T , unsigned D, typename F1 , typename F2 , typename F3 , typename F4 >

void olb::particles::contact::processContactViaVolume	(	const Vector< T, D > &	min,
		const Vector< T, D > &	max,
		T	physDeltaX,
		unsigned	contactBoxResolutionPerDirection,
		F1	resetContactBox,
		F2	processCellWrapped,
		F3	calculateIndentation,
		F4	applyForce )

Definition at line 588 of file particleContactForceFunctions.h.

{
  Vector<T, D> contactPhysDeltaX =
      evalContactDeltaX(min, max, physDeltaX, contactBoxResolutionPerDirection);
 
  // direction of force
  Vector<T, D> contactNormal(T(0.));
  // contact point
  Vector<T, D> center(T(0.));
  // number of cells in contact volume
  unsigned volumeCellCount = 0;
  // number of cells in contact surface
  T surfaceCellCount = 0.;
  // volume of cell describing the contact
  T infinitesimalVolume = contactPhysDeltaX[0] * contactPhysDeltaX[1];
  if constexpr (D == 3) {
    infinitesimalVolume *= contactPhysDeltaX[2];
  }
 
  Vector<long, D> startPos;
  Vector<long, D> endPos;
  evalStartAndEndPos(startPos, endPos, max, min, contactPhysDeltaX);
  // Increase the box 2*deltaX in each direction
  startPos -= Vector<long, D>(1);
  endPos += Vector<long, D>(1);
 
  resetContactBox();
 
  forEachPosition(startPos, endPos, [&](const Vector<long, D>& pos) {
    Vector<T, D> physPos;
    for (unsigned iD = 0; iD < D; ++iD) {
      physPos[iD] = pos[iD] * contactPhysDeltaX[iD];
    }
    processCellWrapped(contactNormal, center, volumeCellCount, surfaceCellCount,
                       physPos, contactPhysDeltaX);
  });
 
  if (!util::nearZero(norm(contactNormal))) {
    if (volumeCellCount > 0 && surfaceCellCount > 0) {
      // finishing center calculation
      center /= volumeCellCount;
      // finishing normal calculation
      contactNormal = normalize(contactNormal);
      // calculate precision depending on the size of the overlap
      T precision {0};
      for (unsigned iD = 0; iD < D; ++iD) {
        precision += util::pow(contactNormal[iD] * contactPhysDeltaX[iD], 2);
      }
      precision = T {1e-2} * util::sqrt(precision / D);
      const T indentation =
          calculateIndentation(center, contactNormal, precision);
      // currently, the normal points outside the particle (because we use SDF), but the force acts in the opposite way
      contactNormal *= -1;
      // Avoid indentation < 0
      if (indentation > 0) {
        applyForce(center, contactNormal, util::max(indentation, T {0}),
                   volumeCellCount * infinitesimalVolume);
      }
    }
  }
#ifdef OLB_DEBUG
  else {
    OstreamManager clout(std::cout, "contactProcessing");
    clout << "WARNING: The contact normal's norm is zero. Therefore, the "
             "detected contact is ignored."
          << std::endl;
  }
#endif
}

References evalContactDeltaX(), evalStartAndEndPos(), forEachPosition(), olb::util::max(), olb::util::nearZero(), olb::norm(), olb::normalize(), olb::util::pow(), and olb::util::sqrt().

Here is the call graph for this function:

Here is the caller graph for this function:

receiveContact() [1/2]

template<typename T , typename PARTICLETYPE , typename CONTACTTYPE , bool CONVEX, bool SANITYCHECK = false>

void olb::particles::contact::receiveContact	(	ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	newcontact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		std::vector< CONTACTTYPE > &	contacts,
		int	rankOrig )

Definition at line 546 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
 
  bool addNewContact = true;
 
  // Optional sanity check
  if constexpr (SANITYCHECK) {
    if (newcontact.getResponsibleRank() == singleton::mpi().getRank()) {
      bool isParticleDataAvailable[2] = {false, false};
      // Iterate over all particles to check for errors
      communication::forParticlesInSuperParticleSystem(
          sParticleSystem,
          [&](Particle<T, PARTICLETYPE>&       particle,
              ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
            for (unsigned short i = 0; i < 2; ++i) {
              if (particle.template getField<PARALLELIZATION, ID>() ==
                  newcontact.getIDs()[i]) {
                isParticleDataAvailable[i] = true;
              }
            }
          });
      if (!isParticleDataAvailable[0] || !isParticleDataAvailable[1]) {
        addNewContact = false;
        std::cout << "WARNING: rank " << singleton::mpi().getRank()
                  << " received a contact of particles with the global ids "
                  << newcontact.getIDs()[0] << " and " << newcontact.getIDs()[1]
                  << " from rank " << rankOrig
                  << " for which no particle data are available." << std::endl;
      }
    }
  }
 
  if (addNewContact) {
    const std::array<std::size_t, 2> particleIDs = newcontact.getIDs();
 
    const std::function<bool(const CONTACTTYPE&)> condition =
        [&particleIDs](const CONTACTTYPE& contact) {
          return particleContactConsistsOfIDs(contact, particleIDs);
        };
 
    auto contactIt = getContactIterator(contacts, condition);
 
    if (contactIt != std::end(contacts)) {
      contactIt->combineWith(newcontact);
    }
    else {
      contacts.push_back(newcontact);
    }
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
    std::cout << "Rank " << singleton::mpi().getRank()
              << " received contact of ids " << newcontact.getIDs()[0] << " & "
              << newcontact.getIDs()[1] << ". Rank "
              << newcontact.getResponsibleRank()
              << " is responsible for the calculation." << std::endl;
#endif
  }
}

References olb::particles::communication::forParticlesInSuperParticleSystem(), getContactIterator(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getIDs(), olb::singleton::MpiManager::getRank(), olb::particles::contact::ParticleContactArbitraryFromOverlapVolume< T, D, CONVEX >::getResponsibleRank(), olb::singleton::mpi(), and particleContactConsistsOfIDs().

Here is the call graph for this function:

Here is the caller graph for this function:

receiveContact() [2/2]

template<typename T , typename PARTICLETYPE , typename CONTACTTYPE , bool CONVEX, bool SANITYCHECK = false>

void olb::particles::contact::receiveContact	(	WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	newcontact,
		SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		std::vector< CONTACTTYPE > &	contacts,
		int	rankOrig )

Definition at line 612 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
 
  bool addNewContact = true;
 
  if constexpr (SANITYCHECK) {
    if (newcontact.getResponsibleRank() == singleton::mpi().getRank()) {
      bool isParticleDataAvailable = false;
      // Iterate over all particles to check for errors
      communication::forParticlesInSuperParticleSystem(
          sParticleSystem,
          [&](Particle<T, PARTICLETYPE>&       particle,
              ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
            if (particle.template getField<PARALLELIZATION, ID>() ==
                newcontact.getParticleID()) {
              isParticleDataAvailable = true;
            }
          });
      if (!isParticleDataAvailable) {
        addNewContact = false;
        std::cout << "WARNING: rank " << singleton::mpi().getRank()
                  << " received a contact of particle with the global id "
                  << newcontact.getParticleID() << " from rank " << rankOrig
                  << " for which no particle data is available." << std::endl;
      }
    }
  }
 
  if (addNewContact) {
    const std::size_t particleID            = newcontact.getParticleID();
    const std::size_t solidBoundaryMaterial = newcontact.getWallID();
 
    const std::function<bool(const CONTACTTYPE&)> condition =
        [&particleID, &solidBoundaryMaterial](const CONTACTTYPE& contact) {
          return particleID == contact.getParticleID() &&
                 solidBoundaryMaterial == contact.getWallID();
        };
 
    auto contactIt = getContactIterator(contacts, condition);
 
    if (contactIt != std::end(contacts)) {
      contactIt->combineWith(newcontact);
    }
    else {
      contacts.push_back(newcontact);
    }
  }
}

References olb::particles::communication::forParticlesInSuperParticleSystem(), getContactIterator(), olb::singleton::MpiManager::getRank(), and olb::singleton::mpi().

Here is the call graph for this function:

receiveContactTreatmentResults()

template<typename T , typename PARTICLETYPE >

void olb::particles::contact::receiveContactTreatmentResults	(	SuperParticleSystem< T, PARTICLETYPE > &	sParticleSystem,
		MPI_Comm	contactTreatmentComm,
		singleton::MpiNonBlockingHelper &	mpiNbHelper )

Definition at line 952 of file particleContactCommunicationFunctions.h.

{
  using namespace descriptors;
 
  // Retrieve dimensions
  constexpr unsigned D    = PARTICLETYPE::d;
  constexpr unsigned Drot = utilities::dimensions::convert<D>::rotation;
 
  // Retrieve rank of direct neighbours
  auto& listNeighbourRanks = sParticleSystem.getNeighbourRanks();
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
  std::cout << "Rank " << singleton::mpi().getRank() << " has the neighbours: ";
  for (auto& neighbourRank : listNeighbourRanks) {
    std::cout << neighbourRank << " ";
  }
  std::cout << std::endl;
#endif
 
  // Setting dimension dependent indices
  const std::vector<unsigned> indicesID {0};
  std::vector<unsigned>       indicesForce;
  if constexpr (D == 3) {
    indicesForce = std::vector<unsigned> {0, 1, 2};
  }
  else {
    indicesForce = std::vector<unsigned> {0, 1};
  }
  std::vector<unsigned> indicesTorque;
  if constexpr (D == 3) {
    indicesTorque = indicesForce;
  }
  else {
    indicesTorque = std::vector<unsigned> {0};
  }
 
  std::array<std::size_t, 1> commGlobalParticleID;
  Vector<T, D>               contactForce;
  Vector<T, Drot>            contactTorque;
 
  // Create communicatables
  ConcreteCommunicatable communicatableID =
      ConcreteCommunicatable(commGlobalParticleID);
  ConcreteCommunicatable communicatableForce =
      ConcreteCommunicatable(contactForce);
  ConcreteCommunicatable communicatableTorque =
      ConcreteCommunicatable(contactTorque);
 
  // Retrieve serial size
  const std::size_t serialSize = communicatableID.size(indicesID) +
                                 communicatableForce.size(indicesForce) +
                                 communicatableTorque.size(indicesTorque);
 
  // Receive data and iterate buffer
  // TODO: Improve processing of data, do not iterate through the all particles every time (maybe also add main iC to the communicated data to reduce processed particles)
  communication::receiveAndExecuteForData(
      listNeighbourRanks, serialSize, contactTreatmentComm, mpiNbHelper,
      [&](int rankOrig, std::uint8_t* buffer) {
        std::size_t serialIdx = communicatableID.deserialize(indicesID, buffer);
        serialIdx +=
            communicatableForce.deserialize(indicesForce, &buffer[serialIdx]);
        serialIdx +=
            communicatableTorque.deserialize(indicesTorque, &buffer[serialIdx]);
 
#ifdef VERBOSE_CONTACT_COMMUNICATION
        std::cout << "Rank " << singleton::mpi().getRank()
                  << " received contact treatment results for particle "
                  << commGlobalParticleID[0] << " (force = " << contactForce
                  << "; torque = " << contactTorque << ") from rank "
                  << rankOrig << std::endl;
#endif
 
        communication::forParticlesInSuperParticleSystem<
            T, PARTICLETYPE, conditions::valid_particles>(
            sParticleSystem,
            [&](Particle<T, PARTICLETYPE>&       particle,
                ParticleSystem<T, PARTICLETYPE>& particleSystem, int globiC) {
              const std::size_t currentGlobalParticleID =
                  particle.template getField<PARALLELIZATION, ID>();
              if (commGlobalParticleID[0] == currentGlobalParticleID) {
                // Process contact force
                const Vector<T, D> force =
                    particle.template getField<FORCING, FORCE>();
                particle.template setField<FORCING, FORCE>(force +
                                                           contactForce);
                // Process torque from contact
                const Vector<T, Drot> torque(
                    particle.template getField<FORCING, TORQUE>());
                particle.template setField<FORCING, TORQUE>(
                    utilities::dimensions::convert<D>::serialize_rotation(
                        torque + contactTorque));
              }
            });
      });
}

References olb::ConcreteCommunicatable< COMMUNICATEE >::deserialize(), olb::particles::communication::forParticlesInSuperParticleSystem(), olb::particles::SuperParticleSystem< T, PARTICLETYPE >::getNeighbourRanks(), olb::singleton::MpiManager::getRank(), olb::singleton::mpi(), olb::particles::communication::receiveAndExecuteForData(), and olb::ConcreteCommunicatable< COMMUNICATEE >::size().

Here is the call graph for this function:

Here is the caller graph for this function:

resetResponsibleRank()

template<typename CONTACTTYPE >

void olb::particles::contact::resetResponsibleRank

(

CONTACTTYPE &

contact

)

Definition at line 117 of file contactFunctions.h.

{
  contact.setResponsibleRank(std::numeric_limits<int>::max());
}

Here is the caller graph for this function:

sortParticleIDs()

std::array< std::size_t, 2 > olb::particles::contact::sortParticleIDs

(

const std::array< std::size_t, 2 > &

ids

)

Definition at line 49 of file particleContactDetectionFunctions.h.

{
  return std::array<std::size_t, 2>(
      {util::min(ids[0], ids[1]), util::max(ids[0], ids[1])});
}

References olb::util::max(), and olb::util::min().

Here is the call graph for this function:

unifyPosition()

template<typename T , typename PARTICLETYPE , typename F >

PhysR< T, PARTICLETYPE::d > olb::particles::contact::unifyPosition	(	Particle< T, PARTICLETYPE > &	particle,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		T	deltaX )

Definition at line 124 of file particleContactDetectionFunctions.h.

{
  using namespace descriptors;
  const PhysR<T, PARTICLETYPE::d> pos =
      particles::access::getPosition(particle);
  PhysR<T, PARTICLETYPE::d> unifiedPosition(pos);
 
  if constexpr (isPeriodic(getSetupPeriodicity())) {
    constexpr Vector<bool, PARTICLETYPE::d> isPeriodic(getSetupPeriodicity());
 
    auto    sIndicator   = particle.template getField<SURFACE, SINDICATOR>();
    const T circumRadius = sIndicator->getCircumRadius() -
                           sIndicator->getEpsilon() +
                           util::max(sIndicator->getEpsilon(), 0.5 * deltaX);
 
    for (unsigned i = 0; i < PARTICLETYPE::d; ++i) {
      T const particleMax = pos[i] + circumRadius;
      //T const particleMin = pos[i] - circumRadius;
 
      if (particleMax > cellMax[i] && isPeriodic[i]) {
        unifiedPosition[i] = communication::movePositionToStart(
            pos[i], cellMax[i], cellMin[i]);
      }
      /*
      else if (particleMin < cellMin[i] && isPeriodic[i]) {
        unifiedPosition[i] = cellMax[i] - (cellMin[i] - pos[i]);
      }
      */
    }
  }
  return unifiedPosition;
}

References olb::particles::access::getPosition(), olb::particles::isPeriodic(), olb::util::max(), and olb::particles::communication::movePositionToStart().

Here is the call graph for this function:

Here is the caller graph for this function:

unifyPositions()

template<typename T , typename PARTICLETYPE , typename F >

void olb::particles::contact::unifyPositions	(	Particle< T, PARTICLETYPE > &	particle1,
		Particle< T, PARTICLETYPE > &	particle2,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		std::array< PhysR< T, PARTICLETYPE::d >, 2 > &	positions,
		T	deltaX )

Definition at line 56 of file particleContactDetectionFunctions.h.

{
  using namespace descriptors;
 
  const PhysR<T, PARTICLETYPE::d> pos1 =
      particles::access::getPosition(particle1);
  const PhysR<T, PARTICLETYPE::d> pos2 =
      particles::access::getPosition(particle2);
  positions[0] = pos1;
  positions[1] = pos2;
 
  if constexpr (isPeriodic(getSetupPeriodicity())) {
    constexpr Vector<bool, PARTICLETYPE::d> isPeriodic(getSetupPeriodicity());
 
    auto    sIndicator1   = particle1.template getField<SURFACE, SINDICATOR>();
    const T circumRadius1 = sIndicator1->getCircumRadius() -
                            sIndicator1->getEpsilon() +
                            util::max(sIndicator1->getEpsilon(), 0.5 * deltaX);
    auto    sIndicator2   = particle2.template getField<SURFACE, SINDICATOR>();
    const T circumRadius2 = sIndicator2->getCircumRadius() -
                            sIndicator2->getEpsilon() +
                            util::max(sIndicator2->getEpsilon(), 0.5 * deltaX);
 
    const PhysR<T, PARTICLETYPE::d> middle = 0.5 * (cellMin + cellMax);
 
    for (unsigned i = 0; i < PARTICLETYPE::d; ++i) {
      T const particleMin1 = pos1[i] - circumRadius1;
      T const particleMax1 = pos1[i] + circumRadius1;
      T const particleMin2 = pos2[i] - circumRadius2;
      T const particleMax2 = pos2[i] + circumRadius2;
 
      // Always solely "move" the second particle
      if (particleMin1 < cellMin[i] && isPeriodic[i]) {
        if (pos2[i] > middle[i]) {
          positions[1][i] = communication::movePositionToStart(
              pos2[i], cellMax[i], cellMin[i]);
        }
      }
      else if (particleMax1 > cellMax[i] && isPeriodic[i]) {
        if (pos2[i] < middle[i]) {
          positions[1][i] = communication::movePositionToEnd(
              pos2[i], cellMax[i], cellMin[i]);
        }
      }
      else if (particleMin2 < cellMin[i] && isPeriodic[i]) {
        if (pos1[i] > middle[i]) {
          positions[1][i] = communication::movePositionToEnd(
              pos2[i], cellMax[i], cellMin[i]);
        }
      }
      else if (particleMax2 > cellMax[i] && isPeriodic[i]) {
        if (pos1[i] < middle[i]) {
          positions[1][i] = communication::movePositionToStart(
              pos2[i], cellMax[i], cellMin[i]);
        }
      }
    }
  }
}

References olb::particles::access::getPosition(), olb::particles::isPeriodic(), olb::util::max(), olb::particles::communication::movePositionToEnd(), and olb::particles::communication::movePositionToStart().

Here is the call graph for this function:

Here is the caller graph for this function:

updateContact() [1/2]

template<typename T , typename PARTICLETYPE , bool CONVEX, typename F >

void olb::particles::contact::updateContact	(	ParticleContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		Particle< T, PARTICLETYPE > &	particle1,
		Particle< T, PARTICLETYPE > &	particle2,
		const PhysR< T, PARTICLETYPE::d > &	contactPos,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		T	deltaX )

Definition at line 215 of file particleContactDetectionFunctions.h.

{
  if (!contact.isParticlePositionUpdated()) {
    std::array<PhysR<T, PARTICLETYPE::d>, 2> positions;
    unifyPositions(particle1, particle2, cellMin, cellMax, getSetupPeriodicity,
                   positions, deltaX);
    contact.setParticlePositions(positions);
  }
  // It doesn't matter which particle is used, since the contact must be in both
  contact.updateMinMax(evalContactPosition(
      particle1, contact.getParticlePosition(contact.getIDs()[0]),
      contactPos, cellMin, cellMax, getSetupPeriodicity, deltaX));
}

References evalContactPosition(), and unifyPositions().

Here is the call graph for this function:

Here is the caller graph for this function:

updateContact() [2/2]

template<typename T , typename PARTICLETYPE , bool CONVEX, typename F >

void olb::particles::contact::updateContact	(	WallContactArbitraryFromOverlapVolume< T, PARTICLETYPE::d, CONVEX > &	contact,
		Particle< T, PARTICLETYPE > &	particle,
		const PhysR< T, PARTICLETYPE::d > &	contactPos,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		T	deltaX )

Definition at line 237 of file particleContactDetectionFunctions.h.

{
  if (!contact.isParticlePositionUpdated()) {
    contact.setParticlePosition(
        unifyPosition(particle, cellMin, cellMax, getSetupPeriodicity, deltaX));
  }
  contact.updateMinMax(
      evalContactPosition(particle, contact.getParticlePosition(), contactPos,
                          cellMin, cellMax, getSetupPeriodicity, deltaX));
}

References evalContactPosition(), and unifyPosition().

Here is the call graph for this function:

updateContacts() [1/2]

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >

void olb::particles::contact::updateContacts	(	ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE > &	contactContainer,
		size_t	particleID,
		unsigned	wallID,
		const PhysR< T, PARTICLETYPE::d > &	pos,
		Particle< T, PARTICLETYPE > &	particle,
		const PhysR< T, PARTICLETYPE::d > &	cellMin,
		const PhysR< T, PARTICLETYPE::d > &	cellMax,
		F	getSetupPeriodicity,
		T	deltaX )

Definition at line 304 of file particleContactDetectionFunctions.h.

{
  // find the contact by checking ids
  const std::function<bool(const WALLCONTACTTYPE&)> condition =
      [&particleID, &wallID](const WALLCONTACTTYPE& contact) {
        return particleID == contact.getParticleID() &&
               wallID == contact.getWallID();
      };
  auto contactIt =
      getContactIterator(contactContainer.wallContacts, condition);
 
  if (contactIt != std::end(contactContainer.wallContacts)) {
    if (contactIt->isNew()) {
      updateContact(*contactIt, particle, pos, cellMin, cellMax,
                    getSetupPeriodicity, deltaX);
    }
  }
  else {
    contactContainer.wallContacts.push_back(
        WALLCONTACTTYPE(particleID, wallID));
    updateContact(*(contactContainer.wallContacts.end() - 1), particle, pos,
                  cellMin, cellMax, getSetupPeriodicity, deltaX);
  }
}

References getContactIterator(), updateContact(), and olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::wallContacts.

Here is the call graph for this function:

updateContacts() [2/2]

template<typename T , typename PARTICLETYPE , typename PARTICLECONTACTTYPE , typename WALLCONTACTTYPE , typename F >

void olb::particles::contact::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 )

Definition at line 256 of file particleContactDetectionFunctions.h.

{
  // find the contact by checking particle ids
  const std::function<bool(const PARTICLECONTACTTYPE&)> condition =
      [&ids](const PARTICLECONTACTTYPE& contact) {
        // we assume that ids aren't sorted
        return particleContactConsistsOfIDs(contact, ids);
      };
  auto contactIt =
      getContactIterator(contactContainer.particleContacts, condition);
 
  if (contactIt != std::end(contactContainer.particleContacts)) {
    if (contactIt->isNew()) {
      if (particleContactConsistsOfIDs<PARTICLECONTACTTYPE, true>(
              *contactIt, ids)) {
        updateContact(*contactIt, particle1, particle2, pos, cellMin, cellMax,
                      getSetupPeriodicity, deltaX);
      }
      else {
        updateContact(*contactIt, particle2, particle1, pos, cellMin, cellMax,
                      getSetupPeriodicity, deltaX);
      }
    }
  }
  else {
    contactContainer.particleContacts.push_back(PARTICLECONTACTTYPE(ids));
    if (particleContactConsistsOfIDs<PARTICLECONTACTTYPE, true>(
            *(contactContainer.particleContacts.end() - 1), ids)) {
      updateContact(*(contactContainer.particleContacts.end() - 1), particle1,
                    particle2, pos, cellMin, cellMax, getSetupPeriodicity,
                    deltaX);
    }
    else {
      updateContact(*(contactContainer.particleContacts.end() - 1), particle2,
                    particle1, pos, cellMin, cellMax, getSetupPeriodicity,
                    deltaX);
    }
  }
}

References getContactIterator(), particleContactConsistsOfIDs(), olb::particles::contact::ContactContainer< T, PARTICLECONTACTTYPE, WALLCONTACTTYPE >::particleContacts, and updateContact().

Here is the call graph for this function:

Here is the caller graph for this function:

updateMinMax()

template<typename T , unsigned D>

void olb::particles::contact::updateMinMax	(	PhysR< T, D > &	min,
		PhysR< T, D > &	max,
		const PhysR< T, D > &	pos )

Definition at line 34 of file contactHelpers.h.

{
  for (unsigned iD = 0; iD < D; ++iD) {
    min[iD] = util::min(min[iD], pos[iD]);
    max[iD] = util::max(max[iD], pos[iD]);
  }
}

References olb::util::max(), and olb::util::min().

Here is the call graph for this function:

Here is the caller graph for this function: