diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc index 7b25049dc2..852608c4e3 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.cc @@ -20,12 +20,16 @@ namespace Spheral { template CircularPlaneSolidBoundary:: -CircularPlaneSolidBoundary(const Vector& point, const Vector& normal, const Scalar& extent): +CircularPlaneSolidBoundary(const Vector& point, + const Vector& normal, + const Scalar extent, + const RotationType& angularVelocity): SolidBoundaryBase(), mPoint(point), mNormal(normal), mExtent(extent), - mVelocity(Vector::zero()){ + mVelocity(Vector::zero()), + mAngularVelocity(angularVelocity){ } template @@ -51,7 +55,15 @@ template typename Dimension::Vector CircularPlaneSolidBoundary:: localVelocity(const Vector& position) const { - return mVelocity; + // Linear velocity component + Vector linearVelocityComponent = mVelocity; + + // Angular velocity component: v = w × r + Vector r = position - mPoint; // Position vector from the point of rotation + Vector angularVelocityComponent = DEMDimension::cross(mAngularVelocity,r); + + // Total velocity is the sum of linear and angular components + return linearVelocityComponent + angularVelocityComponent; } template @@ -63,9 +75,11 @@ registerState(DataBase& dataBase, const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; const auto normalKey = boundaryKey +"_normal"; + const auto angularVelocityKey = boundaryKey + "_angularVelocity"; // New key for angular velocity state.enroll(pointKey,mPoint); - state.enroll(pointKey,mVelocity); - state.enroll(pointKey,mNormal); + state.enroll(velocityKey,mVelocity); + state.enroll(normalKey,mNormal); + state.enroll(angularVelocityKey, mAngularVelocity); // Enroll angular velocity } template @@ -73,6 +87,10 @@ void CircularPlaneSolidBoundary:: update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; + // Update the orientation of the boundary based on angular velocity + // Assuming a simple Euler for rotation + mNormal += multiplier * DEMDimension::cross(mAngularVelocity,mNormal); + mNormal = mNormal.unitVector(); // Normalize to maintain unit length } @@ -87,6 +105,7 @@ dumpState(FileIO& file, const string& pathName) const { file.write(mNormal, pathName + "/normal"); file.write(mExtent, pathName + "/extent"); file.write(mVelocity, pathName + "/velocity"); + file.write(mAngularVelocity, pathName + "/omega"); // Write angular velocity } @@ -98,6 +117,7 @@ restoreState(const FileIO& file, const string& pathName) { file.read(mNormal, pathName + "/normal"); file.read(mExtent, pathName + "/extent"); file.read(mVelocity, pathName + "/velocity"); + file.read(mAngularVelocity, pathName + "/omega"); // Read angular velocity } -} \ No newline at end of file +} diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh index 9b9a674f4d..2f469ba312 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundary.hh @@ -8,6 +8,7 @@ #ifndef __Spheral_CircularPlaneSolidBoundary_hh__ #define __Spheral_CircularPlaneSolidBoundary_hh__ +#include "DEM/DEMDimension.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" namespace Spheral { @@ -21,13 +22,15 @@ class CircularPlaneSolidBoundary : public SolidBoundaryBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; + typedef typename DEMDimension::AngularVector RotationType; public: //--------------------------- Public Interface ---------------------------// CircularPlaneSolidBoundary(const Vector& point, - const Vector& normal, - const Scalar& exent); + const Vector& normal, + const Scalar exent, + const RotationType& angularVelocity); ~CircularPlaneSolidBoundary(); @@ -44,7 +47,7 @@ public: const Vector& point() const; void point(const Vector& value); - const Vector& normal() const; + const Vector& normal() const; void normal(const Vector& value); Scalar extent() const; @@ -53,6 +56,9 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + const RotationType& angularVelocity() const; + void angularVelocity(const RotationType& value); + virtual std::string label() const override { return "CircularPlaneSolidBoundary" ; } virtual void dumpState(FileIO& file, const std::string& pathName) const override; virtual void restoreState(const FileIO& file, const std::string& pathName) override; @@ -63,6 +69,7 @@ protected: Vector mNormal; Scalar mExtent; Vector mVelocity; + RotationType mAngularVelocity; private: //--------------------------- Private Interface ---------------------------// diff --git a/src/DEM/SolidBoundary/CircularPlaneSolidBoundaryInline.hh b/src/DEM/SolidBoundary/CircularPlaneSolidBoundaryInline.hh index e961da3244..ca4ea2d350 100644 --- a/src/DEM/SolidBoundary/CircularPlaneSolidBoundaryInline.hh +++ b/src/DEM/SolidBoundary/CircularPlaneSolidBoundaryInline.hh @@ -33,6 +33,7 @@ normal(const typename Dimension::Vector& value) { mNormal = value; } + template inline typename Dimension::Scalar @@ -49,6 +50,7 @@ extent(typename Dimension::Scalar value) { mExtent=value; } + template inline const typename Dimension::Vector& @@ -65,4 +67,21 @@ velocity(const typename Dimension::Vector& value) { mVelocity=value; } + +template +inline +const typename DEMDimension::AngularVector& +CircularPlaneSolidBoundary:: +angularVelocity() const { + return mAngularVelocity; +} + +template +inline +void +CircularPlaneSolidBoundary:: +angularVelocity(const typename DEMDimension::AngularVector& value) { + mAngularVelocity=value; +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc index 81116c88cc..9252de5e63 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.cc +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.cc @@ -20,15 +20,17 @@ namespace Spheral { template CylinderSolidBoundary:: CylinderSolidBoundary(const Vector& point, - const Vector& axis, - const Scalar radius, - const Scalar length): + const Vector& axis, + const Scalar radius, + const Scalar length, + const RotationType& angularVelocity): SolidBoundaryBase(), mPoint(point), mAxis(axis), mRadius(radius), mLength(length), - mVelocity(Vector::zero()){ + mVelocity(Vector::zero()), + mAngularVelocity(angularVelocity){ } template @@ -52,7 +54,15 @@ template typename Dimension::Vector CylinderSolidBoundary:: localVelocity(const Vector& position) const { - return mVelocity; + // Calculate the vector from the axis of rotation to the position + const auto p = position - mPoint; + const auto pnMag = p.dot(mAxis); + const auto pn = pnMag * mAxis; + const auto r = p - pn; // Radial vector from the axis to the position + + // Calculate the tangential velocity due to angular velocity + const auto tangentialVelocity = DEMDimension::cross(r,mAngularVelocity); + return mVelocity + tangentialVelocity; } template @@ -62,11 +72,17 @@ registerState(DataBase& dataBase, State& state) { const auto boundaryKey = "CylinderSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; + const auto axisKey = boundaryKey +"_axis"; + const auto radiusKey = boundaryKey +"_radius"; + const auto lengthKey = boundaryKey +"_length"; const auto velocityKey = boundaryKey +"_velocity"; - //const auto normalKey = boundaryKey +"_normal"; + const auto angularVelocityKey = boundaryKey + "_angularVelocity"; // Key for angular velocity state.enroll(pointKey,mPoint); - state.enroll(pointKey,mVelocity); - //state.enroll(pointKey,mNormal); + state.enroll(axisKey,mAxis); + state.enroll(radiusKey,mRadius); + state.enroll(lengthKey,mLength); + state.enroll(velocityKey,mVelocity); + state.enroll(angularVelocityKey, mAngularVelocity); // Enroll angular velocity } template @@ -74,6 +90,7 @@ void CylinderSolidBoundary:: update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; + // If we want more complex rotation we'll need more complex logic here. } //------------------------------------------------------------------------------ @@ -88,6 +105,7 @@ dumpState(FileIO& file, const string& pathName) const { file.write(mRadius, pathName + "/radius"); file.write(mLength, pathName + "/length"); file.write(mVelocity, pathName + "/velocity"); + file.write(mAngularVelocity, pathName + "/omega"); // Write angular velocity } @@ -100,6 +118,7 @@ restoreState(const FileIO& file, const string& pathName) { file.read(mRadius, pathName + "/radius"); file.read(mLength, pathName + "/length"); file.read(mVelocity, pathName + "/velocity"); + file.read(mAngularVelocity, pathName + "/omega"); // Read angular velocity } diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh index 14ab303ffb..e00ff14421 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundary.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundary.hh @@ -7,6 +7,7 @@ #ifndef __Spheral_CylinderSolidBoundary_hh__ #define __Spheral_CylinderSolidBoundary_hh__ +#include "DEM/DEMDimension.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" namespace Spheral { @@ -21,14 +22,16 @@ class CylinderSolidBoundary : public SolidBoundaryBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; typedef typename Dimension::Tensor Tensor; + typedef typename DEMDimension::AngularVector RotationType; public: //--------------------------- Public Interface ---------------------------// CylinderSolidBoundary(const Vector& point, - const Vector& axis, - const Scalar radius, - const Scalar length); + const Vector& axis, + const Scalar radius, + const Scalar length, + const RotationType& angularVelocity); ~CylinderSolidBoundary(); @@ -57,6 +60,9 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + const RotationType& angularVelocity() const; + void angularVelocity(const RotationType& value); + virtual std::string label() const override { return "CylinderSolidBoundary" ; } virtual void dumpState(FileIO& file, const std::string& pathName) const override; virtual void restoreState(const FileIO& file, const std::string& pathName) override; @@ -69,6 +75,7 @@ protected: Scalar mLength; Vector mVelocity; + RotationType mAngularVelocity; private: //--------------------------- Private Interface ---------------------------// diff --git a/src/DEM/SolidBoundary/CylinderSolidBoundaryInline.hh b/src/DEM/SolidBoundary/CylinderSolidBoundaryInline.hh index bd7a305b25..cb4772a457 100644 --- a/src/DEM/SolidBoundary/CylinderSolidBoundaryInline.hh +++ b/src/DEM/SolidBoundary/CylinderSolidBoundaryInline.hh @@ -16,6 +16,7 @@ point(const typename Dimension::Vector& value) { mPoint=value; } + template inline const typename Dimension::Vector& @@ -32,6 +33,7 @@ axis(const typename Dimension::Vector& value) { mAxis=value; } + template inline typename Dimension::Scalar @@ -48,6 +50,7 @@ length(typename Dimension::Scalar value) { mLength=value; } + template inline typename Dimension::Scalar @@ -64,6 +67,7 @@ radius(typename Dimension::Scalar value) { mRadius=value; } + template inline const typename Dimension::Vector& @@ -80,4 +84,20 @@ velocity(const typename Dimension::Vector& value) { mVelocity=value; } + +template +inline +const typename DEMDimension::AngularVector& +CylinderSolidBoundary:: +angularVelocity() const { + return mAngularVelocity; +} + +template +inline +void +CylinderSolidBoundary:: +angularVelocity(const typename DEMDimension::AngularVector& value) { + mAngularVelocity=value; +} } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc index 15f3ef5174..33e862dae8 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.cc @@ -19,11 +19,14 @@ namespace Spheral { template InfinitePlaneSolidBoundary:: -InfinitePlaneSolidBoundary(const Vector& point, const Vector& normal): +InfinitePlaneSolidBoundary(const Vector& point, + const Vector& normal, + const RotationType& angularVelocity): SolidBoundaryBase(), mPoint(point), mNormal(normal), - mVelocity(Vector::zero()){ + mVelocity(Vector::zero()), + mAngularVelocity(angularVelocity){ } template @@ -42,7 +45,7 @@ template typename Dimension::Vector InfinitePlaneSolidBoundary:: localVelocity(const Vector& position) const { - return mVelocity; + return mVelocity + DEMDimension::cross((position-mPoint),mAngularVelocity); // Include angular velocity effect } template @@ -54,9 +57,11 @@ registerState(DataBase& dataBase, const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; const auto normalKey = boundaryKey +"_normal"; + const auto angularVelocityKey = boundaryKey +"_angularVelocity"; state.enroll(pointKey,mPoint); state.enroll(velocityKey,mVelocity); state.enroll(normalKey,mNormal); + state.enroll(angularVelocityKey, mAngularVelocity); // Enroll angular velocity } template @@ -64,6 +69,10 @@ void InfinitePlaneSolidBoundary:: update(const double multiplier, const double t, const double dt) { mPoint += multiplier*mVelocity; + // Update the normal vector based on angular velocity + // Using a small angle approximation: newNormal = oldNormal + dt * (angularVelocity x oldNormal) + mNormal += multiplier * DEMDimension::cross(mNormal,mAngularVelocity); + mNormal = mNormal.unitVector(); // Normalize to maintain it as a unit vector } //------------------------------------------------------------------------------ @@ -76,6 +85,7 @@ dumpState(FileIO& file, const string& pathName) const { file.write(mPoint, pathName + "/point"); file.write(mNormal, pathName + "/normal"); file.write(mVelocity, pathName + "/velocity"); + file.write(mAngularVelocity, pathName + "/omega"); // Write angular velocity } @@ -86,6 +96,7 @@ restoreState(const FileIO& file, const string& pathName) { file.read(mPoint, pathName + "/point"); file.read(mNormal, pathName + "/normal"); file.read(mVelocity, pathName + "/velocity"); + file.read(mAngularVelocity, pathName + "/omega"); // Read angular velocity } } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh index 406a0c75e2..d9189345c0 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundary.hh @@ -6,6 +6,7 @@ #ifndef __Spheral_InfinitePlaneSolidBoundary_hh__ #define __Spheral_InfinitePlaneSolidBoundary_hh__ +#include "DEM/DEMDimension.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" namespace Spheral { @@ -20,12 +21,14 @@ class InfinitePlaneSolidBoundary : public SolidBoundaryBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; + typedef typename DEMDimension::AngularVector RotationType; public: //--------------------------- Public Interface ---------------------------// InfinitePlaneSolidBoundary(const Vector& point, - const Vector& normal); + const Vector& normal, + const RotationType& angularVelocity); ~InfinitePlaneSolidBoundary(); @@ -48,6 +51,9 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + const RotationType& angularVelocity() const; + void angularVelocity(const RotationType& value); + virtual std::string label() const override { return "InfinitePlaneSolidBoundary" ; } virtual void dumpState(FileIO& file, const std::string& pathName) const override; virtual void restoreState(const FileIO& file, const std::string& pathName) override; @@ -57,6 +63,7 @@ protected: Vector mPoint; Vector mNormal; Vector mVelocity; + RotationType mAngularVelocity; private: //--------------------------- Private Interface ---------------------------// diff --git a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundaryInline.hh b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundaryInline.hh index 319aa60708..a169fc5617 100644 --- a/src/DEM/SolidBoundary/InfinitePlaneSolidBoundaryInline.hh +++ b/src/DEM/SolidBoundary/InfinitePlaneSolidBoundaryInline.hh @@ -16,6 +16,7 @@ point(const typename Dimension::Vector& value) { mPoint=value; } + template inline const typename Dimension::Vector& @@ -32,6 +33,7 @@ normal(const typename Dimension::Vector& value) { mNormal = value; } + template inline const typename Dimension::Vector& @@ -48,4 +50,21 @@ velocity(const typename Dimension::Vector& value) { mVelocity=value; } + +template +inline +const typename DEMDimension::AngularVector& +InfinitePlaneSolidBoundary:: +angularVelocity() const { + return mAngularVelocity; +} + +template +inline +void +InfinitePlaneSolidBoundary:: +angularVelocity(const typename DEMDimension::AngularVector& value) { + mAngularVelocity=value; +} + } \ No newline at end of file diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc index c7c5bc2806..bc97c2e524 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.cc @@ -19,12 +19,16 @@ namespace Spheral { template RectangularPlaneSolidBoundary:: -RectangularPlaneSolidBoundary(const Vector& point, const Vector& extent, const Tensor& basis): +RectangularPlaneSolidBoundary(const Vector& point, + const Vector& extent, + const Tensor& basis, + const RotationType& angularVelocity): SolidBoundaryBase(), mPoint(point), mBasis(basis), mExtent(extent), - mVelocity(Vector::zero()){ + mVelocity(Vector::zero()), + mAngularVelocity(angularVelocity){ } template @@ -32,7 +36,6 @@ RectangularPlaneSolidBoundary:: ~RectangularPlaneSolidBoundary(){ } - template typename Dimension::Vector RectangularPlaneSolidBoundary:: @@ -46,7 +49,10 @@ template typename Dimension::Vector RectangularPlaneSolidBoundary:: localVelocity(const Vector& position) const { - return mVelocity; + // Calculate the velocity due to angular motion + const auto r = position - mPoint; + const auto angularVelocityContribution = DEMDimension::cross(mAngularVelocity,r); + return mVelocity + angularVelocityContribution; } template @@ -57,8 +63,10 @@ registerState(DataBase& dataBase, const auto boundaryKey = "RectangularPlaneSolidBoundary_" + std::to_string(std::abs(this->uniqueIndex())); const auto pointKey = boundaryKey +"_point"; const auto velocityKey = boundaryKey +"_velocity"; + const auto angularVelocityKey = boundaryKey + "_angularVelocity"; // New key for angular velocity state.enroll(pointKey,mPoint); state.enroll(velocityKey,mVelocity); + state.enroll(angularVelocityKey, mAngularVelocity); // Enroll angular velocity } template void @@ -78,6 +86,7 @@ dumpState(FileIO& file, const string& pathName) const { file.write(mBasis, pathName + "/basis"); file.write(mExtent, pathName + "/extent"); file.write(mVelocity, pathName + "/velocity"); + file.write(mAngularVelocity, pathName + "/omega"); // Write angular velocity } @@ -89,6 +98,7 @@ restoreState(const FileIO& file, const string& pathName) { file.read(mBasis, pathName + "/basis"); file.read(mExtent, pathName + "/extent"); file.read(mVelocity, pathName + "/velocity"); + file.read(mAngularVelocity, pathName + "/omega"); // Read angular velocity } diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh index 9a8f5d4fc6..3ca500162d 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundary.hh @@ -8,6 +8,7 @@ #ifndef __Spheral_RectangularPlaneSolidBoundary_hh__ #define __Spheral_RectangularPlaneSolidBoundary_hh__ +#include "DEM/DEMDimension.hh" #include "DEM/SolidBoundary/SolidBoundaryBase.hh" namespace Spheral { @@ -22,13 +23,15 @@ class RectangularPlaneSolidBoundary : public SolidBoundaryBase { typedef typename Dimension::Scalar Scalar; typedef typename Dimension::Vector Vector; typedef typename Dimension::Tensor Tensor; + typedef typename DEMDimension::AngularVector RotationType; public: //--------------------------- Public Interface ---------------------------// RectangularPlaneSolidBoundary(const Vector& point, - const Vector& exent, - const Tensor& basis); + const Vector& exent, + const Tensor& basis, + const RotationType& angularVelocity); ~RectangularPlaneSolidBoundary(); @@ -54,6 +57,9 @@ public: const Vector& velocity() const; void velocity(const Vector& value); + const RotationType& angularVelocity() const; + void angularVelocity(const RotationType& value); + virtual std::string label() const override { return "RectangularPlaneSolidBoundary" ; } virtual void dumpState(FileIO& file, const std::string& pathName) const override; virtual void restoreState(const FileIO& file, const std::string& pathName) override; @@ -64,6 +70,7 @@ protected: Tensor mBasis; Vector mExtent; Vector mVelocity; + RotationType mAngularVelocity; private: //--------------------------- Private Interface ---------------------------// diff --git a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundaryInline.hh b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundaryInline.hh index eceff05002..671fac8785 100644 --- a/src/DEM/SolidBoundary/RectangularPlaneSolidBoundaryInline.hh +++ b/src/DEM/SolidBoundary/RectangularPlaneSolidBoundaryInline.hh @@ -16,6 +16,7 @@ point(const typename Dimension::Vector& value) { mPoint=value; } + template inline const typename Dimension::Vector& @@ -32,6 +33,7 @@ extent(const typename Dimension::Vector& value) { mExtent=value; } + template inline const typename Dimension::Tensor& @@ -65,4 +67,21 @@ velocity(const typename Dimension::Vector& value) { mVelocity=value; } + +template +inline +const typename DEMDimension::AngularVector& +RectangularPlaneSolidBoundary:: +angularVelocity() const { + return mAngularVelocity; +} + +template +inline +void +RectangularPlaneSolidBoundary:: +angularVelocity(const typename DEMDimension::AngularVector& value) { + mAngularVelocity=value; +} + } \ No newline at end of file diff --git a/src/NodeGenerators/GenerateDEMfromSPHGenerator.py b/src/NodeGenerators/GenerateDEMfromSPHGenerator.py index dcf37960cd..4141fb7588 100644 --- a/src/NodeGenerators/GenerateDEMfromSPHGenerator.py +++ b/src/NodeGenerators/GenerateDEMfromSPHGenerator.py @@ -32,7 +32,15 @@ def constantRadiusFunc(position): # set up our initial radius #-------------------------------------------------- - if type(particleRadius) in [float,int]: + if isinstance(particleRadius, list): + radii = list(particleRadius) + counter = {"i" : 0} + def radiusFunc(position): + r = radii[counter["i"]] + counter["i"] += 1 + return r + self.particleRadiusFunc = radiusFunc + elif type(particleRadius) in [float,int]: self.particleRadiusFunc = constantRadiusFunc elif hasattr(particleRadius,"__call__"): self.particleRadiusFunc = particleRadius diff --git a/src/PYB11/DEM/SolidBoundaries.py b/src/PYB11/DEM/SolidBoundaries.py index 23fda6223e..d814391409 100644 --- a/src/PYB11/DEM/SolidBoundaries.py +++ b/src/PYB11/DEM/SolidBoundaries.py @@ -29,11 +29,13 @@ class InfinitePlaneSolidBoundary(SolidBoundaryBase): PYB11typedefs = """ typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; + typedef typename DEMDimension<%(Dimension)s>::AngularVector RotationType; """ def pyinit(self, point = "const Vector&", - normal = "const Vector&"): + normal = "const Vector&", + angularVelocity = "const RotationType&"): "solid planar boundary" @PYB11virtual @@ -54,7 +56,7 @@ def update(self, @PYB11virtual @PYB11const def localVelocity(self, - position = "const Vector&"): + position = "const Vector&"): "velocity of bc." return "Vector" @@ -68,6 +70,7 @@ def distance(self, velocity = PYB11property("const Vector&", "velocity", "velocity", returnpolicy="reference_internal", doc="velocity of plane") point = PYB11property("const Vector&", "point", "point", returnpolicy="reference_internal", doc="point in plane definition") normal = PYB11property("const Vector&", "normal", "normal",returnpolicy="reference_internal", doc="normal in plane definition") + angularVelocity = PYB11property("const RotationType&", "angularVelocity", "angularVelocity",returnpolicy="reference_internal",doc="angular velocity of plane") #------------------------------------------------------------------------------- # Finite rectangular solid boundary @@ -80,12 +83,14 @@ class RectangularPlaneSolidBoundary(SolidBoundaryBase): typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; typedef typename %(Dimension)s::Tensor Tensor; + typedef typename DEMDimension<%(Dimension)s>::AngularVector RotationType; """ def pyinit(self, point = "const Vector&", extent = "const Vector&", - basis = "const Tensor&"): + basis = "const Tensor&", + angularVelocity = "const RotationType&"): "solid planar boundary" @PYB11virtual @@ -121,6 +126,7 @@ def distance(self, point = PYB11property("const Vector&", "point", "point", returnpolicy="reference_internal", doc="point in plane definition") extent = PYB11property("const Vector&", "extent", "extent", returnpolicy="reference_internal", doc="extent of rectangle") basis = PYB11property("const Tensor&", "basis", "basis",returnpolicy="reference_internal", doc="basis vectors for rectangle") + angularVelocity = PYB11property("const RotationType&", "angularVelocity", "angularVelocity",returnpolicy="reference_internal",doc="angular velocity of plane") #------------------------------------------------------------------------------- # Disk shaped solid boundary @@ -132,12 +138,14 @@ class CircularPlaneSolidBoundary(SolidBoundaryBase): PYB11typedefs = """ typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; + typedef typename DEMDimension<%(Dimension)s>::AngularVector RotationType; """ def pyinit(self, point = "const Vector&", normal = "const Vector&", - extent = "const Scalar"): + extent = "const Scalar", + angularVelocity = "const RotationType&"): "solid planar boundary" @PYB11virtual @@ -173,6 +181,7 @@ def distance(self, point = PYB11property("const Vector&", "point", "point", returnpolicy="reference_internal", doc="point in plane definition") extent = PYB11property("Scalar", "extent", "extent", returnpolicy="reference_internal", doc="extent of rectangle") normal = PYB11property("const Vector&", "normal", "normal",returnpolicy="reference_internal", doc="normal in plane definition") + angularVelocity = PYB11property("const RotationType&", "angularVelocity", "angularVelocity", returnpolicy="reference_internal", doc="angular velocity of plane") #------------------------------------------------------------------------------- # Cylinder solid boundary. In 2d this would be two planes. @@ -184,13 +193,15 @@ class CylinderSolidBoundary(SolidBoundaryBase): PYB11typedefs = """ typedef typename %(Dimension)s::Scalar Scalar; typedef typename %(Dimension)s::Vector Vector; + typedef typename DEMDimension<%(Dimension)s>::AngularVector RotationType; """ def pyinit(self, point = "const Vector&", axis = "const Vector&", radius = "const Scalar", - length = "const Scalar"): + length = "const Scalar", + angularVelocity = "const RotationType&"): "solid planar boundary" @PYB11virtual @@ -227,6 +238,7 @@ def distance(self, axis = PYB11property("const Vector&", "axis", "axis", returnpolicy="reference_internal", doc="extent of rectangle") radius = PYB11property("Scalar", "radius", "radius", doc="normal in plane definition") length = PYB11property("Scalar", "length", "length", doc="normal in plane definition") + angularVelocity = PYB11property("const RotationType&", "angularVelocity", "angularVelocity",returnpolicy="reference_internal",doc="angular velocity of plane") #------------------------------------------------------------------------------- diff --git a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py index 22a4e0f39e..33169bc534 100644 --- a/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/MovingSolidBoundaries/singleParticleMovingBoundaryCollision-3d.py @@ -183,7 +183,9 @@ packages = [dem] # set the solid bcs and add to dem package -solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + Vector(0.0, 0.0, 1.0), + Vector(0.0, 0.0, 0.0)) solidWall.velocity = Vector(0.0,0.0,vbc) dem.appendSolidBoundary(solidWall) diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py index dbdbb6ff7e..f9fedcbcdb 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/finiteSolidPlanes-3d.py @@ -193,15 +193,22 @@ packages = [dem] if planeType == "infinite": - solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) + solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + Vector(0.0, 0.0, 1.0), + Vector(0.0, 0.0, 0.0)) elif planeType == "circular": - solidWall = CircularPlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0),0.25) + solidWall = CircularPlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + Vector(0.0, 0.0, 1.0), + 0.25, + Vector(0.0, 0.0, 0.0)) elif planeType == "rectangular": basis = Tensor(0.0,0.0,1.0, 1.0,0.0,0.0, 0.0,1.0,0.0,) extent = Vector(0.0,0.25,0.25) - solidWall = RectangularPlaneSolidBoundary(Vector(0.0, 0.0, 0.0), extent, basis) + solidWall = RectangularPlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + extent, basis, + Vector(0.0, 0.0, 0.0)) dem.appendSolidBoundary(solidWall) diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py index e34acc3020..8aa379ee48 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/pureTorsion-3d.py @@ -195,7 +195,9 @@ def DEMParticleGenerator(xi,yi,zi,Hi,mi,Ri): # implement boundary condition using the DEM packages solid wall feature if useSolidBoundary: - solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) + solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + Vector(0.0, 0.0, 1.0), + Vector(0.0, 0.0, 0.0)) dem.appendSolidBoundary(solidWall) # implement boundary condition using Spheral's ghost particle reflection diff --git a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py index 2d40345772..5e90d97023 100644 --- a/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py +++ b/tests/functional/DEM/LinearSpringDEM/SolidBoundaryCondition/singleParticleBoundaryCollision-3d.py @@ -201,7 +201,9 @@ packages = [dem] -solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), Vector( 0.0, 0.0, 1.0)) +solidWall = InfinitePlaneSolidBoundary(Vector(0.0, 0.0, 0.0), + Vector(0.0, 0.0, 1.0), + Vector(0.0, 0.0, 0.0)) dem.appendSolidBoundary(solidWall) #-------------------------------------------------------------------------------