fmm_operations.c

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

#define RCSID "$Id: FMM_Operations.c,v 1.15 2006/02/26 00:42:53 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):
 *   Ruth Sabariego
 */

#include "GetDP.h"
#include "Numeric.h"
#include "Data_Passive.h"
#include "Data_DefineE.h"
#include "CurrentData.h"
#include "Cal_Value.h"
#include "Data_FMM.h"

/* ------------------------------------------------------------------------- */
/*  FMM matrix vector product                                                */
/* ------------------------------------------------------------------------- */

void FMM_MatVectorProd(double *x, double *y){ 

  int  NbrFMMEqu, iFMMEqu ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;

  GetDP_Begin("FMM_MatVectorProd");

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;    
    ((void (*)(struct FMMmat*,double*,double*))FMMmat_P->FctProd)(FMMmat_P, x, y) ;
  }

  GetDP_End;
}



void FMM_NoRenumbering( ){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, *NumEqu_A ;
  int  NbrHar, NumDofi, NumEqui, Flag_Iter = 0 ;

  List_T *NumDof_L, *NumEqu_L ;
  List_T *NumDofi_L, *NumEqui_L ;


  GetDP_Begin("FMM_NoRenumbering");
  
  NbrHar = Current.NbrHar;
 
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;

    if(FMMmat_P->NumDofr == NULL) 
      FMMmat_P->NumDofr = FMMmat_P->NumDof ;
    
    if(FMMmat_P->NumEqur == NULL) 
      FMMmat_P->NumEqur = FMMmat_P->NumEqu ;

    if(NbrHar == 2){
      NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );

      if(FMMmat_P->NumDofi == NULL) 
	FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
      else Flag_Iter = 1 ;
      
      if(!Flag_Iter){
	for (iG = 0 ; iG < NbrGroupSrc ; iG++){
	  List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
	  NbrDof = List_Nbr( NumDof_L);   
	  NumDof_A = (int*)(NumDof_L->array) ;
	  
	  NumDofi_L = List_Create(NbrDof, 1, sizeof(int)) ;
	  
	  for (iDof = 0 ; iDof < NbrDof ; iDof++){
	    NumDof  = NumDof_A[iDof] ;
	    NumDofi =  NumDof+1 ;/* Imaginary part */	  
	    List_Add(NumDofi_L, &NumDofi) ;
	  }
	
	  List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
	}

	if(FMMmat_P->NumEqu !=  FMMmat_P->NumDof){
	  NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
	  FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
      
	  for (iG = 0 ; iG < NbrGroupObs ; iG++){
	    List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
	    NbrEqu = List_Nbr( NumEqu_L);   
	    NumEqu_A = (int*)(NumEqu_L->array) ;

	    NumEqui_L = List_Create(NbrEqu,  1, sizeof(int)) ;

	    for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	      NumEqu = NumEqu_A[iEqu] ;
	      NumEqui = NumEqu+1 ;/* Imaginary part */ 
	    
	      List_Add(NumEqui_L, &NumEqui) ;
	    }
	    List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
	  }
	}
	else
	  FMMmat_P->NumEqui = FMMmat_P->NumDofi ;

      }/* Flag_iter */
    }/* NbrHar */
  }/* iFMMEqu */

 GetDP_End;

}


