Dynamic Allocation of 2D Array

[AliMoradi]

I wrote this entire code before I found out that 'new' doesn't work with multidimensional arrays, and now I'm trying to salvage the code without having to redo all my loops and if statements. My code basically has 2D arrays with one fixed dimension and another dimension that will have to be expanded as the algorithm progresses. I know that I can just do U [(M+1)*MAX)] instead of U [M+1][MAX], but that would require me to change my for-loops. Is there any way to get this code to work without having to do that?

If it helps, I'm using CodeBlocks as my compiler and I'm getting the following errors:

error: 'MAX' cannot appear in a constant-expression.
error: invalid types 'double[int]' for array subscript.

I appreciate any help!

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <cmath>
#include <algorithm>

#define pi 3.14159265

using namespace std;

const int M = 40; //X-Grid Intervals
double L0 = -1;
double LF = 1;
double tmax = 0.15; //Final Time
int MAX = M; //Maximum t-vector Size

//Dynamic Memory Allocation
double * U;
double * E;
double * t = new double [MAX];

int SGN(double aP)
{
    if (aP > 0)
    return 1;
    else if (aP == 0)
    return 0;
    else if (aP < 0)
    return -1;
}


int main()
{
    //Create X Vector and Initial Conditions
    double x [M+1];
    U = new double [M+1][MAX];
    E = new double [M+1][MAX];
    x[0] = L0;
    x[M] = LF;
    for (int i = 1; i < M; i++)
    {
        x[i] = -1 + (i - 0.5)*2/M;
        U[i][0] = 0.25 + 0.5*sin(pi*x[i]);
        E[i][0] = 0.5*U[i][0]*U[i][0];
    }
    U[0][0] = U[M-1][0];
    E[0][0] = 0.5*U[0][0]*U[0][0];
    U[M][0] = U[1][0];
    E[M][0] = 0.5*U[M][0]*U[M][0];
    double dx = x[2] - x[1];

    //Allocate Variables and Vectors
    double DUP;
    double DUM;
    double aP [M+1];
    double aM [M+1];
    double dtt [M+1];
    double dt;
    int G;
    int iv;
    double xx;

    t[0] = 0;

    int s = 0;
    int ERR = 0;
    while (ERR == 0)
    {
        //Calculate Alpha and Time Step
        for (int i = 1; i < M; i++)
        {
            DUP = U[i+1][s] - U[i][s];
            DUM = U[i][s] - U[i-1][s];
            if (DUP == 0)
            {
                aP[i] = U[i][s];
            }
            else
            {
                aP[i] = (E[i+1][s] - E[i][s])/(U[i+1][s] - U[i][s]);
            }
            if (DUM == 0)
            {
                aM[i] = U[i][s];
            }
            else
            {
                aM[i] = (E[i][s] - E[i-1][s])/(U[i][s] - U[i-1][s]);
            }

            dtt[i] = dx/abs(aP[i]);

        }

        //Find Minimum Time Step and Round To Low Value
        dt = dtt[1];
        for (int i = 2; i < M; i++)
        {
            if (dtt[i] < dt)
            {
                dt = dtt[i];
            }
        }

        //Even Out Time Step
        G = 0;
        iv = 1;
        while (G == 0)
        {
            xx = tmax/iv;
            if (xx < y)
            {
                dt = xx;
                G = 1;
            }
            else
            {
                iv = iv + 1;
            }
        }

        //Calculate Solution For New Time Step
        for (int i = 1; i < M; i++)
        {
            U[i][s+1] = U[i][s] - 0.5*dt/dx*(1 - SGN(aP[i]))*(E[i+1][s] - E[i][s]) - 0.5*dt/dx*(1 + SGN(aM[i]))*(E[i][s] - E[i-1][s]);
            E[i][s+1] = 0.5*U[i][s+1]*U[i][s+1];
        }
        //Period Boundary Condition
        U[0][s+1] = U[M-1][s+1];
        U[M][s+1] = U[1][s+1];
        E[0][s+1] = 0.5*U[0][s+1]*U[0][s+1];
        E[M][s+1] = 0.5*U[M][s+1]*U[M][s+1];

        //Calculate New Time
        t[s+1] = t[s] + dt;
        s = s + 1;

        //Expand Arrays
        if (s >= MAX)
        {
            MAX = MAX*2;
            double * temp = new double [M+1][MAX];
            double * temp2 = new double [M+1][MAX];
            double * temp3 = new double [MAX];
            for (int j = 0; j < s; j++)
            {
                for (int i = 0; i < M+1; i++)
                {
                    temp[i][j] = U[i][j];
                    temp2[i][j] = E[i][j];
                }

                temp3[j] = t[j];
            }

            delete [] U;
            delete [] E;
            delete [] t;
            U = temp;
            E = temp2;
            t = temp3;

        }

        //Check To See If T-max Has Been Reached
        if (t[s] >= tmax)
        {
            ERR = 1;
        }

    }

    //Output Data Into Text File
    ofstream outputdata("timevector.txt");
    for (int i = 0; i<MAX; i++)
    {
        outputdata << t[i] << endl;
    }
    outputdata.close();

    ofstream outputdata3("xvector.txt");
    for (int j = 0; j<M+1; j++)
    {
        outputdata3 << x[j] << endl;
    }
    outputdata3.close();

    ofstream outputdata2("U.txt");
    for (int i = 0; i<MAX; i++)
    {
        for (int j = 0; j<M+1; j++)
        {
            outputdata2 << U[j][i] << endl;
        }
    }
    outputdata2.close();


    return 0;


}

