Code:
/*Laplace equation 2D*/
/*Build the Mesh*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define max 100
/*====================================================================
Defining function VOID
======================================================================*/
/*====================================================================
Function Contour
======================================================================*/
void contour(float f[max][max],float x[max],float y[max],int jmax, int imax, float ui,
int i, int j, float il,float it)
{
for(j=1;j<=jmax;j++)
{
f[1][j]=ui*x[1];
f[imax][j]=ui*x[imax];
}
for(i=1;i<=imax;i++)
{
f[i][jmax]=ui*x[i];
if((i>=il) && (i<=it))
f[i][1]=f[i][2]-0.05*(y[2]-y[1])*(1.0-2*x[i]);
else
f[i][1]=f[i][2];
}
}// end contour ()
/*===================================================================
end VOID
=====================================================================*/
/*====================================================================
Function Res
======================================================================*/
void resi(float res[max][max],float x[max],float y[max],float f[max][max],int i, int j,
int imax, int jmax, float remax)
{
for (i=2;i<=imax-1;i++)
{
for(j=2;j<=jmax-1;j++)
res[i][j]=(2/(x[i+1]-x[i-1]))*((f[i+1][j]-f[i][j])/(x[i+1]-x[i])-(f[i][j]-f[i-1][j])/(x[i]-x[i-1]))
+(2/(y[j+1]-y[j-1]))*((f[i][j+1]-f[i][j])/(y[j+1]-y[j])-(f[i][j]-f[i][j-1])/(y[j]-y[j-1]));
}
remax=0;
for(i=2;i<=imax-1;i++)
{
for(j=2;j<=jmax-1;j++)
if (fabs(res[i][j])>=remax)
remax=abs(res[i][j]);
}
}
/*================================================================
End resi
==================================================================*/
int
main()
{
FILE *mesh;
float dx,x[max],y[max],f[max][max],u[max][max],v[max][max],c[max][max],df[max][max],res[max][max];
float ys,xs,il,it,ui,dy,remax;
int imax,jmax,nmax,n,i,j,s;
mesh=fopen("mesh.data","w");
il=11;
it=31;
imax=41;
jmax=12;
ys=1.24;
xs=1.24;
ui=1.0;
nmax=10000;
/*=========================================================
Generating the Mesh
===========================================================*/
dx=-1.0/(il-it);
for (i=il;i<=it;i++)
{
x[i]=(i-il)*dx;
}
for (i=it+1;i<=imax;i++)
{
x[i]=x[i-1]+(x[i-1]-x[i-2])*xs;
}
for(i=il-1;i>=1;i--)
{
x[i]=x[i+1]+(x[i+1]-x[i+2])*xs;
}
y[1]=-dx/2;
y[2]=dx/2;
for(j=3;j<=jmax;j++)
{
y[j]=y[j-1]+(y[j-1]-y[j-2])*ys;
}
/*===============================================================*/
/*===============================================================
Initial Conditions
=================================================================*/
for(i=1;i<=imax;i++)
{
for(j=1;j<=jmax;j++)
f[i][j]=ui*x[i];
}
/*============================================================*/
for(n=1;n<=nmax;n++)
{
contour( f, x, y,jmax,imax, ui,i, j, il,it);
resi( res,x,y,f,i,j,imax,jmax,remax);
for(i=2;i<=imax-1;i++)
{
for(j=2;j<=jmax-1;j++)
dx=(x[i+1]-x[i-1])/2;
dy=(y[j+1]-y[j-1])/2;
df[i][j]=-res[i][j]/(-2/(dx*dx)-2/(dy*dy));
}
for(i=2;i<=imax-1;i++)
{
for(j=2;j<=jmax-1;j++)
f[i][j]=f[i][j]+df[i][j];
}
}
/*=============================================================*/
/*===============================================================
Calculus of Velocities
=================================================================*/
for(i=2;i<=imax-1;i++)
{
for(j=2;j<=jmax-1;j++)
u[i][j]=(f[i+1][j]-f[i-1][j])/(x[i+1]-x[i-1]);
v[i][j]=(f[i][j+1]-f[i][j-1])/(y[j+1]-y[j-1]);
c[i][j]=1.0-u[i][j]*u[i][j]+v[i][j]*v[i][j];
}
/*===============================================================
Print the Results
==================================================================*/
for(j=1;j<=jmax;j++)
{
for(i=1;i<=imax;i++)
/*fprintf(mesh,"%1.15f %1.15f\n",u[i][j],v[i][j]);*/
printf("%lf\n",u[i][j]);
printf("\n");
}
fclose(mesh);
exit(0);
}