f_misc.c

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

  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 ) )
    Msg(GERROR, "Function 'TransformTensor' requires 2 Tensors on input (NOT %s %s )",
               Get_StringForDefine(Field_Type,(A+0)->Type),
               Get_StringForDefine(Field_Type,(A+1)->Type) );

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

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

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

  /* begin of transformation ! */

  T_S=dmatrix(0,N-1,0,N-1);
  T_eps=dmatrix(0,N/2-1,0,N/2-1);
  a=dmatrix(0,N/2-1,0,N-1);
  y=dmatrix(0,N/2-1,0,N/2-1);    
  
  T_Strain(T_S, N, alpha, beta, Gamma);   
  T_Epsilon(T_eps, alpha, beta, Gamma);  
  for(i=0;i<N/2;i++){
    for(j=0;j<N/2;j++){
      if (i!=j)  y[i][j]=T_eps[j][i];
      else    y[i][j]=T_eps[i][j];
    }
  }
  for(i=0;i<N/2;i++){
    for (j=0;j<N;j++){
      a[i][j]=0.0;
      for(k=0;k<N/2;k++) 
	a[i][j]=a[i][j]+y[i][k]*e[k][j];  
    }
  }  
  for(i=0;i<N/2;i++){
    for(j=0;j<N;j++){
      e_trformed[i][j]=0.0;
      for(k=0;k<N;k++) 
	e_trformed[i][j]=e_trformed[i][j]+a[i][k]*T_S[k][j];
    }
  }

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

  case 2 :
  V->Val[0] = e_trformed[0][3]; V->Val[1] = e_trformed[0][4]; V->Val[2] = e_trformed[0][5]; 
  V->Val[3] = e_trformed[1][3]; V->Val[4] = e_trformed[1][4]; V->Val[5] = e_trformed[1][5]; 
  V->Val[6] = e_trformed[2][3]; V->Val[7] = e_trformed[2][4]; V->Val[8] = e_trformed[2][5]; 
/*
  fprintf(stderr,"(e 12 0) %.16g %.16g %.16g\n", e_trformed[0][3], e_trformed[0][4], e_trformed[0][5]);
  fprintf(stderr,"(e 12 1) %.16g %.16g %.16g\n", e_trformed[1][3], e_trformed[1][4], e_trformed[1][5]);
  fprintf(stderr,"(e 12 2) %.16g %.16g %.16g\n", e_trformed[2][3], e_trformed[2][4], e_trformed[2][5]); 
*/
  break ;
  }   
  
  free_dmatrix(e, 0,N/2-1,0,N-1);
  free_dmatrix(e_trformed, 0,N/2-1,0,N-1);
  free_dmatrix(T_S, 0,N-1,0,N-1);
  free_dmatrix(T_eps, 0,N/2-1,0,N/2-1);
  free_dmatrix(a, 0,N/2-1,0,N-1);
  free_dmatrix(y, 0,N/2-1,0,N/2-1);

  V->Type = TENSOR; 

/* end of transformation */
  GetDP_End ;
}


void  F_TransformPiezoT (F_ARG) {

  int    ident, i, j, k,*indx, ii = 0, jj = 0, N=6 ;
  double **T_eps, **T_T, **e, **eT, **eT_trformed, d, *col, **a, **y, alpha, beta, Gamma;

  GetDP_Begin("F_TransformCij");

  e = dmatrix(0,N/2-1,0,N-1);
  eT = dmatrix(0,N-1,0,N/2-1);
  eT_trformed = dmatrix(0,N-1,0,N/2-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+3)->Type != TENSOR_SYM && (A+3)->Type != TENSOR_DIAG ) )
    Msg(GERROR, "Function 'TransformTensor' requires 2 Tensors on input (NOT %s %s )",
               Get_StringForDefine(Field_Type,(A+0)->Type),
               Get_StringForDefine(Field_Type,(A+1)->Type) );

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

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

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

  /* begin of transformation ! */

  T_T=dmatrix(0,N-1,0,N-1);
  T_eps=dmatrix(0,N/2-1,0,N/2-1);
  a=dmatrix(0,N-1,0,N/2-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);     
  T_Epsilon(T_eps, alpha, beta, Gamma);  
  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/2;j++){
      if (i!=j)  eT[i][j]=e[j][i];
      else    eT[i][j]=e[i][j];
    }
  }

  for (i=0;i<N;i++){
    for (j=0;j<N/2;j++){
      a[i][j]=0.0;
      for(k=0;k<N;k++) a[i][j]=a[i][j]+y[i][k]*eT[k][j];  
    }
  }
  for (i=0;i<N;i++){
    for (j=0;j<N/2;j++){
      eT_trformed[i][j]=0.0;
      for(k=0;k<N/2;k++) eT_trformed[i][j]=eT_trformed[i][j]+a[i][k]*T_eps[k][j]; 
    }
  }

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

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

  }   
  
  free_dmatrix(e, 0,N/2-1,0,N-1);
  free_dmatrix(eT, 0,N/2,0,N/2-1);
  free_dmatrix(eT_trformed, 0,N-1,0,N/2-1);
  free_dmatrix(T_T, 0,N-1,0,N-1);
  free_dmatrix(T_eps, 0,N/2-1,0,N/2-1);
  free_dmatrix(a, 0,N-1,0,N/2-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; 

/* end of transformation */
  GetDP_End ;
}

