CODEWITHSUNDEEP
c++c++11
Gilfoyle
Gilfoyle

Reputation: 3616

Nested loop over std::vector

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?

Upvotes: 2

Views: 119

Answers (1)

Ton van den Heuvel
Ton van den Heuvel

Reputation: 10528

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.

Upvotes: 3

Related Questions