Assigning values to objects in C++

Question

I can be considered a beginner in programming. I am trying to write a class in C++ that can hold multi-dimensional data (such as an MxN matrix). I don't want to do it via the way of vector of a vector. I wrote the following piece of code, but when I compile it with g++, I get the segmentation fault: 11 error. When I try to compile it with Xcode on the other hand, it gives "Thread 1: EXC_BAD_ACCESS (code = 1, address = 0x0)" at the line

particles[i].x = x0 + R*cos(theta*i);

So, I guess I'm getting a sort of memory error related to assigning objects a value.

Is this kind of syntax even allowed in C++:

particles_old[i].x = particles[i].x;

Or can I something like this:

// say A is a class with x a vector this time, instead of a double data. 
// in nested for loops
vector<A> B;
B[i].x[j] = some value;

I know it's a bit vague, but is it at least correct from a syntax point of view?

#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
#include <vector>
#include <random>
using namespace std;

class Particle{
public:
    double x; // x position
    double y; // y position
    double vx; // velocity in the x direction
    double vy; // velocity in the y direction
    double Fx; // force in the x direction
    double Fy; // force in the y direction

    // Default constructor
    Particle()
    : x(0.0),y(0.0),vx(0.0),vy(0.0),Fx(0.0),Fy(0.0){
    }
};

int main() {

    const float pi = 3.14;
    int N = 30; // Number of 'particles' that make up the cell
    float theta = 2*pi/N; // Angle between two particles in radians
    float x0 = 0; // Center of the cell [x]
    float y0 = 0; // Center of the cell [y]
    float R = 5e-6; // Radius of the cell
    vector<Particle> particles; // particles

    // Assigning the initial points onto the circle
    for(int i = 0; i < N; i++) {
        particles[i].x = x0 + R*cos(theta*i);
        particles[i].y = y0 + R*sin(theta*i);
    }

    float k = 4.3e-7; // Spring constant connecting the particles
    float m = 2e-8; // Mass of the particles

    // Calculating the initial spring force between the particles on the cell
    particles[0].Fx = -k*(particles[1].x - particles[N].x);
    particles[0].Fy = -k*(particles[1].y - particles[N].y);
    for(int i = 1; i < N-1; i++) {
        particles[i].Fx = -k*(particles[i+1].x - particles[i-1].x);
        particles[i].Fy = -k*(particles[i+1].y - particles[i-1].y);
    }
    particles[N].Fx = -k*(particles[0].x - particles[N-1].x);
    particles[N].Fy = -k*(particles[0].y - particles[N-1].y);

    // Initial velocities are given to each particle randomly from a Gaussian distribution
    random_device rdx; // Seed
    default_random_engine generatorx(rdx()); // Default random number generator
    random_device rdy; // Seed
    default_random_engine generatory(rdy()); // Default random number generator
    normal_distribution<float> distributionx(0,1); // Gaussian distribution with 0 mean and 1 variance
    normal_distribution<float> distributiony(0,1); // Gaussian distribution with 0 mean and 1 variance
    for(int i = 0; i < N; i++) {
        float xnumber = distributionx(generatorx);
        float ynumber = distributiony(generatory);
        particles[i].vx = xnumber;
        particles[i].vy = ynumber;
    }


    // Molecular dynamics simulation with velocity Verlet algorithm

    // 'Old' variables
    vector<Particle> particles_old;

    for(int i = 0; i < N; i++) {
        particles_old[i].x = particles[i].x;
        particles_old[i].y = particles[i].y;
        particles_old[i].vx = particles[i].vx;
        particles_old[i].vy = particles[i].vy;
        particles_old[i].Fx = particles[i].Fx;
        particles_old[i].Fy = particles[i].Fy;
    }

    // Sampling variables
    int sampleFreq = 2;
    int sampleCounter = 0;

    // MD variables
    float dt = 1e-7;
    float dt2 = dt*dt;
    float m2 = 2*m;
    int MdS = 1e+4; // Molecular dynamics step number

    // MD
    for(int j = 0; j < MdS; j++) {

        // Update x
        for(int i = 0; i < N; i++) {
            particles[i].x = particles_old[i].x + dt*particles_old[i].vx + dt2*particles_old[i].Fx/m2;
            particles[i].y = particles_old[i].y + dt*particles_old[i].vy + dt2*particles_old[i].Fy/m2;
        }

        // Update F
        particles[0].Fx = -k*(particles[1].x - particles[N].x);
        particles[0].Fy = -k*(particles[1].y - particles[N].y);
        for(int i = 1; i < N-1; i++) {
            particles[i].Fx = -k*(particles[i+1].x - particles[i-1].x);
            particles[i].Fy = -k*(particles[i+1].y - particles[i-1].y);
        }
        particles[N].Fx = -k*(particles[0].x - particles[N-1].x);
        particles[N].Fy = -k*(particles[0].y - particles[N-1].y);

        // Update v
        for(int i = 0; i < N; i++) {
            particles[i].vx = particles_old[i].vx + dt*(particles_old[i].Fx + particles[i].Fx)/m2;
            particles[i].vy = particles_old[i].vy + dt*(particles_old[i].Fy + particles[i].Fy)/m2;
        }

        // Copy new variables to old variables
        for(int i = 0; i < N; i++) {
            particles_old[i].x = particles[i].x;
            particles_old[i].y = particles[i].y;
            particles_old[i].vx = particles[i].vx;
            particles_old[i].vy = particles[i].vy;
            particles_old[i].Fx = particles[i].Fx;
            particles_old[i].Fy = particles[i].Fy;
        }

    }

}

