C++ matrix boundary condition problems

[rolotomassi]

Main Question or Discussion Point

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 im 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 cant 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.