geodata.c

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

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

#include "GetDP.h"
#include "GeoData.h"
#include "Data_Passive.h"
#include "Data_Element.h"
#include "Numeric.h"
#include "Magic.h"

extern double Flag_ORDER ;

FILE  * File_GEO ;

struct GeoData  * CurrentGeoData ;


/* ------------------------------------------------------------------------ */
/*  G e o _ A d d G e o D a t a                                             */
/* ------------------------------------------------------------------------ */

int  Geo_AddGeoData(List_T * GeoData_L,
		    char * Name_MshFile, char * Name_DefaultMshFile,
		    char * Name_AdaptFile, char * Name_DefaultAdaptFile) {

  struct GeoData  GeoData_S ;
  int  i ;

  GetDP_Begin("Geo_AddGeoData");

  if (!Name_MshFile)  Name_MshFile = Name_DefaultMshFile ;

  if ((i = List_ISearchSeq(GeoData_L, Name_MshFile, fcmp_GeoData_Name)) < 0) {
    Msg(LOADING,"Geometric data '%s'", Name_MshFile) ;
    i = List_Nbr(GeoData_L) ;
    Geo_InitGeoData(&GeoData_S, i, Name_MshFile) ;
    Geo_OpenFile(Name_MshFile, "r") ;
    Geo_ReadFile(&GeoData_S) ;
    Geo_CloseFile() ;

    if (!Name_AdaptFile) Name_AdaptFile = Name_DefaultAdaptFile ;
    if (Name_AdaptFile) {
      Msg(LOADING,"Adaptation data '%s'", Name_AdaptFile) ;
      Geo_OpenFile(Name_AdaptFile, "r") ;
      Geo_SetCurrentGeoData(&GeoData_S) ;
      Geo_ReadFileAdapt(&GeoData_S) ;
      Geo_CloseFile() ;
    }
    List_Add(GeoData_L, &GeoData_S) ;
  }

  GetDP_Return(i) ;
}

