field.c

来自「波浪数值模拟」· C语言 代码 · 共 2,306 行 · 第 1/4 页

    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  
    - A, B, and C must have the same sizes                     */

void field2D_self_add_two(field2D *T, const field2D *A, const field2D *B)
{

  int N = T->N;
  int M = T->M;
  int i;
  int NUM;

  double *td = T->data;
  const double *ad;
  const double *bd;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M) ) { 
    field_error_message("Bad M field2D sizes to void field2D_self_add_two()");
  }

  if (A->N != B->N) {
    field_error_message("Bad N field2D sizes to void field2D_self_add_two()");
  }

  switch (A->N - N) {
  case 0:       /* if A & T have same sizes, OK, start at same location */
    ad = A->data;
    bd = B->data;
    td = T->data;
    NUM = T->num;
    break;
  case 2: 
    ad = &(A->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    bd = &(B->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    td = T->data;
    NUM = T->num;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ad = A->data;
    bd = B->data;
    td = &(T->data[M]);
    NUM = A->num;
    break;
  default:   field_error_message("Bad N field2D sizes to void field2D_self_add_two()");
  }


  /* Then go through and do the addition */

  for (i=0;i<NUM;i++) {
    *td++ += (*ad++) + (*bd++); 
  }
}


/* Field2D math routinues :   T += A+B+C  
    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  
    - A, B, and C must have the same sizes                     */

void field2D_self_add_three(field2D *T, const field2D *A, const field2D *B, const field2D *C)
{

  int N = T->N;
  int M = T->M;
  int i;
  int NUM;

  double *td = T->data;
  const double *ad;
  const double *bd;
  const double *cd;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M) ) { 
    field_error_message("Bad M field2D sizes to void field2D_self_add_three()");
  }

  if ((A->N != B->N)||(A->N != C->N)) {
    field_error_message("Bad N field2D sizes to void field2D_self_add_three()");
  }

  switch (A->N - N) {
  case 0:       /* if A & T have same sizes, OK, start at same location */
    ad = A->data;
    bd = B->data;
    cd = C->data;
    td = T->data;
    NUM = T->num;
    break;
  case 2: 
    ad = &(A->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    bd = &(B->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    cd = &(C->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    td = T->data;
    NUM = T->num;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ad = A->data;
    bd = B->data;
    cd = C->data;
    td = &(T->data[M]);
    NUM = A->num;
    break;
  default:   field_error_message("Bad N field2D sizes to void field2D_self_add_three()");
  }


  /* Then go through and do the addition */

  for (i=0;i<NUM;i++) {
    *td++ += (*ad++) + (*bd++) + (*cd++); 
  }
}


void field2D_add_const(field2D *T, const double a)
{
  int i;

  for (i=0;i<T->num;i++)
    T->data[i] += a;

}


/* Field2D math routinues :   T *= A  
    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  
    - if A is (N-2,M) then the multiplyis done from T(2:N-1,M) 
*/

void field2D_self_multiply(field2D *T, const field2D *A)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M) ) {
    field_error_message("Bad M field2D sizes to void field2D_self_multiply()");
  }

  switch (A->N - N) {
  case 0:       /* if A & T have same sizes, OK, start at same location */
    ad = A->data;
    td = T->data;
    NUM = T->num;
    break;
  case 2: 
    ad = &(A->data[M]); /* if A(N+2,M) & T(N,M) then move ad pointer up */
    td = T->data;
    NUM = T->num;
    break;
  case -2:             /* if A(N,M) & T(N+2,M) then move td pointer up */
    ad = A->data;
    td = &(T->data[M]);
    NUM = A->num;
    break;
  default:   field_error_message("Bad N field2D sizes to void field2D_self_multiply()");
  }

  if (A->N == N+2) 
    ad = &(A->data[M]);
  else 
    ad = &(A->data[0]);

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++) {
    //tmp = *ad++;
    *td++ *= (*ad++);
  }

}



/* Field2D math routinues :   T *= (A*B)  
    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  
    - A and B must have the same sizes                        */
