cal_derhamtermoffemequation.c

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

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

#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "Get_DofOfElement.h"
#include "Get_Geometry.h"
#include "Get_Cells.h"
#include "Pos_Search.h" /* xyz2uvw */
#include "Numeric.h"
#include "Data_DefineE.h"
#include "GeoData.h"
#include "CurrentData.h"
#include "Tools.h"

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

void  Cal_InitdeRhamTermOfFemEquation(struct EquationTerm     * EquationTerm_P,
				      struct QuantityStorage  * QuantityStorage_P0,
				      struct QuantityStorage  * QuantityStorageNoDof,
				      struct Dof              * DofForNoDof_P) {
  struct FemLocalTermActive  * FI ;

  GetDP_Begin("Cal_InitdeRhamTermOfFemEquation");

  FI = EquationTerm_P->Case.LocalTerm.Active ;

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

  if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu != NOOP)
    Msg(GERROR, "An aperator cannot act on the \"test function\" in a de Rham Map");

  FI->Type_FormEqu = -1 ;

  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 ;
  }


  FI->SymmetricalMatrix = 0 ;

  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, "A differential operator cannot 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   */
  /*  --------------------------------------------  */
  
  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   */
  /*  --------------------------------------  */
  
  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 = 1 ;
  FI->Flag_InvJac = 
    (FI->Type_FormDof == FORM1) ||
    (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 = NULL ;

  FI->xChangeOfCoordinatesDof =
    (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormDof) ;

  FI->Cal_Productx = NULL ;

  /*  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;
  default      : Msg(GERROR, "Unknown type of operator for de Rham (%d)", 
		     EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative);
  }

  GetDP_End ;
}



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


void  Cal_deRhamTermOfFemEquation(struct Element          * Element,
				  struct EquationTerm     * EquationTerm_P,
				  struct QuantityStorage  * QuantityStorage_P0) {

  struct Element              Cells[NBR_MAX_DERHAM_CELLS];
  struct FemLocalTermActive * FI ;
  struct QuantityStorage    * QuantityStorageEqu_P, * QuantityStorageDof_P ;
  struct IntegrationCase    * IntegrationCase_P ;
  struct Quadrature         * Quadrature_P ;
  struct Value                vBFxDof[NBR_MAX_BASISFUNCTIONS] ;
  struct Value                Cell_Vector[NBR_MAX_DERHAM_CELLS];
  struct Value                Cell_Value[NBR_MAX_BASISFUNCTIONS][NBR_MAX_DERHAM_CELLS];
  struct Group              * Group_P ;

  double  u, v, w, weight, Factor ;
  double  u2, v2, w2, x2, y2, z2 ;
  double  Ek [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_HARMONIC] ;   

  int     Nbr_Equ, Nbr_Dof = 0 ;
  int     i, j, k, n, ii, Nbr_IntPoints, i_IntPoint, i_Cell ;
  int     Type_Dimension, Cell_Type_Dimension ;
  int     Nbr_Cells, Cell_RelativeJacobianType ;
  int     Nbr_Ent1 = 0, Nbr_Ent2 = 0, *Dk = NULL ;
  int     Product ; 

  void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
    (struct Element * Element, int NumEntity, 
     double u, double v, double w, double Value[] ) ;
  double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
  double (*Cell_Get_Jacobian)(struct Element*, MATRIX3x3*) = 0;
  void (*Get_IntPoint)(int,int,double*,double*,double*,double*);

  GetDP_Begin("Cal_deRhamTermOfFemEquation");

  FI = EquationTerm_P->Case.LocalTerm.Active ;

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


  /*  ----------------------------------------------------------------------  */
  /*  E q u a t i o n s   t o   b u i l d ?                                   */
  /*  ----------------------------------------------------------------------  */
  
  /* Msg(INFO, "de Rham for Element %d", Element->Num); */

  if (!(Nbr_Equ = QuantityStorageEqu_P->NbrElementaryBasisFunction)) {
    if(Flag_VERBOSE > 2) 
      Msg(WARNING, "No equation to build in Element %d (de Rham)", Element->Num);
    GetDP_End ;
  }

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

  switch (FI->Type_DefineQuantityDof) {
  case LOCALQUANTITY :
    Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ;
    Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
		      EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
		      QuantityStorageDof_P, &FI->Type_FormDof) ;
    break ;
  case INTEGRALQUANTITY :
    Get_InitElementSource(Element,
			  QuantityStorageDof_P->DefineQuantity->IntegralQuantity.InIndex) ;  
    break ;
  case NODOF : 
    Nbr_Dof = 1 ;  
    break ;
  }
  
  /*  -------------------------------------------------------------------  */
  /*  G e t   I n t e g r a t i o n   M e t h o d                          */
  /*  -------------------------------------------------------------------  */

  IntegrationCase_P = Get_IntegrationCase(Element, 
					  FI->IntegrationCase_L, 
					  FI->CriterionIndex);
 
  /*  -------------------------------------------------------------------  */
  /*  G e t   D e R h a m   S u p p o r t                                  */
  /*  -------------------------------------------------------------------  */

  Get_NodesCoordinatesOfElement(Element) ;

  Group_P = (struct Group*) List_Pointer(Problem_S.Group, 
					 EquationTerm_P->Case.LocalTerm.InIndex) ;

  Get_deRhamCells(Element, QuantityStorageEqu_P, Group_P, &Nbr_Cells,
		  Cells, Cell_Vector, &Cell_RelativeJacobianType);

  /* All the cells created from one element are supposed to be of the same type */
  
  /*  -------------------------------------------------------------------  */
  /*  G e t   J a c o b i a n   M e t h o d   f o r   E l e m e n t        */
  /*  -------------------------------------------------------------------  */
  
  i = 0 ;
  while ((i < FI->NbrJacobianCase) &&
	 ((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) &&
	 !List_Search
	 (((struct Group *)List_Pointer(Problem_S.Group, j))
	  ->InitialList, &Element->Region, fcmp_int) )  i++ ;
  
  if (i == FI->NbrJacobianCase)
    Msg(GERROR, "Undefined Jacobian in Region %d", Element->Region);
  
  Element->JacobianCase = FI->JacobianCase_P0 + i ;

  Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
    Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
			 Element->Type, &Type_Dimension) ;

  /*  -------------------------------------------------------------------  */
  /*  G e t   J a c o b i a n   M e t h o d   f o r   C e l l s            */
  /*  -------------------------------------------------------------------  */

  if(Cells[0].Type != POINT){
    Cell_Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))