f_type.c

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

#define RCSID "$Id: F_Type.c,v 1.17 2006/02/26 00:42:53 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):
 *   Johan Gyselinck
 *   Ruth Sabariego
 */

#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"

/* ------------------------------------------------------------------------ */
/*  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


/* ------------------------------------------------------------------------ */
/*  Create a Complex Value from k Real Values (of same type!)               */
/* ------------------------------------------------------------------------ */

void  F_Complex (F_ARG) {

  /* Warning: this function takes a variable number of arguments (depending on 
     Current.NbrHar). There is no test to check if this number is correct (it 
     just has to be a multiple of 2). */

  int k ;

  GetDP_Begin("F_Complex");

  switch(A->Type){

  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if((A+k)->Type != A->Type) 
	Msg(GERROR, "Mixed type of arguments in function 'Complex'");
      V->Val[MAX_DIM*k] = (A+k)->Val[0] ;
    }
    break;

  case VECTOR :
  case TENSOR_DIAG :
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if((A+k)->Type != A->Type)
	Msg(GERROR, "Mixed type of arguments in function 'Complex'");
      V->Val[MAX_DIM*k  ] = (A+k)->Val[0] ; 
      V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
      V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
    }
    break;

  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if((A+k)->Type != A->Type)
	Msg(GERROR, "Mixed type of arguments in function 'Complex'");
      V->Val[MAX_DIM*k  ] = (A+k)->Val[0] ; 
      V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
      V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
      V->Val[MAX_DIM*k+3] = (A+k)->Val[3] ; 
      V->Val[MAX_DIM*k+4] = (A+k)->Val[4] ;
      V->Val[MAX_DIM*k+5] = (A+k)->Val[5] ;
    }
    break;

  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k++) {
      if((A+k)->Type != A->Type)
	Msg(GERROR, "Mixed type of arguments in function 'Complex'");
      V->Val[MAX_DIM*k  ] = (A+k)->Val[0] ; 
      V->Val[MAX_DIM*k+1] = (A+k)->Val[1] ;
      V->Val[MAX_DIM*k+2] = (A+k)->Val[2] ;
      V->Val[MAX_DIM*k+3] = (A+k)->Val[3] ; 
      V->Val[MAX_DIM*k+4] = (A+k)->Val[4] ;
      V->Val[MAX_DIM*k+5] = (A+k)->Val[5] ;
      V->Val[MAX_DIM*k+6] = (A+k)->Val[6] ; 
      V->Val[MAX_DIM*k+7] = (A+k)->Val[7] ;
      V->Val[MAX_DIM*k+8] = (A+k)->Val[8] ;
    }
    break;

  default :
    Msg(GERROR, "Unknown type of arguments in function 'Complex'");
    break;
  }

  V->Type = A->Type ;

  GetDP_End ;
}


/* ----------------------------------------------------------------------- */
/*  Get the Real Part of a Value                                            */
/* ------------------------------------------------------------------------ */

void  F_Re (F_ARG) {
  int k ;

  GetDP_Begin("F_Re");

  switch (A->Type) {
  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k]     = A->Val[MAX_DIM*k] ;
      V->Val[MAX_DIM*(k+1)] = 0. ;
    }
    break;

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

  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     = A->Val[MAX_DIM*k  ] ;
      V->Val[MAX_DIM*k+1]     = A->Val[MAX_DIM*k+1] ;
      V->Val[MAX_DIM*k+2]     = A->Val[MAX_DIM*k+2] ;
      V->Val[MAX_DIM*k+3]     = A->Val[MAX_DIM*k+3] ;
      V->Val[MAX_DIM*k+4]     = A->Val[MAX_DIM*k+4] ;
      V->Val[MAX_DIM*k+5]     = A->Val[MAX_DIM*k+5] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
      V->Val[MAX_DIM*(k+1)+3] = 0. ;
      V->Val[MAX_DIM*(k+1)+4] = 0. ;
      V->Val[MAX_DIM*(k+1)+5] = 0. ;
    }
    break;

  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     = A->Val[MAX_DIM*k  ] ;
      V->Val[MAX_DIM*k+1]     = A->Val[MAX_DIM*k+1] ;
      V->Val[MAX_DIM*k+2]     = A->Val[MAX_DIM*k+2] ;
      V->Val[MAX_DIM*k+3]     = A->Val[MAX_DIM*k+3] ;
      V->Val[MAX_DIM*k+4]     = A->Val[MAX_DIM*k+4] ;
      V->Val[MAX_DIM*k+5]     = A->Val[MAX_DIM*k+5] ;
      V->Val[MAX_DIM*k+6]     = A->Val[MAX_DIM*k+6] ;
      V->Val[MAX_DIM*k+7]     = A->Val[MAX_DIM*k+7] ;
      V->Val[MAX_DIM*k+8]     = A->Val[MAX_DIM*k+8] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
      V->Val[MAX_DIM*(k+1)+3] = 0. ;
      V->Val[MAX_DIM*(k+1)+4] = 0. ;
      V->Val[MAX_DIM*(k+1)+5] = 0. ;
      V->Val[MAX_DIM*(k+1)+6] = 0. ;
      V->Val[MAX_DIM*(k+1)+7] = 0. ;
      V->Val[MAX_DIM*(k+1)+8] = 0. ;
    }
    break;

  default :
    Msg(GERROR, "Unknown type of arguments in function 'Re'");
    break;
  }
  
  V->Type = A->Type ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Get the Imaginary Part of a Value                                       */
