f_misc.c

cfd求解器使用与gmsh网格的求解

  T_y[4][3] = sin(beta);
  T_y[4][4] = cos(beta);
  T_y[5][0] = cos(beta)*sin(beta);
  T_y[5][2] =-cos(beta)*sin(beta);
  T_y[5][5] = SQU(cos(beta))-SQU(sin(beta));

  T_x[0][0] = 1.0;
  T_x[1][1] = SQU(cos(gamma));
  T_x[1][2] = SQU(sin(gamma));
  T_x[1][4] = 2*cos(gamma)*sin(gamma);
  T_x[2][1] = SQU(sin(gamma));    
  T_x[2][2] = SQU(cos(gamma));
  T_x[2][4] =-2*cos(gamma)*sin(gamma);
  T_x[3][3] = cos(gamma);
  T_x[3][5] = sin(gamma);
  T_x[4][1] =-cos(gamma)*sin(gamma);
  T_x[4][2] = cos(gamma)*sin(gamma);
  T_x[4][4] = SQU(cos(gamma))-SQU(sin(gamma));
  T_x[5][3] =-sin(gamma);
  T_x[5][5] = cos(gamma);
 
  for (i=0;i<n;i++){
    for (j=0;j<n;j++){
      sum=0.0;
      for(k=0;k<n;k++)
	sum=sum+T_z[i][k]*T_y[k][j];
      T_zy[i][j]=sum;
     }
   }
  for (i=0;i<n;i++){
    for (j=0;j<n;j++){
      sum=0.0;
      for(k=0;k<n;k++)
	sum=sum+T_zy[i][k]*T_x[k][j];
      T[i][j]=sum;
     }
   }

  free_dmatrix(T_z, 0,n-1,0,n-1);
  free_dmatrix(T_y, 0,n-1,0,n-1);
  free_dmatrix(T_x, 0,n-1,0,n-1);
  free_dmatrix(T_zy,0,n-1,0,n-1);
}

/* Calculate the transformation strain tensor taking into account only one rotation 
   alpha/beta/gamma around the axis z/y/x respectively */

void T_Strain(double **T, int n, double alpha, double beta, double gamma)
{
  int i,j,k;
  double **T_z, **T_y, **T_x, **T_zy, sum;
  T_z=dmatrix(0,n-1,0,n-1);
  T_y=dmatrix(0,n-1,0,n-1);
  T_x=dmatrix(0,n-1,0,n-1);
  T_zy=dmatrix(0,n-1,0,n-1);
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
       T_z[i][j] = 0;
       T_y[i][j] = 0;
       T_x[i][j] = 0;
    }    
  }
  T_z[0][0] = SQU(cos(alpha));
  T_z[0][1] = SQU(sin(alpha));
  T_z[0][3] = cos(alpha)*sin(alpha);
  T_z[1][0] = SQU(sin(alpha));    
  T_z[1][1] = SQU(cos(alpha));
  T_z[1][3] =-cos(alpha)*sin(alpha);
  T_z[2][2] = 1.0;
  T_z[3][0] =-2*cos(alpha)*sin(alpha);
  T_z[3][1] = 2*cos(alpha)*sin(alpha);
  T_z[3][3] = SQU(cos(alpha))-SQU(sin(alpha));
  T_z[4][4] = cos(alpha);
  T_z[4][5] =-sin(alpha);
  T_z[5][4] = sin(alpha);
  T_z[5][5] = cos(alpha);

  T_y[0][0] = SQU(cos(beta));
  T_y[0][2] = SQU(sin(beta));
  T_y[0][5] =-cos(beta)*sin(beta);
  T_y[1][1] = 1.0; 
  T_y[2][0] = SQU(sin(beta));    
  T_y[2][2] = SQU(cos(beta));
  T_y[2][5] = cos(beta)*sin(beta);
  T_y[3][3] = cos(beta);
  T_y[3][4] =-sin(beta);
  T_y[4][3] = sin(beta);
  T_y[4][4] = cos(beta);
  T_y[5][0] = 2*cos(beta)*sin(beta);
  T_y[5][2] =-2*cos(beta)*sin(beta);
  T_y[5][5] = SQU(cos(beta))-SQU(sin(beta));

  T_x[0][0] = 1.0;
  T_x[1][1] = SQU(cos(gamma));
  T_x[1][2] = SQU(sin(gamma));
  T_x[1][4] = cos(gamma)*sin(gamma);
  T_x[2][1] = SQU(sin(gamma));    
  T_x[2][2] = SQU(cos(gamma));
  T_x[2][4] =-cos(gamma)*sin(gamma);
  T_x[3][3] = cos(gamma);
  T_x[3][5] = sin(gamma);
  T_x[4][1] =-2*cos(gamma)*sin(gamma);
  T_x[4][2] = 2*cos(gamma)*sin(gamma);
  T_x[4][4] = SQU(cos(gamma))-SQU(sin(gamma));
  T_x[5][3] =-sin(gamma);
  T_x[5][5] = cos(gamma);
 
  for (i=0;i<n;i++){
    for (j=0;j<n;j++){
      sum=0.0;
      for(k=0;k<n;k++)
	sum=sum+T_z[i][k]*T_y[k][j];
      T_zy[i][j]=sum;
    }
  }
  for (i=0;i<n;i++){
    for (j=0;j<n;j++){
      sum=0.0;
      for(k=0;k<n;k++)
	sum=sum+T_zy[i][k]*T_x[k][j];
      T[i][j]=sum;
    }
  }
  free_dmatrix(T_z, 0,n-1,0,n-1);
  free_dmatrix(T_y, 0,n-1,0,n-1);
  free_dmatrix(T_x, 0,n-1,0,n-1);
  free_dmatrix(T_zy,0,n-1,0,n-1);

}

