get_geometry.c

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

    case 1: R = (X-Cx) ; dRdx =1.0 ; dRdy =0.0 ; dRdz =0.0 ; break;
    case 2: R = (Y-Cy) ; dRdx =0.0 ; dRdy =1.0 ; dRdz =0.0 ; break;
    case 3: R = (Z-Cz) ; dRdx =0.0 ; dRdy =0.0 ; dRdz =1.0 ; break;
    default: Msg(GERROR, "Bad axis specification: 1 for X, 2 for Y and 3 for Z");
    }
  }

  if ( (fabs(R) > fabs(B) + 1.e-12*fabs(B)) || 
       (fabs(R) < fabs(A) - 1.e-2*fabs(A)) )
    Msg(GERROR, "Bad parameters for transformation Jacobian: %g not in [%g,%g]", R, A, B) ;

  if (B == R) {
    Jac->c11 = 1. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
    Jac->c21 = 0. ; Jac->c22 = 1. ; Jac->c23 = 0. ;
    Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
    GetDP_Return(1.) ;
  }

  f     = power((A*(B-A))/(R*(B-R)), p) ;
  theta = p * (B-2.*R) / (B-R) ;
  XR    = (X-Cx)/R ;
  YR    = (Y-Cy)/R ;
  ZR    = (Z-Cz)/R ;

  Jac->c11  = f * (1.0 - theta * XR * dRdx) ; 
  Jac->c12  = f * (    - theta * XR * dRdy) ;
  Jac->c13  = f * (    - theta * XR * dRdz) ;
  Jac->c21  = f * (    - theta * YR * dRdx) ;
  Jac->c22  = f * (1.0 - theta * YR * dRdy) ;
  Jac->c23  = f * (    - theta * YR * dRdz) ;
  Jac->c31  = f * (    - theta * ZR * dRdx) ;
  Jac->c32  = f * (    - theta * ZR * dRdy) ;
  Jac->c33  = f * (1.0 - theta * ZR * dRdz) ;

  DetJac = f * f * (1.0 - theta) ;

  GetDP_Return(DetJac) ;

}

/* ------------------------------------------------------------------------ */
/*  J a c o b i a n V o l                                                   */
/* ------------------------------------------------------------------------ */

/* 0D */

double  JacobianVol0D (struct Element * Element, MATRIX3x3 * Jac) {
  
  GetDP_Begin("JacobianVol0D");

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

  GetDP_Return(1.) ;
}

/* 1D */

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

  GetDP_Begin("JacobianVol1D");

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

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

  DetJac = Jac->c11 ;

  GetDP_Return(DetJac) ;
}

/* 2D */

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

  GetDP_Begin("JacobianVol2D");

  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 = 1. ;

  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] ;
  }

  DetJac = Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21 ;

  GetDP_Return(DetJac) ;
}

