ClumpParticle.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef MultiParticle_H
#define MultiParticle_H
 
#include "Mercury3D.h"
#include "Math/SmallVector.h"
#include "BaseParticle.h"
#include "BaseInteractable.h"
#include "NonSphericalParticle.h"
#include "Species/ParticleSpecies.h"
#include "DPMBase.h"
class ClumpParticle final: public NonSphericalParticle
{
public:
    ClumpParticle();
 
    ClumpParticle(const ClumpParticle& p);
    
    ~ClumpParticle() override;
 
    ClumpParticle* copy() const override;
 
 
    void write(std::ostream& os) const override
    {
        BaseParticle::write(os);
        os << " nPebble_ " << nPebble_;
        os << " clumpMass_ " << clumpMass_;
        os << " viscousDamping_ " << viscousDamping_;
        os << " pebblePos_ ";
        for (const auto & pebbleP : pebblePos_) {os << pebbleP;}
        os << " pebbleRadius_ ";
        for (double pebbleRad : pebbleRadius_) {os << pebbleRad;}
        //os << " pebbleParticles_ ";
        //for (auto pebblePart : pebbleParticles_) {os << pebblePart;}
        os << " principalDirections_ ";
        os << principalDirections_;
        os << " initPrincipalDirections_ ";
        os << initPrincipalDirections_;
        os << " clumpInitInertia_ " << clumpInitInertia_;
        os << " clumpInertia_ " << clumpInertia_;
        os << " rotationMatrix_ ";
        os << rotationMatrix_;
        os << " isPebble_ " << isPebble_;
        //os << " clumpParticle_ " << clumpParticle_;
        os << " isPebble_ " << isPebble_;
        os << " isDzhanibekovParticle_ " << isDzhanibekovParticle_;
        os << " isVerticallyOriented_ " << isVerticallyOriented_;
    }
 
    void read(std::istream& is) override
    {
        BaseParticle::read(is);
        std::string dummy;
        is >> dummy >> nPebble_;
        is >> dummy >> nPebble_;
        is >> dummy >> clumpMass_;
        is >> dummy >> viscousDamping_;
        is >> dummy;
        for (int i = 0; i<nPebble_; i++) {is >> pebblePos_[i];}
        is >> dummy;
        for (int i = 0; i<nPebble_; i++) {is >> pebbleRadius_[i];}
        //is >> dummy;
        //for (int i = 0; i<nPebble_; i++) {is >> pebbleParticles_[i];}
        is >> dummy;
        is >> principalDirections_;
        is >> dummy;
        is >> initPrincipalDirections_;
        is >> dummy >> clumpInitInertia_;
        is >> dummy >> clumpInertia_;
        is >> dummy;
        is >> rotationMatrix_;
        is >> dummy >> isPebble_;
        //is >> dummy >> clumpParticle_;
        is >> dummy >> isPebble_;
        is >> dummy >> isDzhanibekovParticle_;
        is >> dummy >> isVerticallyOriented_;
 
 
    }
 
    std::string getName() const override;
 
    bool isSphericalParticle() const override
    {
        return true;
    }
 
    void computeMass(const ParticleSpecies& s) override;
 
    void setInitInertia(MatrixSymmetric3D inertia);
 
    void rotateTensorOfInertia();
 
    void addPebble(Vec3D position, Mdouble radius);
 
    void setClump();
 
    void setPrincipalDirections(Matrix3D directions);
 
    void setInitPrincipalDirections(Matrix3D directions);
 
    Vec3D getPrincipalDirections_e1() const
    {
 
        return Vec3D(principalDirections_.XX, principalDirections_.YX, principalDirections_.ZX);
    }
    Vec3D getPrincipalDirections_e2() const
    {
        return Vec3D(principalDirections_.XY, principalDirections_.YY, principalDirections_.ZY);
    }
    Vec3D getPrincipalDirections_e3() const
    {
        return Vec3D(principalDirections_.XZ, principalDirections_.YZ, principalDirections_.ZZ);
    }
 
    // Methods to obtain initial principal directions
    Vec3D getInitPrincipalDirections_e1() const
    {
        return Vec3D(initPrincipalDirections_.XX, initPrincipalDirections_.YX, initPrincipalDirections_.ZX);
    }
    Vec3D getInitPrincipalDirections_e2() const
    {
        return Vec3D(initPrincipalDirections_.XY, initPrincipalDirections_.YY, initPrincipalDirections_.ZY);
    }
    Vec3D getInitPrincipalDirections_e3() const
    {
        return Vec3D(initPrincipalDirections_.XZ, initPrincipalDirections_.YZ, initPrincipalDirections_.ZZ);
    }
 
    Matrix3D getRotationMatrix() const
    {
        return rotationMatrix_;
    }
 
    int NPebble() const; // Number of pebbles (for a clump particle)
 
    void actionsAfterAddObject() override; // The function that updates clump quantities after adding a pebble
 
    void updatePebblesVelPos();
 
    void integrateBeforeForceComputation(double time, double timeStep) override;
 
    void integrateAfterForceComputation(double time, double timeStep) override;
 
    void angularAccelerateClumpIterative(double timeStep); // Specific method of time integration for a clump particle
 
    void rotatePrincipalDirections(Vec3D rotation);
 
    // Principal direction
    void setPrincipalDirections_e1(Vec3D e)
    {
        principalDirections_.XX = e.X;
        principalDirections_.YX = e.Y;
        principalDirections_.ZX = e.Z;
    }
    void setPrincipalDirections_e2(Vec3D e)
    {
        principalDirections_.XY = e.X;
        principalDirections_.YY = e.Y;
        principalDirections_.ZY = e.Z;
    }
    void setPrincipalDirections_e3(Vec3D e)
    {
        principalDirections_.XZ = e.X;
        principalDirections_.YZ = e.Y;
        principalDirections_.ZZ = e.Z;
    }
 
