April 3, 2019 1 / 34
7.1 ( ) 2 Poisson 2 / 34
7.2 femfp.c [1] main( ) input( ) assem( ) ecm( ) f( ) solve( ) gs { solve( ) output( ) 3 / 34
7.3 fopen() #include <stdio.h> FILE *fopen(char *fname, char *mode); fname mode ( ) NULL mode mode r rb r+ rb+ w wb w+ wb+ a ab a ab+ 4 / 34
fclose() #include <stdio.h> int fclose(file *fp); fopen fp 0 EOF fprintf() include <stdio.h> int fprintf(file *fp, char *fmt,...); fp fmt 5 / 34
fputs() #include <stdio.h> int fputs(char *string, FILE *fp); fp string 0 EOF fgets() #include <stdio.h> char *fgets(char *string, int n, FILE *fp); fp string 0 EOF 6 / 34
7.4 main() /* The Finite Element Method for the Poisson Equation */ /* Full Matrix Version */ /* The Gauss Elimination Method */ #include <stdio.h> #include <stdlib.h> #include <math.h> #define NDIM1 200 #define NDIM2 400 7 / 34
void input(int *nnode, int *nelmt, int *nbc, int ielmt[][3], double x[], double y[], int ibc[]); void assem(int nnode, int nelmt, int nbc, int ielmt[][3], double x[], double y[], int ibc[], double am[][ndim1], double fm[]); void ecm(int ielmt[][3], double x[], double y[], int ie, double ae[][3], double fe[]); void solve(int m, double a[][ndim1], double f[]); void gs_solve(int m, double a[][ndim1], double f[]); void output(double fm[], int nnode, double x[], double y[]); double f(double x, double y); 8 / 34
/* Main Program */ int main() { int nnode, nelmt, nbc, ielmt[ndim2][3], ibc[ndim1]; double am[ndim1][ndim1], fm[ndim1], x[ndim1], y[ndim1]; input(&nnode, &nelmt, &nbc, ielmt, x, y, ibc); assem(nnode, nelmt, nbc, ielmt, x, y, ibc, am, fm); solve(nnode, am, fm); //gs_solve(nnode, am, fm); output(fm, nnode, x, y); return(exit_success); [1] : void output(double fm[], int nnode, double x[], double y[]); output(fm, nnode, x, y); 9 / 34
input() /* Input */ void input(int *nnode, int *nelmt, int *nbc, int ielmt[][3], double x[], double y[], int ibc[]) { int j; FILE *fpin, *fpout; if((fpin=fopen("indata.txt","r")) == NULL){ printf(" \n"); exit(exit_failure); if((fpout=fopen("outdata.txt","w")) == NULL){ fclose(fpin); printf(" \n"); exit(exit_failure); 10 / 34
/* Input of Data */ printf("input (nnode, nelmt, nbc): "); fscanf(fpin, "%d%d%d", nnode, nelmt, nbc); printf("\ninput(x,y) of node i: \n"); for(j=0; j<*nnode; j++){ printf(" for i=%d: \n", j+1); fscanf(fpin, "%lg%lg", &(x[j]), &(y[j])); printf("input(1st, 2nd, 3rd) nodes of element i: \n"); for(j=0; j<*nelmt; j++){ printf(" for i=%d: \n", j+1); fscanf(fpin, "%d%d%d", &(ielmt[j][0]), &(ielmt[j][1]), &(ielmt[j][2])); 11 / 34
if(*nbc>0){ printf("input nodes with zero Dirichlet data"); printf(" for i=1 to %d: ", *nbc); for(j=0; j<*nbc; j++) fscanf(fpin, "%d", &(ibc[j])); printf("\n"); printf("nodes of elements \n"); for(j=0; j<2; j++) printf(" i i1 i2 i3 "); printf("\n"); for(j=0; j<*nelmt; j+=2){ printf("*%4d*%4d%4d%4d ", j+1, ielmt[j][0], ielmt[j][1], ielmt[j][2]); if(j<*nelmt-1) printf("*%4d*%4d%4d%4d", j+2, ielmt[j+1][0], ielmt[j+1][1], ielmt[j+1][2]); printf("\n"); 12 / 34
if(*nbc>0){ printf("nodes with zero Dirichlet data \n"); for(j=0; j<*nbc; j++) printf("%4d", ibc[j]); printf("\n"); [1] : FILE *fpin, *fpout;... : indata.txt, : outdata.txt fscanf(fpin,... 13 / 34
assem() /* The Direct Stiffness Method */ void assem(int nnode, int nelmt, int nbc, int ielmt[][3], double x[], double y[], int ibc[], double am[ndim1][ndim1], double fm[]) { int i, j, ie, ii, jj; double ae[3][3], fe[3]; /* Initial Clear */ for(i=0; i<nnode; i++){ fm[i]=0.0; for(j=0; j<nnode; j++) am[i][j]=0.0; 14 / 34
/* Assemblage of Total Matrix and Vector */ for(ie=0; ie<nelmt; ie++){ ecm(ielmt, x, y, ie, ae, fe); for(i=0; i<3; i++){ ii=ielmt[ie][i]-1; fm[ii]+=fe[i]; for(j=0; j<3; j++){ jj=ielmt[ie][j]-1; am[ii][jj]+=ae[i][j]; 15 / 34
/* The Homogeneous Dirichlet Condition */ if (nbc!=0){ for(i=0;i<nbc; i++){ ii=ibc[i]-1; fm[ii]=0.0; for(j=0; j<nnode;j++){ am[ii][j]=0.0; am[j][ii]=0.0; am[ii][ii]=1.0; 16 / 34
ecm() /* Element Coeffcient Matrix and Free Vector */ void ecm(int ielmt[][3], double x[], double y[], int ie, double ae[][3], double fe[]) { int i, j, k; double d, s, xe[3], ye[3], b[3], c[3]; for(i=0; i<3;i++){ j=ielmt[ie][i]-1; xe[i]=x[j]; ye[i]=y[j]; d=xe[0]*(ye[1]-ye[2])+xe[1]*(ye[2]-ye[0]) +xe[2]*(ye[0]-ye[1]); s=fabs(d)/2.0; 17 / 34
/* Calculation of Element Coefficient Matrix */ for(i=0; i<3; i++){ j=i+1; if(i==2) j=0; k=i-1; if(i==0) k=2; b[i]=(ye[j]-ye[k])/d; c[i]=(xe[k]-xe[j])/d; for(i=0; i<3; i++){ for(j=0; j<3; j++) ae[i][j]=s*(b[j]*b[i]+c[j]*c[i]); /* Calculation of Element Free Vector */ for(i=0; i<3; i++) fe[i]=s*f(xe[i], ye[i])/3.0; 18 / 34
solve() /* The Gauss Elimination Method */ void solve(int m, double a[][ndim1], double f[]) /* m=number of unknowns */ { int i, j, k; double aa; /* Forward Elimination */ for(i=0; i<m-1; i++){ for(j=i+1; j<m; j++){ aa=a[j][i]/a[i][i]; f[j]-=aa*f[i]; for(k=i+1; k<m; k++) a[j][k]-=aa*a[i][k]; 19 / 34
/* Backward Substitution */ f[m-1]/=a[m-1][m-1]; for(i=m-2; i>=0; i--){ for(j=i+1; j<m; j++) f[i]-=a[i][j]*f[j]; f[i]/=a[i][i]; 20 / 34
[1] : gs solve() gs solve() /* The Gauss Seidel Method */ void gs_solve(int m, double a[][ndim1], double f[]) { int i, j, k; double aa = 0.0; double u[ndim1],u_old[ndim1]; double max_delta=0.0, max_delta_old; double max_u=0.0; int count = 0; for(i=0; i<m; i++){ u[i] = 0.0; 21 / 34
do{ max_delta_old = max_delta; for(i=0; i<m; i++){ u_old[i] = u[i]; for(i=0; i<m; i++){ aa = 0.0; for(j=0; j<m; j++){ if(i!= j){ aa += a[i][j]*u[j]; u[i]=(f[i] - aa)/a[i][i]; 22 / 34
max_delta = fabs(u[0] - u_old[0]); for(i=1; i<m; i++){ if(max_delta < fabs(u[i] - u_old[i])){ max_delta = fabs(u[i] - u_old[i]); max_u = fabs(u[0]); for(i=1; i<m; i++){ if(max_u < fabs(u[i])){ max_u = fabs(u[i]); 23 / 34
count++; while(fabs((max_delta - max_delta_old)/max_u) >= 1.0E-5); printf("count = %d\n",count); for(i=0; i<m; i++){ f[i] = u[i]; 24 / 34
output() /* 0utput of Approximate Nodal Values of u */ void output(double fm[], int nnode, double x[], double y[]) { int j,m; FILE *fpout; if((fpout=fopen("outdata.txt","a")) == NULL){ printf(" \n"); exit(exit_failure); 25 / 34
fprintf(fpout, "\nnodal values of u\n"); fprintf(fpout, " i u \n"); for(j=0; j<nnode; j++) fprintf(fpout, "%4d%12.3e\n", j+1, fm[j]); fprintf(fpout, "\nmathematica 6 Data\n"); fprintf(fpout, "ListPlot3D[{\n"); for(j=0; j<nnode-1; j++){ fprintf(fpout, "{%f, %f, %f,\n", x[j], y[j], fm[j]); fprintf(fpout, "{%f, %f, %f\n, Mesh -> All]\n", x[j], y[j], fm[j]); fclose(fpout); [1] : FILE *fpout;... 26 / 34
f() /* Function for Free Term */ double f(double x, double y) { return(1.0); 27 / 34
7.5 indata.txt 9 8 5 0 0 0 0.5 0 1 0.5 0 0.5 0.5 0.5 1 1 0 1 0.5 1 1 1 4 5 1 5 2 2 5 6 2 6 3 4 7 8 4 8 5 5 8 9 5 9 6 1 2 3 4 7 outdata.txt Basic data nnode = 9 nelmt = 8 nbc = 5 x,y-coordinates of nodes i x(i) y(i) 1 0.0000 0.0000 2 0.0000 0.5000 3 0.0000 1.0000 4 0.5000 0.0000 5 0.5000 0.5000 6 0.5000 1.0000 7 1.0000 0.0000 8 1.0000 0.5000 9 1.0000 1.0000 Nodes of elements i i1 i2 i3 1 1 4 5 2 1 5 2 3 2 5 6 4 2 6 3 5 4 7 8 6 4 8 5 7 5 8 9 8 5 9 6 28 / 34
Nodes with zero Dirichlet data 1 2 3 4 7 Nodal values of u i u 1 0.000e+000 2 0.000e+000 3 0.000e+000 4 0.000e+000 5 1.771e-001 6 2.292e-001 7 0.000e+000 8 2.292e-001 9 3.125e-001 Mathematica 6 Data ListPlot3D[{ {0.000000, 0.000000, 0.000000, {0.000000, 0.500000, 0.000000, {0.000000, 1.000000, 0.000000, {0.500000, 0.000000, 0.000000, {0.500000, 0.500000, 0.177083, {0.500000, 1.000000, 0.229167, {1.000000, 0.000000, 0.000000, {1.000000, 0.500000, 0.229167, {1.000000, 1.000000, 0.312500, Mesh -> All] 29 / 34
Mathematica outdata.txt ListPlot3D[...] Shift + Enter 30 / 34
7.6 1 2 Poisson Mathematica 31 / 34
2 indata 105.txt femfp.c femfp.c 32 / 34
7.7 2 Poisson 1 C 2 C 3 1 4 33 / 34
[1]. :.,, 1980. 34 / 34