geodata.c

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

	  GeoData_P->Ymax = MAX(GeoData_P->Ymax, Geo_Node.y);
	  GeoData_P->Zmin = MIN(GeoData_P->Zmin, Geo_Node.z); 
	  GeoData_P->Zmax = MAX(GeoData_P->Zmax, Geo_Node.z); 
	}
      }

      if(GeoData_P->Xmin != GeoData_P->Xmax && 
	 GeoData_P->Ymin != GeoData_P->Ymax && 
	 GeoData_P->Zmin != GeoData_P->Zmax)
	GeoData_P->Dimension = _3D;
      else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Ymin != GeoData_P->Ymax)
	GeoData_P->Dimension = _2D;
      else if(GeoData_P->Xmin != GeoData_P->Xmax && GeoData_P->Zmin != GeoData_P->Zmax)
	GeoData_P->Dimension = _2D;
      else if(GeoData_P->Ymin != GeoData_P->Ymax && GeoData_P->Zmin != GeoData_P->Zmax)
	GeoData_P->Dimension = _2D;
      else if(GeoData_P->Xmin != GeoData_P->Xmax)
	GeoData_P->Dimension = _1D;
      else if(GeoData_P->Ymin != GeoData_P->Ymax)
	GeoData_P->Dimension = _1D;
      else if(GeoData_P->Zmin != GeoData_P->Zmax)
	GeoData_P->Dimension = _1D;
      else
	GeoData_P->Dimension = _0D;
      
      GeoData_P->CharacteristicLength = 
	sqrt(SQU(GeoData_P->Xmax - GeoData_P->Xmin) + 
	     SQU(GeoData_P->Ymax - GeoData_P->Ymin) + 
	     SQU(GeoData_P->Zmax - GeoData_P->Zmin));
      if(!GeoData_P->CharacteristicLength) 
	GeoData_P->CharacteristicLength = 1.;
    }

    /*  E L E M E N T S  */

    else if (!strncmp(&String[1], "ELM", 3) ||
	     !strncmp(&String[1], "Elements", 8)) {

      fscanf(File_GEO, "%d", &Nbr) ;
      if (GeoData_P->Elements == NULL)
	GeoData_P->Elements =
	  List_Create(Nbr, 1000, sizeof(struct Geo_Element)) ;

      Geo_Element.NbrEdges = Geo_Element.NbrFacets = 0 ;
      Geo_Element.NumEdges = Geo_Element.NumFacets = NULL ;

      for (i = 0 ; i < Nbr ; i++) {
	if(Version == 1.0){
	  fscanf(File_GEO, "%d %d %d %d %d",
		 &Geo_Element.Num, &Type, &Geo_Element.Region,
		 &iDummy, &Geo_Element.NbrNodes) ;
	  Geo_Element.Type = Geo_GetElementType(FORMAT_GMSH, Type) ;
	}
	else{
	  fscanf(File_GEO, "%d %d %d", &Geo_Element.Num, &Type, &NbTags);
	  Geo_Element.Region = 1;
	  for(j = 0; j < NbTags; j++){
	    fscanf(File_GEO, "%d", &iDummy);
	    if(j == 0)
	      Geo_Element.Region = iDummy;
	    /* ignore any other tags for now */
	  }
	  Geo_Element.Type = Geo_GetElementType(FORMAT_GMSH, Type) ;
	  Geo_Element.NbrNodes = Geo_GetNbNodesPerElement(Geo_Element.Type);
	}
	Geo_Element.FMMGroup = Geo_Element.Region ;
	Geo_Element.NumNodes = (int *)Malloc(Geo_Element.NbrNodes * sizeof(int)) ;
	for (j = 0 ; j < Geo_Element.NbrNodes ; j++)
	  fscanf(File_GEO, "%d", &Geo_Element.NumNodes[j]) ;

	List_Add(GeoData_P->Elements, &Geo_Element) ;
      }

      List_Sort(GeoData_P->Elements, fcmp_Elm) ;
      
    }

    do {
      fgets(String, MAX_STRING_LENGTH, File_GEO) ;
      if (feof(File_GEO)) Msg(GERROR, "Prematured end of file");
    } while (String[0] != '$') ;

  }   /* while 1 ... */

  GetDP_End ;
}


