getsrc.c

二维热能方程的求解

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

void getSRC(nx,ny,level,interpol,n1,n2,
            ax,bx,ay,by,tn,tsbgn,tsend,HS,src,ierr)

int   nx,ny,level,interpol,n1,n2;
VTYPE *ax,*bx,*ay,*by,*tn,*tsbgn,*tsend;
DTYPE HS[][n2];
ATYPE src[][nx];
int *ierr;

{
   int i,ix,iy,nearest;
   VTYPE x,y,dist,tmp;
   VTYPE hx,hy,px,py,xmax,xmin,ymax,ymin;

   if(level>=1) printf("GETSRC: interpol=%d\n",interpol);
   if(level>=3){
      printf(" nx=%d ny=%d",nx,ny);
      printf(" n1=%d n2=%d\n",n1,n2);
      printf(" ax=%g bx=%g ay=%g by=%g\n",*ax,*bx,*ay,*by);
      printf(" tsbgn=%g tsend=%g\n",*tsbgn,*tsend);
      }

   for(iy=0;iy<ny;iy++)
   for(ix=0;ix<nx;ix++)
      src[iy][ix]=0.e0;

   if(*tn<*tsbgn || *tn>*tsend){
      if(level>=1) printf(" tn is out of [tsbgn,tsend].\n");
      return; }

   hx=(*bx-*ax)/(nx-1);
   hy=(*by-*ay)/(ny-1);

   if(interpol){

      xmin=*bx; xmax=*ax;
      ymin=*by; ymax=*ay;
      for(i=0;i<n1;i++){
         if(xmin>HS[i][0]) xmin=HS[i][0];
         if(xmax<HS[i][0]) xmax=HS[i][0];
         if(ymin>HS[i][1]) ymin=HS[i][1];
         if(ymax<HS[i][1]) ymax=HS[i][1];
       }
      for(iy=0;iy<ny;iy++){
         py=hy*iy+(*ay);
         if(py>=ymin && py<=ymax){
            for(ix=0;ix<nx;ix++){
               px=hx*ix+(*ax);
               if(px>=xmin && px<=xmax){
                  dist=2.*(*bx-*ax);
                  for(i=0;i<n1;i++){
                     x=HS[i][0]; y=HS[i][1];
                     tmp=sqrt(x*x+y*y);
                     if(tmp<dist){nearest=i;dist=tmp;}
                  }
                 src[iy][ix]=HS[nearest][2];
               }
             }
          }
       }
    }

   else{
      for(i=0;i<n1;i++){
         x=HS[i][0]; y=HS[i][1];
         ix=(int)((x-*ax)/hx); if(ix<0) ix=0;
         iy=(int)((y-*ay)/hy); if(iy<0) iy=0;
         if((x-(*ax+hx*ix))>(0.5*hx)) ix++; if(ix>(nx-1)) ix=nx-1;
         if((y-(*ay+hy*iy))>(0.5*hy)) iy++; if(iy>(ny-1)) iy=ny-1;
         src[iy][ix]=HS[i][2];
       }
    }

}