cal_value.c

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

#define RCSID "$Id: Cal_Value.c,v 1.26 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):
 *   Johan Gyselinck
 */

#include "GetDP.h"
#include "Data_Passive.h"
#include "Data_DefineE.h"
#include "Cal_Value.h"
#include "CurrentData.h"
#include "Numeric.h"


/* ------------------------------------------------------------------------ */
/*  O p e r a t o r s   o n   V a l u e s                                   */
/* ------------------------------------------------------------------------ */
/* Warning: the pointers V1 and R or V2 and R can be identical. You must    */
/*          use temporary variables in your computations: you can only      */
/*          affect to R at the very last time (when you're sure you will    */
/*          not use V1 or V2 any more).                                     */
/* ------------------------------------------------------------------------ */

static double  a0;
static double  a1     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  a2     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  b1     [NBR_MAX_HARMONIC * MAX_DIM] ; 
static double  b2     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  b3     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  c1     [NBR_MAX_HARMONIC * MAX_DIM] ; 
static double  c2     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  c3     [NBR_MAX_HARMONIC * MAX_DIM] ;
static double  tmp[27][NBR_MAX_HARMONIC * MAX_DIM] ;

static int TENSOR_SYM_MAP[]  = {0,1,2,1,3,4,2,4,5};
static int TENSOR_DIAG_MAP[] = {0,-1,-1,-1,1,-1,-1,-1,2};

void Cal_ComplexProduct(double V1[], double V2[], double P[]) {

  GetDP_Begin("Cal_ComplexProduct");

  P[0]       = V1[0] * V2[0]       - V1[MAX_DIM] * V2[MAX_DIM] ;
  P[MAX_DIM] = V1[0] * V2[MAX_DIM] + V1[MAX_DIM] * V2[0] ;

  GetDP_End ;
}

void Cal_ComplexDivision(double V1[], double V2[], double P[]) {
  double Mod2 ;
  
  GetDP_Begin("Cal_ComplexDivision");

  Mod2       = DSQU(V2[0]) + DSQU(V2[MAX_DIM]) ;
  if(!Mod2) Msg(GERROR, "Division by zero in 'Cal_ComplexDivision'");
  P[0]       = (  V1[0] * V2[0]       + V1[MAX_DIM] * V2[MAX_DIM]) / Mod2 ;
  P[MAX_DIM] = (- V1[0] * V2[MAX_DIM] + V1[MAX_DIM] * V2[0])       / Mod2 ;

  GetDP_End ;
}

void Cal_ComplexInvert(double V1[], double P[]) {
  double Mod2 ;

  GetDP_Begin("Cal_ComplexInvert");

  Mod2       = DSQU(V1[0]) + DSQU(V1[MAX_DIM]) ;
  if(!Mod2) Msg(GERROR, "Division by zero in 'Cal_ComplexInvert'");
  P[0]       =   V1[0]       / Mod2 ;
  P[MAX_DIM] = - V1[MAX_DIM] / Mod2 ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- V1 
   ------------------------------------------------------------------------ */

void  Cal_CopyValue(struct Value * V1, struct Value * R) {
  int  k ;

  GetDP_Begin("Cal_CopyValue");
  
  if (V1->Type == SCALAR) {
    R->Type = SCALAR ;
    for (k = 0 ; k < Current.NbrHar ; k++)
      R->Val[MAX_DIM*k  ] = V1->Val[MAX_DIM*k  ] ;
  }
  else if (V1->Type == VECTOR || V1->Type == TENSOR_DIAG){
    R->Type = V1->Type ;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      R->Val[MAX_DIM*k  ] = V1->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
    }
  }
  else if (V1->Type == TENSOR_SYM){
    R->Type = TENSOR_SYM ;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      R->Val[MAX_DIM*k  ] = V1->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
      R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] ;
      R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] ;
      R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] ;
    }
  }
  else if (V1->Type == TENSOR){
    R->Type = TENSOR ;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      R->Val[MAX_DIM*k  ] = V1->Val[MAX_DIM*k  ] ;
      R->Val[MAX_DIM*k+1] = V1->Val[MAX_DIM*k+1] ;
      R->Val[MAX_DIM*k+2] = V1->Val[MAX_DIM*k+2] ;
      R->Val[MAX_DIM*k+3] = V1->Val[MAX_DIM*k+3] ;
      R->Val[MAX_DIM*k+4] = V1->Val[MAX_DIM*k+4] ;
      R->Val[MAX_DIM*k+5] = V1->Val[MAX_DIM*k+5] ;
      R->Val[MAX_DIM*k+6] = V1->Val[MAX_DIM*k+6] ;
      R->Val[MAX_DIM*k+7] = V1->Val[MAX_DIM*k+7] ;
      R->Val[MAX_DIM*k+8] = V1->Val[MAX_DIM*k+8] ;
    }
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ 
   R <- V1 
   ------------------------------------------------------------------------ */