double  JacobianVolSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolSphShell2D");

  DetJac1 = JacobianVol2D (Element, &Jac1) ;
  DetJac2 = Transformation(_2D, JACOBIAN_SPH, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = 1. ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolRectShell2D");

  DetJac1 = JacobianVol2D (Element, &Jac1) ;
  DetJac2 = Transformation (_2D, JACOBIAN_RECT, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = 1. ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolPlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolPlpdX2D");

  DetJac1 = JacobianVol2D (Element, &Jac1) ;
  DetJac2 = PlpdX2D       (Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = 1. ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

/* 2D Axi (Attention, l'axe doit etre x=z=0) */

double  JacobianVolAxi2D (struct Element * Element, MATRIX3x3 * Jac) {
  int  i ;
  double s = 0., DetJac ;

  GetDP_Begin("JacobianVolAxi2D");

  DetJac = JacobianVol2D(Element, Jac) ;

  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
    s += Element->x[i] * Element->n[i] ;

  /* Warning! For evaluations on the symmetry axis */
  if (s==0.0) {
    for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) 
      s += Element->x[i] ;
    s /= (double)Element->GeoElement->NbrNodes  ;
  }

  Jac->c33 = s ;

  GetDP_Return(DetJac * Jac->c33) ;
}

double  JacobianVolAxiSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolAxiSphShell2D");

  DetJac1 = JacobianVolAxi2D  (Element, &Jac1) ;
  DetJac2 = Transformation    (_2D, JACOBIAN_SPH, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = Jac1.c33 ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolAxiRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolAxiRectShell2D");

  DetJac1 = JacobianVolAxi2D  (Element, &Jac1) ;
  DetJac2 = Transformation    (_2D, JACOBIAN_RECT, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = Jac1.c33 ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolAxiPlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolAxiPlpdX2D");

  DetJac1 = JacobianVolAxi2D (Element, &Jac1) ;
  DetJac2 = PlpdX2D          (Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = Jac1.c33 ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

/* 2D Axi avec transformation quadratique (Attention, l'axe doit etre x=z=0) */

double  JacobianVolAxiSqu2D (struct Element * Element, MATRIX3x3 * Jac) {
  int    i ;
  double s = 0., r, DetJac ;

  GetDP_Begin("JacobianVolAxiSqu2D");

  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 = 1. ;

  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
    s += DSQU(Element->x[i]) * Element->n[i] ;


  /* Warning! For evaluations on the symmetry axis */
  if (s==0.0) {
    for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) 
      s += Element->x[i] * Element->x[i] ;
    s /= (double)Element->GeoElement->NbrNodes  ;
  }

  r = sqrt(s);

  for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
    Jac->c11 += 0.5/r * DSQU(Element->x[i]) * Element->dndu[i][0] ;
    Jac->c21 += 0.5/r * DSQU(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->c33 = r ;

  DetJac = (Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21) * Jac->c33 ;

  GetDP_Return(DetJac) ;
}

double  JacobianVolAxiSquSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolAxiSquSphShell2D");

  DetJac1 = JacobianVolAxiSqu2D(Element, &Jac1) ;
  DetJac2 = Transformation     (_2D, JACOBIAN_SPH, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = Jac1.c33 ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolAxiSquRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolAxiSquRectShell2D");

  DetJac1 = JacobianVolAxiSqu2D(Element, &Jac1) ;
  DetJac2 = Transformation     (_2D, JACOBIAN_RECT, Element, &Jac2) ;

  Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;

                                    Jac->c13 = 0. ;
                                    Jac->c23 = 0. ;
  Jac->c31 = 0. ;  Jac->c32 = 0. ;  Jac->c33 = Jac1.c33 ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

/* 3D */

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

  GetDP_Begin("JacobianVol3D");

  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->c21 += Element->x[i] * Element->dndu[i][1] ;
    Jac->c31 += Element->x[i] * Element->dndu[i][2] ;

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

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

  DetJac = Jac->c11 * Jac->c22 * Jac->c33 + Jac->c13 * Jac->c21 * Jac->c32
         + Jac->c12 * Jac->c23 * Jac->c31 - Jac->c13 * Jac->c22 * Jac->c31
         - Jac->c11 * Jac->c23 * Jac->c32 - Jac->c12 * Jac->c21 * Jac->c33 ;

  GetDP_Return(DetJac) ;
}

double  JacobianVolSphShell3D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolShell3D");

  DetJac1 = JacobianVol3D (Element, &Jac1) ;
  DetJac2 = Transformation(_3D, JACOBIAN_SPH, Element, &Jac2) ;

  Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianVolRectShell3D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianVolRectShell3D");

  DetJac1 = JacobianVol3D (Element, &Jac1) ;
  DetJac2 = Transformation (_3D, JACOBIAN_RECT, Element, &Jac2) ;

  Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

/* ------------------------------------------------------------------------ */
/*  J a c o b i a n S u r                                                   */
/* ------------------------------------------------------------------------ */

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

  GetDP_Begin("JacobianSur2D");

  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 = 1. ;

  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] ;
  }

  DetJac = HYPOT(Jac->c11, Jac->c12) ;

  GetDP_Return(DetJac) ;
}

double  JacobianSurSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianSurSphShell2D");

  DetJac1 = JacobianSur2D (Element, &Jac1) ;
  DetJac2 = Transformation(_2D, JACOBIAN_SPH, Element, &Jac2) ;

  Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

double  JacobianSurRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
  MATRIX3x3  Jac1, Jac2 ;
  double     DetJac1, DetJac2 ;

  GetDP_Begin("JacobianSurRectShell2D");

  DetJac1 = JacobianSur2D (Element, &Jac1) ;
  DetJac2 = Transformation(_2D, JACOBIAN_RECT, Element, &Jac2) ;

  Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;

  GetDP_Return(DetJac1 * DetJac2) ;
}

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

  GetDP_Begin("JacobianSurAxi2D");

  DetJac = JacobianSur2D(Element, Jac) ;

  Jac->c33 = 0. ;
  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
    Jac->c33 += Element->x[i] * Element->n[i] ;

  GetDP_Return(DetJac * Jac->c33) ;
}

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

  GetDP_Begin("JacobianSur3D");

  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. ;