From CFD-Wiki
#include <stdio.h>
#include <stdlib.h>
#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));
}