0

Hi
I wrote a code for Galerkin Finite Element .The program runs with no error but have some problem breaking the loop also the results are all wrong .Can anyone please help me? This is my first time writting aprogram and I am so confused.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>



#define Limit  20

int main()
{
    double eps = 10E-6;              /* Tolerance                   */
    double Sep  = 1.0;               /*                   */
    int K    ;                    /* Counter for iterations      */
    int Max  = 20;                   /* Maximum number of iterat.   */
    int Cond = 1;                    /*         */
    int j, L, i,T,r,P;                     /* Loop counters               */
    double X[Limit];
    double CON,CONX;                      /* B in  AX = B ,  INPUT       */
    double DX;
    int N;                           /* No.of Nodes      */   
    int NE;
    double A[Limit][Limit];
                                /* No. Of elements      */
    double SUM, M;                                                       
    int  Pivot;
    int L1;
    double satisfied = 0.0;
    int Row[Limit];                      /* Variable in dominance check */
    double C[Limit];  
    double DCNEW[Limit];
    double Cnew[Limit];
    double Error;
    double SumF[Limit],SumJ[Limit][Limit];
    double PHI[2];
    double PHIX[2];
    double GP[3]={0.1127,0.5,0.8877};
    double w[3]={0.2778,0.4444,0.2778};
    #define Gamma 0.1


    {
      printf("Please enter number of elements [Not more than %d]\n",Limit);
      scanf_s("%d", &NE);
      N =   NE+1;

     /* Define mesh size*/
      for (i = 1; i<= N; i++)
          X[i] = (i-1)/NE;

      /* Initial guess for N.R */
      for(i=1;i<N;i++) C[i]=.5;//inital guess vector
      C[N] = 1.0;


      while (satisfied == 0)
      {
      for(i=1;i<N;i++)
      {
          SumF[i] = 0.0;
      for(j=1;j<N;j++)
          SumJ[i][j] = 0.0;
      }
      DX = X[N] - X[1]/ NE;

      CON = 0.0;
      CONX = 0.0;

      for(j=1;j<=3;j++)
      {
      PHI[1] = 1-GP[j];
      PHI[2] = GP[j];
      PHIX[1] = -1/DX;
      PHIX[2] = 1/DX;
      for(L=1;L<=2;L++)
      {
          L1 = j+L-1;
      CON = CON + PHI[L]*C[L1];
      CONX = CONX + PHIX[L]*C[L1];
      SumF[L1] = SumF[L1] - DX* w[j]*(-CONX*PHIX[L] - Gamma* pow(CON,3)*PHI[L]);
      SumJ [L1][L] = SumJ[L1][L1] + DX*w[j] *(-PHIX[L]*PHIX[L1] - 3.0*Gamma*pow(CON,2)*PHI[L]*PHI[L1]);
      }

      }
      /* Apply Boundry Condition */

    for(i=1;i<N;i++) SumJ[N][i]=0.0;//inital guess vector
      SumJ[N][N] = 1.0;
      SumF[N] = 0.0;

      for (i = 1; i<=N; i++)
      {
          for (j = 1; j<=N; j++)
              A[i][j] = SumJ[i][j];
          A[i][N+1] = SumF[j];
      }

      //Gaussian Elimination Method
      for (r = 1; r<= N; r++) Row[r-1] = r - 1;

    /* Start upper-triangularization */

    for (P = 1; P <= N - 1; P++)
    {

     /* Pivoting */ 

      for (K = P + 1; K <= N; K++)
      {
        if ( fabs(A[Row[K-1]][P-1]) > fabs(A[Row[P-1]][P-1]) )
        {

           Pivot      = Row[P-1];
           Row[P-1] = Row[K-1];
           Row[K-1] = Pivot;
        }
      } 

      /* Process: a21/a11, a'32/a'22, ... */

     for (K = P + 1; K <= N; K++)
     { 
       M = A[Row[K-1]][P-1] / A[Row[P-1]][P-1];


       for (T = P + 1; T <= N + 1; T++)
       {
          A[Row[K-1]][T-1] -= M * A[Row[P-1]][T-1];
       }
     } 

  } 

    if( A[Row[N-1]][N-1] == 0)
    {
      printf("The matrix is SINGULAR as the determinant is zero !\n");
      printf(" Gaussian Elimination Algorithm failed --> exit\n");
      exit(1);
    }

    /* Back substitution */

    DCNEW[N-1] = A[Row[N-1]][N] / A[Row[N-1]][N-1];

    for (K = N - 1; K >= 1; K--)
    {
      SUM = 0;
      for (T = K + 1; T <= N; T++)
      {
           SUM += A[Row[K-1]][T-1] * DCNEW[T-1];   
      }

      DCNEW[K-1] = ( A[Row[K-1]][N] - SUM) / A[Row[K-1]][K-1];

    }  /* End of back substitution */


    for (i = 1; i <= N; i++)
        {

            Cnew[i] = C[i] + DCNEW[i];

    /* Calculating Error */
    Error = abs ((Cnew[i]-C[i]/Cnew[i])*100);
    if (Error < eps)  
    {
        satisfied=1;

        break;

    }//if end

    for (i = 1; i <= N; i++)
    printf("C[%d] = %1f\n",i,Cnew[i]);


 } //while end

 }
}
}

Edited by mike_2000_17: Fixed formatting

2
Contributors
1
Reply
2
Views
6 Years
Discussion Span
Last Post by frogboy77
0

1. use code tags
2.

This is my first time writting aprogram

what happened to the hello world one?

Either this is your first program(doubtful) or it's not your's. If it's the latter then i doubt you'll get much help here.

Tip: try being honest and specific (works a treat).

This topic has been dead for over six months. Start a new discussion instead.
Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.