TerrainContact.h

Go to the documentation of this file.

// SPDX-License-Identifier: MIT
// Copyright (c) 2021-2025 MPM-Geomechanics Development Team
 
#ifndef TERRAIN_CONTACT_H
#define TERRAIN_CONTACT_H
 
#include "Mesh/STLReader.h"
#include "Mesh/Mesh.h"
 
class TerrainContact {
 
private:
    
    STLReader* stlMesh; 
    
    double frictionCoefficient; 
 
    std::vector<double> densityLevelSet; 
 
    std::vector<std::pair<Particle*, Triangle*>> contactPairs; 
 
    double scalingFactor = 2.0; 
 
    bool usePenaltyContact = false; 
 
    double penaltyStiffness = 0.0; 
 
public:
 
    inline std::vector<std::pair<Particle*, Triangle*>> getContactPairs() const { return contactPairs; }
 
    TerrainContact(STLReader* mesh, double friction)
        : stlMesh(mesh), frictionCoefficient(friction) {
            // initialize density level set
            size_t npoints = mesh->getTriangles().size();
            densityLevelSet.resize(npoints, 0.0);
        }
 
    void computeDistanceLevelSetFunction(Mesh* mesh);
 
    void particleDistanceLevelSet(Mesh* mesh, vector< Particle* >* particles);
 
    void nodalDensityLevelSet(Mesh* mesh, vector< Particle* >* particles);
 
    void trianglesDensityLevelSet(Mesh* mesh);
 
    void determineContactPotentialPairs(Mesh* mesh, std::vector< Particle* >* particles);
 
    STLReader* getSTLMesh() { return stlMesh; }
 
    const std::vector<double>& getTriangleDensityLevelSet() const { return densityLevelSet; }
 
    // \brief Compute the contact forces
    // \f$ f_n = -m_p \frac{v_p^n}{dt} e_n \f$
    // \f$ f_t = -m_p \frac{v_p - v_p^n e_n}{dt} \f$
    void computeContactForces(double dt);
 
    // \brief Project particles that have penetrated the STL mesh
    // If a particle is inside the STL mesh, move it to the surface
    // \f$ x_projected = x_p - d * e_n \f$
    void projectParticles(Mesh *mesh, std::vector<Particle *> *particles);
 
    void setDistanceThreshold(double threshold) { scalingFactor = threshold > 0.0 ? threshold : 2.0; };
 
    void apply(Mesh* mesh, std::vector<Particle*>* particles, double dt);
 
    void enablePenaltyContact(bool enable);
    
    void setPenaltyStiffness(double k);
};
 
#endif // TERRAINCONTACT_H