int  fcmp_GeoData_Name(const void * a, const void * b) {
  return ( strcmp((char *)a, ((struct GeoData *)b)->Name ) ) ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ I n i t G e o D a t a                                           */
/* ------------------------------------------------------------------------ */

void  Geo_InitGeoData(struct GeoData * GeoData_P, int Num, char * Name) {

  GetDP_Begin("Geo_InitGeoData");

  GeoData_P->Num  = Num ;
  GeoData_P->Name = Name ;

  GeoData_P->Nodes = NULL ;
  GeoData_P->Elements = NULL ;

  GeoData_P->NbrElementsWithEdges = GeoData_P->NbrElementsWithFacets = 0 ;
  GeoData_P->NumCurrentEdge = GeoData_P->NumCurrentFacet = 0 ;
  GeoData_P->EdgesXNodes =
    Tree_Create(sizeof(struct Entity2XEntity1), fcmp_E2XE1) ;
  GeoData_P->FacetsXEdges =
    Tree_Create(sizeof(struct Entity2XEntity1), fcmp_E2XE1) ;

  GeoData_P->NodesXElements = NULL ;
  GeoData_P->Normals =
    Tree_Create(sizeof(struct EntityXVector), fcmp_EXVector) ;

  GeoData_P->GroupForPRE = NULL ;

  GeoData_P->Grid.Init = 0 ;

  GeoData_P->H = GeoData_P->P = NULL ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ S e t C u r r e n t G e o D a t a B a s e                       */
/* ------------------------------------------------------------------------ */

void  Geo_SetCurrentGeoData(struct GeoData * GeoData_P) {

  GetDP_Begin("Geo_SetCurrentGeoData");

  CurrentGeoData = GeoData_P ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ O p e n F i l e                                                 */
/* ------------------------------------------------------------------------ */

void  Geo_OpenFile(char * Name, char * Mode) {

  GetDP_Begin("Geo_OpenFile");

  File_GEO = fopen(Name, Mode) ;

  if (!File_GEO) Msg(GERROR, "Unable to open file '%s'", Name);

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ C l o s e F i l e                                               */
/* ------------------------------------------------------------------------ */

void  Geo_CloseFile(void) {

  GetDP_Begin("Geo_CloseFile");
  
  fclose(File_GEO) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ R e a d F i l e                                                 */
/* ------------------------------------------------------------------------ */

int Geo_GetElementType(int Format, int Type){

  GetDP_Begin("Geo_GetElementType");

  switch(Format){
  case FORMAT_GMSH :
    switch(Type){
    case 15 : GetDP_Return(POINT) ;
    case 1  : GetDP_Return(LINE) ;
    case 2  : GetDP_Return(TRIANGLE) ;
    case 3  : GetDP_Return(QUADRANGLE) ;
    case 4  : GetDP_Return(TETRAHEDRON) ;      
    case 5  : GetDP_Return(HEXAHEDRON) ;
    case 6  : GetDP_Return(PRISM) ;
    case 7  : GetDP_Return(PYRAMID) ;
    case 8  : GetDP_Return(LINE_2) ;
    case 9  : GetDP_Return(TRIANGLE_2) ;
    case 10 : GetDP_Return(QUADRANGLE_2) ;
    case 11 : GetDP_Return(TETRAHEDRON_2) ;      
    case 12 : GetDP_Return(HEXAHEDRON_2) ;
    case 13 : GetDP_Return(PRISM_2) ;
    case 14 : GetDP_Return(PYRAMID_2) ;
    default : Msg(GERROR, "Unkown type of Element in Gmsh format") ; GetDP_Return(-1) ;
    }
    break ;
  default :
    Msg(GERROR, "Unkown mesh format") ;
    GetDP_Return(-1) ;
  }

}


int Geo_GetElementTypeInv(int Format, int Type){

  GetDP_Begin("Geo_GetElementTypeInv");

  switch(Format){
  case FORMAT_GMSH :
    switch(Type){
    case POINT         : GetDP_Return(15) ;
    case LINE          : GetDP_Return(1) ;
    case TRIANGLE      : GetDP_Return(2) ;
    case QUADRANGLE    : GetDP_Return(3) ;
    case TETRAHEDRON   : GetDP_Return(4) ;      
    case HEXAHEDRON    : GetDP_Return(5) ;
    case PRISM         : GetDP_Return(6) ;
    case PYRAMID       : GetDP_Return(7) ;
    case LINE_2        : GetDP_Return(8) ;
    case TRIANGLE_2    : GetDP_Return(9) ;
    case QUADRANGLE_2  : GetDP_Return(10) ;
    case TETRAHEDRON_2 : GetDP_Return(11) ;      
    case HEXAHEDRON_2  : GetDP_Return(12) ;
    case PRISM_2       : GetDP_Return(13) ;
    case PYRAMID_2     : GetDP_Return(14) ;
    default : Msg(GERROR, "Unkown type of Element in Gmsh format") ; GetDP_Return(-1) ;
    }
    break ;
  default :
    Msg(GERROR, "Unkown mesh format") ;
    GetDP_Return(-1) ;
  }

}

int Geo_GetNbNodesPerElement(int Type){

  GetDP_Begin("Geo_GetElementNumNodes");

  switch(Type){
  case POINT         : GetDP_Return(1) ;
  case LINE          : GetDP_Return(2) ;
  case TRIANGLE      : GetDP_Return(3) ;
  case QUADRANGLE    : GetDP_Return(4) ;
  case TETRAHEDRON   : GetDP_Return(4) ;      
  case HEXAHEDRON    : GetDP_Return(8) ;
  case PRISM         : GetDP_Return(6) ;
  case PYRAMID       : GetDP_Return(5) ;
  case LINE_2        : GetDP_Return(3) ;
  case TRIANGLE_2    : GetDP_Return(6) ;
  case QUADRANGLE_2  : GetDP_Return(8) ;
  case TETRAHEDRON_2 : GetDP_Return(10) ;      
  case HEXAHEDRON_2  : GetDP_Return(20) ;
  case PRISM_2       : GetDP_Return(15) ;
  case PYRAMID_2     : GetDP_Return(13) ;
  default : Msg(GERROR, "Unkown type of Element") ; GetDP_Return(-1) ;
  }

}

void  Geo_SaveMesh(struct GeoData * GeoData_P, List_T * InitialList, char * FileName) {

  FILE * file;
  struct Geo_Node     Geo_Node ;
  struct Geo_Node   * Geo_Node_P ;
  struct Geo_Element  Geo_Element ;
  struct GeoData GeoData ;
  int i, j, Type, iDummy=0;
  int  fcmp_int(const void * a, const void * b);

  GetDP_Begin("Geo_SaveMesh");

  GeoData.Nodes    = List_Create(1000, 1000, sizeof(struct Geo_Node)) ;
  GeoData.Elements = List_Create(1000, 1000, sizeof(struct Geo_Node)) ;

  for (i = 0 ; i < List_Nbr(GeoData_P->Elements) ; i++) {
    List_Read(GeoData_P->Elements, i, &Geo_Element) ;
    if (List_Search(InitialList, &Geo_Element.Region, fcmp_int) ) {
      List_Add(GeoData.Elements, &Geo_Element) ;
      for (j = 0 ; j < Geo_Element.NbrNodes ; j++)
	if (!List_Search(GeoData.Nodes, 
	 		 Geo_Node_P = Geo_GetGeoNodeOfNum(Geo_Element.NumNodes[j]), fcmp_Nod) ) 
	  List_Add(GeoData.Nodes, Geo_Node_P) ;
    }
  }
 
  file = fopen(FileName,"w");
  Msg(INFO,"Saving mesh in file \"%s\" (%d nodes, %d elements)", 
      FileName, List_Nbr(GeoData.Nodes), List_Nbr(GeoData.Elements));

  fprintf(file, "$NOD\n%d\n", List_Nbr(GeoData.Nodes));

  for (i = 0 ; i < List_Nbr(GeoData.Nodes) ; i++) {
    List_Read(GeoData.Nodes, i, &Geo_Node) ;
    fprintf(file, "%d %.16g %.16g %.16g\n",
	   Geo_Node.Num, Geo_Node.x, Geo_Node.y, Geo_Node.z) ;
  }

  fprintf(file, "$ENDNOD\n$ELM\n%d\n", List_Nbr(GeoData.Elements));

  for (i = 0 ; i < List_Nbr(GeoData.Elements) ; i++) {
    List_Read(GeoData.Elements, i, &Geo_Element) ;
    Type = Geo_GetElementTypeInv(FORMAT_GMSH, Geo_Element.Type) ;
    fprintf(file, "%d %d %d %d %d ",
	    Geo_Element.Num, Type, Geo_Element.Region,
	    iDummy, Geo_Element.NbrNodes) ;
    for (j = 0 ; j < Geo_Element.NbrNodes ; j++)
      fprintf(file, "%d ", Geo_Element.NumNodes[j]) ;
    fprintf(file, "\n") ;
  }
    
  fprintf(file, "$ENDELM\n");

  fclose(file);
  List_Delete(GeoData.Nodes);
  List_Delete(GeoData.Elements);
  GetDP_End ;
}


void  Geo_ReadFile(struct GeoData * GeoData_P) {

  int                 Nbr, i, j, Type, iDummy, Format, Size, NbTags ;
  double              Version = 1.0 ;
  struct Geo_Node     Geo_Node ;
  struct Geo_Element  Geo_Element ;
  char                String[MAX_STRING_LENGTH] ;

  GetDP_Begin("Geo_ReadFile");

  while (1) {

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

    /*  F O R M A T  */

    if(!strncmp(&String[1], "MeshFormat", 10)) {

      fscanf(File_GEO, "%lf %d %d\n", &Version, &Format, &Size);
      if(Version != 2.0){
	Msg(GERROR, "Unknown mesh file version (%g)", Version);
	return;
      }

    }

    /*  N O D E S  */

    else if (!strncmp(&String[1], "NOE", 3) ||
	     !strncmp(&String[1], "NOD", 3) ||
	     !strncmp(&String[1], "Nodes", 5)) {

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

      for (i = 0 ; i < Nbr ; i++) {
	fscanf(File_GEO, "%d %lf %lf %lf",
	       &Geo_Node.Num, &Geo_Node.x, &Geo_Node.y, &Geo_Node.z) ;
	List_Add(GeoData_P->Nodes, &Geo_Node) ;
	if(!i){
	  GeoData_P->Xmin = GeoData_P->Xmax = Geo_Node.x;
	  GeoData_P->Ymin = GeoData_P->Ymax = Geo_Node.y;
	  GeoData_P->Zmin = GeoData_P->Zmax = Geo_Node.z; 
	}
	else{
	  GeoData_P->Xmin = MIN(GeoData_P->Xmin, Geo_Node.x);
	  GeoData_P->Xmax = MAX(GeoData_P->Xmax, Geo_Node.x);
	  GeoData_P->Ymin = MIN(GeoData_P->Ymin, Geo_Node.y);