Hey guys,
In the following code to find the determinant, the program crashes when it attempts to delete the matrix and gives an error about heap corruption.
I'm having trouble determining why I'm getting this error. The code has worked fine until I realised my derivative definitions were wrong and added the brackets to N1E1, N2E2, N3E3 and N4E4 at the top.
I've double checked, by returning them to the previous values and the code worked again.
Could someone possibly shed some light on why this is happening?
Many Thanks,
Dillinger
//headers etc..
#define N1E1 4*(E1-1)
#define N5E1 4*E2
#define N7E1 4*E3
#define N8E1 4*E4
#define N2E2 4*(E2-1)
#define N5E2 4*E1
#define N6E2 4*E3
#define N9E2 4*E4
#define N3E3 4*(E3-1)
#define N6E3 4*E2
#define N7E3 4*E1
#define N10E3 4*E4
#define N4E4 4*(E4-1)
#define N8E4 4*E1
#define N9E4 4*E2
#define N10E4 4*E3
#define C elMod/((1 + poisRat)*(1 - 2*poisRat))
#define G elMod/(2*(1 + poisRat))
int N2E1, N3E1, N4E1, N6E1, N9E1, N10E1, N1E2, N3E2, N4E2, N7E2, N8E2, N10E2, N1E3, N2E3, N4E3, N5E3, N8E3, N9E3, N1E4, N2E4, N3E4, N5E4, N6E4, N7E4 = 0;
Mesh Mesh1;
float force[30] = {0};
double* gauss(int N, // number of unknowns
double A [30] [31], // coefficients and constants
double result[30],
bool& err);
double matrix_det(double **in_matrix, int n);
double matrix_cofact(double **in_matrix, int row, int col, int order);
double** matrix_multi(double **matrix_1, int n1, int m1, double **matrix_2, int n2, int m2);
using namespace std;
int main(int argc, char** argv)
{
double A[30][31];
bool err;
double **E;
int elMod = 60;
float poisRat = 0.25;
vector<float> row;
int i, currentElement;
float Ke[30][30] = {0};
double **J, detJ;
float K[12][12] = {0};
double **B;
float Nd[10][4] = {0};
double E1, E2, E3, E4;
E = new double *[6];
for (i = 0; i < 6; i++)
{
E[i] = new double [6];
for (int j = 0; j < 6; j++)
E[i][j] = 0;
}
//Build stress-straincompliance matrix
E[0][0] = E[1][1] = E[2][2] = (1 - poisRat)*C;
E[0][1] = E[0][2] = E[1][2] = E[1][0] = E[2][0] = E[2][1] = poisRat*C;
E[3][3] = E[4][4] = E[5][5] = G;
//print to check
cout << "E = \n";
for(i = 0; i < 6; i++)
{
for(int j = 0; j < 6; j++)
cout << E[i][j] << " ";
cout << endl;
}
//declare temp element
Element Element1;
//Nodes
Mesh1.createNode(0.0f, 0.0f, 0.0f);
Mesh1.createNode(2.0f, 0.0f, 0.0f);
Mesh1.createNode(2.0f, 0.0f, 2.0f);
Mesh1.createNode(2.0f, 2.0f, 0.0f);
Mesh1.createNode(1.0f, 0.0f, 0.0f);
Mesh1.createNode(2.0f, 0.0f, 1.0f);
Mesh1.createNode(1.0f, 0.0f, 1.0f);
Mesh1.createNode(1.0f, 1.0f, 0.0f);
Mesh1.createNode(2.0f, 1.0f, 0.0f);
Mesh1.createNode(2.0f, 1.0f, 1.0f);
//point to mesh elements
Element1.elementNodes[0] = &Mesh1.meshNodes[0];
Element1.elementNodes[1] = &Mesh1.meshNodes[1];
Element1.elementNodes[2] = &Mesh1.meshNodes[2];
Element1.elementNodes[3] = &Mesh1.meshNodes[3];
Element1.elementNodes[4] = &Mesh1.meshNodes[4];
Element1.elementNodes[5] = &Mesh1.meshNodes[5];
Element1.elementNodes[6] = &Mesh1.meshNodes[6];
Element1.elementNodes[7] = &Mesh1.meshNodes[7];
Element1.elementNodes[8] = &Mesh1.meshNodes[8];
Element1.elementNodes[9] = &Mesh1.meshNodes[9];
//mesh node integer references
Element1.elementNodeIndex[0] = 0;
Element1.elementNodeIndex[1] = 1;
Element1.elementNodeIndex[2] = 2;
Element1.elementNodeIndex[3] = 3;
Element1.elementNodeIndex[4] = 4;
Element1.elementNodeIndex[5] = 5;
Element1.elementNodeIndex[6] = 6;
Element1.elementNodeIndex[7] = 7;
Element1.elementNodeIndex[8] = 8;
Element1.elementNodeIndex[9] = 9;
// store element in mesh
Mesh1.meshElements.push_back(Element1);
//declare memory for 4x4 matrix
J = new double *[4];
for(int j = 0; j < 4; j++)
J[j] = new double [4];
currentElement = 0;
for(i = 0; i < 4; i++)
{
//determine local coords for gauss point i.
E1 = E2 = E3 = E4 = 0.13819660;
switch(i)
{
case 0: E1 = 0.58541020;
break;
case 1: E2 = 0.58541020;
break;
case 2: E3 = 0.58541020;
break;
case 3: E4 = 0.58541020;
break;
}
//first row of Jacobian is 1
J[0][0] = J[0][1] = J[0][2] = J[0][3] = 1;
for(int j = 1; j < 4; j++)
{
J[j][0] = Mesh1.meshElements[currentElement].elementNodes[0]->Coords[j-1]*N1E1 +
Mesh1.meshElements[currentElement].elementNodes[4]->Coords[j-1]*N5E1 +
Mesh1.meshElements[currentElement].elementNodes[6]->Coords[j-1]*N7E1 +
Mesh1.meshElements[currentElement].elementNodes[7]->Coords[j-1]*N8E1;
J[j][1] = Mesh1.meshElements[currentElement].elementNodes[1]->Coords[j-1]*N2E2 +
Mesh1.meshElements[currentElement].elementNodes[5]->Coords[j-1]*N6E2 +
Mesh1.meshElements[currentElement].elementNodes[4]->Coords[j-1]*N5E2 +
Mesh1.meshElements[currentElement].elementNodes[8]->Coords[j-1]*N9E2;
J[j][2] = Mesh1.meshElements[currentElement].elementNodes[2]->Coords[j-1]*N2E3 +
Mesh1.meshElements[currentElement].elementNodes[6]->Coords[j-1]*N7E3 +
Mesh1.meshElements[currentElement].elementNodes[5]->Coords[j-1]*N6E3 +
Mesh1.meshElements[currentElement].elementNodes[9]->Coords[j-1]*N10E3;
J[j][3] = Mesh1.meshElements[currentElement].elementNodes[3]->Coords[j-1]*N4E4 +
Mesh1.meshElements[currentElement].elementNodes[7]->Coords[j-1]*N8E4 +
Mesh1.meshElements[currentElement].elementNodes[8]->Coords[j-1]*N9E4 +
Mesh1.meshElements[currentElement].elementNodes[9]->Coords[j-1]*N10E4;
}
//store dervatives for easier access later
Nd[0][0] = N1E1;
Nd[4][0] = N5E1;
Nd[6][0] = N7E1;
Nd[7][0] = N8E1;
Nd[1][1] = N2E2;
Nd[4][1] = N5E2;
Nd[5][1] = N6E2;
Nd[8][1] = N9E2;
Nd[2][2] = N3E3;
Nd[5][2] = N6E3;
Nd[6][2] = N7E3;
Nd[9][2] = N10E3;
Nd[3][3] = N4E4;
Nd[7][3] = N8E4;
Nd[8][3] = N9E4;
Nd[9][3] = N10E4;
//print Jacobian
for(int j = 0; j < 4; j++)
{
for(int k = 0; k < 4; k++)
{
cout << J[j][k] << "\t";
}
cout << endl;
}
//find determinant of Jacobian
detJ = matrix_det(J, 4);
/* REST OF MAIN SNIPPED DUE TO IRRELEVANCE */
//borrowed determinant code
double matrix_det(double **in_matrix, int n)
{
int i, j, k;
double **matrix;
double det = 1;
matrix = new double *[n];
for ( i = 0; i < n; i++ )
matrix[i] = new double[n];
for ( i = 0; i < n; i++ ) {
for ( j = 0; j < n; j++ )
matrix[i][j] = in_matrix[i][j];
}
for ( k = 0; k < n; k++ ) {
if ( matrix[k][k] == 0 ) {
bool ok = false;
for ( j = k; j < n; j++ ) {
if ( matrix[j][k] != 0 )
ok = true;
}
if (!ok)
return 0;
for ( i = k; i < n; i++ )
swap ( matrix[i][j], matrix[i][k] );
det = -det;
}
det *= matrix[k][k];
if ( k + 1 < n ) {
for ( i = k + 1; i < n; i++ ) {
for ( j = k + 1; j < n; j++ )
matrix[i][j] = matrix[i][j] - matrix[i][k] *
matrix[k][j] / matrix[k][k];
}
}
}
for ( i = 0; i < n; i++ )
delete [] matrix[i];
delete [] matrix;
return det;
}