cal_quantity.c

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

#define RCSID "$Id: Cal_Quantity.c,v 1.44 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 "Cal_Quantity.h"
#include "Treatment_Formulation.h"
#include "Get_Geometry.h"
#include "Pos_Search.h"
#include "F_Function.h"
#include "CurrentData.h"
#include "Numeric.h"
#include "Tools.h"

void  Cal_SolidAngle(int Source, struct Element *Element,
		     struct QuantityStorage * QuantityStorage,
		     int Nbr_Dof, int Index, struct Value **Stack);

/* ------------------------------------------------------------------------ */
/*  G e t _ V a l u e O f E x p r e s s i o n                               */
/* ------------------------------------------------------------------------ */

void  Get_ValueOfExpression(struct Expression * Expression_P, 
			    struct QuantityStorage * QuantityStorage_P0, 
			    double u, double v, double w,
			    struct Value * Value) {
  int k ;
  struct ExpressionPerRegion  * ExpressionPerRegion_P ;

  GetDP_Begin("Get_ValueOfExpression");

  switch (Expression_P->Type) {

  case CONSTANT :
    if (Current.NbrHar == 1) {
      Value->Val[0] = Expression_P->Case.Constant ;
    }
    else {
      for (k = 0 ; k < Current.NbrHar ; k += 2) {
	Value->Val[MAX_DIM* k   ] = Expression_P->Case.Constant ;
	Value->Val[MAX_DIM*(k+1)] = 0. ;
      }
    }
    Value->Type = SCALAR ;
    break ;

  case WHOLEQUANTITY :
    Cal_WholeQuantity(Current.Element, QuantityStorage_P0, 
		      Expression_P->Case.WholeQuantity,
		      u,v,w, -1, 0, Value) ;
    break ;

  case PIECEWISEFUNCTION :
    if (Current.Region != Expression_P->Case.PieceWiseFunction.NumLastRegion) {
      if ((ExpressionPerRegion_P = (struct ExpressionPerRegion*)
	   List_PQuery(Expression_P->Case.PieceWiseFunction.ExpressionPerRegion,
		       &Current.Region, fcmp_int))) {
	Expression_P->Case.PieceWiseFunction.NumLastRegion = Current.Region ;
	Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion =
	  (struct Expression*)List_Pointer(Problem_S.Expression,
					   ExpressionPerRegion_P->ExpressionIndex) ;
      }
      else {
	if(Current.Region == NO_REGION)
	  Msg(GERROR, "Function '%s' undefined in expressions without support",
	      Expression_P->Name);
	else
	  Msg(GERROR, "Function '%s' undefined in Region %d",
	      Expression_P->Name, Current.Region);
      }
    }
    Get_ValueOfExpression
      (Expression_P->Case.PieceWiseFunction.ExpressionForLastRegion, 
       QuantityStorage_P0, u, v, w, Value) ;
    break ;

  default :
    Msg(GERROR, "Unknown type (%d) of Expression (%s)",
	Expression_P->Type, Expression_P->Name) ;  
    break;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  G e t _ V a l u e O f E x p r e s s i o n B y I n d e x                 */
/* ------------------------------------------------------------------------ */

void  Get_ValueOfExpressionByIndex(int Index_Expression,
				   struct QuantityStorage * QuantityStorage_P0, 
				   double u, double v, double w,
				   struct Value * Value) {

  GetDP_Begin("Get_ValueOfExpressionByIndex");

  Get_ValueOfExpression
    ((struct Expression *)List_Pointer(Problem_S.Expression, Index_Expression),
     QuantityStorage_P0,  u, v, w, Value) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  C a l _ W h o l e Q u a n t i t y                                       */
/* ------------------------------------------------------------------------ */

#define MAX_REGISTER_SIZE   100

#define CAST3V    void(*)(struct Value*, struct Value*, struct Value*)
#define CAST1V    void(*)(struct Value*)
#define CASTF2V   void(*)(struct Function*, struct Value*, struct Value*)

/* There can be at max one "Dof{op qty}" per WholeQuantity, but as
   many {op qty} as you want. 

   Note that Stack[8][MAX_STACK_SIZE] should actually be
   Stack[NBR_MAX_BASISFUNCTION][MAX_STACK_SIZE]. But this tends to
   overflow the stack when we don't use USE_GLOBAL_STACK. Also, 8 is
   OK since the 'multi' feature is only used for SolidAngle
   computations at the moment.

   A better solution would be to build a single stack (instead of 8
   stacks), where a Value could be multiple. But this requires to
   change the way we deal with function arguments. */

static struct Value ValueSaved[MAX_REGISTER_SIZE] ;  

void Cal_WholeQuantity(struct Element * Element,
		       struct QuantityStorage * QuantityStorage_P0,
		       List_T * WholeQuantity_L,
		       double u, double v, double w, 
		       int DofIndexInWholeQuantity, 
		       int Nbr_Dof, struct Value DofValue[]) {

  static int Flag_WarningMissSolForDt = 0 ;
  static int Flag_WarningMissSolForTime_ntime = 0 ;
  static int Flag_InfoForTime_ntime = 0 ;

  int     i_WQ, j, k, Flag_True, Index, DofIndex, Multi[MAX_STACK_SIZE] ;
  int     Save_NbrHar, Save_Region, Type_Dimension, ntime ;
  double  Save_Time, Save_TimeImag, Save_TimeStep, X, Y, Z, Order ;

  struct WholeQuantity   *WholeQuantity_P0, *WholeQuantity_P ;
  struct DofData         *Save_DofData ;
  struct Solution        *Solution_P0, *Solution_PN ;

#define USE_GLOBAL_STACK

#if defined(USE_GLOBAL_STACK)

  /* Use a single global (static) stack for all quantity evaluations.

     Stack size 
       = MAX_RECURSION * MAX_STACK_SIZE * 8 * sizeof(struct Value)
       = 50 * 40 * 8 * (MAX_DIM * NBR_MAX_HARMONIC * sizeof(double))
       = 50 * 40 * 8 * (9 * 2 * 8)
      ~= 2 Mb

     Beware that for NBR_MAX_HARMONIC=40, the size would grow to
     40Mb... So let's define MAX_RECURSION as follows : */

#if NBR_MAX_HARMONIC <= 10
#define MAX_RECURSION 50
#else
#define MAX_RECURSION 10
#endif

  /* We need MAX_RECURSION sufficiently large for expressions like
     (a?b:(c?d:(e?...))) with all n<MAX_RECURSION first tests
     evaluating to false. This case happens quite often when
     specifying piecewise defined physical characteristics.

     If this poses problem in the future, we could still think about
     reallocating ths stack as it grows. The same holds for
     MAX_STACK_SIZE, of course. */

  static struct Value ***StaticStack;
  static int RecursionIndex = -1, first = 1;
  struct Value **Stack;

#else

  /* Use a local stack: this can quickly overflow the stack on Windows
     and Mac OS X */
  struct Value Stack[8][MAX_STACK_SIZE] ;

#endif

  GetDP_Begin("Cal_WholeQuantity");

#if defined(USE_GLOBAL_STACK)
  if(first){
    StaticStack = (struct Value***)Malloc(MAX_RECURSION*sizeof(struct Value**));
    for(j = 0; j < MAX_RECURSION; j++){
      StaticStack[j] = (struct Value**)Malloc(8*sizeof(struct Value*));
      for(k = 0; k < 8; k++){
	StaticStack[j][k] = (struct Value*)Malloc(MAX_STACK_SIZE*sizeof(struct Value));
      }
    }
    first = 0;
  }
  RecursionIndex++;
  if(RecursionIndex < 0 || RecursionIndex >= MAX_RECURSION) 
    Msg(GERROR, "Recursion problem in Cal_WholeQuantity (%d outside [0,%d])", 
	RecursionIndex, MAX_RECURSION);
  Stack = StaticStack[RecursionIndex];
#endif

  WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ;

  Index     = 0 ;
  DofIndex  = -1 ;

  for (i_WQ = 0 ; i_WQ < List_Nbr(WholeQuantity_L) ; i_WQ++) {

    if(Index >= MAX_STACK_SIZE) Msg(GERROR, "Stack size exceeded (%d)", MAX_STACK_SIZE);

    WholeQuantity_P = WholeQuantity_P0 + i_WQ ;

    switch (WholeQuantity_P->Type) {

    case WQ_OPERATORANDQUANTITY : /* {op qty}  Dof{op qty}  BF{op qty} */
      if (i_WQ != DofIndexInWholeQuantity){ /* Attention!!! || TreatmentStatus == _POS){  */
	Pos_FemInterpolation
	  (Element,
	   QuantityStorage_P0,
	   QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
	   WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
	   WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0,
	   u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
	Multi[Index] = 0 ;
      }
      else {
	DofIndex = Index ;
      }
      Index++ ;  
      break ;

    case WQ_ORDER : /* Order[{qty}] */
      Order = Cal_InterpolationOrder
	(Element, QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index) ;
      for (k = 0 ; k < Current.NbrHar ; k += 2) {
	Stack[0][Index].Val[MAX_DIM* k   ] = Order ;
	Stack[0][Index].Val[MAX_DIM*(k+1)] = 0. ;
      }
      Stack[0][Index].Type = SCALAR ;
      Multi[Index] = 0 ;
      Index++ ;  
      break ;

    case WQ_OPERATORANDQUANTITYEVAL : 
      /* {op qty}[x,y,z], {op qty}[x,y,z,dimension] or {op qty}[ntime] */
      if (i_WQ != DofIndexInWholeQuantity || TreatmentStatus == _POS){
	j = WholeQuantity_P->Case.OperatorAndQuantity.NbrArguments;
	if (j == 3 || j == 4) { 
	  Index -= j ;
	  X = Stack[0][Index  ].Val[0] ;
	  Y = Stack[0][Index+1].Val[0] ;
	  Z = Stack[0][Index+2].Val[0] ;
	  if(j == 4) 
	    Type_Dimension = (int)Stack[0][Index+3].Val[0] ;
	  else
	    Type_Dimension = -1 ;
	  Pos_FemInterpolation
	    (Element,
	     QuantityStorage_P0,
	     QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
	     WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
	     WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, Type_Dimension, 1,
	     u, v, w, X, Y, Z, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
	  Multi[Index] = 0 ;
	  Index++ ;
	} 
	else if (j == 1) { 
	  Index -= j ;
	  ntime = (int)Stack[0][Index].Val[0] ;

	  for (k = 0 ; k < Current.NbrSystem ; k++){
	    Solution_P0 = (struct Solution*)List_Pointer((Current.DofData_P0+k)->Solutions, 0);
	    if(((Current.DofData_P0+k)->CurrentSolution - Solution_P0) >= ntime){ 
	      ((Current.DofData_P0+k)->CurrentSolution) -= ntime ;
	      if (Flag_InfoForTime_ntime != List_Nbr((Current.DofData_P0+k)->Solutions)) {
		Msg(INFO, "Accessing solution from %d time steps ago", ntime);
		Msg(INFO, "  -> System %d/%d: TimeStep = %d, Time = %g + i * %g",
		    k+1, Current.NbrSystem, 
		    (Current.DofData_P0+k)->CurrentSolution->TimeStep,
		    (Current.DofData_P0+k)->CurrentSolution->Time,
		    (Current.DofData_P0+k)->CurrentSolution->TimeImag);
		Flag_InfoForTime_ntime = List_Nbr((Current.DofData_P0+k)->Solutions);
	      }
	    }
	    else {
	      if (!Flag_WarningMissSolForTime_ntime) {
		Msg(WARNING,
		    "Missing solution for time step -%d computation (System #%d/%d)",
		    ntime, k+1, Current.NbrSystem);
		Flag_WarningMissSolForTime_ntime = 1 ;
	      }
	    }
	  }

	  Pos_FemInterpolation
	    (Element,
	     QuantityStorage_P0,
	     QuantityStorage_P0 + WholeQuantity_P->Case.OperatorAndQuantity.Index,
	     WholeQuantity_P->Case.OperatorAndQuantity.TypeQuantity,
	     WholeQuantity_P->Case.OperatorAndQuantity.TypeOperator, -1, 0,
	     u, v, w, 0, 0, 0, Stack[0][Index].Val, &Stack[0][Index].Type, 1) ;
	  
	  Multi[Index] = 0 ;
	  Index++ ;

	  for (k = 0 ; k < Current.NbrSystem ; k++){
	    Solution_PN = (struct Solution*)
	      List_Pointer((Current.DofData_P0+k)->Solutions, 
			   List_Nbr((Current.DofData_P0+k)->Solutions)-1);
	    if((Solution_PN - (Current.DofData_P0+k)->CurrentSolution) >= ntime)
	      ((Current.DofData_P0+k)->CurrentSolution) += ntime ;
	  }
	  
	} 
	else
	  Msg(GERROR, "Explicit (x,y,z,time) evaluation not implemented");
      }
      else{
	Msg(GERROR, "Explicit Dof{} evaluation out of context");
      }
      break ;

    case WQ_TRACE : /* Trace[WholeQuantity, Group] */
      Save_Region = Current.Region ; 

      if(!Element->ElementTrace)
	Msg(GERROR, "The trace operator should act on a discrete Quantity");
      
      Current.Region = Element->ElementTrace->Region ;
      
      if(WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity >= 0){
	Cal_WholeQuantity(Element->ElementTrace, QuantityStorage_P0,
			  WholeQuantity_P->Case.Trace.WholeQuantity,
			  Current.ut, Current.vt, Current.wt, 
			  WholeQuantity_P->Case.Trace.DofIndexInWholeQuantity,
			  Nbr_Dof, DofValue) ;
	DofIndexInWholeQuantity = DofIndex = Index ;
      }
      else{
	Current.x = Current.y = Current.z = 0. ;
	for (j = 0 ; j < Element->GeoElement->NbrNodes ; j++) {
	  Current.x += Element->x[j] * Element->n[j] ;
	  Current.y += Element->y[j] * Element->n[j] ;
	  Current.z += Element->z[j] * Element->n[j] ;
	}