void field2D_self_multiply_two(field2D *T, const field2D *A, const field2D *B)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M)  || (B->M != M) || (A->N != B->N) || !( ( A->N == N) || (A->N == N+2)) ) {
    field_error_message("Bad field2D sizes to void field2D_multiply()");
  }

  if (A->N == N+2) { 
    ad = &(A->data[M]);
    bd = &(B->data[M]);
  }
  else {
    ad = &(A->data[0]);
    bd = &(B->data[0]);
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++) {
    *td++ *= (*ad++) * (*bd++);
  }

}


void field2D_multiply_add(field2D *T, const field2D *A, const field2D *B)   /* T  += A * B */
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;
  const double *bd = B->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->N != N) || (B->N != N) || (A->M != M) || (B->M != M) ) {
    field_error_message("Bad field2D sizes to void field2D_multiply()");
  }

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ += (*ad++) * (*bd++);

}




/* Field2D math routinues :   T /= A  
    If T has size (N,M), then A must also be size (N,M) OR
    - if A is (N+2,M) then the multiply is done from (2:N-1,M)  */

void field2D_self_divide(field2D *T, const field2D *A)
{

  int NUM = T->num;
  int N = T->N;
  int M = T->M;
  int i;

  double *td = T->data;
  const double *ad = A->data;

  /* Check that the dimensions are ok with an #ifdef */

  if ( (A->M != M)  || !( ( A->N == N) || (A->N == N+2)) ) {
    field_error_message("Bad field2D sizes to void field2D_self_divide()");
  }

  if (A->N == N+2) 
    ad = REF2(A, 1);
  else 
    ad = REF2(A, 0); 

  /* Then go through and do the multiplication */

  for (i=0;i<NUM;i++)
    *td++ /= (*ad++);

}

/* File handling functions */


FILE*  open_file_ptr_for_writing(char *s)
{
  FILE *fp;

  if ((fp=fopen(s,"w"))==NULL) {
    fprintf(stderr,"Cannot open file for writing: %s\n",s);
    exit(-1);
  }

  return fp;
}

FILE*  open_file_ptr_for_reading(char *s)
{
  FILE *fp;

  if ((fp=fopen(s,"r"))==NULL) {
    fprintf(stderr,"Cannot open file for writing: %s\n",s);
    exit(-1);
  }

  return fp;
}

/* field1d Input functions */


int load_field1D_ascii(char *fname, field1D *T)    /* returns 1 if successfull */
{
   FILE *fp;
   int i,numc,rval=1;
   double a;

   fp = open_file_ptr_for_reading(fname);

   for (i=0;i<T->N;i++) {
     numc = fscanf(fp,"%lf",&a);
     if ((numc == 0)||(numc==EOF))
       rval = 0;
     DR1(T,i) = a;
   }

   fclose(fp);
   return rval;
}

int load_field1D_binary(char *fname, field1D *T)    /* returns 1 if successfull */
{
   FILE *fp;
   int rval=1;
   size_t sized,numt,count2;

   fp = open_file_ptr_for_reading(fname);

   numt = (size_t ) T->N;

   sized = sizeof(double);
   count2 = fread(T->data,sized,numt,fp);

   if (numt != count2)
     rval=0;

   fclose(fp);
   return rval; 
}

/* field1d Output functions */


void save_field1D_ascii(char *fname, field1D *T)
{
  FILE *fp;
  int i;

  fp = open_file_ptr_for_writing(fname);

  for (i=0;i<T->N;i++) {
    fprintf(fp,"%g\n",DR1(T,i));
  }
  fprintf(fp,"\n");


  fclose(fp);

}


void save_field1D_binary(char *fname, field1D *T)
{
  FILE *fp;
  size_t size;

  fp = open_file_ptr_for_writing(fname);

  size = sizeof(double);

  fwrite(T->data,size,T->N,fp);

  fclose(fp);

}


/* field2d Input functions */


