f_misc.c

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

#define RCSID "$Id: F_Misc.c,v 1.32 2006/03/02 22:04:12 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):
 *   Tuan Ledinh
 */

#include "GetDP.h"
#include "Data_DefineE.h"
#include "F_Function.h"
#include "GeoData.h"
#include "Get_Geometry.h"
#include "Cal_Value.h" 
#include "CurrentData.h"
#include "Numeric.h"
#include "Tools.h"

#if !defined(HAVE_GSL)
#define NRANSI
#include "nrutil.h"    /* pour Tuan */
#endif

/* ------------------------------------------------------------------------ */
/*  Warning: the pointers A and V can be identical. You must                */
/*           use temporary variables in your computations: you can only     */
/*           affect to V at the very last time (when you're sure you will   */
/*           not use A any more).                                           */
/* ------------------------------------------------------------------------ */

#define F_ARG    struct Function * Fct, struct Value * A, struct Value * V

/* ------------------------------------------------------------------------ */
/*  Printf                                                                  */
/* ------------------------------------------------------------------------ */

void  F_Printf (F_ARG) {

  GetDP_Begin("F_Printf");

  V = A ;  /* Attention: ne sert a rien!?! */
  Print_Value(A) ;
  printf("\n") ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Computes the normal to an element                                       */
/* ------------------------------------------------------------------------ */

void  F_Normal(F_ARG) {
  int  k ;

  GetDP_Begin("F_Normal");

  if(!Current.Element || Current.Element->Num == NO_ELEMENT)
    Msg(GERROR, "No element on which to compute 'F_Normal'");

  Geo_CreateNormal(Current.Element->Type, 
		   Current.Element->x, 
		   Current.Element->y, 
		   Current.Element->z, 
		   V->Val);

  if (Current.NbrHar != 1) {
    V->Val[MAX_DIM] = 0. ;
    V->Val[MAX_DIM+1] = 0. ;
    V->Val[MAX_DIM+2] = 0. ;
    for (k = 2 ; k < Current.NbrHar ; k += 2) {
      V->Val[MAX_DIM* k   ] = V->Val[0] ;
      V->Val[MAX_DIM* k +1] = V->Val[1] ;
      V->Val[MAX_DIM* k +2] = V->Val[2] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
    }
  }
  V->Type = VECTOR ;
 
  GetDP_End ;
}

void  F_NormalSource(F_ARG) {
  int  k ;

  GetDP_Begin("F_NormalSource");

  if(!Current.ElementSource || Current.ElementSource->Num == NO_ELEMENT)
    Msg(GERROR, "No element on which to compute 'F_NormalSource'");

  Geo_CreateNormal(Current.ElementSource->Type, 
		   Current.ElementSource->x, 
		   Current.ElementSource->y, 
		   Current.ElementSource->z, 
		   V->Val);

  if (Current.NbrHar != 1) {
    V->Val[MAX_DIM] = 0. ;
    V->Val[MAX_DIM+1] = 0. ;
    V->Val[MAX_DIM+2] = 0. ;
    for (k = 2 ; k < Current.NbrHar ; k += 2) {
      V->Val[MAX_DIM* k   ] = V->Val[0] ;
      V->Val[MAX_DIM* k +1] = V->Val[1] ;
      V->Val[MAX_DIM* k +2] = V->Val[2] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
    }
  }
  V->Type = VECTOR ;

  GetDP_End ;
}


void  F_Tangent(F_ARG) {
  int  k ;
  double  tx, ty, tz, norm ;

  GetDP_Begin("F_Tangent");

  if(!Current.Element || Current.Element->Num == NO_ELEMENT)
    Msg(GERROR, "No element on which to compute 'F_Tangent'");

  switch (Current.Element->Type) {
   
  case LINE :
    tx = Current.Element->x[1] - Current.Element->x[0] ;
    ty = Current.Element->y[1] - Current.Element->y[0] ;
    tz = Current.Element->z[1] - Current.Element->z[0] ;
    norm = sqrt(DSQU(tx)+DSQU(ty)+DSQU(tz)) ;      
    V->Val[0] = tx/norm ;
    V->Val[1] = ty/norm ;
    V->Val[2] = tz/norm ;
    break ;

  default :
    Msg(GERROR, "Function 'Tangent' only valid for Line Elements");
  }

  if (Current.NbrHar != 1) {
    V->Val[MAX_DIM] = 0. ;
    V->Val[MAX_DIM+1] = 0. ;
    V->Val[MAX_DIM+2] = 0. ;
    for (k = 2 ; k < Current.NbrHar ; k += 2) {
      V->Val[MAX_DIM* k   ] = V->Val[0] ;
      V->Val[MAX_DIM* k +1] = V->Val[1] ;
      V->Val[MAX_DIM* k +2] = V->Val[2] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
    }
  }
  V->Type = VECTOR ;

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  Comparison                                                              */
/* ------------------------------------------------------------------------ */

void  F_CompElementNum (F_ARG) {

  GetDP_Begin("F_CompElementNum");

  if(!Current.Element || !Current.ElementSource)
    Msg(GERROR, "Uninitialized Element in 'F_CompElementNum'");

  V->Type = SCALAR ;
  V->Val[0] = (Current.Element->Num == Current.ElementSource->Num) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Volume                                                                  */
/* ------------------------------------------------------------------------ */

void F_ElementVol (F_ARG) {
  int k;
  double Vol = 0.;
  MATRIX3x3 Jac;

  /* It would be more efficient to compute the volumes directly from
     the node coordinates, but I'm lazy... */

  Get_NodesCoordinatesOfElement(Current.Element) ;
  Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w) ;

  /* we use the most general case (3D embedding) */

  switch(Current.Element->Type){
  case LINE:
    Vol = 2. * JacobianLin3D(Current.Element,&Jac);
    break;
  case TRIANGLE:
    Vol = 0.5 * JacobianSur3D(Current.Element,&Jac) ;
    break;
  case QUADRANGLE:
    Vol = 4. * JacobianSur3D(Current.Element,&Jac) ;
    break;
  case TETRAHEDRON:
    Vol = 1./6. * JacobianVol3D(Current.Element,&Jac) ;
    break;
  case HEXAHEDRON:
    Vol = 8. * JacobianVol3D(Current.Element,&Jac) ;
    break;
  case PRISM:
    Vol = JacobianVol3D(Current.Element,&Jac) ;
    break;
  default :
    Msg(GERROR, "F_ElementVol not implemented for %s",
	Get_StringForDefine(Element_Type, Current.Element->Type));
  }

  V->Type = SCALAR ;
  V->Val[0] = fabs(Vol);

  for (k = 2 ; k < Current.NbrHar ; k += 2) 
    V->Val[MAX_DIM* k] = V->Val[0] ;
  
}

/* ------------------------------------------------------------------------ */
/*  SurfaceArea                                                             */
/* ------------------------------------------------------------------------ */

void  F_SurfaceArea (F_ARG) {

  struct Element  Element ;
  List_T  * InitialList_L;

  int     Index_Region, Nbr_Element, i_Element ;
  double  Val_Surface ;
  double  c11, c21, c12, c22, DetJac ;
  int     i, k ;

  GetDP_Begin("F_SurfaceArea");

  if (!Fct->Active) {
    Fct->Active = (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;

    if (Fct->NbrParameters == 1) {
      Index_Region = (int)(Fct->Para[0]) ;

      InitialList_L = List_Create(1,1,sizeof(int));
      List_Reset(InitialList_L);
      List_Add(InitialList_L,&Index_Region);


      /*
      InitialList_L = ((struct Group *)
		       List_Pointer(Problem_S.Group, Index_Region))->InitialList ;
      */
    }
    else {
      Index_Region = -1 ;
      InitialList_L = NULL ;
    }

    Val_Surface = 0. ;
    Nbr_Element = Geo_GetNbrGeoElements() ;
    for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
      Element.GeoElement = Geo_GetGeoElement(i_Element) ;
      if ((InitialList_L && 
	   List_Search(InitialList_L, &(Element.GeoElement->Region), fcmp_int)) ||
	  (!InitialList_L && Element.GeoElement->Region == Current.Region)) {
	Element.Num    = Element.GeoElement->Num ;
	Element.Type   = Element.GeoElement->Type ;

	if (Element.Type == TRIANGLE || Element.Type == QUADRANGLE) {

	  Get_NodesCoordinatesOfElement(&Element) ;
	  Get_BFGeoElement(&Element, 0., 0., 0.) ;

	  c11 = c21 = c12 = c22 = 0. ;
	  for ( i = 0 ; i < Element.GeoElement->NbrNodes ; i++ ) {
	    c11 += Element.x[i] * Element.dndu[i][0] ;
	    c21 += Element.x[i] * Element.dndu[i][1] ;
	    c12 += Element.y[i] * Element.dndu[i][0] ;
	    c22 += Element.y[i] * Element.dndu[i][1] ;
	  }
	  DetJac = c11 * c22 - c12 * c21 ;

	  if (Element.Type == TRIANGLE)
	    Val_Surface += fabs(DetJac) * 0.5 ;
	  else if (Element.Type == QUADRANGLE)
	    Val_Surface += fabs(DetJac) * 4. ;

	}
	else {
	  Msg(GERROR, "Function 'SurfaceArea' only valid for Triangle or Quandrangle Elements");
	}
      }
    }
    Fct->Active->Case.SurfaceArea.Value = Val_Surface ;
  }

  V->Type = SCALAR ;
  V->Val[0] = Fct->Active->Case.SurfaceArea.Value ;
  V->Val[MAX_DIM] = 0;

  for (k = 2 ; k < Current.NbrHar ; k += 2) { 
    V->Val[MAX_DIM* k] = V->Val[0] ;
    V->Val[MAX_DIM* (k+1)] = 0 ;
  }

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  Transformation of a stiffness matrix (6x6) with 3 given angles          */
/* ------------------------------------------------------------------------ */

#if defined(HAVE_GSL)

/* All the following should really be rewritten and generalized... */

void  F_TransformTensor (F_ARG) {
  Msg(GERROR, "Tuan's routines are not available with the GSL");
}

void  F_TransformPerm (F_ARG) {
  Msg(GERROR, "Tuan's routines are not available with the GSL");
}

void  F_TransformPiezo (F_ARG) {
  Msg(GERROR, "Tuan's routines are not available with the GSL");
}

void  F_TransformPiezoT (F_ARG) {
  Msg(GERROR, "Tuan's routines are not available with the GSL");
}

#else

void ludcmp1(double **a, int n, int *indx, double *d);
void lubksb1(double **a, int n, int *indx, double b[]);

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

void T_Stress(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] = 2*cos(alpha)*sin(alpha);
  T_z[1][0] = SQU(sin(alpha));    
  T_z[1][1] = SQU(cos(alpha));
  T_z[1][3] =-2*cos(alpha)*sin(alpha);
  T_z[2][2] = 1.0;
  T_z[3][0] =-cos(alpha)*sin(alpha);
  T_z[3][1] = 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] =-2*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] = 2*cos(beta)*sin(beta);
  T_y[3][3] = cos(beta);
  T_y[3][4] =-sin(beta);