movingband2d.c

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

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

#include "GetDP.h"
#include "ExtendedGroup.h"
#include "Data_DefineE.h"
#include "CurrentData.h"
#include "GeoData.h"
#include "Numeric.h"
#include "Tools.h"

/* ------------------------------------------------------------------------ */
/*              M o v i n g B a n d 2 D                                     */
/* ------------------------------------------------------------------------ */

/*
struct MovingBand2D {
  List_T  * InitialList1, * ExtendedList1, * InitialList2, * ExtendedList2;
  int NbrNodes1, *NumNodes1, NbrNodes2, *NumNodes2 ;
  double *x1, *y1, *z1, *x2, *y2, *z2, Area;
  int  Period2, ntr1, ntr2, Closed1, Closed2;
  int PhysNum, StartNumTr, StartIndexTr ;
  int *b1_p1, *b1_p2, *b1_p3, *b2_p1, *b2_p2, *b2_p3;
} ;
*/

struct ThreeInt { int  Int1, Int2, Int3 ; } ;  

int  fcmp_int3(const void * a, const void * b) {
  return  ((struct ThreeInt *)a)->Int1 - ((struct ThreeInt *)b)->Int1 ;
}

int  fcmp_int32(const void * a, const void * b) {
  return  ((struct ThreeInt *)a)->Int2 - ((struct ThreeInt *)b)->Int2 ;
}

void  Contour_MovingBand2D(List_T * InitialList, List_T ** ExtendedList,
			   int * NbrNodes, int ** NumNodes) {
  Tree_T  * Element_Tr ;
  struct  Geo_Element  * GeoElement;
  struct ThreeInt *ThreeInt, ThreeInt1 ;
  int *Nodes, LeftNode, RightNode, LeftInt, RightInt, i_El, i;

  GetDP_Begin("Contour_MovingBand2D");

  Element_Tr = Tree_Create(sizeof(struct ThreeInt), fcmp_int3) ;

  for (i_El = 0 ; i_El < Geo_GetNbrGeoElements() ; i_El++) {
    GeoElement = Geo_GetGeoElement(i_El) ;
    if (List_Search(InitialList, &GeoElement->Region, fcmp_int)) {
      if (GeoElement->Type != LINE) 
	Msg(GERROR, "MovingBand2D contour contains no-LINE elements") ;
      ThreeInt1.Int1 = i_El; ThreeInt1.Int2 = 0; ThreeInt1.Int3 = 0;
      Tree_Add(Element_Tr, &ThreeInt1) ;
    }
  }
  *ExtendedList = Tree2List(Element_Tr) ; Tree_Delete(Element_Tr) ;

  ThreeInt = (struct ThreeInt *)List_Pointer(*ExtendedList, 0) ;
  Nodes = Geo_GetGeoElement(ThreeInt->Int1)->NumNodes;
  RightNode = Nodes[1] ; RightInt = 0;
  LeftNode = Nodes[0] ; LeftInt = 0;
  ThreeInt->Int2 = 0; ThreeInt->Int3 = 1;

  /* printf("i %d   Int2 %d Int3 %d  node1/2 %d %d \n", 
     i, ThreeInt->Int2, ThreeInt->Int3, Nodes[0], Nodes[1]); */

  for (i = 1 ; i < List_Nbr(*ExtendedList) ; i++) {
    for (i_El = 1 ; i_El < List_Nbr(*ExtendedList) ; i_El++) {
      if ( ((ThreeInt = (struct ThreeInt *)List_Pointer(*ExtendedList, i_El))->Int3) == 0) {	
	Nodes = Geo_GetGeoElement(ThreeInt->Int1)->NumNodes;
	if (Nodes[0] == RightNode) {
	  ThreeInt->Int2 = ++RightInt ; ThreeInt->Int3 =  1 ; RightNode = Nodes[1] ;
	}
	else if (Nodes[1] == RightNode) {	      
	  ThreeInt->Int2 = ++RightInt ; ThreeInt->Int3 = -1 ; RightNode = Nodes[0] ;
	}
	else if (Nodes[0] == LeftNode) {	      
	  ThreeInt->Int2 = --LeftInt  ; ThreeInt->Int3 = -1 ; LeftNode  = Nodes[1] ;
	}
	else if (Nodes[1] == LeftNode) {	      
	  ThreeInt->Int2 = --LeftInt  ; ThreeInt->Int3 =  1 ; LeftNode  = Nodes[0] ;
	}
	if (ThreeInt->Int3) break;
      }
    }
    if (!ThreeInt->Int3) Msg(GERROR, "Moving Band contour is not connected !!") ;
  }

  List_Sort(*ExtendedList, fcmp_int32) ;
  *NbrNodes = List_Nbr(*ExtendedList)+1 ;
  *NumNodes = (int *)Malloc(*NbrNodes*sizeof(int)) ;  

  ThreeInt = (struct ThreeInt *)List_Pointer(*ExtendedList, 0) ;
  if (ThreeInt->Int3 == 1) {
    (*NumNodes)[0] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[0];
    (*NumNodes)[1] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[1];
  }
  else {
    (*NumNodes)[0] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[1];
    (*NumNodes)[1] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[0];
  }
  for (i_El = 1 ; i_El < List_Nbr(*ExtendedList) ; i_El++) {
    ThreeInt = (struct ThreeInt *)List_Pointer(*ExtendedList, i_El) ;
    if (ThreeInt->Int3 == 1) {
      (*NumNodes)[i_El+1] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[1];
    }
    else {
      (*NumNodes)[i_El+1] = (Geo_GetGeoElement(ThreeInt->Int1)->NumNodes)[0];
    }
  }    

  GetDP_End ;
} 




