f_multihar.c

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

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

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

struct MH_InitData{
  int Case ;
  int NbrPoints, NbrPointsX;  /* number of samples per smallest and fundamental period resp. */
  struct DofData *DofData ;
  double **H, ***HH ;
  double *t, *w;
};

List_T * MH_InitData_L = NULL;

int fcmp_MH_InitData(const void * a, const void * b) {
  int Result ;
  if ((Result = ((struct MH_InitData *)a)->DofData  - ((struct MH_InitData *)b)->DofData) != 0)
    return Result ;
  if ((Result = ((struct MH_InitData *)a)->Case  - ((struct MH_InitData *)b)->Case) != 0)
    return Result ;
  if (((struct MH_InitData *)a)->Case != 3)
    return ((struct MH_InitData *)a)->NbrPoints  - ((struct MH_InitData *)b)->NbrPoints ;
  else 
    return ((struct MH_InitData *)a)->NbrPointsX - ((struct MH_InitData *)b)->NbrPointsX ;
  
}


int NbrValues_Type (int Type){
  switch (Type){
  case SCALAR : return 1 ;
  case VECTOR : case TENSOR_DIAG : return 3 ; 
  case TENSOR_SYM : return 6 ;
  case TENSOR : return 9 ;
  default :
    Msg(GERROR, "Unknown type in NbrValues_Type");
    return 0;
  }
}

double Product_SCALARxSCALARxSCALAR (double *V1, double *V2, double *V3){
  return V1[0] * V2[0] * V3[0] ; 
}
double Product_VECTORxTENSOR_SYMxVECTOR (double *V1, double *V2, double *V3){
  return V3[0] * (V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2]) + 
         V3[1] * (V1[0] * V2[1] + V1[1] * V2[3] + V1[2] * V2[4]) + 
         V3[2] * (V1[0] * V2[2] + V1[1] * V2[4] + V1[2] * V2[5]) ; 
}
double Product_VECTORxTENSOR_DIAGxVECTOR (double *V1, double *V2, double *V3){
  return V1[0] * V2[0] * V3[0] + V1[1] * V2[1] * V3[1] + V1[2] * V2[2] * V3[2] ;
}
double Product_VECTORxSCALARxVECTOR (double *V1, double *V2, double *V3){
  return V2[0] * (V1[0] * V3[0] + V1[1] * V3[1] + V1[2] * V3[2]) ;
}



void  * Get_RealProductFunction_Type1xType2xType1 (int Type1, int Type2) {

/*  Type1 * Type2 * Type1 */
 
  GetDP_Begin("Get_RealProductFunction_Type1xType2xType1");

  if (Type1 == SCALAR && Type2 == SCALAR) {
    GetDP_Return((void *)Product_SCALARxSCALARxSCALAR) ;
  } else if (Type1 == VECTOR && Type2 == TENSOR_SYM) {
    GetDP_Return((void *)Product_VECTORxTENSOR_SYMxVECTOR) ;
  } else if (Type1 == VECTOR && Type2 == TENSOR_DIAG) {
    GetDP_Return((void *)Product_VECTORxTENSOR_DIAGxVECTOR) ;
  } else if (Type1 == VECTOR && Type1 == SCALAR) {
    GetDP_Return((void *)Product_VECTORxSCALARxVECTOR) ;
  } else {
    Msg(GERROR, "Not allowed types in Get_RealProductFunction_Type1xType2xType1");
    GetDP_Return(NULL) ;
  }

}



/* ------------------------------------------------------------------------ */
/*   MH_Get_InitData                                                        */
/* ------------------------------------------------------------------------ */

/*
  Case = 1 : MHTransform       NbrPoints (samples per smallest period) is given, 
                               NbrPointsX (samples per fundamental period) is derived
  Case = 2 : MHJacNL           NbrPoints given,  NbrPointsX derived
  Case = 3 : HarmonicToTime    NbrPointsX given,  NbrPoints derived
*/

