Code: One dimensional steady state conduction with heat generation

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