/* ------------------------------------------------------------------------ */

void  F_Im (F_ARG) {
  int k ;

  GetDP_Begin("F_Im");

  switch (A->Type) {
  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k]     = A->Val[MAX_DIM*(k+1)] ;
      V->Val[MAX_DIM*(k+1)] = 0. ;
    }
    break;

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

  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     = A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*k+1]     = A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*k+2]     = A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*k+3]     = A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*k+4]     = A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*k+5]     = A->Val[MAX_DIM*(k+1)+5] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
      V->Val[MAX_DIM*(k+1)+3] = 0. ;
      V->Val[MAX_DIM*(k+1)+4] = 0. ;
      V->Val[MAX_DIM*(k+1)+5] = 0. ;
    }
    break;

  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     = A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*k+1]     = A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*k+2]     = A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*k+3]     = A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*k+4]     = A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*k+5]     = A->Val[MAX_DIM*(k+1)+5] ;
      V->Val[MAX_DIM*k+6]     = A->Val[MAX_DIM*(k+1)+6] ;
      V->Val[MAX_DIM*k+7]     = A->Val[MAX_DIM*(k+1)+7] ;
      V->Val[MAX_DIM*k+8]     = A->Val[MAX_DIM*(k+1)+8] ;
      V->Val[MAX_DIM*(k+1)  ] = 0. ;
      V->Val[MAX_DIM*(k+1)+1] = 0. ;
      V->Val[MAX_DIM*(k+1)+2] = 0. ;
      V->Val[MAX_DIM*(k+1)+3] = 0. ;
      V->Val[MAX_DIM*(k+1)+4] = 0. ;
      V->Val[MAX_DIM*(k+1)+5] = 0. ;
      V->Val[MAX_DIM*(k+1)+6] = 0. ;
      V->Val[MAX_DIM*(k+1)+7] = 0. ;
      V->Val[MAX_DIM*(k+1)+8] = 0. ;
    }
    break;

  default :
    Msg(GERROR, "Unknown type of arguments in function 'Re'");
    break;
  }
  
  V->Type = A->Type ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  Conjugate                                                               */
/* ------------------------------------------------------------------------ */

