cal_galerkintermoffemequation.c

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

#define RCSID "$Id: Cal_GalerkinTermOfFemEquation.c,v 1.25 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
 *   Ruth Sabariego
 */

#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "Get_DofOfElement.h"
#include "Get_Geometry.h"
#include "Data_DefineE.h"
#include "GeoData.h"
#include "CurrentData.h"
#include "Tools.h"
#include "F_FMM.h"


void  Cal_InitGalerkinTermOfFemEquation_MHJacNL(struct EquationTerm  * EquationTerm_P) ;
void  Cal_GalerkinTermOfFemEquation_MHJacNL(struct Element          * Element,
					    struct EquationTerm     * EquationTerm_P,
					    struct QuantityStorage  * QuantityStorage_P0) ;
  

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

void  Cal_InitGalerkinTermOfFemEquation(struct EquationTerm     * EquationTerm_P,
					struct QuantityStorage  * QuantityStorage_P0,
					struct QuantityStorage  * QuantityStorageNoDof,
					struct Dof              * DofForNoDof_P) {
  struct FemLocalTermActive  * FI ;
  extern int MH_Moving_Matrix_simple, MH_Moving_Matrix_probe, MH_Moving_Matrix_separate; 

  GetDP_Begin("Cal_InitGalerkinTermOfFemEquation");

  FI = EquationTerm_P->Case.LocalTerm.Active ;

  FI->QuantityStorageEqu_P = QuantityStorage_P0 +
    EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;

  Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu,
			FI->QuantityStorageEqu_P, &FI->Type_FormEqu) ;

  if (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof >= 0) {
    FI->QuantityStorageDof_P = QuantityStorage_P0 +
      EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof ;
    FI->Type_DefineQuantityDof = FI->QuantityStorageDof_P->DefineQuantity->Type ;
  }
  else {
    FI->QuantityStorageDof_P = QuantityStorageNoDof ;
    FI->Type_DefineQuantityDof = NODOF ;
    FI->DofForNoDof_P = DofForNoDof_P ;
    Dof_InitDofForNoDof(DofForNoDof_P, Current.NbrHar) ;
    QuantityStorageNoDof->BasisFunction[0].Dof = DofForNoDof_P ;
  }

  /* Warning: not correct if nonsymmetrical tensor in expression */
  FI->SymmetricalMatrix =
    (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ==
     EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof) &&
    (EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu ==
     EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof) ;

  if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
    FI->SymmetricalMatrix = 0 ;

  if (FI->SymmetricalMatrix) {
    FI->Type_FormDof = FI->Type_FormEqu ;
  }
  else {
    switch (FI->Type_DefineQuantityDof) {
    case LOCALQUANTITY :
      Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
			    FI->QuantityStorageDof_P, &FI->Type_FormDof) ;
      break ;
    case INTEGRALQUANTITY :
      if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof != NOOP){
	Msg(GERROR, "No operator can act on an Integral Quantity");
      }
      FI->Type_FormDof = VECTOR ; /* we don't know the type a priori */
      FI->IntegralQuantityActive.IntegrationCase_L =
	((struct IntegrationMethod *)
	 List_Pointer(Problem_S.IntegrationMethod,
		      FI->QuantityStorageDof_P->DefineQuantity->
		      IntegralQuantity.IntegrationMethodIndex)) ->IntegrationCase ;
      FI->IntegralQuantityActive.CriterionIndex =
	((struct IntegrationMethod *)
	 List_Pointer(Problem_S.IntegrationMethod,
		      FI->QuantityStorageDof_P->DefineQuantity->
		      IntegralQuantity.IntegrationMethodIndex)) ->CriterionIndex ;
      FI->IntegralQuantityActive.JacobianCase_L =
	((struct JacobianMethod *)
	 List_Pointer(Problem_S.JacobianMethod,
		      FI->QuantityStorageDof_P->DefineQuantity->
		      IntegralQuantity.JacobianMethodIndex)) ->JacobianCase ;
      break ;
    case NODOF :
      FI->Type_FormDof = SCALAR ;
      break ;
    }
  }

  FI->Type_ValueDof = Get_ValueFromForm(FI->Type_FormDof);

  /*  G e t   I n t e g r a t i o n   M e t h o d   */
  /*  --------------------------------------------  */
  
  if(EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex < 0)
    Msg(GERROR, "Integration method missing in equation term");

  FI->IntegrationCase_L = 
    ((struct IntegrationMethod *)
     List_Pointer(Problem_S.IntegrationMethod,
		  EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))
    ->IntegrationCase ;

  FI->CriterionIndex = 
    ((struct IntegrationMethod *)
     List_Pointer(Problem_S.IntegrationMethod,
		  EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))
    ->CriterionIndex ;


  /*  G e t   J a c o b i a n   M e t h o d   */
  /*  --------------------------------------  */

  if(EquationTerm_P->Case.LocalTerm.JacobianMethodIndex < 0)
    Msg(GERROR, "Jacobian method missing in equation term");
  
  FI->JacobianCase_L =
    ((struct JacobianMethod *)
     List_Pointer(Problem_S.JacobianMethod,
		  EquationTerm_P->Case.LocalTerm.JacobianMethodIndex))
    ->JacobianCase ;

  FI->JacobianCase_P0 =
    (FI->NbrJacobianCase = List_Nbr(FI->JacobianCase_L)) ?
    (struct JacobianCase*)List_Pointer(FI->JacobianCase_L, 0) : NULL ;

  FI->Flag_ChangeCoord = 
    ( FI->SymmetricalMatrix ||
      !( 
	( (FI->Type_FormEqu == FORM0 || FI->Type_FormEqu == FORM0P) &&
	  (FI->Type_FormDof == FORM3 || FI->Type_FormDof == FORM3P) ) ||
	/*
        ( (FI->Type_FormEqu == FORM1 || FI->Type_FormEqu == FORM1P)  &&
	  (FI->Type_FormDof == FORM2 || FI->Type_FormDof == FORM2P) ) ||
	( (FI->Type_FormEqu == FORM2 || FI->Type_FormEqu == FORM2P)  &&
	  (FI->Type_FormDof == FORM1 || FI->Type_FormDof == FORM1P) ) || 
	 */
	( (FI->Type_FormEqu == FORM3 || FI->Type_FormEqu == FORM3P) &&
	  (FI->Type_FormDof == FORM0 || FI->Type_FormDof == FORM0P) ) 
	)
      )
    ||  /* For operators on VECTOR's (To be extended) */
    (FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR)
    ||
    (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) ;

  if (FI->Flag_ChangeCoord){
    FI->Flag_InvJac = ( (FI->Type_FormEqu == FORM1) ||
			(!FI->SymmetricalMatrix && (FI->Type_FormDof == FORM1)) ||
                        /* For operators on VECTOR's (To be extended) */
			(FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR) ||
			(EquationTerm_P->Case.LocalTerm.Term.QuantityIndexPost == 1) ) ;
  }
  
  /*  G e t   C h a n g e O f C o o r d i n a t e s   */
  /*  ----------------------------------------------  */

  FI->xChangeOfCoordinatesEqu = 
    (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormEqu) ;

  if (!FI->SymmetricalMatrix)
    FI->xChangeOfCoordinatesDof =
      (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormDof) ;
  else
    FI->xChangeOfCoordinatesDof =
      (void (*)())Get_ChangeOfCoordinates(0, FI->Type_FormDof) ; /* Used only for transfer */


  /*  G e t   C a l _ P r o d u c t x  */
  /*  -------------------------------- */
  
  switch (FI->Type_FormEqu) {
  case FORM1   : case FORM1S :
  case FORM2   : case FORM2P : case FORM2S :
  case VECTOR  :
    FI->Cal_Productx = (double (*)())Cal_Product123 ; break ;
  case FORM1P  :  
  case VECTORP :
    FI->Cal_Productx = (double (*)())Cal_Product3 ; break ;
  case FORM0   : 
  case FORM3   :  case FORM3P :  
  case SCALAR  :
    FI->Cal_Productx = (double (*)())Cal_Product1 ; break ;
  default      : 
    Msg(GERROR, "Unknown type of Form (%d)", FI->Type_FormEqu);
  }

  /*  G e t   F u n c t i o n _ A s s e m b l e T e r m  */
  /*  -------------------------------------------------  */

  switch (EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative) {
  case NODT_   : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NoDt   ; break;
  case DT_     : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_Dt     ; break;
  case DTDOF_  : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDof  ; break;
  case DTDT_   : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDt   ; break;
  case DTDTDOF_: FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDof; break;
  case JACNL_  : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_JacNL  ; break;
  case NEVERDT_: FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NeverDt; break;
  case DTNL_   : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtNL   ; break;
  default      : Msg(GERROR, "Unknown type of Operator for Galerkin term (%d)", 
		     EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative);
  }


  if (MH_Moving_Matrix_simple) {
    /* Msg (INFO, "AssembleTerm_MH_Moving") ; */
    FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_simple ; 
  }
  if (MH_Moving_Matrix_probe) {
    /* Msg (INFO, "AssembleTerm_MH_Moving") ; */
    FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_probe ; 
  }
  if (MH_Moving_Matrix_separate) {
    /* Msg (INFO, "AssembleTerm_MH_Moving") ; */
    FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_separate ; 
  }


  /*  initialisation of MHJacNL-term (nonlinear multi-harmonics) if necessary */
  Cal_InitGalerkinTermOfFemEquation_MHJacNL(EquationTerm_P);


  /* Full_Matrix */

  if (EquationTerm_P->Case.LocalTerm.Full_Matrix) {
    FI->Full_Matrix = 1;
    FI->FirstElements = List_Create(20, 10, sizeof(struct FirstElement)) ;    
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  C a l _ T e r m O f F e m E q u a t i o n                               */
/* ------------------------------------------------------------------------ */


void  Cal_GalerkinTermOfFemEquation(struct Element          * Element,
				    struct EquationTerm     * EquationTerm_P,
				    struct QuantityStorage  * QuantityStorage_P0) {
  
  struct FemLocalTermActive * FI ;
  struct QuantityStorage    * QuantityStorageEqu_P, * QuantityStorageDof_P ;
  struct IntegrationCase    * IntegrationCase_P ;
  struct Quadrature         * Quadrature_P ;
  struct Value                vBFxDof [NBR_MAX_BASISFUNCTIONS], CoefPhys ;
  struct Value                CanonicExpression_Equ, V1, V2;

  int     Nbr_Equ, Nbr_Dof = 0;
  int     i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
  int     NextElement ;
  
  double  weight, Factor = 1. ;
  double  vBFuEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ;
  double  vBFxEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ;
  double  Ek [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_HARMONIC] ;

  void (*xFunctionBFEqu[NBR_MAX_BASISFUNCTIONS])
    (struct Element * Element, int NumEntity, 
     double u, double v, double w, double Value[] ) ;
  void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
    (struct Element * Element, int NumEntity, 
     double u, double v, double w, double Value[] ) ;
  double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
  void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
  
  extern int Flag_RHS;  
 

  GetDP_Begin("Cal_GalerkinTermOfFemEquation");

  FI = EquationTerm_P->Case.LocalTerm.Active ;


  /* treatment of MHJacNL-term in separate routine */ 
  if (FI->MHJacNL) {  
    /* if only the RHS of the system is to be calculated 
       (in case of adaptive relaxation of the Newton-Raphson scheme)
       the (expensive and redundant) calculation of this term can be skipped */ 
    if (!Flag_RHS)  
      Cal_GalerkinTermOfFemEquation_MHJacNL(Element, EquationTerm_P, QuantityStorage_P0) ;
    GetDP_End ; 
  }

  QuantityStorageEqu_P = FI->QuantityStorageEqu_P ;
  QuantityStorageDof_P = FI->QuantityStorageDof_P ;

  /*  ----------------------------------------------------------------------  */
  /*  G e t   F u n c t i o n V a l u e  f o r  t e s t  f u n c t i o n s    */