Thanks in advance.

Answer 1

// Assigning the initial points onto the circle
for(int i = 0; i < N; i++) {
    particles[i].x = x0 + R*cos(theta*i);
    particles[i].y = y0 + R*sin(theta*i);
}

At the point you enter that loop, your vector is empty. You then attempt to access elements in it using operator[]. If you look at information on std::vector operator[] it states the following:

"Portable programs should never call this function with an argument n that is out of range, since this causes undefined behavior."

You are doing precisely this, accessing an element that is out of range, and in this case the undefined behavior is a crash.

Change that loop to:

// Assigning the initial points onto the circle
for(int i = 0; i < N; i++) {
    // Create and initialize a new particle
    Particle newParticle;
    newParticle.x = x0 + R*cos(theta*i);
    newParticle.y = y0 + R*sin(theta*i);
    // add it to the vector
    particles.push_back(newParticle);
}

That explicitly creates a new Particle for you, and then adds it to your std::vector with the .push_back() method.

Note that these two lines:

particles[N].Fx = -k*(particles[0].x - particles[N-1].x);
particles[N].Fy = -k*(particles[0].y - particles[N-1].y);

will probably crash for the exact same reason, since at that point you've populated exactly [0 .. N-1], you have not put an element in slot N. I'll leave you to fix that, since the solution is exactly the same as in the for loop.

Answer 2

Hey you need to initialize your particles when adding them to your list like below for instance

for (int i = 0; i < N+1; i++) { // You are accessing [N] directly in your code..
    Particle p;
    p.x = x0 + R*cos(theta*i);
    p.y = y0 + R*sin(theta*i);
    particles.push_back(p);
}

I simplified the code above a bit, since you use [N] directly the +1 is there, you should really do a [N-1] instead in your code, or you will go outside your vector.

You push your particle to the vector, at the end of the vector, remember your vector is empty! I'm sure this can be done better, to be honest I'm a bit rusty here my self :) Picking up my "baby language" and all the cool features that have been implemented there.

I'm not sure if I followed your questions, but if you want to have a look into simple 2 dimentional arrays there is a rather usefull post here at SO.

Well anyhow, I hope it helps you.

Answer 3

You have already two excellent answers. Here is another variant:

If you know in advance the size of your vector, you may as well define it as follows:

vector<Particle> particles(N);  // create a vector with N particles default initialized. 

The rest of your code should then work, because particles[i] will refer to a particle that already exists.

And yes, old_particles = particles; will copy the whole vector. No need to copyelement by element in atedious loop.

Answer 4

In addition to other answers, I would do either

particles.reserve(N); (or N+1) before creating particles, so the vector doesn't have to reallocate memory during the loop

or even declare particles as vector<Particle> particles(N); (or again N+1) so the vector is already filled with default-initialized particles, and you just assign their values like you did before.

Choosing N or N+1 depends on the line

particles[0].Fx = -k*(particles[1].x - particles[N].x);

you seem to access Nth element here, and if it's not a typo, the vector size should be N + 1.