/* Tensor Transformation */

void  F_TransformTensor (F_ARG) {

  int    ident, i, j, k, ii = 0, jj = 0, *indx, N = 6 ;
  double **T_T, **T_S, **C, **C_trformed, **a, **y, d, *col, alpha, beta, Gamma;

  GetDP_Begin("F_TransformTensor");

  C = dmatrix(0,N-1,0,N-1);
  C_trformed = dmatrix(0,N-1,0,N-1);

  ident = (int)Fct->Para[0];
  alpha = Fct->Para[1];
  beta  = Fct->Para[2];
  Gamma = Fct->Para[3];

  if( ( (A+0)->Type != TENSOR && (A+0)->Type != TENSOR_SYM && (A+0)->Type != TENSOR_DIAG ) ||
      ( (A+1)->Type != TENSOR && (A+1)->Type != TENSOR_SYM && (A+1)->Type != TENSOR_DIAG ) ||
      ( (A+2)->Type != TENSOR && (A+2)->Type != TENSOR_SYM && (A+2)->Type != TENSOR_DIAG ) ||
      ( (A+3)->Type != TENSOR && (A+3)->Type != TENSOR_SYM && (A+3)->Type != TENSOR_DIAG ) )
    Msg(GERROR, "Function 'TransformTensor' requires 4 Tensors on input (NOT %s %s %s %s)",
               Get_StringForDefine(Field_Type,(A+0)->Type),
               Get_StringForDefine(Field_Type,(A+1)->Type),
               Get_StringForDefine(Field_Type,(A+2)->Type),
               Get_StringForDefine(Field_Type,(A+3)->Type) );

  for(i=0;i<N;i++)
    for(j=0;j<N;j++)
      C[i][j] = 0.0;   /* Reset C[i][j] to zeros */

  for(i=0;i<4;i++) { 
    switch(i){
    case 0 : ii=0; jj=0; break;
    case 1 : ii=0; jj=3; break;
    case 2 : ii=3; jj=0; break;
    case 3 : ii=3; jj=3; break;
    }      

    switch((A+i)->Type){
    case TENSOR :
      C[ii+0][jj+0] = (A+i)->Val[0]; C[ii+0][jj+1] = (A+i)->Val[1];  C[ii+0][jj+2] = (A+i)->Val[2]; 
      C[ii+1][jj+0] = (A+i)->Val[3]; C[ii+1][jj+1] = (A+i)->Val[4];  C[ii+1][jj+2] = (A+i)->Val[5];  
      C[ii+2][jj+0] = (A+i)->Val[6]; C[ii+2][jj+1] = (A+i)->Val[7];  C[ii+2][jj+2] = (A+i)->Val[8];  
      break;
  
    case TENSOR_DIAG :
      C[ii+0][jj+0]=(A+i)->Val[0];  C[ii+1][jj+1]=(A+i)->Val[1];  C[ii+2][jj+2]=(A+i)->Val[2]; 
      break;
      
    case TENSOR_SYM :
      C[ii+0][jj+0]=(A+i)->Val[0];  C[ii+0][jj+1]=(A+i)->Val[1];  C[ii+0][jj+2]=(A+i)->Val[2]; 
      C[ii+1][jj+0]=C[ii+0][jj+1];  C[ii+1][jj+1]=(A+i)->Val[3];  C[ii+1][jj+2]=(A+i)->Val[4];
      C[ii+2][jj+0]=C[ii+0][jj+2];  C[ii+2][jj+1]=C[ii+1][jj+2];  C[ii+2][jj+2]=(A+i)->Val[5];
      break;
    }
  }

  /* begin of transformation ! */

  T_T=dmatrix(0,N-1,0,N-1);
  T_S=dmatrix(0,N-1,0,N-1);
  a=dmatrix(0,N-1,0,N-1);
  y=dmatrix(0,N-1,0,N-1);    
  indx=ivector(0,N-1);
  col=dvector(0,N-1);
  
  T_Stress(T_T, N, alpha, beta, Gamma);     /* Create the T_sigma */
  T_Strain(T_S, N, alpha, beta, Gamma);   /* Create the T_epsilon */

  ludcmp1(T_T,N,indx,&d);                  /* LU decomposition of T_sigma */     
  for(j=0;j<N;j++) {
    for(i=0;i<N;i++) col[i]=0.0;
    col[j]=1.0;
    lubksb1(T_T,N,indx,col);               /* Forward and back substitution */
    for(i=0;i<N;i++) y[i][j]= col[i];        /* inverse( T_sigma ) = y */
  }
  for (i=0;i<N;i++){
    for (j=0;j<N;j++){
      a[i][j]=0.0;
      for(k=0;k<N;k++) a[i][j]=a[i][j]+y[i][k]*C[k][j];   /* inverse(T_sigma)*C[i][j] */
    }
  }
  for (i=0;i<N;i++){
    for (j=0;j<N;j++){
      C_trformed[i][j]=0.0;
      for(k=0;k<N;k++) C_trformed[i][j]=C_trformed[i][j]+a[i][k]*T_S[k][j]; 
      /* C[i][j]' = inverse(T_sigma)*C[i][j]*T_epsilon */
    }
  }

  switch (ident) {   
  case 11 : 
  (V+0)->Val[0] = C_trformed[0][0]; (V+0)->Val[1] = C_trformed[0][1]; (V+0)->Val[2] = C_trformed[0][2]; 
  (V+0)->Val[3] = C_trformed[1][0]; (V+0)->Val[4] = C_trformed[1][1]; (V+0)->Val[5] = C_trformed[1][2]; 
  (V+0)->Val[6] = C_trformed[2][0]; (V+0)->Val[7] = C_trformed[2][1]; (V+0)->Val[8] = C_trformed[2][2]; 
/* 
  fprintf(stderr,"(C11 0) %.16g %.16g %.16g\n", C_trformed[0][0], C_trformed[0][1], C_trformed[0][2]);
  fprintf(stderr,"(C11 1) %.16g %.16g %.16g\n", C_trformed[1][0], C_trformed[1][1], C_trformed[1][2]);
  fprintf(stderr,"(C11 2) %.16g %.16g %.16g\n", C_trformed[2][0], C_trformed[2][1], C_trformed[2][2]); 
*/
  break ;

  case 12 :
  (V+0)->Val[0] = C_trformed[0][3]; (V+0)->Val[1] = C_trformed[0][4]; (V+0)->Val[2] = C_trformed[0][5]; 
  (V+0)->Val[3] = C_trformed[1][3]; (V+0)->Val[4] = C_trformed[1][4]; (V+0)->Val[5] = C_trformed[1][5]; 
  (V+0)->Val[6] = C_trformed[2][3]; (V+0)->Val[7] = C_trformed[2][4]; (V+0)->Val[8] = C_trformed[2][5]; 
/*
  fprintf(stderr,"(C12 0) %.16g %.16g %.16g\n", C_trformed[0][3], C_trformed[0][4], C_trformed[0][5]);
  fprintf(stderr,"(C12 1) %.16g %.16g %.16g\n", C_trformed[1][3], C_trformed[1][4], C_trformed[1][5]);
  fprintf(stderr,"(C12 2) %.16g %.16g %.16g\n", C_trformed[2][3], C_trformed[2][4], C_trformed[2][5]); 
*/
  break ;

  case 21 :
  (V+0)->Val[0] = C_trformed[3][0]; (V+0)->Val[1] = C_trformed[3][1]; (V+0)->Val[2] = C_trformed[3][2]; 
  (V+0)->Val[3] = C_trformed[4][0]; (V+0)->Val[4] = C_trformed[4][1]; (V+0)->Val[5] = C_trformed[4][2]; 
  (V+0)->Val[6] = C_trformed[5][0]; (V+0)->Val[7] = C_trformed[5][1]; (V+0)->Val[8] = C_trformed[5][2]; 
/*
  fprintf(stderr,"(C21 0) %.16g %.16g %.16g\n", C_trformed[3][0], C_trformed[3][1], C_trformed[3][2]);
  fprintf(stderr,"(C21 1) %.16g %.16g %.16g\n", C_trformed[4][0], C_trformed[4][1], C_trformed[4][2]);
  fprintf(stderr,"(C21 2) %.16g %.16g %.16g\n", C_trformed[5][0], C_trformed[5][1], C_trformed[5][2]); 
*/
  break ;

  case 22 :
  (V+0)->Val[0] = C_trformed[3][3]; (V+0)->Val[1] = C_trformed[3][4]; (V+0)->Val[2] = C_trformed[3][5]; 
  (V+0)->Val[3] = C_trformed[4][3]; (V+0)->Val[4] = C_trformed[4][4]; (V+0)->Val[5] = C_trformed[4][5]; 
  (V+0)->Val[6] = C_trformed[5][3]; (V+0)->Val[7] = C_trformed[5][4]; (V+0)->Val[8] = C_trformed[5][5]; 
/*
  fprintf(stderr,"(C22 0) %.16g %.16g %.16g\n", C_trformed[3][3], C_trformed[3][4], C_trformed[3][5]);
  fprintf(stderr,"(C22 1) %.16g %.16g %.16g\n", C_trformed[4][3], C_trformed[4][4], C_trformed[4][5]);
  fprintf(stderr,"(C22 2) %.16g %.16g %.16g\n", C_trformed[5][3], C_trformed[5][4], C_trformed[5][5]);
*/
  break ;
  }   
  
  free_dmatrix(C, 0,N-1,0,N-1);
  free_dmatrix(C_trformed, 0,N-1,0,N-1);
  free_dmatrix(T_T, 0,N-1,0,N-1);
  free_dmatrix(T_S, 0,N-1,0,N-1);
  free_dmatrix(a, 0,N-1,0,N-1);
  free_dmatrix(y, 0,N-1,0,N-1);
  free_ivector(indx, 0,N-1);
  free_dvector(col, 0,N-1);

  V->Type = TENSOR; /* We do not know exactly the type of tensor, 
		       so TENSOR is a most general choice */

/* end of transformation */
  GetDP_End ;
}