    std::vector<Mdouble> getPebbleRadius() const {
        return pebbleRadius_;
    }
 
    // add pointer to pebble pointers list
    void setPebble(int kPebble, ClumpParticle* pPebble) {
        pebbleParticles_[kPebble] = pPebble;
    }
 
    // Sets the particle to be a pebble of a given clump
    void setClump(ClumpParticle* master) {
        isClump_ = false;
        isPebble_ = true;
        clumpParticle_ = master;
    }
 
    void setClumpMass(Mdouble mass)
    {
        clumpMass_ = mass;
        invMass_ = 1 / clumpMass_;
    }
 
    // Extra viscous damping on a clump
    void setDamping(Mdouble damp)
    {
        viscousDamping_ = damp;
    }
 
    Mdouble getKineticEnergy() const override{
        Mdouble res = 0;
        if (isClump_) {
            Vec3D v = getVelocity();
            res = 0.5 * clumpMass_ * (v.X * v.X + v.Y * v.Y + v.Z * v.Z );
        }
        return res;
    }
 
    Mdouble getRotationalEnergy() const override{
        Mdouble res = 0;
        if (isClump_) {
            Vec3D nn = getAngularVelocity();
            Mdouble nl = nn.getLength();
            Mdouble tol = 1e-10;
            if (nl > tol) {
                nn /= nl;
                Mdouble ww = getAngularVelocity().getLengthSquared();
                Mdouble II = Vec3D::dot(nn, (getInertia() * nn));
                res = 0.5 * II * ww;
            }
        }
        return res;
 
    }
 
    std::vector <Vec3D> getPebblePositions(){
        std::vector <Vec3D> globalPos;
        Vec3D e1 = getPrincipalDirections_e1();
        Vec3D e2 = getPrincipalDirections_e2();
        Vec3D e3 = getPrincipalDirections_e3();
        for (int i = 1; i <= NPebble(); i++){
        globalPos.push_back(getPosition() + e1 * pebblePos_[i - 1].X + e2 * pebblePos_[i - 1].Y + e3 * pebblePos_[i - 1].Z);
        }
        return globalPos;
    }
 
    std::vector <Mdouble> getPebbleRadii(){ return pebbleRadius_; }
 
    // Methods setting and getting some extra Boolean properties
 
    // check if particle is in "Dzhanibekov" state
    bool getDzhanibekovParticle()
    {
        return isDzhanibekovParticle_;
    }
 
    // check if particle is "vertically oriented"
    bool getVerticallyOriented()
    {
        return isVerticallyOriented_;
    }
 
    // set the "Dzhanibekov" state
    void setDzhanibekovParticle( bool d)
    {
        isDzhanibekovParticle_ = d;
    }
 
    // set "vertically oriented" state
    void setVerticallyOriented( bool d)
    {
        isVerticallyOriented_ = d;
    }
 
    unsigned getNumberOfFieldsVTK() const override
    {
        return 2;
    }
 
    std::string getTypeVTK(unsigned i) const override
    {
        return "Int8";
    }
 
 
    std::string getNameVTK(unsigned i) const override
    {
        if (i==0)
            return "DzhanibekovParticle";
        else
            return "VerticallyOriented";
    }
 
 
    std::vector<Mdouble> getFieldVTK(unsigned i) const override
    {
        if (i==0)
            return std::vector<Mdouble>(1, isDzhanibekovParticle_);
        else
            return std::vector<Mdouble>(1, isVerticallyOriented_);
    }
 
    void updateExtraQuantities();
        
 
private:
 
    int nPebble_;               // Number of pebbles
    
    bool isDzhanibekovParticle_;  // This property is needed to quantify Dzhanibekov gas properties
 
    bool isVerticallyOriented_;   // This property is useful for mechnical stability simulations (Gomboc, Dominos)
 
    Vec3D angularAcceleration_; // Clump angular acceleration
 
    Mdouble clumpMass_;         // Clump mass
 
    Mdouble viscousDamping_;    // Viscous damping parameter - for extra clump damping
 
    MatrixSymmetric3D clumpInertia_;        // Clump tensor of inertia
 
    MatrixSymmetric3D clumpInitInertia_;    // Clump initial tensor of inertia
 
    std::vector<Vec3D> pebblePos_;          // positions of pebbles in a local coordinate system
 
    std::vector<Mdouble> pebbleRadius_;     // radii of pebbles
 
    Matrix3D principalDirections_;          // clump's principal directions
 
    Matrix3D initPrincipalDirections_;      // clump's initial principal directions
 
    std::vector<ClumpParticle*> pebbleParticles_; // pointers to pebbles
 
    Matrix3D rotationMatrix_; // Rotation matrix
 
    //Helper functions
 
    // Converts a Matrix3D into a MatrixSymmetric3D
    MatrixSymmetric3D MtoS( Matrix3D M){ return MatrixSymmetric3D(M.XX, M.XY, M.XZ, M.YY, M.YZ, M.ZZ);}
 
    // Converts a MatrixSymmetric3D into a Matrix3D
    Matrix3D StoM( MatrixSymmetric3D M){ return Matrix3D(M.XX, M.XY, M.XZ, M.XY, M.YY, M.YZ, M.XZ, M.YZ, M.ZZ);}
 
    // Transposes the matrix
    Matrix3D transpose(Matrix3D M){ return Matrix3D(M.XX, M.YX, M.ZX, M.XY, M.YY, M.ZY, M.XZ, M.YZ, M.ZZ);}
};
 
#endif