void  Cal_CopyValueArray(struct Value *V1, struct Value *R, int Nbr_Values){
  int  k, i;

  GetDP_Begin("Cal_CopyValueArray");

  
  if (V1[0].Type == SCALAR) {
    for (i = 0 ; i < Nbr_Values ; i++){
      R[i].Type = SCALAR ;
      for (k = 0 ; k < Current.NbrHar ; k++)
	R[i].Val[MAX_DIM*k  ] = V1[i].Val[MAX_DIM*k  ] ;
    }
  }
  else if (V1[0].Type == VECTOR || V1[0].Type == TENSOR_DIAG){
    for (i = 0 ; i < Nbr_Values ; i++){
      R[i].Type = V1[i].Type ;
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R[i].Val[MAX_DIM*k  ] = V1[i].Val[MAX_DIM*k  ] ;
	R[i].Val[MAX_DIM*k+1] = V1[i].Val[MAX_DIM*k+1] ;
	R[i].Val[MAX_DIM*k+2] = V1[i].Val[MAX_DIM*k+2] ;
      }
    }
  }
  else if (V1[0].Type == TENSOR_SYM){
    for (i = 0 ; i < Nbr_Values ; i++){
      R[i].Type = TENSOR_SYM ;
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R[i].Val[MAX_DIM*k  ] = V1[i].Val[MAX_DIM*k  ] ;
	R[i].Val[MAX_DIM*k+1] = V1[i].Val[MAX_DIM*k+1] ;
	R[i].Val[MAX_DIM*k+2] = V1[i].Val[MAX_DIM*k+2] ;
	R[i].Val[MAX_DIM*k+3] = V1[i].Val[MAX_DIM*k+3] ;
	R[i].Val[MAX_DIM*k+4] = V1[i].Val[MAX_DIM*k+4] ;
	R[i].Val[MAX_DIM*k+5] = V1[i].Val[MAX_DIM*k+5] ;
      }
    }
  }
  else if (V1[0].Type == TENSOR){
    for (i = 0 ; i < Nbr_Values ; i++){
      R[i].Type = TENSOR ;
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R[i].Val[MAX_DIM*k  ] = V1[i].Val[MAX_DIM*k  ] ;
	R[i].Val[MAX_DIM*k+1] = V1[i].Val[MAX_DIM*k+1] ;
	R[i].Val[MAX_DIM*k+2] = V1[i].Val[MAX_DIM*k+2] ;
	R[i].Val[MAX_DIM*k+3] = V1[i].Val[MAX_DIM*k+3] ;
	R[i].Val[MAX_DIM*k+4] = V1[i].Val[MAX_DIM*k+4] ;
	R[i].Val[MAX_DIM*k+5] = V1[i].Val[MAX_DIM*k+5] ;
	R[i].Val[MAX_DIM*k+6] = V1[i].Val[MAX_DIM*k+6] ;
	R[i].Val[MAX_DIM*k+7] = V1[i].Val[MAX_DIM*k+7] ;
	R[i].Val[MAX_DIM*k+8] = V1[i].Val[MAX_DIM*k+8] ;
      }
    }
  }
  
  GetDP_End ;
}


