get_geometry.c

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

  for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
    Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
    Jac->c21 += Element->x[i] * Element->dndu[i][1] ;

    Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
    Jac->c22 += Element->y[i] * Element->dndu[i][1] ;

    Jac->c13 += Element->z[i] * Element->dndu[i][0] ;
    Jac->c23 += Element->z[i] * Element->dndu[i][1] ;
  }

  DetJac = sqrt( DSQU(Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21)
		 + DSQU(Jac->c13 * Jac->c21 - Jac->c11 * Jac->c23)
		 + DSQU(Jac->c12 * Jac->c23 - Jac->c13 * Jac->c22) ) ;

  GetDP_Return(DetJac) ;
}

/* ------------------------------------------------------------------------ */
/*  J a c o b i a n L i n                                                   */
/* ------------------------------------------------------------------------ */

double  JacobianLin3D (struct Element * Element, MATRIX3x3 * Jac) {
  int  i ;
  double DetJac ;

  GetDP_Begin("JacobianLin3D");

  Jac->c11 = 0. ;  Jac->c12 = 0. ;  Jac->c13 = 0. ;
  Jac->c21 = 0. ;  Jac->c22 = 0. ;  Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = 0. ;

  for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
    Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
    Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
    Jac->c13 += Element->z[i] * Element->dndu[i][0] ;
  }

  DetJac = sqrt(DSQU(Jac->c11)+DSQU(Jac->c12)+DSQU(Jac->c13)) ;

  GetDP_Return(DetJac) ;
}


/* ------------------------------------------------------------------------ */
/*  G e t _ I n v e r s e M a t r i x                                       */
/* ------------------------------------------------------------------------ */

