alg125.c

c语言版的 数值计算方法

      /* output gamma */
      fprintf(*OUP, "Coefficients of Basis Functions:\n");
      for (I=1; I<=LM; I++) {
	 LLL = LL[I-1];
	 fprintf(*OUP, "%3d %13.8f %d", I, GAMMA[I-1], I);
	 for (J=1; J<=LLL; J++) fprintf(*OUP, " %4d", II[I-1][J-1]);
	 fprintf(*OUP, "\n");
      }
      fclose(*OUP);
      /* STEP 23 */
   }
   return 0;
}

double P(double X, double Y)
{
   double p; 

   p =  1.0;
   return p;
}

double Q(double X, double Y)
{
   double q; 

   q =  1.0;
   return q;
}

double R(double X, double Y)
{
   double r; 

   r =  0.0;
   return r;
}

double F(double X, double Y)
{
   double f; 

   f =  0.0;
   return f;
}

double G(double X, double Y)
{
   double g; 

   g =  4.0;
   return g;
}

double G1(double X, double Y)
{
   double g1; 

   g1 =  0.0;
   return g1;
}

double G2(double X, double Y)
{
   double g2, T, Z1; 

   T = 1.0e-05;
   Z1 = 0.0;
   if ( ( (0.2-T) <= X) && (X <= (0.4+T)) && (absval(Y-0.2) <= T)) 
      Z1 = X;
   if ( ( (0.5-T) <= X) && (X <= (0.6+T)) && (absval(Y-0.1) <= T))
      Z1 = X;
   if ( ( -T <= Y) && (Y <= (0.1+T)) && (absval(X-0.6) <= T)) 
      Z1 = Y;
   if ( ( -T <= X) && (X <= (0.2+T)) && (absval(Y+X-0.4) <= T)) 
      Z1 = (X+Y)/sqrt(2.0);
   if ( (0.4 -T <= X) && (X <= (0.5+T)) && (absval(Y+X-0.6) <= T)) 
      Z1 = (X+Y)/sqrt(2.0);
   g2 = Z1;
   return g2;
}

double RR(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I1, int I2)
{
   double rr; 

   rr = R(X,Y)*(A[I-1][I1-1]+B[I-1][I1-1]*X+C[I-1][I1-1]*Y)*(A[I-1][I2-1]
	+B[I-1][I2-1]*X+C[I-1][I2-1]*Y);
   return rr;
}

double FFF(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I1)
{
   double fff; 

   fff = F(X,Y)*(A[I-1][I1-1]+B[I-1][I1-1]*X+C[I-1][I1-1]*Y);
   return fff;
}

double GG1(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I1, int I2)
{
   double gg1; 

   gg1 = G1(X,Y)*(A[I-1][I1-1]+B[I-1][I1-1]*X+C[I-1][I1-1]*Y)*(A[I-1][I2-1]
	 +B[I-1][I2-1]*X+C[I-1][I2-1]*Y);
   return gg1;
}

double GG2(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I1)
{
   double gg2; 

   gg2 = G2(X,Y)*(A[I-1][I1-1]+B[I-1][I1-1]*X+C[I-1][I1-1]*Y);  
   return gg2;
}

double GG3(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I2)
{
   double gg3; 

   gg3 = G2(X,Y)*(A[I-1][I2-1]+B[I-1][I2-1]*X+C[I-1][I2-1]*Y);
   return gg3;
}

double GG4(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I1)
{
   double gg4; 

   gg4 = G1(X,Y) * (A[I-1][I1-1]+B[I-1][I1-1]*X+C[I-1][I1-1]*Y) * (A[I-1][I1-1]
	 +B[I-1][I1-1]*X+C[I-1][I1-1]*Y);
   return gg4;
}

double GG5(double X, double Y, double A[][3], double B[][3], double C[][3], int I, int I2)
{
   double gg5; 

   gg5 = G1(X,Y) * (A[I-1][I2-1]+B[I-1][I2-1]*X+C[I-1][I2-1]*Y) * 
	 (A[I-1][I2-1]+B[I-1][I2-1]*X+C[I-1][I2-1]*Y);
   return gg5;
}

double QQ(int L, double *XX, double *YY, double DELTA, int J1, int J2, int J3, int I1, int I2, double A[][3], double B[][3], double C[][3])
{
   double X[7], Y[7], S[7];
   double qq,QQQ,T1,T2,T3; 
   int I;

   X[0] = XX[J1-1]; Y[0] = YY[J1-1];
   X[1] = XX[J2-1]; Y[1] = YY[J2-1];
   X[2] = XX[J3-1]; Y[2] = YY[J3-1];
   X[3] = 0.5*(X[0]+X[1]); Y[3] = 0.5*(Y[0]+Y[1]);
   X[4] = 0.5*(X[0]+X[2]); Y[4] = 0.5*(Y[0]+Y[2]);
   X[5] = 0.5*(X[1]+X[2]); Y[5] = 0.5*(Y[1]+Y[2]);
   X[6] = (X[0]+X[1]+X[2])/3.0; Y[6] = (Y[0]+Y[1]+Y[2])/3.0;
   switch (L) {
      case 1:
	 for (I=1; I<=7; I++) S[I-1] = P(X[I-1],Y[I-1]);
	 break;
      case 2:
	 for (I=1; I<=7; I++) S[I-1] = Q(X[I-1],Y[I-1]);
	 break;
      case 3:
	 for (I=1; I<=7; I++) S[I-1] = RR(X[I-1],Y[I-1],A,B,C,I,I1,I2);
	 break;
      case 4:
	 for (I=1; I<=7; I++) S[I-1] = FFF(X[I-1],Y[I-1],A,B,C,I,I1);
	 break;
   }
   T1 = 0.0;
   for (I=1; I<=3; I++)  T1 = T1 + S[I-1];
   T2 = 0.0;
   for (I=4; I<=6; I++) T2 = T2 + S[I-1];
   T3 = S[6];
   QQQ = 0.5*(T1/20.0 + 2.0*T2/15.0 + 9.0*T3/20.0)*absval(DELTA);
   qq = QQQ;
   return qq;
}