void  Cal_ValueArray2DoubleArray(struct Value *V1, double *R, int Nbr_Values){
  int  k, i;

  GetDP_Begin("Cal_ValueArray2DoubleArray");

  
  if (V1[0].Type == SCALAR) {
    for (i = 0 ; i < Nbr_Values ; i++){
      for (k = 0 ; k < Current.NbrHar ; k++)
	R[Current.NbrHar*i+k] = V1[i].Val[MAX_DIM*k  ] ;
    }
  }
  else if (V1[0].Type == VECTOR){
    for (i = 0 ; i < Nbr_Values ; i++){
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R[3*(Current.NbrHar*i+k)  ] = V1[i].Val[MAX_DIM*k  ] ;
	R[3*(Current.NbrHar*i+k)+1] = V1[i].Val[MAX_DIM*k+1] ;
	R[3*(Current.NbrHar*i+k)+2] = V1[i].Val[MAX_DIM*k+2] ;
      }
    }
  }
  else {
    Msg(GERROR, "Wrong type conversion: %s ",
	Get_StringForDefine(Field_Type, V1[0].Type));
  }

  GetDP_End ;
}



void  Cal_AddValueArray2DoubleArray(struct Value *V1, double *R, int Nbr_Values){
  int  k, i;

  GetDP_Begin("Cal_AddValueArray2DoubleArray");

  
  if (V1[0].Type == SCALAR) {
    for (i = 0 ; i < Nbr_Values ; i++){
      for (k = 0 ; k < Current.NbrHar ; k++){
	R[Current.NbrHar*i+k] += V1[i].Val[MAX_DIM*k  ] ;
      } 
    }
  }
  else if (V1[0].Type == VECTOR){
    for (i = 0 ; i < Nbr_Values ; i++){
      for (k = 0 ; k < Current.NbrHar ; k++) {
	R[3*(Current.NbrHar*i+k)  ] += V1[i].Val[MAX_DIM*k  ] ;
	R[3*(Current.NbrHar*i+k)+1] += V1[i].Val[MAX_DIM*k+1] ;
	R[3*(Current.NbrHar*i+k)+2] += V1[i].Val[MAX_DIM*k+2] ;
      }
    }
  }
  else {
    Msg(GERROR, "Wrong type conversion: %s ",
	Get_StringForDefine(Field_Type, V1[0].Type));
  }
  GetDP_End ;
}



/* ------------------------------------------------------------------------ 
   R <- 0 
   ------------------------------------------------------------------------ */

static double VALUE_ZERO [NBR_MAX_HARMONIC * MAX_DIM] = 
{0.,0.,0., 0.,0.,0.,
 0.,0.,0., 0.,0.,0.,
 0.,0.,0., 0.,0.,0.} ;

void  Cal_ZeroValue(struct Value * R) {

  GetDP_Begin("Cal_ZeroValue");

  memcpy(R->Val, VALUE_ZERO, Current.NbrHar*MAX_DIM*sizeof(double));

  GetDP_End ;
}


/* ------------------------------------------------------------------------ 
   R <- V1 + V2 
   ------------------------------------------------------------------------ */

#define ADD(i)   R->Val[i] = V1->Val[i] + V2->Val[i]
#define CADD(i)  R->Val[MAX_DIM*k+i] = V1->Val[MAX_DIM*k+i] + V2->Val[MAX_DIM*k+i]