/* Calculate the transformation permitivity tensor taking into account only one rotation 
   alpha/beta/gamma around the axis z/y/x respectively */

void T_Epsilon(double **T, double alpha, double beta, double gamma)
{
  T[0][0] = cos(alpha)* cos(beta);
  T[0][1] = sin(alpha)*cos(gamma)+cos(alpha)*sin(beta)*sin(gamma);
  T[0][2] = sin(alpha)*sin(gamma)-cos(alpha)*sin(beta)*cos(gamma);
  T[1][0] =-sin(alpha)*cos(beta);
  T[1][1] = cos(alpha)*cos(gamma)-sin(alpha)*sin(beta)*sin(gamma);
  T[1][2] = cos(alpha)*sin(gamma)+sin(alpha)*sin(beta)*cos(gamma);
  T[2][0] = sin(beta);
  T[2][1] =-cos(beta)*sin(gamma);
  T[2][2] = cos(beta)*cos(gamma);
}

void  F_TransformPerm (F_ARG) {
  
  int    i, j, k, N=3 ;
  double **T_eps,**EPSr,**EPSr_trformed,**a,**y,alpha,beta,Gamma;
  
  GetDP_Begin("F_TransformPerm");
  
  EPSr = dmatrix(0,N-1,0,N-1);
  EPSr_trformed = dmatrix(0,N-1,0,N-1);
  
  alpha = Fct->Para[0];
  beta  = Fct->Para[1];
  Gamma = Fct->Para[2];
  
  if ( A->Type != TENSOR && A->Type != TENSOR_SYM && A->Type != TENSOR_DIAG )
    Msg(GERROR, "Wrong type of argument for function 'TransformTensor2' (NOT %s) ",
	Get_StringForDefine(Field_Type,A->Type));
  
  for(i=0;i<N;i++)
    for(j=0;j<N;j++)
      EPSr[i][j] = 0.0;   /* Reset EPSr[i][j] to zeros */
  
  switch(A->Type){
  case TENSOR :
    EPSr[0][0] = A->Val[0]; EPSr[0][1] = A->Val[1];  EPSr[0][2] = A->Val[2]; 
    EPSr[1][0] = A->Val[3]; EPSr[1][1] = A->Val[4];  EPSr[1][2] = A->Val[5]; 
    EPSr[2][0] = A->Val[6]; EPSr[2][1] = A->Val[7];  EPSr[2][2] = A->Val[8]; 
    break;
    
  case TENSOR_DIAG :
    EPSr[0][0] = A->Val[0];  EPSr[1][1] = A->Val[1];  EPSr[2][2]=A->Val[2]; 
    break;
    
  case TENSOR_SYM :
    EPSr[0][0] = A->Val[0];   EPSr[0][1] = A->Val[1];   EPSr[0][2]=A->Val[2];
    EPSr[1][0] = EPSr[0][1];  EPSr[1][1] = A->Val[3];   EPSr[1][2]=A->Val[4];
    EPSr[2][0] = EPSr[0][2];  EPSr[2][1] = EPSr[1][2];  EPSr[2][2]=A->Val[5];
    break;
  }
  
  /* begin of transformation ! */
  
  T_eps=dmatrix(0,N-1,0,N-1);
  a=dmatrix(0,N-1,0,N-1);
  y=dmatrix(0,N-1,0,N-1);    
  
  T_Epsilon(T_eps, alpha, beta, Gamma);  
  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if (i!=j)  y[i][j]=T_eps[j][i];
      else       y[i][j]=T_eps[i][j];
    }
  }
  for(i=0;i<N;i++){
    for (j=0;j<N;j++){
      a[i][j]=0.0;
      for(k=0;k<N;k++) 
	a[i][j]=a[i][j]+y[i][k]*EPSr[k][j];  
    }
  }  
  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      EPSr_trformed[i][j]=0.0;
      for(k=0;k<N;k++) 
	EPSr_trformed[i][j]=EPSr_trformed[i][j]+a[i][k]*T_eps[k][j];
    }
  }
  
  V->Val[0] = EPSr_trformed[0][0]; V->Val[1] = EPSr_trformed[0][1]; V->Val[2] = EPSr_trformed[0][2]; 
  V->Val[3] = EPSr_trformed[1][0]; V->Val[4] = EPSr_trformed[1][1]; V->Val[5] = EPSr_trformed[1][2]; 
  V->Val[6] = EPSr_trformed[2][0]; V->Val[7] = EPSr_trformed[2][1]; V->Val[8] = EPSr_trformed[2][2]; 
