Nested loop over std::vector

Question

In my simulation I have a particle class and a nbody class.

struct Particle{
    double m;           // mass
    double x[DIM];      // position
    double v[DIM];      // velocity
    double F[DIM];      // force 
    double F_old[DIM];  // force past time step
    void update_position(double dt);
    void update_velocity(double dt);
    void update_force(Particle*, Particle*);
};

class Nbody {
    private:
        const double dt;                // step size
        const double t_max;             // max simulation time
        std::vector<Particle> p;
    public:
        Nbody(unsigned int n_, double dt_, double t_max_);
        void comp_force();
        void comp_position();
        void comp_velocity();
};

I update the position of every particle by looping over every particle and calling void Particle::update_position().

void Nbody::comp_position() {
    for (Particle& particle: p) { 
        particle.update_position(dt);
    }
}

In the following code snippet, I show the update rule to compute the new position.

void Particle::update_position(double dt) {
    const double a = dt * 0.5 / m;
    for (unsigned int d=0; d<DIM; ++d) {
        x[d] += dt * (v[d] + a * F[d]);
        F_old[d] = F[d];
    }
}

However, I have problems with the force calculation where I want to compute the interaction of particles and therefore need a nested loop. How can I do that? Here is my approach so far, which gives me an error:

void Particle::update_force(Particle* i, Particle* j) {
    double r=EPS; // smoothing
    for (unsigned int d=0; d<DIM; d++) {
        r += sqr(j->x[d] - i->x[d]);
    }
    const double f = i->m * j->m / (sqrt(r) * r);
    for (unsigned int d=0; d<DIM; d++) {
        i->F[d] += f * (j->x[d] - i->x[d]);
        j->F[d] -= f * (j->x[d] - i->x[d]);
    }
}

void Nbody::comp_force() {
    for (Particle& particle: p) {
        for (unsigned int d=0; d<DIM; ++d) {
            particle.F[d] = 0.;
        }
    }
    for (unsigned int i=0; i<n; i++) {
        for (unsigned int j=i+1; j<n; j++) {
            update_force(&p[i], &p[j]);
        }
    }
}

Here is the error I get:

nbody.cpp: In member function ‘void Nbody::comp_force()’:
nbody.cpp:84:13: error: ‘update_force’ was not declared in this scope
             update_force(&p[i], &p[j]);
             ^~~~~~~~~~~~

Can I write this similar as the update_position method?

Suggestion: Your position, velocity and force could all benefit from a vector class (not std::vector but a math vector). template <size_t Dim> class Vector {... — Ted Lyngmo
– Ted Lyngmo, Commented Oct 25, 2019 at 16:33
Particle::update_force is being called from Nbody::comp_force. Note the differing scopes. — user4581301
– user4581301, Commented Oct 25, 2019 at 16:38
Looks like Particle::update_force could be made static (or a friend free function, but static seems more appealing to me here). — user4581301
– user4581301, Commented Oct 25, 2019 at 16:41
@TedLyngmo Could you elaborate on your suggestion a little? I can't imagine how I could implement that. — Gilfoyle
– Gilfoyle, Commented Oct 25, 2019 at 17:33

Ton van den Heuvel · Accepted Answer · 2019-10-25 19:06:41Z

3

You are calling update_force defined as a member function on a Particle from a member function of NBody. There is no function update_force in the scope of NBody, hence the error. My suggestion would be to not define update_force as a member function on a particle, but as a free function outside Particle:

struct Particle{
    double m;           // mass
    double x[DIM];      // position
    double v[DIM];      // velocity
    double F[DIM];      // force 
    double F_old[DIM];  // force past time step
    void update_position(double dt);
    void update_velocity(double dt);
};

...

void update_particle_force(Particle& i, Particle& j) {
    double r=EPS; // smoothing
    for (unsigned int d=0; d<DIM; d++) {
        r += sqr(j.x[d] - i.x[d]);
    }
    const double f = i.m * j.m / (sqrt(r) * r);
    for (unsigned int d=0; d<DIM; d++) {
        i.F[d] += f * (j.x[d] - i.x[d]);
        j.F[d] -= f * (j.x[d] - i.x[d]);
    }
}

which is then called as:

...
update_particle_force(p[i], p[j]);
...

The advantage of using references over pointers for the parameters of update_particle_force is that to a certain extent you will avoid null pointer dereferences; the potential burden of the null pointer check is now on the caller of update_particle_force.

edited Oct 25, 2019 at 19:06

answered Oct 25, 2019 at 18:57

Ton van den Heuvel

10.6k6 gold badges49 silver badges92 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

user4581301 Over a year ago

Huh. It didn't click that Particle was a struct That eliminates my objections to update_force being a free function.

Collectives™ on Stack Overflow

Nested loop over std::vector

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related