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