eval_ax.c

二维热能方程的求解

#include <stdio.h>
#include <math.h>
typedef  double ATYPE;
typedef  double VTYPE;
typedef  float DTYPE;

void eval_Ax(nx,ny,level,coefI,coefA1,coefA2,A1,A2,x,b,ierr)

int nx,ny,level;
VTYPE coefI,coefA1,coefA2;
ATYPE A1[][nx][3];
ATYPE A2[][ny][3];
ATYPE x[][nx];
ATYPE b[][nx];
int *ierr;

{
   int ix,iy,nxm1,nym1;
   VTYPE eps;

   if(level>=1) printf("EVAL_AX: cI=%g cA1=%g cA2=%g\n",coefI,coefA1,coefA2);
   if(level>=2)
      printf(" nx=%d ny=%d level=%d\n",nx,ny,level);

   nxm1=nx-1;
   nym1=ny-1;
   eps=1.e-12;

   /*** contribution from I  ***/

   if(fabs(coefI)>eps)
      for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            b[iy][ix]=coefI*x[iy][ix];
   else               /** initialization **/
      for(iy=0;iy<ny;iy++)
         for(ix=0;ix<nx;ix++)
            b[iy][ix]=0.0;

   /*** contribution from A1 ***/

   if(fabs(coefA1)>eps){

      for(iy=0;iy<ny;iy++){
         ix=0;
            b[iy][ix]+=coefA1*(A1[iy][ix][1]*x[iy][ix]
                              +A1[iy][ix][2]*x[iy][ix+1]);
         for(ix=1;ix<nxm1;ix++)
            b[iy][ix]+=coefA1*(A1[iy][ix][0]*x[iy][ix-1]
                              +A1[iy][ix][1]*x[iy][ix]
                              +A1[iy][ix][2]*x[iy][ix+1]);
         ix=nxm1;
            b[iy][ix]+=coefA1*(A1[iy][ix][0]*x[iy][ix-1]
                              +A1[iy][ix][1]*x[iy][ix]);
      }
   }

   /*** contribution from A2 ***/

   if(fabs(coefA2)>eps){

      for(ix=0;ix<nx;ix++){
         iy=0;
            b[iy][ix]+=coefA2*(A2[ix][iy][1]*x[iy][ix]
                              +A2[ix][iy][2]*x[iy+1][ix]);
         for(iy=1;iy<nym1;iy++)
            b[iy][ix]+=coefA2*(A2[ix][iy][0]*x[iy-1][ix]
                              +A2[ix][iy][1]*x[iy][ix]
                              +A2[ix][iy][2]*x[iy+1][ix]);
         iy=nym1;
            b[iy][ix]+=coefA2*(A2[ix][iy][0]*x[iy-1][ix]
                              +A2[ix][iy][1]*x[iy][ix]);
      }
   }

}