solvinganalyse.c

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

#define RCSID "$Id: SolvingAnalyse.c,v 1.34 2006/02/28 12:16:29 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 "GeoData.h"
#include "DofData.h"
#include "Init_Problem.h"
#include "Cal_Quantity.h"
#include "Tools.h"
#include "Data_DefineE.h"
#include "Numeric.h"
#include "Get_DofOfElement.h"
#include "CurrentData.h"
#include "Graph.h"
#include "Parser.h"

/* ------------------------------------------------------------------------ */
/*  S o l v i n g A n a l y s e                                             */
/* ------------------------------------------------------------------------ */
/*! Global analyse of a problem */

void  SolvingAnalyse (void) {

  struct Resolution     * Resolution_P     , * Resolution2_P ;
  struct DefineSystem   * DefineSystem_P0  , * DefineSystem2_P0, * DefineSystem_P ;
  struct Solution       * Solution_P       , Solution_S ;
  struct GeoData        * GeoData_P0 ;
  struct DofData        * DofData_P0       , * DofData2_P0 ;
  List_T                * DofData_L        , * DofData2_L ;

  int                     Num_Resolution   ,   Num_Resolution2 ;
  int                     Nbr_DefineSystem ,   Nbr_DefineSystem2 ;
  int                     Nbr_Solution     ;

  struct DofData        * DofData_P ;
  struct Dof            * Dof_P ;
  struct PostOperation  * PostOperation_P[NBR_MAX_POS] ;
  struct PostProcessing * PostProcessing_P[NBR_MAX_POS] ;
  struct PreResolutionInfo  PreResolutionInfo_S ;

  double  d;
  int  i, j ;
  int  Num, Nbr_GeoData = 0;
  int  Nbr_PreResolution, Nbr_OtherSystem ;


  GetDP_Begin("SolvingAnalyse");

  GeoData_L = List_Create( 1, 5, sizeof(struct GeoData)) ;

  /* -------------------- */
  /* Treatment Resolution */
  /* -------------------- */


  if (Flag_PRE) {
    Nbr_OtherSystem = 0 ;
    if (Name_Resolution)
      Num_Resolution = List_ISearchSeq(Problem_S.Resolution, Name_Resolution,
                                       fcmp_Resolution_Name) ;
    else
      Msg(GERROR, "Missing Resolution");
  }
  else if (Flag_PAR || Flag_CAL || Flag_POS) {
    Dof_OpenFile(DOF_PRE, Name_Generic, "r") ;
    Dof_ReadFilePRE0(&Num_Resolution, &Nbr_DefineSystem) ;
    Nbr_OtherSystem = Nbr_DefineSystem ;
  }


  if (Num_Resolution < 0 ||
      Num_Resolution + 1 > List_Nbr(Problem_S.Resolution))
    Msg(GERROR, "Unknown Resolution (%s)", Name_Resolution);

  Treatment_Resolution(Num_Resolution, &Nbr_DefineSystem, &Nbr_OtherSystem,
                       &Resolution_P, &DefineSystem_P0, &DofData_P0,
		       &DofData_L, GeoData_L, &GeoData_P0) ;

  /* -------------- */
  /* Pre-processing */
  /* -------------- */

  TreatmentStatus = _PRE ;

  Msg(DIRECT, "P r e - P r o c e s s i n g . . .") ;

  if (Flag_PRE) {
    
    PreResolutionIndex_L = List_Create(10, 10, sizeof(struct PreResolutionInfo)) ;

    Treatment_Preprocessing(Nbr_DefineSystem, DofData_P0,
			    DefineSystem_P0, GeoData_P0) ;

    Nbr_PreResolution = List_Nbr(PreResolutionIndex_L) ;

    for (i = 0 ; i < Nbr_PreResolution ; i++) {

      Msg(DIRECT, "P r e - R e s o l u t i o n  (%d/%d) . . .",
          i+1, Nbr_PreResolution) ;
      
      List_Read(PreResolutionIndex_L, i, &PreResolutionInfo_S) ;
      Num_Resolution2 = PreResolutionInfo_S.Index ;

      Nbr_OtherSystem = 0 ;
      Treatment_Resolution(Num_Resolution2, &Nbr_DefineSystem2, &Nbr_OtherSystem,
                           &Resolution2_P, &DefineSystem2_P0,
                           &DofData2_P0, &DofData2_L, GeoData_L, &GeoData_P0) ;

      TreatmentStatus = _PRE ;
      Treatment_Preprocessing(Nbr_DefineSystem2, DofData2_P0,
			      DefineSystem2_P0, GeoData_P0) ;

      for (j = 0 ; j < Nbr_DefineSystem2 ; j++)
        Dof_TransferDofTreeToList(DofData2_P0 + j) ;

      Init_DofDataInFunctionSpace(Nbr_DefineSystem2, DofData2_P0) ;

      Current.TypeTime = TIME_STATIC  ; 
      Current.Time = 0. ; Current.TimeImag = 0. ; Current.TimeStep = 0. ;
      Current.RelativeDifference = 0. ; Current.RelaxationFactor = 1. ;

      TreatmentStatus = _CAL ;

      Current.NbrSystem  = Nbr_DefineSystem2 ;  /* Attention: init for Dt[] */
      Current.DofData_P0 = DofData2_P0 ;

      Treatment_Operation(Resolution2_P, Resolution2_P->Operation, 
                          DofData2_P0, GeoData_P0, Resolution_P, DofData_P0) ;
      
      if (PreResolutionInfo_S.Type == PR_GLOBALBASISFUNCTION) {
        for (j = 0 ; j < Nbr_DefineSystem2 ; j++) {
          DofData_P = DofData2_P0 + j ;
          Dof_TransferSolutionToConstraint(DofData_P) ;
          DofData_P->Num += Nbr_DefineSystem ;
          List_Add(DofData_L, DofData_P) ;
        }
        Nbr_DefineSystem = List_Nbr(DofData_L) ;  /* New Value ... */
        DofData_P0 = (struct DofData*)List_Pointer(DofData_L, 0) ;  /* New Value ... */
      }

      Msg(DIRECT, "E n d   P r e - R e s o l u t i o n  (%d/%d)",
          i+1, Nbr_PreResolution) ;
    }
    
    Dof_OpenFile(DOF_PRE, Name_Generic, "w+") ;
    Dof_WriteFilePRE0(Num_Resolution, Resolution_P->Name, Nbr_DefineSystem) ;

    for (i = 0 ; i < Nbr_DefineSystem ; i++){
      Init_PartInDofData(DofData_P0 + i, Flag_PAR) ;
      Dof_WriteFilePRE(DofData_P0 + i) ;
    }
    Nbr_GeoData = List_Nbr(GeoData_L) ;
    for (i = 0 ; i < Nbr_GeoData ; i++)
      Geo_WriteFilePRE(GeoData_P0 + i, Problem_S.Group) ;

    Dof_CloseFile(DOF_PRE) ;

    if (Flag_PAR || Flag_CAL || Flag_POS)
      for (i = 0 ; i < Nbr_DefineSystem ; i++)
        Dof_TransferDofTreeToList(DofData_P0 + i) ;

  }  

  else if (Flag_PAR || Flag_CAL || Flag_POS) {
    
    if(Flag_PAR && (Flag_CAL || Flag_POS))
      Msg(GERROR, "Please partition independently of -cal or -pos");
      
    if(Flag_PAR > NBR_MAX_PARTITION)
      Msg(GERROR, "Too many partitions");

    Msg(LOADING,"Pre-Processing data '%s.pre'", Name_Generic) ;

    for(i = 0 ; i < Nbr_DefineSystem ; i++)
      Dof_ReadFilePRE(DofData_P0 + i) ;
    for(i = 0 ; i < Nbr_OtherSystem ; i++) {
      DofData_P = DofData_P0 + Nbr_DefineSystem + i ;
      Dof_ReadFilePRE(DofData_P) ;
      DefineSystem_P = (struct DefineSystem*)
        List_Pointer((((struct Resolution*)List_Pointer(Problem_S.Resolution,
                                                        DofData_P->ResolutionIndex))
                      ->DefineSystem), DofData_P->SystemIndex) ;
      DofData_P->GeoDataIndex =
	Geo_AddGeoData(GeoData_L, DefineSystem_P->MeshName, Name_MshFile,
		       DefineSystem_P->AdaptName, Name_AdaptFile) ;
      Init_HarInDofData(DefineSystem_P, DofData_P) ;
    }
    Nbr_DefineSystem = List_Nbr(DofData_L) ;  /* New Value ... */

    Nbr_GeoData = List_Nbr(GeoData_L) ;
    
    Geo_ReadFilePRE(GeoData_P0, Nbr_GeoData, Problem_S.Group) ;

    Dof_CloseFile(DOF_PRE) ;
    
  }  

  Msg(SUMMARY, "");
  Msg(DIRECT, "E n d   P r e - P r o c e s s i n g");

  /* ------------- */
  /* Partitionning */
  /* ------------- */

  if (Flag_PAR > 1) {
    TreatmentStatus = _PAR ;
    Msg(DIRECT, "P a r t i t i o n n i n g . . .") ;

    Init_DofDataInFunctionSpace(Nbr_DefineSystem, DofData_P0) ;

    for (i = 0 ; i < Nbr_DefineSystem ; i++) {
      DofData_P = DofData_P0 + i ;
      InitGraph(DofData_P->NbrDof, &DofData_P->Graph);
    }

    Treatment_Preprocessing(Nbr_DefineSystem, DofData_P0, DefineSystem_P0, GeoData_P0) ;

    Dof_OpenFile(DOF_PRE, Name_Generic, "w") ;

    Dof_WriteFilePRE0(Num_Resolution, Resolution_P->Name, Nbr_DefineSystem) ;

    for (i = 0 ; i < Nbr_DefineSystem ; i++){
      DofData_P = DofData_P0 + i ;
      PartitionGraph(DofData_P, Flag_PAR);
      DofData_P->DofTree = NULL ;
      Dof_WriteFilePRE(DofData_P) ;
    }

    for (i = 0 ; i < Nbr_GeoData ; i++) 
      Geo_WriteFilePRE(GeoData_P0 + i, Problem_S.Group) ;

    Dof_CloseFile(DOF_PRE) ;

    Msg(SUMMARY, "");
    Msg(DIRECT, "E n d   P a r t i t i o n n i n g");
  }


  /* ---------- */
  /* Processing */
  /* ---------- */

  if (Flag_CAL) {
    TreatmentStatus = _CAL ;
    Msg(DIRECT, "P r o c e s s i n g . . .") ;

    Init_DofDataInFunctionSpace(Nbr_DefineSystem, DofData_P0) ;

    if(Flag_RESTART) {
      i = 0 ;
      while(Name_ResFile[i]){
	Msg(LOADING, "Processing data '%s'", Name_ResFile[i]) ;
	Dof_OpenFile(DOF_RES, Name_ResFile[i], "rb");
	Dof_ReadFileRES(DofData_L, NULL, -1, &Current.Time, &Current.TimeImag, 
			&Current.TimeStep) ;
	Dof_CloseFile(DOF_RES);
	i++ ;
      }
      Msg(BIGINFO, "Restarting computation (time = %g) s (TimeStep %g)", 
	  Current.Time, Current.TimeStep) ;
    }
    else{
      Current.Time = Current.TimeImag = Current.TimeStep = 0. ;
    }

    Current.NbrHar = 1 ; /* Bug : peut ne pas etre initialise si -cal sans 
			    -pre et evaluation d'expression sans init de systeme avant */
    Current.TypeTime = TIME_STATIC ;
    Current.RelativeDifference = 0. ; 
    Current.RelaxationFactor = 1. ;
    Current.NbrSystem  = Nbr_DefineSystem ;  /* Attention: init for Dt[] */
    Current.DofData_P0 = DofData_P0 ;

    Treatment_Operation(Resolution_P, Resolution_P->Operation, 
                        DofData_P0, GeoData_P0, NULL, NULL) ;

    Msg(SUMMARY, "");
    Msg(DIRECT, "E n d   P r o c e s s i n g");
  }

  /* --------------- */
  /* Post-processing */
  /* --------------- */

  if (Flag_POS) {
    TreatmentStatus = _POS ;

    Flag_FMM = 0 ; Current.FMM.Flag_GF = 0 ;

    Msg(DIRECT, "P o s t - P r o c e s s i n g . . .") ;

    i = 0 ;
    if(Flag_IPOS){
      while(Name_PostProcessing[i]){
	if((Num = List_ISearchSeq(Problem_S.PostProcessing, Name_PostProcessing[i],
				  fcmp_PostProcessing_Name)) < 0)
	  Msg(GERROR, "Unknown PostProcessing (%s)", Name_PostProcessing[i]) ;
	PostProcessing_P[i] = (struct PostProcessing *)
	  List_Pointer(Problem_S.PostProcessing, Num) ;
	PostOperation_P[i] = NULL ;
	i++ ;
      }
    }
    else{
      while(Name_PostOperation[i]){
	if((Num = List_ISearchSeq(Problem_S.PostOperation, Name_PostOperation[i],
				  fcmp_PostOperation_Name)) < 0)
	  Msg(GERROR, "Unknown PostOperation (%s)", Name_PostOperation[i]) ;
	PostOperation_P[i] = (struct PostOperation*)
	  List_Pointer(Problem_S.PostOperation, Num) ;
	PostProcessing_P[i] = (struct PostProcessing *)
	  List_Pointer(Problem_S.PostProcessing, PostOperation_P[i]->PostProcessingIndex) ;
	i++ ;
      }
    }
    PostProcessing_P[i] = NULL ;

    if (!Flag_CAL) {
      i = 0 ;
      while(Name_ResFile[i]){
	Msg(LOADING, "Processing data '%s'", Name_ResFile[i]) ;
	Dof_OpenFile(DOF_RES, Name_ResFile[i], "rb");
	Dof_ReadFileRES(DofData_L, NULL, -1, &d, &d, &d) ;
	Dof_CloseFile(DOF_RES) ;
	i++ ;
      }
    }

    for (i = 0 ; i < Nbr_DefineSystem ; i++) {
      Current.DofData = DofData_P = DofData_P0 + i ;

      for(j=0 ; j<DofData_P->NbrAnyDof ; j++){
        Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, j) ;
        if(Dof_P->Type == DOF_UNKNOWN_INIT){
          Dof_P->Type = DOF_UNKNOWN ;
          LinAlg_ZeroScalar(&Dof_P->Val) ;
        }
      }

      Current.NbrHar = Current.DofData->NbrHar ;
      Nbr_Solution = List_Nbr(DofData_P->Solutions) ;
      if (Nbr_Solution == 0) { /* en cas de pos sans cal, apres calcul de
                                  function de base globale... a reorganiser */
        if (DofData_P->Solutions == NULL)
          DofData_P->Solutions = List_Create( 1, 1, sizeof(struct Solution)) ;
        Solution_S.Time = 0. ;  
	Solution_S.SolutionExist = 0 ;
        List_Add(DofData_P->Solutions, &Solution_S) ;
        Nbr_Solution = 1 ;
      }
      if (!Flag_CAL) { /* Pas necessaire si Flag_CAL */
	for (j = 0 ; j < Nbr_Solution ; j++) {
	  Solution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, j) ;
	  Current.Time = Solution_P->Time ;
	  Current.TimeImag = Solution_P->TimeImag ;
	  Current.TimeStep = 0.;
	  Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ;
	}
      }
      DofData_P->CurrentSolution = (Nbr_Solution)?
	(struct Solution*)List_Pointer(DofData_P->Solutions, 0) : NULL ;
      /* La solution courante est la 1ere. A mieux gerer ? */
    }
    Init_DofDataInFunctionSpace(Nbr_DefineSystem, DofData_P0) ;