void  Get_InverseSingularMatrix(MATRIX3x3 * Mat, MATRIX3x3 * InvMat) {
  double T[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
  int i, j, k;
#if defined(HAVE_GSL)
  gsl_matrix *M, *V;
  gsl_vector *W, *TMPVEC;
  double ww;
#else
  double **M, **V, *W;
#endif

  GetDP_Begin("Get_InverseSingularMatrix");

  InvMat->c11 = InvMat->c12 = InvMat->c13 = 0.0;
  InvMat->c21 = InvMat->c22 = InvMat->c23 = 0.0;
  InvMat->c31 = InvMat->c32 = InvMat->c33 = 0.0;

#if defined(HAVE_GSL)
  M = gsl_matrix_alloc(3, 3);
  V = gsl_matrix_alloc(3, 3);
  W = gsl_vector_alloc(3);
  TMPVEC = gsl_vector_alloc(3);

  gsl_matrix_set(M, 0, 0, Mat->c11);
  gsl_matrix_set(M, 0, 1, Mat->c12);
  gsl_matrix_set(M, 0, 2, Mat->c13);
  gsl_matrix_set(M, 1, 0, Mat->c21);
  gsl_matrix_set(M, 1, 1, Mat->c22);
  gsl_matrix_set(M, 1, 2, Mat->c23);
  gsl_matrix_set(M, 2, 0, Mat->c31);
  gsl_matrix_set(M, 2, 1, Mat->c32);
  gsl_matrix_set(M, 2, 2, Mat->c33);

  gsl_linalg_SV_decomp(M, V, W, TMPVEC);
  for(i = 0; i < 3; i++) {
    for(j = 0; j < 3; j++) {
      ww = gsl_vector_get(W, i);
      if(fabs(ww) > 1.e-16) {   /* singular value precision */
        T[i][j] += gsl_matrix_get(M, j, i) / ww;
      }
    }
  }
  for(k=0 ; k<3 ; k++){
    InvMat->c11 += gsl_matrix_get(V, 0, k) * T[k][0] ;
    InvMat->c12 += gsl_matrix_get(V, 0, k) * T[k][1] ;
    InvMat->c13 += gsl_matrix_get(V, 0, k) * T[k][2] ;
    InvMat->c21 += gsl_matrix_get(V, 1, k) * T[k][0] ;
    InvMat->c22 += gsl_matrix_get(V, 1, k) * T[k][1] ;
    InvMat->c23 += gsl_matrix_get(V, 1, k) * T[k][2] ;
    InvMat->c31 += gsl_matrix_get(V, 2, k) * T[k][0] ;
    InvMat->c32 += gsl_matrix_get(V, 2, k) * T[k][1] ;
    InvMat->c33 += gsl_matrix_get(V, 2, k) * T[k][2] ;
  }
  gsl_matrix_free(M);
  gsl_matrix_free(V);
  gsl_vector_free(W);
  gsl_vector_free(TMPVEC);
#else
  M = dmatrix(1,3,1,3);
  V = dmatrix(1,3,1,3);
  W = dvector(1,3);

  M[1][1] = Mat->c11 ; M[1][2] = Mat->c12 ; M[1][3] = Mat->c13 ;
  M[2][1] = Mat->c21 ; M[2][2] = Mat->c22 ; M[2][3] = Mat->c23 ;
  M[3][1] = Mat->c31 ; M[3][2] = Mat->c32 ; M[3][3] = Mat->c33 ;

  dsvdcmp(M, 3, 3, W, V);

  /* cf. Numerical Recipes in C, p. 62 */

  for(i=1 ; i<=3 ; i++)
    for(j=1 ; j<=3 ; j++)
      if(fabs(W[i]) > 1.e-16) /* precision */
	T[i-1][j-1] += M[j][i] / W[i] ;
  
  for(k=1 ; k<=3 ; k++){
    InvMat->c11 += V[1][k] * T[k-1][0] ;
    InvMat->c12 += V[1][k] * T[k-1][1] ;
    InvMat->c13 += V[1][k] * T[k-1][2] ;
    InvMat->c21 += V[2][k] * T[k-1][0] ;
    InvMat->c22 += V[2][k] * T[k-1][1] ;
    InvMat->c23 += V[2][k] * T[k-1][2] ;
    InvMat->c31 += V[3][k] * T[k-1][0] ;
    InvMat->c32 += V[3][k] * T[k-1][1] ;
    InvMat->c33 += V[3][k] * T[k-1][2] ;
  }

  free_dmatrix(M,1,3,1,3);
  free_dmatrix(V,1,3,1,3);
  free_dvector(W,1,3);
#endif

  GetDP_End ;
}

void  Get_InverseMatrix(int Type_Dimension, int Type_Element, double DetMat,
			MATRIX3x3 * Mat, MATRIX3x3 * InvMat) {

  GetDP_Begin("Get_InverseMatrix");

  switch (Type_Dimension) {

  case _0D :
    InvMat->c11 = InvMat->c22 = InvMat->c33 = 1. ;
    InvMat->c12 = InvMat->c21 = 0. ;
    InvMat->c13 = InvMat->c31 = 0. ;
    InvMat->c23 = InvMat->c32 = 0. ;
    break ;

  case _1D :
    InvMat->c11 = 1. / Mat->c11 ;
    InvMat->c22 = 1. / Mat->c22 ;
    InvMat->c33 = 1. / Mat->c33 ;
    InvMat->c12 = InvMat->c21 = 0. ;
    InvMat->c13 = InvMat->c31 = 0. ;
    InvMat->c23 = InvMat->c32 = 0. ;
    break ;

  case _2D :
    
    switch(Type_Element){

    case TRIANGLE   : case TRIANGLE_2   : 
    case QUADRANGLE : case QUADRANGLE_2 :
      if(!DetMat) Msg(GERROR, "Null determinant in 'Get_InverseMatrix'");
      InvMat->c11 =   Mat->c22 * Mat->c33 / DetMat ;
      InvMat->c21 = - Mat->c21 * Mat->c33 / DetMat ;
      InvMat->c12 = - Mat->c12 * Mat->c33 / DetMat ;
      InvMat->c22 =   Mat->c11 * Mat->c33 / DetMat ;
      InvMat->c13 = InvMat->c23 = InvMat->c31 = InvMat->c32 = 0. ;
      InvMat->c33 =   1. / Mat->c33 ;
      break ;

    default : 
      Get_InverseSingularMatrix(Mat, InvMat);
      break;

    }
    break;

  case _3D :

    switch(Type_Element){

    case TETRAHEDRON : case TETRAHEDRON_2 : 
    case HEXAHEDRON  : case HEXAHEDRON_2  : 
    case PRISM       : case PRISM_2       :
    case PYRAMID     : case PYRAMID_2     :
      if(!DetMat) Msg(GERROR, "Null determinant in 'Get_InverseMatrix'");
      InvMat->c11 =  ( Mat->c22 * Mat->c33 - Mat->c23 * Mat->c32 ) / DetMat ;
      InvMat->c21 = -( Mat->c21 * Mat->c33 - Mat->c23 * Mat->c31 ) / DetMat ;
      InvMat->c31 =  ( Mat->c21 * Mat->c32 - Mat->c22 * Mat->c31 ) / DetMat ;
      InvMat->c12 = -( Mat->c12 * Mat->c33 - Mat->c13 * Mat->c32 ) / DetMat ;
      InvMat->c22 =  ( Mat->c11 * Mat->c33 - Mat->c13 * Mat->c31 ) / DetMat ;
      InvMat->c32 = -( Mat->c11 * Mat->c32 - Mat->c12 * Mat->c31 ) / DetMat ;
      InvMat->c13 =  ( Mat->c12 * Mat->c23 - Mat->c13 * Mat->c22 ) / DetMat ;
      InvMat->c23 = -( Mat->c11 * Mat->c23 - Mat->c13 * Mat->c21 ) / DetMat ;
      InvMat->c33 =  ( Mat->c11 * Mat->c22 - Mat->c12 * Mat->c21 ) / DetMat ;
      break;

    default :
      Get_InverseSingularMatrix(Mat, InvMat);
      break;

    }
    break;

  default :
    Msg(GERROR, "Wrong dimension in 'Get_InverseMatrix'");
    break ;

  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e t _ P r o d u c t M a t r i x                                       */
/* ------------------------------------------------------------------------ */

void  Get_ProductMatrix(int Type_Dimension,
			MATRIX3x3 * A, MATRIX3x3 * B, MATRIX3x3 * AB) {

  GetDP_Begin("Get_ProductMatrix");

  switch (Type_Dimension) {

  case _2D :
    AB->c11 = A->c11 * B->c11 + A->c12 * B->c21 ;
    AB->c12 = A->c11 * B->c12 + A->c12 * B->c22 ;
    AB->c21 = A->c21 * B->c11 + A->c22 * B->c21 ;
    AB->c22 = A->c21 * B->c12 + A->c22 * B->c22 ;
    break ;

  case _3D :
    AB->c11 = A->c11 * B->c11 + A->c12 * B->c21 + A->c13 * B->c31 ;
    AB->c12 = A->c11 * B->c12 + A->c12 * B->c22 + A->c13 * B->c32 ;
    AB->c13 = A->c11 * B->c13 + A->c12 * B->c23 + A->c13 * B->c33 ;
    AB->c21 = A->c21 * B->c11 + A->c22 * B->c21 + A->c23 * B->c31 ;
    AB->c22 = A->c21 * B->c12 + A->c22 * B->c22 + A->c23 * B->c32 ;
    AB->c23 = A->c21 * B->c13 + A->c22 * B->c23 + A->c23 * B->c33 ;
    AB->c31 = A->c31 * B->c11 + A->c32 * B->c21 + A->c33 * B->c31 ;
    AB->c32 = A->c31 * B->c12 + A->c32 * B->c22 + A->c33 * B->c32 ;
    AB->c33 = A->c31 * B->c13 + A->c32 * B->c23 + A->c33 * B->c33 ;
    break ;

  }

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*   G e t _ C h a n g e O f C o o r d i n a t e s                          */
/* ------------------------------------------------------------------------ */


void *Get_ChangeOfCoordinates(int Flag_ChangeCoord, int Type_Form) {

  GetDP_Begin("Get_ChangeOfCoordinates");

  switch (Type_Form) {
    
  case SCALAR :
  case FORM0  :  
    GetDP_Return((void *)ChangeOfCoord_No1) ;

  case FORM1 :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1  :
		 (void *)ChangeOfCoord_No123) ;

  case FORM2 :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form2 :
		 (void *)ChangeOfCoord_No123) ;
    
  case FORM3 :  case FORM3P :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form3 :
		 (void *)ChangeOfCoord_No1) ;
    
  case FORM1P :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1P :
		 (void *)ChangeOfCoord_No123) ;

  case FORM2P :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form2P :
		 (void *)ChangeOfCoord_No123) ;
    
  case VECTOR :
    GetDP_Return((void *)ChangeOfCoord_No123) ;

  case FORM1S :
    GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1S :
		 (void *)ChangeOfCoord_No123) ;

  default :
    Msg(GERROR, "Unknown type of field (%s)",
	Get_StringForDefine(Field_Type, Type_Form)) ; 
    GetDP_Return(NULL) ;
  }
}