#endif /* #if !defined(HAVE_GSL) */



/* ------------------------------------------------------------------------ */
/*          Start:            V I R T U A L   W O R K                       */
/* ------------------------------------------------------------------------ */

void  JacobianVol_dx (struct Element * Element, 
			MATRIX3x3 * Jac, double DetJac,
			MATRIX3x3 * Jac_dx, double * DetJac_dx, int i) {

  GetDP_Begin("JacobianVol2D_dx");

  Jac_dx->c11 = Element->dndu[i][0] ;  Jac_dx->c12 = Element->dndu[i][0] ;  Jac_dx->c13 = Element->dndu[i][0] ;
  Jac_dx->c21 = Element->dndu[i][1] ;  Jac_dx->c22 = Element->dndu[i][1] ;  Jac_dx->c23 = Element->dndu[i][1] ;
  Jac_dx->c31 = Element->dndu[i][2] ;  Jac_dx->c32 = Element->dndu[i][2] ;  Jac_dx->c33 = Element->dndu[i][2] ;

  DetJac_dx[0] = Jac_dx->c11 * ( Jac->c22 * Jac->c33 - Jac->c23 * Jac->c32 )
               - Jac_dx->c21 * ( Jac->c12 * Jac->c33 - Jac->c13 * Jac->c32 )
               + Jac_dx->c31 * ( Jac->c12 * Jac->c23 - Jac->c22 * Jac->c13 );

  DetJac_dx[1] = - Jac_dx->c12 * ( Jac->c21 * Jac->c33 - Jac->c23 * Jac->c31 )
                 + Jac_dx->c22 * ( Jac->c11 * Jac->c33 - Jac->c13 * Jac->c31 )
                 - Jac_dx->c32 * ( Jac->c11 * Jac->c23 - Jac->c13 * Jac->c21 );

  DetJac_dx[2] = Jac_dx->c11 * ( Jac->c21 * Jac->c32 - Jac->c22 * Jac->c31 )
               - Jac_dx->c21 * ( Jac->c11 * Jac->c32 - Jac->c12 * Jac->c31 )
               + Jac_dx->c31 * ( Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21 );


  /* DetJac_dx[0] = Jac_dx->c11 * Jac->c22 - Jac_dx->c21 * Jac->c12 ; */
  /* DetJac_dx[1] = Jac_dx->c22 * Jac->c11 - Jac_dx->c12 * Jac->c21 ; */
  /* DetJac_dx[2] = 0. ; */

  if (DetJac < 0) {
    DetJac_dx[0] *= -1.;
    DetJac_dx[1] *= -1.;
    DetJac_dx[2] *= -1.;
  }

  GetDP_End ;
}


/* F_VirtualWork */

