Member Avatar for dillinger88

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;
}
  • 3 Contributors
  • 2 Replies
  • 174 Views
  • 8 Hours Discussion Span
  • Latest Post Latest Post by Salem

Recommended Answers

All 2 Replies

Member Avatar for Narue

>I'm having trouble determining why I'm getting this error.
Obviously you're writing to an out of bounds index. Step through your code with a test set that reproduces the error and verify all of your indices.

Member Avatar for Salem

What about the declarations of mesh and element?
And the associated functions?

Plus, remove all the unused variables from main()

Edit:
Which OS/Compiler are you using?

Edited by Salem because: n/a
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.