CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Source code archive - educational

Source code archive - educational

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
m (Reverted edit of Engware, changed back to last version by Jola)
m (One dimensional stedy state conduction with Heat generation code using C)
Line 1: Line 1:
-
Codes and scripts used in education. This is the place to put all your model-problem codes in Matlab etc.
+
#include <stdio.h>
 +
#include <stdlib.h>
-
{{stub}}
+
#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));
 +
}

Revision as of 07:29, 28 March 2006

  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)); }

My wiki