c++ - 气体粒子模拟碰撞计算C++

Question

这个我觉得是个比较复杂的问题，我希望能把它放到足够小的篇幅里，让大家看得懂。我目前正在编写代码

simulate Ideal gas particles inside a box. I'm calculating if two particles will collide having calculated the time taken for them to reach their closest point. (using an example where they have head on collision).

在这部分代码中，我需要先确定两个粒子是否会发生碰撞，然后再计算碰撞时间和碰撞方式等。因此对于我的两个粒子:

Main.cpp

Vector vp1(0,0,0);
Vector vv1(1,0,0);
Vector vp2(12,0,0);
Vector vv2(-1,0,0);
Particle Particle1(1, vp1, vv1);
Particle Particle2(1, vp2, vv2);
Particle1.timeToCollision(Particle2);

在我的程序中，我将粒子定义为:

头文件

class Particle {
private:

    Vector p;   //position
    Vector v;   //velocity
    double radius;   //radius

public:

    Particle();

    Particle(double r, const Vector Vecp, const Vector Vecv);

    void setPosition(Vector);
    void setVelocity(Vector);

    Vector getPosition() const;
    Vector getVelocity() const;

    double getRadius() const;

    void move(double t);

    double timeToCollision(const Particle particle);

    void collideParticles(Particle);

    ~Particle();
};

Vector是另一个类，简而言之 x , y , z值。它还包含用于操作这些的多个函数。

---

以及我需要帮助的部分，在 .cpp 中(忽略 cout 开始和字母等，它们是简单的检查我的代码退出的位置以进行测试。)

给定方程式:

我已经编写了代码来为我做点积和模数，并且:

在哪里

s是时间旅行的距离 tac .

double Particle::timeToCollision(const Particle particle){
    Vector r2 = particle.getPosition();
    Vector r1 = p;
    Vector v2 = particle.getVelocity();
    Vector v1 = v;
    Vector r0 = r2 - r1;
    Vector v = v2 - v1;

    double modv;
    double tca;
    double result = 0;
    double bsqr;

    modv = getVelocity().modulus();

    cout << "start" << endl;

    if(modv < 0.0000001){
        cout << "a" << endl;
        result = FLT_MAX;
    }else{
        cout << "b" << endl;
        tca = ((--r0).dot(v)) / v.modulusSqr();

// -- is an overridden operator that gives the negation ( eg (2, 3, 4) to (-2, -3, -4) )
        if (tca < 0) {
            cout << "c" << endl;
            result = FLT_MAX;
        }else{
            cout << "d" << endl;
            Vector s(v.GetX(), v.GetY(), v.GetZ());
            s.Scale(tca);
            cout << getVelocity().GetX() << endl;
            cout << getVelocity().GetY() << endl;
            cout << getVelocity().GetZ() << endl;

            double radsqr = radius * radius;
            double bx = (r0.GetX() * r0.GetX() - (((r0).dot(v)) *((r0).dot(v)) / v.modulusSqr()));
            double by = (r0.GetY() * r0.GetY() - (((r0).dot(v)) *((r0).dot(v)) / v.modulusSqr()));
            double bz=(r0.GetZ() * r0.GetZ() - (((r0).dot(v)) * ((r0).dot(v)) / v.modulusSqr()));

            if (bsqr < 4 * radsqr) {
                cout << "e" << endl;
                result = FLT_MAX;
            } else {
            }
            cout << "tca: " << tca << endl;
        }
    }
    cout << "fin" << endl;
    return result;
}

我有计算几个方面的方程式，tca指最接近时间。

如代码中所写，我需要检查是否 b > 4 r^2 ，我做了一些尝试并编写了X , Y和 Z b 的组件出去。但我得到的是垃圾答案。我只需要帮助来确定我是否已经犯了错误或我应该朝着哪个方向前进。

---

在此之前我的所有代码都按预期工作，并且我为每个代码编写了多个测试以进行检查。

如果您觉得我遗漏了任何信息，请在评论中告诉我。

非常感谢任何帮助。

Answer 1

您的代码中有几个错误。你永远不会设置 result到不同于 0 的值或 FLT_MAX .你也永远不会计算bsqr .我想如果 bsqr < 4r^2 就会发生碰撞而不是相反。 (我不明白为什么 4r^2 而不是 r^2 但没关系)。因为你隐藏了你的 vector 实现，所以我使用了一个通用的 vector 库。我还建议无论如何不要使用手工制作的东西。查看 Armadillo 或 Eigen。

在这里，您可以在 Eigen 中尝试一下。

#include <iostream>
#include <limits>
#include <type_traits>
#include "Eigen/Dense"

struct Particle {
    double radius;
    Eigen::Vector3d p;
    Eigen::Vector3d v;
};

template <class FloatingPoint>
    std::enable_if_t<std::is_floating_point<FloatingPoint>::value, bool>
        almost_equal(FloatingPoint x, FloatingPoint y, unsigned ulp=1)
        {
            FloatingPoint max = std::max(std::abs(x), std::abs(y));
            return std::abs(x-y) <= std::numeric_limits<FloatingPoint>::epsilon()*max*ulp;
        }

double timeToCollision(const Particle& left, const Particle& right){
    Eigen::Vector3d r0 = right.p - left.p;
    Eigen::Vector3d v  = right.v - left.v;

    double result = std::numeric_limits<double>::infinity();

    double vv = v.dot(v);
    if (!almost_equal(vv, 0.)) {
        double tca = (-r0).dot(v) / vv;
        if (tca >= 0) {
            Eigen::Vector3d s = tca*v;
            double bb = r0.dot(r0) - s.dot(s);
            double radius = std::max(left.radius, right.radius);
            if (bb < 4*radius*radius)
                result = tca;
        }
    }
    return result;
}

int main()
{
    Eigen::Vector3d vp1 {0,0,0};
    Eigen::Vector3d vv1 {1,0,0};
    Eigen::Vector3d vp2 {12,0,0};
    Eigen::Vector3d vv2 {-1,0,0};
    Particle p1 {1, vp1, vv1};
    Particle p2 {1, vp2, vv2};
    std::cout << timeToCollision(p1, p2) << '\n';
}