Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

C++ matrix boundary condition problems

  1. Apr 30, 2016 #1
    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.

    Code (Text):

    #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
     
  2. jcsd
  3. Apr 30, 2016 #2

    Svein

    User Avatar
    Science Advisor

    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.
     
  4. Apr 30, 2016 #3
    I freaking love you! I just assumed negative modulo was fair game. Thanks a lot.
     
  5. May 2, 2016 #4
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: C++ matrix boundary condition problems
  1. C++ Matrix Solver (Replies: 0)

Loading...