void  Cal_AddValue (struct Value * V1, struct Value * V2, struct Value * R) {
  int           i, k;
  int           i1,i2;
  struct Value  A;

  GetDP_Begin("Cal_AddValue");

  if (V1->Type == SCALAR && V2->Type == SCALAR) {
    if (Current.NbrHar == 1) {
      ADD(0);
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	CADD(0);
      }
    }
    R->Type = SCALAR ;
  }

  else if ((V1->Type == VECTOR      && V2->Type == VECTOR) ||
	   (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_DIAG)) {
    if (Current.NbrHar == 1) {
      ADD(0); ADD(1); ADD(2); 
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	CADD(0); CADD(1); CADD(2);
      }
    }
    R->Type = V1->Type;
  }

  else if (V1->Type == TENSOR_SYM && V2->Type == TENSOR_SYM) {
    if (Current.NbrHar == 1) {
      ADD(0); ADD(1); ADD(2); ADD(3); ADD(4); ADD(5);
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	CADD(0); CADD(1); CADD(2); CADD(3); CADD(4); CADD(5);
      }
    }
    R->Type = TENSOR_SYM;
  }

  else if (V1->Type == TENSOR && V2->Type == TENSOR) {
    if (Current.NbrHar == 1) {
      ADD(0); ADD(1); ADD(2); ADD(3); ADD(4); ADD(5); ADD(6); ADD(7); ADD(8);
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k++) {
	CADD(0); CADD(1); CADD(2); CADD(3); CADD(4); CADD(5); CADD(6); CADD(7); CADD(8);
      }
    }
    R->Type = TENSOR;
  }

  else if ((V1->Type == TENSOR     && V2->Type == TENSOR_SYM) ||
	   (V1->Type == TENSOR     && V2->Type == TENSOR_DIAG)||
	   (V1->Type == TENSOR_SYM && V2->Type == TENSOR_DIAG)){
    A.Type = V1->Type;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      for(i=0 ; i<9 ; i++){
	i1 = (V1->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
	i2 = (V2->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
	A.Val[MAX_DIM*k+i1] = 
	  V1->Val[MAX_DIM*k+i1] + ((i2<0)?0.0:V2->Val[MAX_DIM*k+i2]);
      }
    }
    Cal_CopyValue(&A,R);
  }
  
  else if ((V1->Type == TENSOR_SYM  && V2->Type == TENSOR)      ||
	   (V1->Type == TENSOR_DIAG && V2->Type == TENSOR)      ||
	   (V1->Type == TENSOR_DIAG && V2->Type == TENSOR_SYM)){
    A.Type = V2->Type;
    for (k = 0 ; k < Current.NbrHar ; k++) {
      for(i=0 ; i<9 ; i++){
	i1 = (V1->Type==TENSOR_SYM)?TENSOR_SYM_MAP[i]:TENSOR_DIAG_MAP[i];
	i2 = (V2->Type==TENSOR)?i:TENSOR_SYM_MAP[i];
	A.Val[MAX_DIM*k+i2] = 
	  ((i1<0)?0.0:V1->Val[MAX_DIM*k+i1]) + V2->Val[MAX_DIM*k+i2];
      }
    }
    Cal_CopyValue(&A,R);
  }
  
  else {
    Msg(GERROR, "Addition of different quantities: %s + %s",
	Get_StringForDefine(Field_Type, V1->Type),
	Get_StringForDefine(Field_Type, V2->Type));
  }

  GetDP_End ;
}

#undef ADD
#undef CADD



/* ------------------------------------------------------------------------ 
   R <- V1 + V2 
   ------------------------------------------------------------------------ */

#define ADD(i,j)   R[i].Val[j] = V1[i].Val[j] + V2[i].Val[j]
#define CADD(i,j)  R[i].Val[MAX_DIM*k+j] = V1[i].Val[MAX_DIM*k+j] + V2[i].Val[MAX_DIM*k+j]


void  Cal_AddValueArray (struct Value *V1, struct Value *V2, struct Value *R, int Nbr_Values) {
  int           i, ii, k, i1,i2;
  struct Value  A;


  GetDP_Begin("Cal_AddValueArray");


  if (V1[0].Type == SCALAR && V2[0].Type == SCALAR){
    if (Current.NbrHar == 1)
      for(i = 0 ; i < Nbr_Values; i++){
	R[i].Type = SCALAR ;
	ADD(i,0);      
      }
    else
      for(i = 0 ; i < Nbr_Values; i++){
	R[i].Type = SCALAR ;
	for (k = 0 ; k < Current.NbrHar ; k++)
	  CADD(i,0);
      }
  }
  else if ((V1[0].Type == VECTOR      && V2[0].Type == VECTOR) ||
	   (V1[0].Type == TENSOR_DIAG && V2[0].Type == TENSOR_DIAG)){
    if (Current.NbrHar == 1)
      for(i = 0 ; i < Nbr_Values; i++){ 
	R[i].Type = V1[i].Type;
	ADD(i,0); ADD(i,1); ADD(i,2); 
      }
    else
      for(i = 0 ; i < Nbr_Values; i++){