Source code archive - educational

From CFD-Wiki

1. include <stdio.h>
2. include <stdlib.h>

1. define MAXROWS 10000

double **GetMatrix(double, double, double, double, double, int); void squareoutput(double **, int); void nonsquareoutput(double **, int, int); double **triangularize(double *[MAXROWS], int, int); double *returnsolvector(double **c, int nrows); void writeoutputvector(double *p, int nrows);

int main(int argc, char **argv) {

/* Variables chosen: GAMMA: Thermal Conductivity L  : Length of the domain Dx  : Node spacing

• /

double GAMMA, L, Dx, Ta , Tb, Area; int N, row; double **temp; double **aug; double *solvector;

printf("\nThis program works for only simple 1-D conduction"); printf("Is thermal Conductivity constant?");

printf("\nEnter the Value for Heat conductivity"); scanf("%lf", &GAMMA); printf("\nEnter the Length of the computational domain"); scanf("%lf", &L); printf("\nEnter the cross-sectional area of the computational domain"); scanf("%lf", &Area); printf("\nEnter the temperature at the ends of the domain"); scanf("%lf %lf", &Ta, &Tb); printf("\nEnter the number of Nodes u want to be chosen for analysis"); scanf("%d", &N);

temp=(double **) malloc(N*sizeof(double *)); for(row=0;row<N;row++) temp[row]=(double *) malloc(N*sizeof(double));

aug= (double **) malloc ((N+1)*sizeof(double)); for(row=0;row<N;row++) aug[row]=(double *) malloc((N+1)*sizeof(double));

solvector= (double *) malloc (N*sizeof(double));

Dx=L/(N);

printf("\nThe Node spacing is: %lf", Dx);

aug=GetMatrix(GAMMA, Dx, Ta, Tb, Area, N); aug=triangularize(aug, N, N+1); nonsquareoutput(aug, N, N+1); printf("\n\n"); solvector=returnsolvector(aug,N); writeoutputvector(solvector, N); printf("\n\n\n"); return 0; }

double **GetMatrix(double GAMMA, double Dx, double Ta, double Tb, double Area, int N) { int row, col; double DxW, DxE; double Aw, Ae, Ap;

/* Memory Allocation for the matrix */

double **matrix; double *solvector; double **aug;

solvector= (double *) malloc (N*sizeof(double));

matrix=(double **) malloc(N*sizeof(double)); for(row=0;row<N;row++) matrix[row]=(double *) malloc(N*sizeof(double));

aug= (double **) malloc ((N+1)*sizeof(double)); for(row=0;row<N;row++) aug[row]=(double *) malloc((N+1)*sizeof(double));

/* Uniform Mesh has been chosen */

DxW=DxE=Dx;

Aw=(GAMMA/DxW)*Area; Ae=(GAMMA/DxE)*Area;

/* Generation of the Matrix */

row=0; for(col=0; col<N; col++) { if(col==row) *(*(matrix+row)+col)= (Ae+2*Aw); else if(col==row+1) *(*(matrix+row)+col)= -Ae ; else *(*(matrix+row)+col)=0; }

for (row=1; row<N-1; row++) { for(col=0; col<N; col++) { if(col==row+1) *(*(matrix+row)+col)= -Aw; else if (col==row-1) *(*(matrix+row)+col)= -Ae; else if (col==row) *(*(matrix+row)+col)= (Aw+Ae); else *(*(matrix+row)+col)=0; } }

row=N-1; for(col=0; col<N; col++) { if(col==row) *(*(matrix+row)+col)= (2*Ae+Aw); else if(col==row-1) *(*(matrix+row)+col)=-Aw; else *(*(matrix+row)+col)=0; }

/* Getting the solution vector */

*(solvector+0)= 2*Aw*Ta;

for (row=1; row<=N-2; row++) *(solvector+row)=0;

*(solvector+N-1)= 2*Ae*Tb;

/* Getting the augmented matrix */

for (row=0; row<N; row++) { for(col=0; col<N; col++) { *(*(aug+row)+col) = *(*(matrix+row)+col); } }

for (row=0; row<N; row++) *(*(aug+row)+N)= *(solvector+row);

return(aug); }

/* Displaying the matrix */

void squareoutput(double **c, int nrows) { int row, col; for(row=0; row<nrows; row++) { printf("\n");printf("\n"); for (col=0; col<nrows; col++) { printf("%8.2lf", *(c[row]+col)); } } return; }

void nonsquareoutput(double **c, int nrows, int ncols)

{ int row, col; for(row=0; row<nrows; row++) { printf("\n");printf("\n"); for (col=0; col<ncols; col++) { printf("%12.2lf", *(c[row]+col)); } } return; }

double **triangularize(double *a[MAXROWS], int nrows, int ncols) { int row, col, k, temp=0; double coef; double **tri; tri=(double **)malloc(ncols*sizeof(double *)); for(row=0;row<nrows;row++) tri[row]=(double*)malloc(ncols*sizeof(double));

{ int row, col; for (row=0; row<nrows; row++) for(col=0; col<ncols; col++) *(tri[row]+col)=*(a[row]+col); }

for(temp=0;temp<nrows-1; temp++) { if ((tri[temp]+temp)==0)continue; else { for(row=temp; row<nrows; row++) { for(k=row+1; k< nrows; k++) { if( *(tri[k]+temp)==0) continue; else coef= (*(tri[k]+temp))/(*(tri[temp]+temp)); for(col=0; col<ncols; col++) *(tri[k]+col)=*(tri[k]+col)-(coef* *(tri[row]+col)); } } } }

for(temp=nrows-1;temp>0; temp--) { if ((tri[temp]+temp)==0)continue; else { for(row=temp; row>0; row--) { for(k=row-1; k>=0; k--) { if( *(tri[k]+temp)==0) continue; else coef= (*(tri[k]+temp))/(*(tri[temp]+temp)); for(col=0; col<ncols; col++) *(tri[k]+col)=*(tri[k]+col)-(coef* *(tri[row]+col)); } } } }

return(tri); }

double *returnsolvector(double **c, int nrows) { int temp; double *solvector; solvector= (double *) malloc (nrows*sizeof(double));

for(temp=0; temp<nrows; temp++) { *(solvector+temp)= *(*(c+temp)+nrows)/ *(*(c+temp)+temp); } return(solvector); }

void writeoutputvector(double *p, int nrows) { int temp; for(temp=0; temp<nrows; temp++) printf("\nTemperature at Node %d = %lf", temp+1, *(p+temp)); }