int load_field2D_ascii(char *fname, field2D *T)    /* returns 1 if successfull */
{
   FILE *fp;
   int i,j,numc,rval=1;
   double a;

   fp = open_file_ptr_for_reading(fname);

   for (i=0;i<T->N;i++) {
     for (j=0;j<T->M;j++) {
        numc = fscanf(fp,"%lf",&a);
        if ((numc == 0)||(numc==EOF)) {
	  rval = 0;
          return rval;
	}
        DR2(T,i,j) = a;
     }
   }
   fclose(fp);
   return rval;
}

int load_field2D_binary(char *fname, field2D *T)    /* returns 1 if successfull */
{
   FILE *fp;
   int rval=1;
   size_t sized,numt,count2;

   fp = open_file_ptr_for_reading(fname);

   numt = (size_t ) T->num;

   sized = sizeof(double);
   count2 = fread(T->data,sized,numt,fp);

   if (numt != count2)
     rval=0;

   fclose(fp);
   return rval; 
}



void print_field1D(field1D *T)
{
  int i;

  for (i=0;i<T->N;i++) {
    fprintf(stdout,"%g ",DR1(T,i));
  }
  fprintf(stdout,"\n");
}


/* field2d Output functions */


void print_field2D(field2D *T)
{

  int i,j;

  for (i=0;i<T->N;i++) {
    for (j=0;j<T->M;j++) 
      fprintf(stdout,"%g ",DR2(T,i,j));
    fprintf(stdout,"\n");
  }

}


void save_field2D_ascii(char *fname, field2D *T)
{
  FILE *fp;
  int i,j;

  fp = open_file_ptr_for_writing(fname);

  for (i=0;i<T->N;i++) {
    for (j=0;j<T->M;j++) 
      fprintf(fp,"%g ",DR2(T,i,j));
    fprintf(fp,"\n");
  }

  fclose(fp);

}


void save_field2D_binary(char *fname, field2D *T)
{
  FILE *fp;
  size_t size;

  fp = open_file_ptr_for_writing(fname);

  size = sizeof(double);

  fwrite(T->data,size,T->num,fp);

  fclose(fp);

}


double max_array(double *arr, int num)
{
  int i,index=-1;
  double tmp,max=-1.0;

  for (i=0;i<num;i++) {
    tmp = fabs(arr[i]);
    if (tmp>max) {
      index = i; 
      max = tmp;
    }
  }

  printf("Max array = %g  at location %d\n",max,index);
  
  return max;
}

/* ---------- C-GRID OPERATIONS FOR FIELD2D ----------- 

   INCLUDES: \delta_x, \delta_y  (transforming from one grid to another)
             \bar^x , \bar^y  and adding routines 
                                       
	     ------------------------------------------------------------ */


/* Calculates a \bar^x F and puts it into FB 
   - requires FB have size (N-1,M) or (N+1,M) relative to F (N,M)        
   - if size=(N+1,M) then the i=0 & N points are left blank          */

void field2D_barx(field2D *FB, const field2D *F)
{
  int NF = F->N;
  int NFB = FB->N;
  const int M = F->M;

  int i, j;

  const double *fdup = &(F->data[M]);
  const double *fddn = &(F->data[0]);
  double *dfd;
  double *dfd0;

  if ( (NFB != NF-1) && (NFB != NF+1) ) {
    field_error_message("Bad N field2D sizes to field2D_barx()");
  }

  if (FB->M != M) {
    field_error_message("Bad M field2D sizes to field2D_barx()");
  }

  if (NFB == NF-1) {
    dfd = &(FB->data[0]);
    for (i=0;i<NF-1;i++) {
      for (j=0;j<M;j++) {
	*dfd++ = 0.5 * ( *fdup++ + *fddn++);
      }
    }
  }
  else {
    dfd = &(FB->data[M]);
    for (i=0;i<NF-1;i++) {
      for (j=0;j<M;j++) {
	*dfd++ = 0.5 * ( *fdup++ + *fddn++);
      }
    }
    /* Now do the near and far boundaries */
    /* if FB has N+2 dimension relative to F */

    dfd0 = REF2( FB, 0 );
    fddn = REF2( F, 0);         /* point to the first row of F */
    fdup = REF2( F, NF-1);      /* point to the last row of F */
    for (j=0;j<M;j++) {
      *dfd++ = *fdup++;           /* first the far boundary */