[Cubbi]

You could either use a proper Matrix class (from a third-party library, such as boost or Eigen), or write your own matrix, or just use vectors of vectors.

Incidentally, all those 1D arrays should have been vectors as well:

int main()
{
    //Create X Vector and Initial Conditions
    vector<double> x(M+1);
    vector<vector<double>> U(M+1, vector<double>(MAX));
    vector<vector<double>> E(M+1, vector<double>(MAX));
...

or

int main()
{
    //Create X Vector and Initial Conditions
    vector<double> x(M+1);
    boost::numeric::ublas::matrix<double> U(M+1, MAX);
    boost::numeric::ublas::matrix<double> E(M+1, MAX);

[vlad from moscow]

These statements are invalid

    U = new double [M+1][MAX];
    E = new double [M+1][MAX];

because MAX (do not name variables with all uppercase letters!) is not a const expression.

Either define MAX as a constant, for example,

const int MAX = M;

and then define U and E as

double ( * U )[MAX];
double ( * E )[MAX];

or you may define U and E as

double ** U;
double ** E;

and then allocate memory for them the following way

U = new double * [M+1];

for ( int i = 0; i < M + 1; i++ ) *U = new double [MAX];

[AliMoradi]

(do not name variables with all uppercase letters!)

Why is that? I didn't know it created any issues.

or you may define U and E as

double ** U;
double ** E;


and then allocate memory for them the following way


U = new double * [M+1];

for ( int i = 0; i < M + 1; i++ ) *U = new double [MAX];

When I did this, I run into a lot of other errors. I'm not sure if they're resulting from the changes I made, or there are different errors that I have. This is the list of errors I'm getting:

|24| error: reference to 'max' is ambiguous|

|19| error: candidates are: int max|

\include\c++\bits\stl_algobase.h |253| error: template<class _Tp, class _Compare> const _Tp& std::max(const _Tp&, const _Tp&, _Compare)|

\include\c++\bits\stl_algobase.h|209|error: template<class _Tp> const _Tp& std::max(const _Tp&, const _Tp&)|

I'm not sure why 'max' reference is ambiguous. The t-vector is only 1D and it shouldn't have been a problem dynamically allocating that with 'new'. And for whatever reason, these errors keep on repeating themselves.

Here is the changed version of the code:

//1st-Order TVD For Inviscid Burgers Equation
//Ali Moradi

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <cmath>
#include <algorithm>

#define pi 3.14159265

using namespace std;

const int M = 40; //X-Grid Intervals
double L0 = -1;
double LF = 1;
double tmax = 0.15; //Final Time
int max = 40; //Maximum t-vector Size

//Dynamic Memory Allocation
double ** U;
double ** E;
double * t = new double [max];

int SGN(double aP)
{
    if (aP > 0)
    return 1;
    else if (aP == 0)
    return 0;
    else if (aP < 0)
    return -1;
}


int main()
{
    //Create X Vector and Initial Conditions
    double x [M+1];
    U = new double * [M+1];
    E = new double * [M+1];
    for ( int i = 0; i < M + 1; i++ ) * U = new double [max];
    for ( int i = 0; i < M + 1; i++ ) * E = new double [max];
    x[0] = L0;
    x[M] = LF;
    for (int i = 1; i < M; i++)
    {
        x[i] = -1 + (i - 0.5)*2/M;
        U[i][0] = 0.25 + 0.5*sin(pi*x[i]);
        E[i][0] = 0.5*U[i][0]*U[i][0];
    }
    U[0][0] = U[M-1][0];
    E[0][0] = 0.5*U[0][0]*U[0][0];
    U[M][0] = U[1][0];
    E[M][0] = 0.5*U[M][0]*U[M][0];
    double dx = x[2] - x[1];

    //Allocate Variables and Vectors
    double DUP;
    double DUM;
    double aP [M+1];
    double aM [M+1];
    double dtt [M+1];
    double dt;
    int G;
    int iv;
    double xx;

    t[0] = 0;

    int s = 0;
    int ERR = 0;
    while (ERR == 0)
    {
        //Calculate Alpha and Time Step
        for (int i = 1; i < M; i++)
        {
            DUP = U[i+1][s] - U[i][s];
            DUM = U[i][s] - U[i-1][s];
            if (DUP == 0)
            {
                aP[i] = U[i][s];
            }
            else
            {
                aP[i] = (E[i+1][s] - E[i][s])/(U[i+1][s] - U[i][s]);
            }
            if (DUM == 0)
            {
                aM[i] = U[i][s];
            }
            else
            {
                aM[i] = (E[i][s] - E[i-1][s])/(U[i][s] - U[i-1][s]);
            }

            dtt[i] = dx/abs(aP[i]);

        }

        //Find Minimum Time Step and Round To Low Value
        dt = dtt[1];
        for (int i = 2; i < M; i++)
        {
            if (dtt[i] < dt)
            {
                dt = dtt[i];
            }
        }

        //Even Out Time Step
        G = 0;
        iv = 1;
        while (G == 0)
        {
            xx = tmax/iv;
            if (xx < y)
            {
                dt = xx;
                G = 1;
            }
            else
            {
                iv = iv + 1;
            }
        }

        //Calculate Solution For New Time Step
        for (int i = 1; i < M; i++)
        {
            U[i][s+1] = U[i][s] - 0.5*dt/dx*(1 - SGN(aP[i]))*(E[i+1][s] - E[i][s]) - 0.5*dt/dx*(1 + SGN(aM[i]))*(E[i][s] - E[i-1][s]);
            E[i][s+1] = 0.5*U[i][s+1]*U[i][s+1];
        }
        //Period Boundary Condition
        U[0][s+1] = U[M-1][s+1];
        U[M][s+1] = U[1][s+1];
        E[0][s+1] = 0.5*U[0][s+1]*U[0][s+1];
        E[M][s+1] = 0.5*U[M][s+1]*U[M][s+1];

        //Calculate New Time
        t[s+1] = t[s] + dt;
        s = s + 1;

        //Expand Arrays
        if (s >= max)
        {
            max = max*2;
            double ** temp;
            double ** temp2;
            double * temp3 = new double [max];
            temp = new double * [M+1];
            temp2 = new double * [M+1];
            for ( int i = 0; i < M + 1; i++ ) * temp = new double [max];
            for ( int i = 0; i < M + 1; i++ ) * temp2 = new double [max];
            for (int j = 0; j < s; j++)
            {
                for (int i = 0; i < M+1; i++)
                {
                    temp[i][j] = U[i][j];
                    temp2[i][j] = E[i][j];
                }

                temp3[j] = t[j];
            }

            delete [] U;
            delete [] E;
            delete [] t;
            U = temp;
            E = temp2;
            t = temp3;

        }

        //Check To See If T-max Has Been Reached
        if (t[s] >= tmax)
        {
            ERR = 1;
        }

    }

    //Output Data Into Text File
    ofstream outputdata("timevector.txt");
    for (int i = 0; i<max; i++)
    {
        outputdata << t[i] << endl;
    }
    outputdata.close();

    ofstream outputdata3("xvector.txt");
    for (int j = 0; j<M+1; j++)
    {
        outputdata3 << x[j] << endl;
    }
    outputdata3.close();

    ofstream outputdata2("U.txt");
    for (int i = 0; i<max; i++)
    {
        for (int j = 0; j<M+1; j++)
        {
            outputdata2 << U[j][i] << endl;
        }
    }
    outputdata2.close();


    return 0;


}

Thanks again for all the help!

[vlad from moscow]

The errors are consequence of that the name max conflicts with the name std::max. If you define a constant you may use name that consists of all uppercase letters.

For whom did I write const int MAX = M;? Can you simplly rewrite it?!
Or use other name than max, for example, max_size or max_dim and so on.