C++ matrix boundary condition problems

[rolotomassi]

I have created a matrix with a class called Lattice. The lattice is filled with objects of type 'Dipole' which is created with another class. The problem I am having is with boundary conditions when I look for a neighbour. e.g If i pick a dipole on the top row, I want its 'above' neighbour to be the site on the bottom row etc...

Below is the .cpp and .h file which contains the functions I am having problems with. Note that the value passed to both functions can be negative. I thought the %xSize should take care of this but still running into problems. I believe the dipole functions are working correctly.

#ifndef lattice
#define lattice
#include "dipole.h"

class Lattice{

    private:
 
    Dipole* dipoleArray;
 
    double E;    // pointing in +ve y-axis V/m
 
    public:
        const int xSize;
        const int ySize;
        Lattice(const int, const int);
        ~Lattice();
        Dipole* GetDipole(int x, int y);  //  return a dipole
        double calcLocalEnergy( int x, int y);
 
 
};
#endif///////////////////////////////////////////////////////////#include "dipole.h"
#include "lattice.h"
#include <cmath>
#include <cstdlib>
#include<iostream>
Lattice::Lattice(const int xSize, const int ySize): xSize(xSize), ySize(ySize){
    dipoleArray =  new Dipole[xSize*ySize];    //store 1D array list
 
    for(int i=0; i<xSize; i++){
        for(int j=0; j<ySize; j++){
            //loop through the entire array and create dipoles
            dipoleArray[ i+j*xSize ] = Dipole();  // i + j*xSize to mimic 2D lattice form
         
        }
    }
}

Dipole* Lattice::GetDipole(int x, int y){

    return &dipoleArray[ (x%xSize) + (y%ySize)*xSize  ];    // modulo is for loop behaviour
}
double Lattice::calcLocalEnergy( int x, int y){
     
    double k = 1/(4*3.14158*8.85E-12);     double p1=0;     double p2 =0;        double r  = 1E-8;     double energy = 0;
    double E =1e4;    double pZero = 1e-27;    double pOne = 1e-29;
 
//    std::cout<<"0"<<std::endl;
    Dipole* central = GetDipole(x,y);  // define central dipole in the lattice at position [x,y] given by the argument
                                                    // take modulo because x,y can potentially  = -1 OR xSize for neighbour dipoles
                                                    // modulo ensures neighbours are found correctly
                                                    // for 0 <=  x,y < xSize modulo has no effect 
                                                     
    Dipole* neighbour[4];
    neighbour[0] = GetDipole(x,y-1);         //up = 0, right = 1, down = 2, left = 3
    neighbour[1] = GetDipole(x+1,y);
    neighbour[2] = GetDipole(x,y+1);
    neighbour[3] = GetDipole(x-1,y);
    int myDirection =  central->GetDirection() ;
    int type = central->GetType();
    if(type==0){p1 = pZero;    }    else{p1 = pOne;    }
 
    for(int i=0; i<4; i++){
     
        int neighbourDirection =  neighbour[i]->GetDirection();
        int Ntype = neighbour[i]->GetType();
        if(Ntype==1){    p2 = pOne;    }    else{    p2 = pZero;    }
    //    std::cout<<"1"<<std::endl;
        double    Energy[4][4]=
               {
                {-2*k*p1*p2/r/r/r - (p1+p2)*E, -p1*E, 2*k*p1*p2/r/r/r - (p1-p2)*E, -p1*E },    // ABOVE OR BELOW DIPOLES
                {-p2*E, 2*k*p1*p2/3/r/r/r, p2*E, -2*k*p1*p2/3/r/r/r    },
                {2*k*p1*p2/r/r/r + (p1-p2)*E, p1*E, -2*k*p1*p2/r/r/r + (p1+p2)*E, p1*E},
                {-p2*E, -2*k*p1*p2/3/r/r/r, p2*E, 2*k*p1*p2/3/r/r/r}
               };
        //       std::cout<<"2"<<std::endl;
       double    Energy2[4][4] =
               {
                { 2*k*p1*p2/3/r/r/r -(p1+p2)*E, -p1*E, -2*k*p1*p2/3/r/r/r -(p1-p2)*E, -p1*E },  // RIGHT OR LEFT ENERGIES
                { -p2*E, -2*k*p1*p2/r/r/r, p2*E, 2*k*p1*p2/r/r/r },
                { -2*k*p1*p2/3/r/r/r +(p1-p2)*E, p1*E, 2*k*p1*p2/3/r/r/r +(p1+p2)*E, p1*E },
                {-p2*E, 2*k*p1*p2/r/r/r, p2*E, -2*k*p1*p2/r/r/r }
               };
    //std::cout<<"3"<<std::endl;
    std::cout<<"i = "<< i << "\t" << type << "\t "<<Ntype <<"\t"<< neighbourDirection<<"\t"<<myDirection<<std::endl;
    std::cout<<"x = " << x << "y = "<< y <<std::endl;
            if(i==0 || i==2){
                energy += Energy[myDirection][neighbourDirection]; 
            }     
         
            else if(i==1 || i==3){          
                energy += Energy2[myDirection][neighbourDirection];         
            }     
    }//std::cout<<"4"<<std::endl;
    std::cout<<"4end clac local E"<<std::endl;
     return energy;
   
}
 
Lattice::~Lattice(){
 
    delete[] dipoleArray;
}

EDIT:: When it crashes it is always at an edge of the matrix. Usually the top row.
If my dipole is top row, it crashes when I find the 'above' neighbour.
When returning properties of this neighbour, "type" and "direction," defined in the dipole class to be 0,1 or 0,1,2,3 it returns large +ve or -ve values for the neighbour

PLEASE help as I can't see any problems :S

 

[Svein]

I am not sure about the value of x%xsize if x = -1. Try using ((x+xsize)%xsize) instead. And, of course, the same for y%ysize.

 

[rolotomassi]

I freaking love you! I just assumed negative modulo was fair game. Thanks a lot.

 

[newjerseyrunner]

Learn how to step through your program line by line using your debugger, you'll be able to see things like that easily next time.