/* ------------------------------------------------------------------------ */
/*   C h a n g e O f C o o r d _ X X X                                      */
/* ------------------------------------------------------------------------ */

void  ChangeOfCoord_No1(struct Element * Element,
			double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_No1");

  vBFx[0] = vBFu[0] ;

  GetDP_End ;
}

void  ChangeOfCoord_No123(struct Element * Element,
			  double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_No123");

  vBFx[0] = vBFu[0] ; vBFx[1] = vBFu[1] ; vBFx[2] = vBFu[2] ;

  GetDP_End ;
}

void  ChangeOfCoord_Form1(struct Element * Element,
			  double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form1");

  vBFx[0] = vBFu[0] * Element->InvJac.c11 + vBFu[1] * Element->InvJac.c12
          + vBFu[2] * Element->InvJac.c13 ;
  vBFx[1] = vBFu[0] * Element->InvJac.c21 + vBFu[1] * Element->InvJac.c22
          + vBFu[2] * Element->InvJac.c23 ;
  vBFx[2] = vBFu[0] * Element->InvJac.c31 + vBFu[1] * Element->InvJac.c32
          + vBFu[2] * Element->InvJac.c33 ;

  GetDP_End ;
}

void  ChangeOfCoord_Form2(struct Element * Element,
			  double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form2");

  if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form2'");

  vBFx[0] = (vBFu[0] * Element->Jac.c11 + vBFu[1] * Element->Jac.c21
	     + vBFu[2] * Element->Jac.c31) / Element->DetJac ;
  vBFx[1] = (vBFu[0] * Element->Jac.c12 + vBFu[1] * Element->Jac.c22
	     + vBFu[2] * Element->Jac.c32) / Element->DetJac ;
  vBFx[2] = (vBFu[0] * Element->Jac.c13 + vBFu[1] * Element->Jac.c23
	     + vBFu[2] * Element->Jac.c33) / Element->DetJac ;

  GetDP_End ;
}