/*
  fprintf(stderr,"(epsr 0) %.16g %.16g %.16g\n", EPSr_trformed[0][0], EPSr_trformed[0][1], EPSr_trformed[0][2]);
  fprintf(stderr,"(epsr 1) %.16g %.16g %.16g\n", EPSr_trformed[1][0], EPSr_trformed[1][1], EPSr_trformed[1][2]);
  fprintf(stderr,"(epsr 2) %.16g %.16g %.16g\n", EPSr_trformed[2][0], EPSr_trformed[2][1], EPSr_trformed[2][2]);
*/  
  free_dmatrix(EPSr, 0,N-1,0,N-1);
  free_dmatrix(EPSr_trformed, 0,N-1,0,N-1);
  free_dmatrix(T_eps, 0,N-1,0,N-1);
  free_dmatrix(a, 0,N-1,0,N-1);
  free_dmatrix(y, 0,N-1,0,N-1);
  
  V->Type = TENSOR; 
  
  /* end of transformation */
  GetDP_End ;
}


void  F_TransformPiezo (F_ARG) {

  int    ident, i, j, k, ii = 0, jj = 0, N = 6 ;
  double **T_eps, **T_S, **e, **e_trformed, **a, **y, alpha, beta, Gamma;

  GetDP_Begin("F_TransformPiezo");

  e = dmatrix(0,N/2-1,0,N-1);
  e_trformed = dmatrix(0,N/2-1,0,N-1);

  ident = (int)Fct->Para[0];
  alpha = Fct->Para[1];
  beta  = Fct->Para[2];
  Gamma = Fct->Para[3];