void MH_Get_InitData(int Case, int NbrPoints, int *NbrPointsX_P,
		     double ***H_P, double ****HH_P, double **t_P, double **w_P){
  int NbrHar, iPul, iTime, iHar, jHar, NbrPointsX ;
  double *Val_Pulsation, MaxPuls, MinPuls ;
  double **H, ***HH = 0, *t, *w ;
  struct MH_InitData MH_InitData_S, *MH_InitData_P ;
  
  GetDP_Begin("MH_Get_InitData(");

  MH_InitData_S.Case       = Case;
  MH_InitData_S.DofData    = Current.DofData; 
  MH_InitData_S.NbrPoints  = NbrPoints;
  MH_InitData_S.NbrPointsX = NbrPointsX = *NbrPointsX_P;

  if (MH_InitData_L == NULL)
    MH_InitData_L = List_Create(1, 1, sizeof(struct MH_InitData)) ; 

  if ((MH_InitData_P = (struct MH_InitData *)
      List_PQuery(MH_InitData_L, &MH_InitData_S, fcmp_MH_InitData))){
    *H_P  = MH_InitData_P->H;
    *HH_P = MH_InitData_P->HH;
    *t_P  = MH_InitData_P->t;
    *w_P  = MH_InitData_P->w;
    *NbrPointsX_P = MH_InitData_P->NbrPointsX;
    GetDP_End;
  }
  
  NbrHar = Current.NbrHar;
  Val_Pulsation = Current.DofData->Val_Pulsation;
  MaxPuls = 0. ; MinPuls = 1.e99 ;
  for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
    if (Val_Pulsation[iPul] &&  Val_Pulsation[iPul] < MinPuls)
      MinPuls = Val_Pulsation[iPul] ;
    if (Val_Pulsation[iPul] &&  Val_Pulsation[iPul] > MaxPuls)
      MaxPuls = Val_Pulsation[iPul] ;
  }

  if (Case != 3)
    NbrPointsX = (int)((MaxPuls/MinPuls*(double)NbrPoints));
  else
    NbrPoints  = (int)((MinPuls/MaxPuls*(double)NbrPointsX));

  Msg(INFO, "MH_Get_InitData => NbrHar = %d  NbrPoints = %d|%d  Case = %d",
      NbrHar, NbrPoints, NbrPointsX, Case);

  t = (double *)Malloc(sizeof(double)*NbrPointsX) ;
  if (Case != 3) 
    for (iTime = 0 ; iTime < NbrPointsX ; iTime++) 
      t[iTime] = (double)iTime/(double)NbrPointsX/(MinPuls/TWO_PI) ;
  else
    for (iTime = 0 ; iTime < NbrPointsX ; iTime++) 
      t[iTime] = (double)iTime/((double)NbrPointsX-1.)/(MinPuls/TWO_PI) ;

  /*   
     w = (double *)Malloc(sizeof(double)*NbrPointsX) ;
       for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
         w[iTime] = 2. / (double)NbrPointsX ;
  */

  w = (double *)Malloc(sizeof(double)*NbrHar) ;
  for (iPul = 0 ; iPul < NbrHar/2 ; iPul++)
    if (Val_Pulsation[iPul]){
      w[2*iPul  ] = 2. / (double)NbrPointsX ;
      w[2*iPul+1] = 2. / (double)NbrPointsX ;
    }else{
      w[2*iPul  ] = 1. / (double)NbrPointsX ;
      w[2*iPul+1] = 1. / (double)NbrPointsX ;
    }  

  H = (double **)Malloc(sizeof(double *)*NbrPointsX) ;
  for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
    H[iTime] = (double *)Malloc(sizeof(double)*NbrHar) ;
    for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
      /* if (Val_Pulsation [iPul]){ */
	H[iTime][2*iPul  ] =   cos(Val_Pulsation[iPul] * t[iTime]) ;
	H[iTime][2*iPul+1] = - sin(Val_Pulsation[iPul] * t[iTime]) ;
    }
    /*
      } 
      else {
        H[iTime][2*iPul    ] = 0.5 ;
        H[iTime][2*iPul + 1] = 0 ;
      }
    */
  }

 /*
   for (iHar = 0 ; iHar < NbrHar ; iHar++)
     for (jHar = iHar ; jHar < NbrHar ; jHar++){
       sum = 0.;
       for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
  	 sum += w[iTime] * H[iTime][iHar] * H[iTime][jHar] ;
         sum -= (iHar==jHar)? 1. : 0. ;
         printf("iHar %d jHar %d sum %e\n", iHar, jHar, sum);
     } 
 */
  

  if (Case == 2) {
    HH = (double ***)Malloc(sizeof(double **)*NbrPointsX) ;
    for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
      HH[iTime] = (double **)Malloc(sizeof(double *)*NbrHar) ;
      for (iHar = 0 ; iHar < NbrHar ; iHar++){
	HH[iTime][iHar] = (double *)Malloc(sizeof(double)*NbrHar) ;
	for (jHar = 0 ; jHar < NbrHar ; jHar++){
	  if (Val_Pulsation [iHar/2] && Val_Pulsation [jHar/2] )
	    HH[iTime][iHar][jHar] = 2. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
	  else
	    HH[iTime][iHar][jHar] = 1. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
	}
      }
    }
  }
     
  *H_P  = MH_InitData_S.H  = H;
  *t_P  = MH_InitData_S.t  = t;
  *w_P  = MH_InitData_S.w  = w;
  *HH_P = MH_InitData_S.HH = HH;
  *NbrPointsX_P = MH_InitData_S.NbrPointsX = NbrPointsX;
  List_Add (MH_InitData_L, &MH_InitData_S);
 
  GetDP_End ;
}