double SQ(int L, double *XX, double *YY, int J1, int J2, int I, int I1, int I2, double A[][3], double B[][3], double C[][3])
{
   double S[101];
   double sq,SSQ,X,T1,T2,X1,X2,Y1,Y2,H,Q1,Q2,Q3,T3; 
   int K;

   X1 = XX[J1-1];
   Y1 = YY[J1-1];
   X2 = XX[J2-1];
   Y2 = YY[J2-1];
   H = 0.01;
   T1 = X2-X1;
   T2 = Y2-Y1;
   T3 = sqrt(T1*T1+T2*T2);
   for (K=1; K<=101; K++) {
       X = (K-1.0)*H;
       switch (L) {
	  case 1:
	     S[K-1] = T3*GG1(T1*X+X1,T2*X+Y1,A,B,C,I,I1,I2);
	     break;
	  case 2:
	     S[K-1] = T3*GG2(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
	     break;
	  case 3:
	     S[K-1] = T3*GG3(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
	     break;
	  case 4:
	     S[K-1] = T3*GG4(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
	     break;
	  case 5:
	     S[K-1] = T3*GG5(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
	     break;
      }      
   }
   Q3 = S[0]+S[100];
   Q1 = 0.0;
   Q2 = S[99];
   for (K=1; K<=49; K++) {
      Q1 = Q1 + S[2*K];
      Q2 = Q2 + S[2*K-1];
   }
   SSQ = (Q3 + 2.0*Q1 + 4.0*Q2)*H/3.0;
   sq = SSQ;
   return sq;
}

void INPUT(int *OK, int *K, int *N, int *M, int *LN1, int *LM, int *NL, double *XX, double *YY, int *LLL, int II[][5], int NV[][5], int *LL, int *N1, int *N2, int NTR[][3], int LINE[][2])
{
   int KK, I, J;
   char AA;
   char NAME[30];
   FILE *INP;

   printf("This is the Finite Element Method.\n");
   *OK = false;
   printf("The input will be from a text file in the following form:\n");
   printf("K     N     M     n     m     nl\n\n");
   printf("Each of the above is an integer -");
   printf("separate with at least one blank.\n\n");
   printf("Follow with the input for each node in the form:\n");
   printf("x-coor., y-coord., number of triangles in which the");
   printf(" node is a vertex.\n\n");
   printf("Continue with the triangle number and vertex number for\n");
   printf("each triangle in which the node is a vertex.\n");
   printf("Separate each entry with at least one blank.\n\n");
   printf("After all nodes have been entered follow with information\n");
   printf("on the lines over which line integrals must be computed.\n");
   printf("The format of this data will be the node number of the\n");
   printf("starting node, followed by the node number of the ending\n");
   printf("node for each line, taken in the positive direction.\n");
   printf("There should be 2 * nl such entries, each an integer\n");
   printf("separated by a blank.\n\n");
   printf("Have the functions P,Q,R,F,G,G1,G2 been created and\n");
   printf("has the input file been created?  Answer Y or N.\n");
   scanf("\n%c", &AA);
   if ((AA == 'Y') || (AA == 'y')) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("for example:   A:DATA.DTA\n");
      scanf("%s", NAME);
      INP = fopen(NAME, "r");
      *OK = true;
      fscanf(INP, "%d %d %d %d %d %d", K, N, M, LN1, LM, NL);
      for (KK=1; KK<=*LM; KK++) {
	 fscanf(INP, "%lf %lf %d", &XX[KK-1], &YY[KK-1], LLL);
	 for (J=1; J<=*LLL; J++) fscanf(INP, "%d %d", &II[KK-1][J-1], &NV[KK-1][J-1]);
	 LL[KK-1] = *LLL;
	 for (J=1; J<=*LLL; J++) {
	    *N1 = II[KK-1][J-1];
	    *N2 = NV[KK-1][J-1];
	    NTR[*N1-1][*N2-1] = KK;
	 }  
      }  
      if (*NL > 0) {
	 for (I=1; I<=*NL; I++)
	    for (J=1; J<=2; J++) fscanf(INP, "%d", &LINE[I-1][J-1]);
      }  
      fclose(INP);
   }
   else printf("The program will end so that the input file can be created.\n");
}

void OUTPUT(FILE **OUP)
{
   int FLAG;
   char NAME[30];

   printf("Choice of output method:\n");
   printf("1. Output to screen\n");
   printf("2. Output to text file\n");
   printf("Please enter 1 or 2.\n");
   scanf("%d", &FLAG);
   if (FLAG == 2) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("for example:   A:OUTPUT.DTA\n");
      scanf("%s", NAME);
      *OUP = fopen(NAME, "w");
   }
   else *OUP = stdout;
   fprintf(*OUP, "FINITE ELEMENT METHOD\n\n\n");
}

/* Absolute Value Function */
double absval(double val)
{
   if (val >= 0) return val;
   else return -val;
}