void  Geo_ReadFileAdapt(struct GeoData * GeoData_P) {

  struct Geo_Element Geo_Element, * Geo_Element_P ;
  int        Nbr, i, Index_GeoElement ;
  double     E, H, P, Max_Order = -1.0 ;
  char       String[MAX_STRING_LENGTH] ;

  GetDP_Begin("Geo_ReadFileAdapt");

  Nbr = List_Nbr(GeoData_P->Elements) ;

  if(!GeoData_P->H){
    GeoData_P->H = (double*)Malloc((Nbr+2)*sizeof(double)) ;
    for (i = 0 ; i < Nbr ; i++) GeoData_P->H[i+1] = -1.0 ;
  }
  if(!GeoData_P->P){
    GeoData_P->P = (double*)Malloc((Nbr+2)*sizeof(double)) ;
    for (i = 0 ; i < Nbr ; i++) GeoData_P->P[i+1] = -1.0 ;
  }

  while (1) {

    do { 
      fgets(String, MAX_STRING_LENGTH, File_GEO) ; 
      if (feof(File_GEO))  break ;
    } while (String[0] != '$') ;  
    
    if (feof(File_GEO))  break ;

    if (!strncmp(&String[1], "Adapt", 5)) {
      fscanf(File_GEO, "%d", &Nbr) ;
      for (i = 0 ; i < Nbr ; i++) {
	fscanf(File_GEO, "%d %lf %lf %lf", &Geo_Element.Num, &E, &H, &P) ;
	if(!(Geo_Element_P = (struct Geo_Element *)
	     List_PQuery(GeoData_P->Elements, &Geo_Element, fcmp_Elm)))
	  Msg(GERROR, "Element %d not found in database", Geo_Element.Num) ;
	Index_GeoElement = Geo_GetGeoElementIndex(Geo_Element_P) ;
	GeoData_P->H[Index_GeoElement+1] = H ;
	GeoData_P->P[Index_GeoElement+1] = P ;
	if(P > Max_Order) Max_Order = P ;
      }
    }

    do {
      fgets(String, MAX_STRING_LENGTH, File_GEO) ;
      if (feof(File_GEO)) Msg(GERROR, "Prematured end of file");
    } while (String[0] != '$') ;

  }   /* while 1 ... */

  if(Flag_ORDER < 0) Flag_ORDER = Max_Order ;
 
  Msg(INFO, "Maximum interpolation order = %g", Flag_ORDER) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  f c m p _ E l m   &   f c m p _ N o d                                   */
/* ------------------------------------------------------------------------ */

int  fcmp_Elm(const void * a, const void * b) {
  return  ((struct Geo_Element *)a)->Num - ((struct Geo_Element *)b)->Num ;
}

int  fcmp_Nod(const void * a, const void * b) {
  return  ((struct Geo_Node *)a)->Num - ((struct Geo_Node *)b)->Num ;
}



/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N b r G e o E l e m e n t s                               */
/* ------------------------------------------------------------------------ */

int  Geo_GetNbrGeoElements(void) {

  GetDP_Begin("Geo_GetNbrGeoElements");

  GetDP_Return(List_Nbr(CurrentGeoData->Elements)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t G e o E l e m e n t                                       */
/* ------------------------------------------------------------------------ */

struct Geo_Element  * Geo_GetGeoElement(int Index_Element) {

  GetDP_Begin("Geo_GetGeoElement");

  GetDP_Return((struct Geo_Element *)List_Pointer(CurrentGeoData->Elements, 
						  Index_Element)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t G e o E l e m e n t I n d e x                             */
/* ------------------------------------------------------------------------ */

int Geo_GetGeoElementIndex(struct Geo_Element * GeoElement) {

  GetDP_Begin("Geo_GetGeoElementIndex");

  GetDP_Return(GeoElement - 
	       (struct Geo_Element*)List_Pointer(CurrentGeoData->Elements, 0)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t G e o E l e m e n t O f N u m                             */
/* ------------------------------------------------------------------------ */

struct Geo_Element  * Geo_GetGeoElementOfNum(int Num_Element) {
  struct Geo_Element  elm ;

  GetDP_Begin("Geo_GetGeoElementOfNum");

  elm.Num = Num_Element ;

  GetDP_Return((struct Geo_Element*)List_PQuery(CurrentGeoData->Elements,
						&elm, fcmp_Elm)) ;
}

/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N b r G e o N o d e s                                     */
/* ------------------------------------------------------------------------ */

int  Geo_GetNbrGeoNodes(void) {

  GetDP_Begin("Geo_GetNbrGeoNodes");

  GetDP_Return(List_Nbr(CurrentGeoData->Nodes)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t G e o N o d e                                             */
/* ------------------------------------------------------------------------ */

struct Geo_Node  * Geo_GetGeoNode(int Index_Node) {

  GetDP_Begin("Geo_GetGeoNode");

  GetDP_Return((struct Geo_Node *)List_Pointer(CurrentGeoData->Nodes, 
					       Index_Node)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t G e o N o d e O f N u m                                   */
/* ------------------------------------------------------------------------ */

struct Geo_Node  * Geo_GetGeoNodeOfNum(int Num_Node) {
  struct Geo_Node  node ;

  GetDP_Begin("Geo_GetGeoNodeOfNum");

  node.Num = Num_Node ;

  GetDP_Return((struct Geo_Node*)List_PQuery(CurrentGeoData->Nodes,
					     &node, fcmp_Nod)) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N o d e s C o o r d i n a t e s                           */
/* ------------------------------------------------------------------------ */

void Geo_GetNodesCoordinates(int Nbr_Node, int * Num_Node,
			     double * x, double * y, double * z) {
  int    i ;
  struct Geo_Node  Geo_Node, * Geo_Node_P ;

  GetDP_Begin("Geo_GetNodesCoordinates");

  for (i = 0 ; i < Nbr_Node ; i++) {
    Geo_Node.Num = abs(Num_Node[i]) ;

    if(!(Geo_Node_P = (struct Geo_Node*)
	 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod)))
      Msg(GERROR, "Node %d does not exist", Geo_Node.Num) ;
    
    x[i] = Geo_Node_P->x ;  y[i] = Geo_Node_P->y ;  z[i] = Geo_Node_P->z ;
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ S e t N o d e s C o o r d i n a t e s                           */
/* ------------------------------------------------------------------------ */

void Geo_SetNodesCoordinates(int Nbr_Node, int * Num_Node,
			     double * x, double * y, double * z) {
  int    i ;
  struct Geo_Node  Geo_Node, * Geo_Node_P ;

  GetDP_Begin("Geo_SetNodesCoordinates");

  for (i = 0 ; i < Nbr_Node ; i++) {
    Geo_Node.Num = abs(Num_Node[i]) ;

    if(!(Geo_Node_P = (struct Geo_Node*)
	 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod)))
      Msg(GERROR, "Node %d does not exist", Geo_Node.Num) ;
    
    Geo_Node_P->x = x[i] ;  Geo_Node_P->y = y[i] ;  Geo_Node_P->z = z[i] ;
  }

  GetDP_End ;
}


void Geo_SetNodesCoordinatesX(int Nbr_Node, int * Num_Node,
			     double *x) {
  int    i ;
  struct Geo_Node  Geo_Node, * Geo_Node_P ;

  GetDP_Begin("Geo_SetNodesCoordinatesX");

  for (i = 0 ; i < Nbr_Node ; i++) {
    Geo_Node.Num = abs(Num_Node[i]) ;

    if(!(Geo_Node_P = (struct Geo_Node*)
	 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod)))
      Msg(GERROR, "Node %d does not exist", Geo_Node.Num) ;
    
    Geo_Node_P->x = x[i] ;
  }

  GetDP_End ;
}


void Geo_SetNodesCoordinatesY(int Nbr_Node, int * Num_Node,
			     double *y) {
  int    i ;
  struct Geo_Node  Geo_Node, * Geo_Node_P ;

  GetDP_Begin("Geo_SetNodesCoordinatesY");

  for (i = 0 ; i < Nbr_Node ; i++) {
    Geo_Node.Num = abs(Num_Node[i]) ;

    if(!(Geo_Node_P = (struct Geo_Node*)
	 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod)))
      Msg(GERROR, "Node %d does not exist", Geo_Node.Num) ;
    
    Geo_Node_P->y = y[i] ;
  }

  GetDP_End ;
}

void Geo_SetNodesCoordinatesZ(int Nbr_Node, int * Num_Node,
			     double *z) {
  int    i ;
  struct Geo_Node  Geo_Node, * Geo_Node_P ;

  GetDP_Begin("Geo_SetNodesCoordinatesZ");

  for (i = 0 ; i < Nbr_Node ; i++) {
    Geo_Node.Num = abs(Num_Node[i]) ;

    if(!(Geo_Node_P = (struct Geo_Node*)
	 List_PQuery(CurrentGeoData->Nodes, &Geo_Node, fcmp_Nod)))
      Msg(GERROR, "Node %d does not exist", Geo_Node.Num) ;
    
    Geo_Node_P->z = z[i] ;
  }

  GetDP_End ;
}