void  F_VirtualWork (F_ARG) {

  int i, numNode, Type_Dimension;
  double u, v, w;
  MATRIX3x3 Jac;
  MATRIX3x3 Jac_dx;
  double DetJac;
  double DetJac_dx [3];
  double valField[3], s[3] ;

  GetDP_Begin("F_VirtualWork");
  
  /*  numNode = (int)((A+1)->Val[0]); */

  numNode = Current.NumEntity;
  
  i = 0;
  while (i < Current.Element->GeoElement->NbrNodes && 
	 Current.Element->GeoElement->NumNodes[i]!=numNode)  i++;
  
  if (i < Current.Element->GeoElement->NbrNodes ) {
    /*
      printf("element : %d %d\n", Current.Element->GeoElement->Num, i+1);
    */
    
    valField[0] = A->Val[0];
    valField[1] = A->Val[1];
    valField[2] = A->Val[2];
    
    u = Current.u; v = Current.v; w = Current.w;
    Get_BFGeoElement(Current.Element, u,v,w);
    
    Type_Dimension = Current.GeoData->Dimension;
    
    switch (Type_Dimension) {
    case _2D :
      DetJac = JacobianVol2D (Current.Element, &Jac) ;
      JacobianVol_dx (Current.Element, &Jac, DetJac, &Jac_dx, DetJac_dx, i) ;
      
      s[0] =   DetJac_dx[0] *  ( - valField[0] * valField[0] + 
				   valField[1] * valField[1] + 
				   valField[2] * valField[2] )
	-  2 * DetJac_dx[1] *  valField[0] * valField[1]
	-  2 * DetJac_dx[2] *  valField[0] * valField[2];
      
      /*  DetJac_dx[0] *  (valField[1] * valField[1] - valField[0] * valField[0]) */
      /*       -  2 * DetJac_dx[1] *  valField[0] * valField[1] ; */
    
      s[1] =   DetJac_dx[1] *  ( valField[0] * valField[0] - 
				 valField[1] * valField[1] + 
				 valField[2] * valField[2] )
	-  2 * DetJac_dx[0] *  valField[1] * valField[0]
	-  2 * DetJac_dx[2] *  valField[1] * valField[2];
      
      /*      DetJac_dx[1] * (valField[0] * valField[0] - valField[1] * valField[1]) */
      /*       -  2 * DetJac_dx[0] *  valField[0] * valField[1] ; */
	
      s[2] =   DetJac_dx[2] *  ( valField[0] * valField[0] + 
				 valField[1] * valField[1] - 
				 valField[2] * valField[2] )
	-  2 * DetJac_dx[0] *  valField[2] * valField[0]
	-  2 * DetJac_dx[1] *  valField[2] * valField[1];
      
      s[0] /= fabs(DetJac);
      s[1] /= fabs(DetJac);
      s[2] /= fabs(DetJac);
      break;

    case _3D :
      DetJac = JacobianVol3D (Current.Element, &Jac) ;
      JacobianVol_dx (Current.Element, &Jac, DetJac, &Jac_dx, DetJac_dx, i) ;
      
      s[0] =   DetJac_dx[0] *  ( - valField[0] * valField[0] + 
				   valField[1] * valField[1] + 
				   valField[2] * valField[2] )
	-  2 * DetJac_dx[1] *  valField[0] * valField[1]
	-  2 * DetJac_dx[2] *  valField[0] * valField[2];
      
      /*  DetJac_dx[0] *  (valField[1] * valField[1] - valField[0] * valField[0]) */
      /*       -  2 * DetJac_dx[1] *  valField[0] * valField[1] ; */
      
      s[1] =   DetJac_dx[1] *  ( valField[0] * valField[0] - 
				 valField[1] * valField[1] + 
				 valField[2] * valField[2] )
	-  2 * DetJac_dx[0] *  valField[1] * valField[0]
	-  2 * DetJac_dx[2] *  valField[1] * valField[2];
      
      /*      DetJac_dx[1] * (valField[0] * valField[0] - valField[1] * valField[1]) */
      /*       -  2 * DetJac_dx[0] *  valField[0] * valField[1] ; */
      
      s[2] =   DetJac_dx[2] *  ( valField[0] * valField[0] +
				 valField[1] * valField[1] - 
				 valField[2] * valField[2] )
	-  2 * DetJac_dx[0] *  valField[2] * valField[0]
	-  2 * DetJac_dx[1] *  valField[2] * valField[1];
      
      s[0] /= fabs(DetJac);
      s[1] /= fabs(DetJac);
      s[2] /= fabs(DetJac);
      break;
      
    default :
      s[0] = 0.; s[1] = 0.; s[2] = 0.;
      break;
      
    }
  }
  else {
    s[0] = 0.; s[1] = 0.; s[2] = 0.;
  }
  
  V->Type = VECTOR ;
  V->Val[0] = s[0] ;
  V->Val[1] = s[1] ;
  V->Val[2] = s[2] ;

  GetDP_End ;
}

#undef F_ARG