void  Init_MovingBand2D (struct Group * Group_P) {
  struct MovingBand2D * MB ;
  int i, Different_Sense, dummy;
  int Different_Sense_MB2D(int nth1, int nth2, int ntr1, int ntr2, int closed1, int closed2,
			   double x1[], double y1[], double x2[], double y2[]);
  GetDP_Begin("Init_MovingBand2D");

  MB = Group_P->MovingBand2D ;
  if (MB->ExtendedList1) {
    Msg(INFO, "Init_MovingBand has already been done ! So nothing to do !?");
    GetDP_End ;
  }

  /* Msg(INFO, "Init_MovingBand!"); */

  Contour_MovingBand2D(MB->InitialList1, &MB->ExtendedList1, &MB->NbrNodes1, &MB->NumNodes1);
  Contour_MovingBand2D(MB->InitialList2, &MB->ExtendedList2, &MB->NbrNodes2, &MB->NumNodes2);

  MB->x1 = (double *)Malloc(MB->NbrNodes1*sizeof(double)) ;
  MB->y1 = (double *)Malloc(MB->NbrNodes1*sizeof(double)) ;
  MB->z1 = (double *)Malloc(MB->NbrNodes1*sizeof(double)) ;
 
  MB->x2 = (double *)Malloc(MB->NbrNodes2*sizeof(double)) ;
  MB->y2 = (double *)Malloc(MB->NbrNodes2*sizeof(double)) ;
  MB->z2 = (double *)Malloc(MB->NbrNodes2*sizeof(double)) ;
 
  MB->Closed1 = MB->Closed2 = 0;
  if (MB->NumNodes1[0] == MB->NumNodes1[MB->NbrNodes1-1]) MB->Closed1 = 1; 
  if (MB->NumNodes2[0] == MB->NumNodes2[MB->NbrNodes2-1]) MB->Closed2 = 1; 
  
  MB->ntr1 = MB->NbrNodes1-1;
  MB->ntr2 = MB->NbrNodes2-1;
  if (MB->Period2 != 1) {
    if ((MB->NbrNodes2-1)%MB->Period2 !=0)
      Msg(WARNING, "Strange periodicity stuff  (%d %d) ! Do you know what you're doing ?", 
	  MB->NbrNodes2-1, MB->Period2);
    MB->ntr2 = (MB->NbrNodes2-1)/MB->Period2;
  }

  Geo_GetNodesCoordinates (MB->NbrNodes1, MB->NumNodes1, MB->x1, MB->y1, MB->z1) ; 
  Geo_GetNodesCoordinates (MB->NbrNodes2, MB->NumNodes2, MB->x2, MB->y2, MB->z2) ; 

  
  Different_Sense = Different_Sense_MB2D(MB->NbrNodes1, MB->NbrNodes2, MB->ntr1, MB->ntr2, 
  		 MB->Closed1, MB->Closed2, MB->x1, MB->y1, MB->x2, MB->y2); 
  if (Different_Sense) { 
    Msg(DEBUG2,"invertinggggggggggggggggggggggggggggggggggg\n");
    for (i=0 ; i<MB->NbrNodes2/2 ; i++) {
      dummy = MB->NumNodes2[i];
      MB->NumNodes2[i] = MB->NumNodes2[MB->NbrNodes2-1-i];
      MB->NumNodes2[MB->NbrNodes2-1-i] = dummy;
    }
  }

  Msg(DEBUG2,"Moving Band Contour 1 has %d nodes :", MB->NbrNodes1);
  for (i=0 ; i<MB->NbrNodes1 ; i++) Msg(DEBUG2," %d ", MB->NumNodes1[i]);
  if (MB->Closed1) Msg(DEBUG2, "  (closed)\n"); else Msg(DEBUG2, "  (open)\n");

  Msg(DEBUG2, "Moving Band Contour 2 has %d nodes :", MB->NbrNodes2);
  for (i=0 ; i<MB->NbrNodes2 ; i++) Msg(DEBUG2, " %d ", MB->NumNodes2[i]);
  if (MB->Closed2) Msg(DEBUG2, "  (closed, "); else Msg(DEBUG2, "  (open, ");
  Msg(DEBUG2, "periodicity 1/%d, ", MB->Period2);
  if (Different_Sense) Msg(DEBUG2, "inversed sense)\n"); else Msg(DEBUG2, "same sense)\n") ;


  MB->b1_p1 = (int *)Malloc((MB->NbrNodes1-1)*sizeof(int)) ;  
  MB->b1_p2 = (int *)Malloc((MB->NbrNodes1-1)*sizeof(int)) ;  
  MB->b1_p3 = (int *)Malloc((MB->NbrNodes1-1)*sizeof(int)) ;  

  MB->b2_p1 = (int *)Malloc((MB->NbrNodes2-1)*sizeof(int)) ;  
  MB->b2_p2 = (int *)Malloc((MB->NbrNodes2-1)*sizeof(int)) ;  
  MB->b2_p3 = (int *)Malloc((MB->NbrNodes2-1)*sizeof(int)) ;  
  
  MB->StartIndexTr = Geo_GetNbrGeoElements() ;
  MB->StartNumTr = Geo_GetGeoElement(0)->Num ;