void  ChangeOfCoord_Form3(struct Element * Element,
			  double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form3");

  if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form3'");

  vBFx[0] = vBFu[0] / Element->DetJac ;

  GetDP_End ;
}


/* Form1P, 2P, 1S : valid in 2D only ! */

void  ChangeOfCoord_Form1P(struct Element * Element,
			   double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form1P");

  vBFx[0] = 0. ;  
  vBFx[1] = 0. ;
  vBFx[2] = vBFu[2] / Element->Jac.c33 ;  /* ... * Element->InvJac.c33 */

  GetDP_End ;
}

void  ChangeOfCoord_Form2P(struct Element * Element,
			   double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form2P");

  if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form2P'");

  vBFx[0] = (vBFu[0] * Element->Jac.c11 + vBFu[1] * Element->Jac.c21)
    / Element->DetJac ;
  vBFx[1] = (vBFu[0] * Element->Jac.c12 + vBFu[1] * Element->Jac.c22)
    / Element->DetJac ;
  vBFx[2] = 0. ;

  GetDP_End ;
}


void  ChangeOfCoord_Form1S(struct Element * Element,
			   double vBFu[], double vBFx[]) {

  GetDP_Begin("ChangeOfCoord_Form1S");

  if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form1S'");

  vBFx[0] = 0. ;  
  vBFx[1] = 0. ;
  vBFx[2] = vBFu[0] / Element->DetJac ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*   C a l _ P r o d u c t X X X                                            */
/* ------------------------------------------------------------------------ */

double  Cal_Product123(double v1[], double v2[]) {
  return  v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] ;
}

double  Cal_Product12 (double v1[], double v2[]) {
  return  v1[0]*v2[0] + v1[1]*v2[1] ;
}

double  Cal_Product3  (double v1[], double v2[]) {
  return  v1[2]*v2[2] ;
}

double  Cal_Product1  (double v1[], double v2[]) {
  return  v1[0]*v2[0] ;
}