/* ------------------------------------------------------------------------ */
/*  F_MHToTime0    (HarmonicToTime in PostOperation)                        */
/* ------------------------------------------------------------------------ */

void  F_MHToTime0 (int init, struct Value * A, struct Value * V, 
		   int iTime, int NbrPointsX, double * TimeMH) {

  static double **H, ***HH, *t, *weight;
  int iVal, nVal, iHar;

  GetDP_Begin("F_MHToTime0");

  if (Current.NbrHar == 1) GetDP_End ;

  if (!init) MH_Get_InitData(3, 0, &NbrPointsX, &H, &HH, &t, &weight);

  *TimeMH = t[iTime] ;
  V->Type = A->Type ;
  nVal = NbrValues_Type (A->Type) ;

  for (iVal = 0 ; iVal < nVal ; iVal++){
    V->Val[iVal] = 0; 
    for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
      V->Val[iVal] += H[iTime][iHar] * A->Val[iHar*MAX_DIM+iVal] ;
  }      

}



/* ---------------------------------------------------------------------- */
/*  F_MHToTime                                                            */
/* ---------------------------------------------------------------------- */

void  F_MHToTime (struct Function * Fct, struct Value * A, struct Value * V) {
  int iHar, iVal, nVal ;
  double time, H[NBR_MAX_HARMONIC];
  struct Value Vtemp;

  GetDP_Begin("F_MHToTime");

  if (Current.NbrHar == 1) Msg(GERROR, "'F_MHToTime' only for Multi-Harmonic stuff") ;
  if((A+1)->Type != SCALAR)			
    Msg(GERROR, "'F_MHToTime' requires second scalar argument (time)");	
  time = (A+1)->Val[0] ;

  for (iHar = 0 ; iHar < Current.NbrHar/2 ; iHar++) {
    /* if (Current.DofData->Val_Pulsation [iHar]){ */
      H[2*iHar  ] =   cos(Current.DofData->Val_Pulsation[iHar] * time) ;
      H[2*iHar+1] = - sin(Current.DofData->Val_Pulsation[iHar] * time) ;
    /*
    }
    else {
      H[2*iHar  ] = 0.5 ;
      H[2*iHar+1] = 0 ;
    }
    */
  }

  nVal = NbrValues_Type (A->Type) ;

  for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
    for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
      Vtemp.Val[iHar*MAX_DIM+iVal] = 0.;
  
  for (iVal = 0 ; iVal < nVal ; iVal++)
    for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
      Vtemp.Val[iVal] += H[iHar] * A->Val[iHar*MAX_DIM+iVal] ;
  
  V->Type = A->Type ;
  for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
    for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
      V->Val[iHar*MAX_DIM+iVal] = Vtemp.Val[iHar*MAX_DIM+iVal] ;

  GetDP_End ;
}




