treatment_formulation.c

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

#define RCSID "$Id: Treatment_Formulation.c,v 1.19 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):
 *   David Colignon
 *   Johan Gyselinck
 *   Ruth Sabariego
 */

#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Get_DofOfElement.h"
#include "ExtendedGroup.h"
#include "GeoData.h"
#include "DofData.h"
#include "CurrentData.h"
#include "Tools.h"

#include "F_FMM.h"
#include "F_FMM_DTA.h"
#include "Pre_GetDofFMMGroup.h"

extern List_T  * GeoData_L ;

int LastElementforAggreg;

/* ------------------------------------------------------------------------ */
/*  T r e a t m e n t _ F o r m u l a t i o n                               */
/* ------------------------------------------------------------------------ */

void  Treatment_Formulation(struct Formulation * Formulation_P) {

  GetDP_Begin("Treatment_Formulation");

  switch (Formulation_P->Type) {

  case FEMEQUATION :
    Treatment_FemFormulation(Formulation_P) ;
    break ;

  case GLOBALEQUATION :
    Treatment_GlobalFormulation(Formulation_P) ;
    break ;

  default :
    Msg(GERROR, "Unknown type for Formulation '%s'", Formulation_P->Name) ;
    break ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  T r e a t m e n t _ F e m F o r m u l a t i o n                         */
/* ------------------------------------------------------------------------ */

void  Treatment_FemFormulation(struct Formulation * Formulation_P) {

  struct Element           Element ;

  struct QuantityStorage   * QuantityStorage_P0, * QuantityStorage_P ;
  struct QuantityStorage   QuantityStorage_S ;
  struct Dof               DofForNoDof_P [NBR_MAX_HARMONIC] ;
  struct EquationTerm      * EquationTerm_P0   , * EquationTerm_P ;
  struct GlobalQuantity    * GlobalQuantity_P ;

  int                      Nbr_DefineQuantity ;
  struct DefineQuantity    * DefineQuantity_P0 , * DefineQuantity_P ;

  List_T                   * FemLocalTermActive_L ;
  struct FemLocalTermActive* FemLocalTermActive_S ;
  List_T                   * QuantityStorage_L ;

  int    i, j, Nbr_Element, i_Element, Nbr_EquationTerm, i_EquTerm ;
  int    Index_DefineQuantity, 	TraceGroupIndex_DefineQuantity ;
  int    Nbr_UnusedEquation ;

  List_T  * InitialListInIndex_L ;
  int     Nbr_Region, i_Region, Num_Region ;
  char    FileFMM[MAX_FILE_NAME_LENGTH];

  extern struct Group * Generate_Group ;
  extern double ** MH_Moving_Matrix ;

  gMatrix  A ;
  gVector  b ;
  
  int Flag_Only ;

  GetDP_Begin("Treatment_FemFormulation");

  /* --------------------------------------------------------------- */
  /* 0.  Initialization of an active zone (QuantityStorage) for each */
  /*     DefineQuantity                                              */
  /* --------------------------------------------------------------- */

  if (!(Nbr_EquationTerm = List_Nbr(Formulation_P->Equation)))
    Msg(GERROR, "No equation in Formulation '%s'", Formulation_P->Name);
  
  if (!(Nbr_DefineQuantity = List_Nbr(Formulation_P->DefineQuantity)))
    Msg(GERROR, "No Quantity in Formulation '%s'", Formulation_P->Name);

  DefineQuantity_P0 = (struct DefineQuantity*)
    List_Pointer(Formulation_P->DefineQuantity, 0) ;

  QuantityStorage_L = List_Create(Nbr_DefineQuantity,  1, 
				  sizeof (struct QuantityStorage) ) ;

  QuantityStorage_S.NumLastElementForFunctionSpace = 
    QuantityStorage_S.NumLastElementForDofDefinition = 
      QuantityStorage_S.NumLastElementForEquDefinition = 0 ;

  for (i = 0 ; i < Nbr_DefineQuantity ; i++) {
    QuantityStorage_S.DefineQuantity = DefineQuantity_P0 + i ;

    if(QuantityStorage_S.DefineQuantity->Type == INTEGRALQUANTITY &&
       QuantityStorage_S.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){
      QuantityStorage_S.FunctionSpace = NULL ; 
      QuantityStorage_S.TypeQuantity = VECTOR ; /* to change */
    }
    else{
      QuantityStorage_S.FunctionSpace = (struct FunctionSpace*)
	List_Pointer(Problem_S.FunctionSpace,
		     QuantityStorage_S.DefineQuantity->FunctionSpaceIndex) ;
      QuantityStorage_S.TypeQuantity = QuantityStorage_S.FunctionSpace->Type ;
    }

    List_Add(QuantityStorage_L, &QuantityStorage_S) ;
  }

  QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L, 0) ;

  Get_InitDofOfElement(&Element) ;


  /* --------------------------------------------------------------- */
  /* 1.  Initialization of equation terms                            */
  /* --------------------------------------------------------------- */

  EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation, 0) ;
  FemLocalTermActive_L = List_Create(Nbr_EquationTerm,  1,
				     sizeof (struct FemLocalTermActive) ) ;
  
  for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
    List_Add(FemLocalTermActive_L, &FemLocalTermActive_S) ;
    EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
    
    switch(EquationTerm_P->Type){

    case GALERKIN :
      EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)
	List_Pointer(FemLocalTermActive_L, i_EquTerm) ;
      switch (TreatmentStatus) {
      case _PRE :
	Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
	break ;
      case _CAL :
	Cal_InitGalerkinTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
					  &QuantityStorage_S, DofForNoDof_P) ;
	break ;
      }
      break;

    case DERHAM :
      EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)
	List_Pointer(FemLocalTermActive_L, i_EquTerm) ;
      switch (TreatmentStatus) {
      case _PRE :
	Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
	break ;
      case _CAL :
	Cal_InitdeRhamTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
					&QuantityStorage_S, DofForNoDof_P) ;
	break ;
      }
      break;
    
    case GLOBALTERM :
      switch (TreatmentStatus) {
      case _PRE :
	Pre_InitGlobalTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
	break ;
      }
      break;

    case GLOBALEQUATION :
      break ;

    default :
      Msg(GERROR, "Unknown type of equation term") ;
      break ;
    }

  }  

  /* ---------------------------------------------------------- */
  /* PREprocessing FMM: Assigning NumDof to Groups              */
  /*  --------------------------------------------------------- */

  if(TreatmentStatus==_CAL && Flag_FMM == 1){
    Msg(DIRECT, "P r e - P r o c e s s i n g F M M. . .") ;

    Pre_GetDofFMMGroup(Formulation_P, QuantityStorage_P0) ;
    if (!Current.FMM.NbrDir){ Msg(INFO, "No FMM far groups -> Flag_FMM = 0 ") ; Flag_FMM = 0 ; }
 
    strcpy(FileFMM, Name_Generic); strcat(FileFMM, ".fmm");   Print_FMMGroupInfo(FileFMM);

    Msg(DIRECT, "E n d  P r e - P r o c e s s i n g F M M. . .") ; 
  }


  /* ---------------------------------------------------------- */
  /* 2.  Loop on geometrical elements :                         */
  /*     Treatment of eventual GALERKIN / DERAHM terms          */
  /*  --------------------------------------------------------- */
  
  Nbr_Element = Geo_GetNbrGeoElements() ;

  for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {

    if (Generate_Group) {
      Element.Region = Geo_GetGeoElement(i_Element)->Region ;
      while (i_Element < Nbr_Element && 
	     !List_Search(Generate_Group->InitialList, 
			  &Element.Region, fcmp_int) ) {
	i_Element++ ;
	if (i_Element < Nbr_Element) Element.Region = Geo_GetGeoElement(i_Element)->Region ;
      }
      if (i_Element == Nbr_Element) break ;
    }


    Progress(i_Element, Nbr_Element, "") ;

    Element.GeoElement = Geo_GetGeoElement(i_Element) ;
    Element.Num    = Element.GeoElement->Num ;
    Element.Type   = Element.GeoElement->Type ;
    Current.Region = Element.Region = Element.GeoElement->Region ;

    if (Flag_FMM) Element.FMMGroup = Element.GeoElement->FMMGroup ;

    /* ---------------------------- */
    /* 2.1.  Loop on equation terms */
    /* ---------------------------- */
    
    for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
      EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
      
      if (Flag_FMM){
	Current.FMM.Src = EquationTerm_P->Case.LocalTerm.FMMSource ;
	Current.FMM.Obs = EquationTerm_P->Case.LocalTerm.FMMObservation ;
      }
     
      
      
      if (EquationTerm_P->Type == GALERKIN ||
	  EquationTerm_P->Type == DERHAM) {

	/* if the element is in the support of integration of the term */
	/* ----------------------------------------------------------- */

	if (List_Search(((struct Group *)
			 List_Pointer(Problem_S.Group,
				      EquationTerm_P->Case.
				      LocalTerm.InIndex))->InitialList, 
			&Element.Region, fcmp_int ) ) {

	  if(Flag_VERBOSE == 10)
	    printf("==> Element #%d, EquationTerm #%d/%d\n",
		   Element.Num, i_EquTerm+1, Nbr_EquationTerm) ;

	  Current.IntegrationSupportIndex = EquationTerm_P->Case.LocalTerm.InIndex ;

	  /* ---------------------------------------------------------- */
	  /* 2.1.1.  Loop on quantities (test fcts and shape functions) */
	  /* ---------------------------------------------------------- */

	  for (i = 0 ; i < EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex ; i++) {

	    Index_DefineQuantity =
	      EquationTerm_P->Case.LocalTerm.Term.QuantityIndexTable[i] ;
	    DefineQuantity_P  = DefineQuantity_P0  + Index_DefineQuantity ;
	    QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;

	    TraceGroupIndex_DefineQuantity =
	      EquationTerm_P->Case.LocalTerm.Term.QuantityTraceGroupIndexTable[i] ;

	    /* Only one analysis for each function space */

	    /*
	     * Attention : l'operateur de trace ne fonctionne que si le champ 
	     * dont on prend la trace n'intervient qu'une seule fois dans le terme.
	     * du a - manque de generalite du code au niveau de la gestion des 
	     *        espaces fonctionnels des fcts tests  pour 'Trace de Dof'
	     *      - et Christophe fatigu