get_geometry.c

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

#define RCSID "$Id: Get_Geometry.c,v 1.36 2006/02/26 00:42:54 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   Patrick Lefevre
 */

#if defined(HAVE_GSL)
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#else
#define NRANSI
#include "nrutil.h"
#endif

#include "GetDP.h"
#include "Get_Geometry.h"
#include "BF_Function.h"
#include "Numeric.h"
#include "Data_DefineE.h"

/* ------------------------------------------------------------------------ */
/*  G e t _ N o d e s C o o r d i n a t e s O f E l e m e n t               */
/* ------------------------------------------------------------------------ */

void  Get_NodesCoordinatesOfElement(struct Element * Element) {
  
  GetDP_Begin("Get_NodesCoordinatesOfElement");

  if (Element->NumLastElementForNodesCoordinates != Element->Num) {
    Element->NumLastElementForNodesCoordinates = Element->Num ;
    Geo_GetNodesCoordinates
      (Element->GeoElement->NbrNodes, Element->GeoElement->NumNodes,
       Element->x, Element->y, Element->z) ;
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e t _ G e o E l e m e n t                                             */
/* ------------------------------------------------------------------------ */

void  Get_BFGeoElement(struct Element * Element, double u, double v, double w) {
  int  i ;

  GetDP_Begin("Get_BFGeoElement");

  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
    BF_Node    (Element, i+1, u, v, w, &(Element->n[i])) ;
    BF_GradNode(Element, i+1, u, v, w, Element->dndu[i]) ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  G e t _ J a c o b i a n F u n c t i o n                                 */
/* ------------------------------------------------------------------------ */

void  * Get_JacobianFunction (int Type_Jacobian, int Type_Element,
			      int * Type_Dimension) {

  GetDP_Begin("Get_JacobianFunction");

  switch (Type_Jacobian) {

  case JACOBIAN_VOL :

    switch (Type_Element) {

    case POINT       : 
      *Type_Dimension = _0D ; GetDP_Return((void *)JacobianVol0D) ;

    case LINE        : case LINE_2 : 
      *Type_Dimension = _1D ; GetDP_Return((void *)JacobianVol1D) ;

    case TRIANGLE    : case TRIANGLE_2   :
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVol2D) ;

    case TETRAHEDRON : case TETRAHEDRON_2 : 
    case HEXAHEDRON  : case HEXAHEDRON_2  :
    case PRISM       : case PRISM_2       :
    case PYRAMID     : case PYRAMID_2     :
      *Type_Dimension = _3D ; GetDP_Return((void *)JacobianVol3D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian Vol for Element Type (%s)", 
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_SPH_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolSphShell2D) ;

    case TETRAHEDRON : case TETRAHEDRON_2 :  
    case HEXAHEDRON  : case HEXAHEDRON_2  :  
    case PRISM       : case PRISM_2       : 
    case PYRAMID     : case PYRAMID_2     : 
      *Type_Dimension = _3D ; GetDP_Return((void *)JacobianVolSphShell3D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolSphShell for Element Type (%s)", 
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_RECT_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolRectShell2D) ;

    case TETRAHEDRON : case TETRAHEDRON_2 :  
    case HEXAHEDRON  : case HEXAHEDRON_2  :  
    case PRISM       : case PRISM_2       : 
    case PYRAMID     : case PYRAMID_2     : 
      *Type_Dimension = _3D ; GetDP_Return((void *)JacobianVolRectShell3D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolRectShell for Element Type (%s)", 
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_PLPD_X :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolPlpdX2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolPlpdX for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_AXI :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxi2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxi for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element)); 
    }
    
  case JACOBIAN_VOL_AXI_SPH_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSphShell2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiSphShell for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_AXI_RECT_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiRectShell2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiRectShell for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_AXI_PLPD_X :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiPlpdX2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiPlpdX for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }


  case JACOBIAN_VOL_AXI_SQU :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSqu2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiSqu for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_AXI_SQU_SPH_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSquSphShell2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiSquSphShell for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_VOL_AXI_SQU_RECT_SHELL :

    switch (Type_Element) {

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSquRectShell2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian VolAxiSquRectShell for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_SUR :

    switch (Type_Element) {

    case POINT :
      *Type_Dimension = _1D ; GetDP_Return((void *)JacobianVol0D) ;

    case LINE        : case LINE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianSur2D) ;

    case TRIANGLE    : case TRIANGLE_2   :  
    case QUADRANGLE  : case QUADRANGLE_2 : 
      *Type_Dimension = _3D ; GetDP_Return((void *)JacobianSur3D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian Sur for element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_SUR_SPH_SHELL :

    switch (Type_Element) {

    case LINE        : case LINE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianSurSphShell2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian SurSphShell for element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_SUR_AXI :

    switch (Type_Element) {

    case LINE        : case LINE_2 : 
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianSurAxi2D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian SurAxi for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  case JACOBIAN_LIN :

    switch (Type_Element) {

    case POINT :
      *Type_Dimension = _2D ; GetDP_Return((void *)JacobianVol0D) ;

    case LINE        : case LINE_2 :
      *Type_Dimension = _3D ; GetDP_Return((void *)JacobianLin3D) ;

    default : 
      Msg(GERROR, "Unknown Jacobian Lin for Element Type (%s)",
	  Get_StringForDefine(Element_Type, Type_Element));
    }

  default :
    Msg(GERROR, "Unknown Jacobian"); GetDP_Return(NULL) ;
  }

}




/* ------------------------------------------------------------------------ */
/*  G e o m e t r i c a l   T r a n s f o r m a t i o n s                   */
/* ------------------------------------------------------------------------ */

double  PlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
  int  i ;
  double  CoorX, CoorY, A, B, R, theta, f ;
  double  DetJac ;

  GetDP_Begin("PlpdX2D");

  CoorX = CoorY = 0. ;
  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
    CoorX += Element->x[i] * Element->n[i] ;
    CoorY += Element->y[i] * Element->n[i] ;
  }

  A = Element->JacobianCase->Para[0] ;  B = Element->JacobianCase->Para[1] ;

  R = CoorX ;

  if ( (R > B+1.e-12*B) || (R < A-1.e-12*A) )
    Msg(GERROR, "Bad parameters for unidirectional transformation Jacobian: "
	       "Rint=%g, Rext=%g, R=%g", A, B, R) ;

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

  f  = (A*(B-A)) / (R*(B-R)) ;
  theta = (B-2.*R) / (B-R) ;

  Jac->c11  = f * (1.- theta) ;  Jac->c12  = 0. ;
  Jac->c21  = 0. ;               Jac->c22  = 1. ;

  DetJac = f*( 1.-theta) ;

  GetDP_Return(DetJac) ;
}

double power(double x, double y) {
  if (y == 1.0) return (x);
  else if (y == 2.0) return (x*x);
  else if (y == 0.5) return (sqrt(x));
  else return (pow(x,y));
}

double  Transformation (int Dim, int Type, struct Element * Element, MATRIX3x3 * Jac) {
  int     i, Axis = 0 ;
  double  X = 0., Y = 0., Z = 0. ;
  double  p = 1., L= 0. ;
  double  Cx = 0., Cy = 0., Cz = 0., A = 0., B = 0., R = 0. ;
  double  theta, XR, YR, ZR, f, dRdx = 0., dRdy = 0., dRdz = 0. ;
  double  DetJac ;

  /*
    A          = interior radius
    B          = exterior radius
    Cx, Cy, Cz = coord of centre
    Axis       = direction of the transformation
    p          = exponant
    1/L        = f(B)
  */

  GetDP_Begin("Transformation");

  for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
    X += Element->x[i] * Element->n[i] ;
    Y += Element->y[i] * Element->n[i] ;
    Z += Element->z[i] * Element->n[i] ;
  }

  if(Element->JacobianCase->NbrParameters >= 2){
    A = Element->JacobianCase->Para[0] ;  
    B = Element->JacobianCase->Para[1] ;
  }
  else
    Msg(GERROR, "Missing interior and/or exterior radius for transformation Jacobian");

  if(Type == JACOBIAN_RECT){
    if(Element->JacobianCase->NbrParameters >= 3)
      Axis = (int)Element->JacobianCase->Para[2] ; 
    if(Element->JacobianCase->NbrParameters >= 6){
      Cx = Element->JacobianCase->Para[3] ;  
      Cy = Element->JacobianCase->Para[4] ; 
      Cz = Element->JacobianCase->Para[5] ;
    }
    if(Element->JacobianCase->NbrParameters >= 7)
      p = Element->JacobianCase->Para[6] ;
    if(Element->JacobianCase->NbrParameters >= 8)
      L = Element->JacobianCase->Para[7] ;
    if(Element->JacobianCase->NbrParameters >= 9){
      Msg(DIRECT, ERROR_STR "Too many parameters for rectangular transformation Jacobian");
      Msg(GERROR, "Valid parameters: Radius1 Radius2 Axis CenterX CenterY CenterZ Power 1/Infinity");
    }
  }
  else if(Type == JACOBIAN_SPH){
    if(Element->JacobianCase->NbrParameters >= 5){
      Cx = Element->JacobianCase->Para[2] ;  
      Cy = Element->JacobianCase->Para[3] ; 
      Cz = Element->JacobianCase->Para[4] ;
    }
    if(Element->JacobianCase->NbrParameters >= 6)
      p = Element->JacobianCase->Para[5] ;
    if(Element->JacobianCase->NbrParameters >= 7)
      L = Element->JacobianCase->Para[6] ;
    if(Element->JacobianCase->NbrParameters >= 8){
      Msg(DIRECT, ERROR_STR "Too many parameters for spherical transformation Jacobian");
      Msg(GERROR, "Valid parameters: Radius1 Radius2 CenterX CenterY CenterZ Power 1/Infinity");
    }
  }
  else
    Msg(GERROR, "Unknown type of transformation Jacobian");

  if(L) B = (B*B-A*A*L)/(B-A*L);

  if(Type == JACOBIAN_SPH){
    R    = sqrt( DSQU(X-Cx) + DSQU(Y-Cy) + DSQU(Z-Cz) ) ;
    dRdx = (X-Cx)/R ; 
    dRdy = (Y-Cy)/R ; 
    dRdz = (Z-Cz)/R ;
  }
  else{
    switch(Axis) {