/* ------------------------------------------------------------------------ */
/*  MHTransform                                                             */
/* ------------------------------------------------------------------------ */
		
void MHTransform(struct Element * Element, struct QuantityStorage * QuantityStorage_P0,
		 double u, double v, double w, struct Value *MH_Value,
		 struct Expression * Expression_P, int NbrPoints){
  double **H, ***HH, *t, *weight ;
  int NbrHar;
  struct Value t_Value, MH_Value_Tr;
  int NbrPointsX, iVal, nVal1, nVal2 = 0, iHar, iTime;

  MH_Get_InitData(1, NbrPoints, &NbrPointsX, &H, &HH, &t, &weight);

  nVal1 = NbrValues_Type (MH_Value->Type) ;
  t_Value.Type = MH_Value_Tr.Type = MH_Value->Type;

  NbrHar = Current.NbrHar;   /* save NbrHar */
  Current.NbrHar = 1;        /* evaluation in time domain ! */

  for (iVal = 0 ; iVal < MAX_DIM ; iVal++)  for (iHar = 0 ; iHar < NbrHar ; iHar++)
    MH_Value_Tr.Val[iHar*MAX_DIM+iVal] = 0. ; 
    
  for (iTime = 0 ; iTime < NbrPointsX ; iTime++) {

    for (iVal = 0 ; iVal < nVal1 ; iVal++){ /* evaluation of MH_Value at iTime-th time point */
      t_Value.Val[iVal] = 0; 
      for (iHar = 0 ; iHar < NbrHar ; iHar++)
	t_Value.Val[iVal] += H[iTime][iHar] * MH_Value->Val[iHar*MAX_DIM+iVal] ;
    }      

    /* evaluation of the function */
    Get_ValueOfExpression(Expression_P, QuantityStorage_P0, u, v, w, &t_Value);

    if (!iTime) nVal2 = NbrValues_Type (t_Value.Type) ;

    for (iVal = 0 ; iVal < nVal2 ; iVal++)  for (iHar = 0 ; iHar < NbrHar ; iHar++)
      MH_Value_Tr.Val[iHar*MAX_DIM+iVal] += 
	weight[iHar] * H[iTime][iHar] * t_Value.Val[iVal] ;
    /* weight[iTime] * H[iTime][iHar] * t_Value.Val[iVal] ; */

  } /*  for iTime  */

  for (iVal = 0 ; iVal < nVal2 ; iVal++)  for (iHar = 0 ; iHar < NbrHar ; iHar++)
    MH_Value->Val[iHar*MAX_DIM+iVal] = MH_Value_Tr.Val[iHar*MAX_DIM+iVal] ; 
  MH_Value->Type = t_Value.Type ;

  Current.NbrHar = NbrHar ;
  GetDP_End ;
} 




/* ----------------------------------------------------------------------------------- */
/*  C a l _ I n i t G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L  */
/* ----------------------------------------------------------------------------------- */

void  Cal_InitGalerkinTermOfFemEquation_MHJacNL(struct EquationTerm  * EquationTerm_P) {