void FMM_Renumbering( int N, int *permr, int *permp){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;
  int  NbrHar, NumDofri, NumEquri, Flag_Renum = 0 ;

  List_T *NumDof_L, *NumDofr_L, *NumEqu_L , *NumEqur_L;
  List_T *NumDofi_L, *NumEqui_L ;


  GetDP_Begin("FMM_Renumbering");
  
  NbrHar = Current.NbrHar;
 

  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );

    if( FMMmat_P->NumDofr == NULL ) 
      FMMmat_P->NumDofr = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
    else Flag_Renum = 1 ;

    if(NbrHar == 2  &&  FMMmat_P->NumDofi == NULL) 
      FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;

    for (iG = 0 ; iG < NbrGroupSrc ; iG++){
      List_Read(FMMmat_P->NumDof, iG, &NumDof_L ) ;
      NbrDof = List_Nbr( NumDof_L);   
      NumDof_A = (int*)(NumDof_L->array) ;
     
      if(!Flag_Renum)
	NumDofr_L = List_Create(NbrDof,  1, sizeof(int)) ;
      else
	List_Read(FMMmat_P->NumDofr, iG, &NumDofr_L ) ;
      

      if(NbrHar == 2){
	if(!Flag_Renum)
	  NumDofi_L = List_Create(NbrDof, 1, sizeof(int)) ;
        else
	  List_Read(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
      }

      for (iDof = 0 ; iDof < NbrDof ; iDof++){
	NumDof  = NumDof_A[iDof] - 1 ;
	NumDofr = permp[(permr[NumDof]-1) + N] ;
	
	if(!Flag_Renum)
	  List_Add(NumDofr_L, &NumDofr) ;
	else
	  List_Write(NumDofr_L, iDof, &NumDofr) ;

	if(NbrHar == 2){
	  NumDof++ ; /* Imaginary part */
	  NumDofri = permp[(permr[NumDof]-1) + N] ;
	  if(!Flag_Renum)
	    List_Add(NumDofi_L, &NumDofri) ;
	  else 
	    List_Write(NumDofi_L, iDof, &NumDofri) ;
	}
      }

      if(!Flag_Renum)
	List_Add(FMMmat_P->NumDofr, &NumDofr_L ) ;
      else 
	List_Write(FMMmat_P->NumDofr, iG, &NumDofr_L ) ;

      if(NbrHar == 2){
	if(!Flag_Renum)
	  List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
	else 
	  List_Write(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
      }
    }

    Flag_Renum = 0 ;

    if(FMMmat_P->NumEqu !=  FMMmat_P->NumDof){
      NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
 
      if( FMMmat_P->NumEqur == NULL ) 
	FMMmat_P->NumEqur = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
      else Flag_Renum = 1 ;

      if (NbrHar == 2 && FMMmat_P->NumEqui == NULL)
	FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
      
      for (iG = 0 ; iG < NbrGroupObs ; iG++){
	List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
	NbrEqu = List_Nbr( NumEqu_L);   
	NumEqu_A = (int*)(NumEqu_L->array) ;

	if(!Flag_Renum)
	  NumEqur_L = List_Create(NbrEqu,  1, sizeof(int)) ;
	else
	  List_Read(FMMmat_P->NumEqur, iG, &NumEqur_L ) ;

	if(NbrHar == 2){
	  if(!Flag_Renum)
	    NumEqui_L = List_Create(NbrEqu,  1, sizeof(int)) ;
	  else
	    List_Read(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
	}

	for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	  NumEqu = NumEqu_A[iEqu] - 1 ;
	  NumEqur  = permp[(permr[NumEqu]-1) + N] ;

	  if(!Flag_Renum)
	    List_Add(NumEqur_L, &NumEqur) ;
	  else
	    List_Write(NumEqur_L, iEqu, &NumEqur) ;

	  if(NbrHar == 2){
	    NumEqu++ ; /* Imaginary part */
	    NumEquri = permp[(permr[NumEqu]-1) + N] ;
	    if(!Flag_Renum)
	      List_Add(NumEqui_L, &NumEquri) ;
	    else 
	      List_Write(NumEqui_L, iEqu, &NumEquri) ;
	  }
	}

	if(!Flag_Renum)
	  List_Add(FMMmat_P->NumEqur, &NumEqur_L ) ;
	else
	  List_Write(FMMmat_P->NumEqur, iG, &NumEqur_L ) ;

	if(NbrHar == 2){
	  if(!Flag_Renum)
	    List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
	  else
	    List_Write(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
	}
      }
    }
    else{ 
      FMMmat_P->NumEqur = FMMmat_P->NumDofr ;
      if (NbrHar == 2) FMMmat_P->NumEqui = FMMmat_P->NumDofi ;
    }

    Flag_Renum = 0;
  }
  
  GetDP_End;
}




void FMM_RenumberingOLD( int N, int *permr, int *permp){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;
  int  NbrHar, NumDofri, NumEquri, Flag_Renum = 0 ;

  List_T *NumDof_L, *NumEqu_L ;
  List_T *NumDofi_L, *NumEqui_L ;


  GetDP_Begin("FMM_Renumbering");

  /* For using together with InverseRenumbering, the NumDof and NumEqu
     are changed and for a second iteration they must be the original
     ones. The Renumbering is done properly but not the
     InverseRenumbering */

  NbrHar = Current.NbrHar;
 
  NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
  FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;
  
  for(iFMMEqu = 0 ; iFMMEqu < NbrFMMEqu ; iFMMEqu ++ ){
    FMMmat_P = FMMmat_P0 + iFMMEqu ;
    NbrGroupSrc = List_Nbr( FMMmat_P->NumDof );
  
    if(NbrHar == 2  &&  FMMmat_P->NumDofi == NULL) 
      FMMmat_P->NumDofi = List_Create(NbrGroupSrc, 1, sizeof(List_T*)) ;
    else Flag_Renum = 1 ;

    for (iG = 0 ; iG < NbrGroupSrc ; iG++){
      List_Read(FMMmat_P->NumDof, iG, &NumDof_L) ;
      NbrDof = List_Nbr( NumDof_L);   
      NumDof_A = (int*)(NumDof_L->array) ;

      if(NbrHar == 2){
	if(!Flag_Renum)
	  NumDofi_L = List_Create(NbrDof,  1, sizeof(int)) ;
        else
	  List_Read(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
      }

      for (iDof = 0 ; iDof < NbrDof ; iDof++){
	NumDof  = NumDof_A[iDof] - 1 ;
	NumDofr = permp[(permr[NumDof]-1) + N] ;
	List_Write(NumDof_L, iDof, &NumDofr) ;

	if(NbrHar == 2){
	  NumDof++ ; /* Imaginary part */
	  NumDofri = permp[(permr[NumDof]-1) + N] ;
	  if(!Flag_Renum)
	    List_Add(NumDofi_L, &NumDofri) ;
	  else 
	    List_Write(NumDofi_L, iDof, &NumDofri) ;
	}
      }
      List_Write(FMMmat_P->NumDof, iG, &NumDof_L ) ;
      if(NbrHar == 2){
	if(!Flag_Renum)
	  List_Add(FMMmat_P->NumDofi, &NumDofi_L ) ;
	else 
	  List_Write(FMMmat_P->NumDofi, iG, &NumDofi_L ) ;
      }
    }
 
    Flag_Renum = 0 ;
    
    if(FMMmat_P->NumEqu !=  FMMmat_P->NumDof){
      NbrGroupObs = List_Nbr( FMMmat_P->NumEqu );
 
      if (NbrHar == 2 && FMMmat_P->NumEqui == NULL)
	FMMmat_P->NumEqui = List_Create(NbrGroupObs, 1, sizeof(List_T*)) ;
      else Flag_Renum = 1 ;

      for (iG = 0 ; iG < NbrGroupObs ; iG++){
	List_Read(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;
	NbrEqu = List_Nbr( NumEqu_L);   
	NumEqu_A = (int*)(NumEqu_L->array) ;

	if(NbrHar == 2){
	  if(!Flag_Renum)
	    NumEqui_L = List_Create(NbrEqu,  1, sizeof(int)) ;
	  else
	    List_Read(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
	}

	for (iEqu = 0 ; iEqu < NbrEqu ; iEqu++){
	  NumEqu = NumEqu_A[iEqu] - 1 ;
	  NumEqur  = permp[(permr[NumEqu]-1) + N] ;
	  List_Write(NumEqu_L, iEqu, &NumEqur) ;

	  if(NbrHar == 2){
	    NumEqu++ ; /* Imaginary part */
	    NumEquri = permp[(permr[NumEqu]-1) + N] ;
	    if(!Flag_Renum)
	      List_Add(NumEqui_L, &NumEquri) ;
	    else 
	      List_Write(NumEqui_L, iEqu, &NumEquri) ;
	  }
	}

	List_Write(FMMmat_P->NumEqu, iG, &NumEqu_L ) ;

	if(NbrHar == 2){
	  if(!Flag_Renum)
	    List_Add(FMMmat_P->NumEqui, &NumEqui_L ) ;
	  else
	    List_Write(FMMmat_P->NumEqui, iG, &NumEqui_L ) ;
	}
      }
    }
    else 
      if(NbrHar == 2) FMMmat_P->NumEqui = FMMmat_P->NumDofi ;
  }
 
  GetDP_End;
}



void FMM_InverseRenumbering( int *rpermr){

  int  NbrFMMEqu, iFMMEqu, iG, NbrGroupSrc, NbrGroupObs ;
  struct FMMmat *FMMmat_P0, *FMMmat_P ;
  int  iDof, NbrDof, NumDof, NumDofr, *NumDof_A ;
  int  iEqu, NbrEqu, NumEqu, NumEqur, *NumEqu_A ;