void  F_Conj (F_ARG) {
  int k ;

  GetDP_Begin("F_Conj");

  switch (A->Type) {
  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k]     =  A->Val[MAX_DIM*k] ;
      V->Val[MAX_DIM*(k+1)] = -A->Val[MAX_DIM*(k+1)] ;
    }
    break;

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

  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     =  A->Val[MAX_DIM*k  ] ;
      V->Val[MAX_DIM*k+1]     =  A->Val[MAX_DIM*k+1] ;
      V->Val[MAX_DIM*k+2]     =  A->Val[MAX_DIM*k+2] ;
      V->Val[MAX_DIM*k+3]     =  A->Val[MAX_DIM*k+3] ;
      V->Val[MAX_DIM*k+4]     =  A->Val[MAX_DIM*k+4] ;
      V->Val[MAX_DIM*k+5]     =  A->Val[MAX_DIM*k+5] ;
      V->Val[MAX_DIM*(k+1)  ] = -A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*(k+1)+1] = -A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*(k+1)+2] = -A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*(k+1)+3] = -A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*(k+1)+4] = -A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*(k+1)+5] = -A->Val[MAX_DIM*(k+1)+5] ;
    }
    break;

  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      V->Val[MAX_DIM*k  ]     =  A->Val[MAX_DIM*k  ] ;
      V->Val[MAX_DIM*k+1]     =  A->Val[MAX_DIM*k+1] ;
      V->Val[MAX_DIM*k+2]     =  A->Val[MAX_DIM*k+2] ;
      V->Val[MAX_DIM*k+3]     =  A->Val[MAX_DIM*k+3] ;
      V->Val[MAX_DIM*k+4]     =  A->Val[MAX_DIM*k+4] ;
      V->Val[MAX_DIM*k+5]     =  A->Val[MAX_DIM*k+5] ;
      V->Val[MAX_DIM*k+6]     =  A->Val[MAX_DIM*k+6] ;
      V->Val[MAX_DIM*k+7]     =  A->Val[MAX_DIM*k+7] ;
      V->Val[MAX_DIM*k+8]     =  A->Val[MAX_DIM*k+8] ;
      V->Val[MAX_DIM*(k+1)  ] = -A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*(k+1)+1] = -A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*(k+1)+2] = -A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*(k+1)+3] = -A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*(k+1)+4] = -A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*(k+1)+5] = -A->Val[MAX_DIM*(k+1)+5] ;
      V->Val[MAX_DIM*(k+1)+6] = -A->Val[MAX_DIM*(k+1)+6] ;
      V->Val[MAX_DIM*(k+1)+7] = -A->Val[MAX_DIM*(k+1)+7] ;
      V->Val[MAX_DIM*(k+1)+8] = -A->Val[MAX_DIM*(k+1)+8] ;
    }
    break;

  default :
    Msg(GERROR, "Unknown type of arguments in function 'Conj'");
    break;
  }
  
  V->Type = A->Type ;

  GetDP_End ;
}


/* -------------------------------------------------------------------------------- */
/*  Cartesian coordinates (Re,Im) to polar coordinates (Amplitude,phase[Radians])   */
/* -------------------------------------------------------------------------------- */

void  F_Cart2Pol (F_ARG) {
  int k ;
  double Re, Im;

  GetDP_Begin("F_Cart2Pol");

  switch (A->Type) {
  case SCALAR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      Re = A->Val[MAX_DIM*k] ; Im = A->Val[MAX_DIM*(k+1)] ;
      V->Val[MAX_DIM*k] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)] = atan2(Im,Re);   
    }
    break;

  case VECTOR :
  case TENSOR_DIAG :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      Re = A->Val[MAX_DIM*k  ] ; Im = A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*k  ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)  ] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);   
    }
    break;

  case TENSOR_SYM :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      Re = A->Val[MAX_DIM*k  ] ; Im = A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*k  ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)  ] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+3] ; Im = A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*k+3] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+3] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+4] ; Im = A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*k+4] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+4] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+5] ; Im = A->Val[MAX_DIM*(k+1)+5] ;
      V->Val[MAX_DIM*k+5] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+5] = atan2(Im,Re);   
    }
    break;

  case TENSOR :
    for (k = 0 ; k < Current.NbrHar ; k+=2) {
      Re = A->Val[MAX_DIM*k  ] ; Im = A->Val[MAX_DIM*(k+1)  ] ;
      V->Val[MAX_DIM*k  ] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)  ] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+1] ; Im = A->Val[MAX_DIM*(k+1)+1] ;
      V->Val[MAX_DIM*k+1] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+1] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+2] ; Im = A->Val[MAX_DIM*(k+1)+2] ;
      V->Val[MAX_DIM*k+2] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+2] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+3] ; Im = A->Val[MAX_DIM*(k+1)+3] ;
      V->Val[MAX_DIM*k+3] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+3] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+4] ; Im = A->Val[MAX_DIM*(k+1)+4] ;
      V->Val[MAX_DIM*k+4] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+4] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+5] ; Im = A->Val[MAX_DIM*(k+1)+5] ;
      V->Val[MAX_DIM*k+5] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+5] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+6] ; Im = A->Val[MAX_DIM*(k+1)+6] ;
      V->Val[MAX_DIM*k+6] = sqrt(SQU(Re)+SQU(Im)) ; V->Val[MAX_DIM*(k+1)+6] = atan2(Im,Re);   
      Re = A->Val[MAX_DIM*k+7] ; Im = A->Val[MAX_DIM*(k+1)+7] ;