geoentity.c

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

#define RCSID "$Id: GeoEntity.c,v 1.13 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 "GeoEntity.h"
#include "Magic.h"
#include "Tools.h"

extern FILE            * File_PRE ;
extern struct GeoData  * CurrentGeoData ;
extern int               Flag_XDATA ;


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N o d e s _ u v w                                         */
/* ------------------------------------------------------------------------ */

double * Geo_GetNodes_uvw(int Type, int *nbn){

  GetDP_Begin("Geo_GetNodes_uvw");

  switch(Type){
  case POINT       : *nbn = NbrNodes_Point ;       GetDP_Return(*Nodes_Point) ;
  case LINE        : *nbn = NbrNodes_Line ;        GetDP_Return(*Nodes_Line) ;
  case TRIANGLE    : *nbn = NbrNodes_Triangle ;    GetDP_Return(*Nodes_Triangle) ;
  case QUADRANGLE  : *nbn = NbrNodes_Quadrangle ;  GetDP_Return(*Nodes_Quadrangle) ;
  case TETRAHEDRON : *nbn = NbrNodes_Tetrahedron ; GetDP_Return(*Nodes_Tetrahedron) ;
  case HEXAHEDRON  : *nbn = NbrNodes_Hexahedron ;  GetDP_Return(*Nodes_Hexahedron) ;
  case PRISM       : *nbn = NbrNodes_Prism ;       GetDP_Return(*Nodes_Prism) ;
  case PYRAMID     : *nbn = NbrNodes_Pyramid ;     GetDP_Return(*Nodes_Pyramid) ;
  default : 
    Msg(GERROR, "Unknown type of Element in Geo_GetNodes_uvw") ; GetDP_Return(NULL) ;
  }
}

/* ------------------------------------------------------------------------ */
/*  G e o _ C r e a t e E d g e s O f E l e m e n t                         */
/* ------------------------------------------------------------------------ */

void  Geo_CreateEdgesOfElement(struct Geo_Element * Geo_Element) {

  int  Nbr_Entities2, * D_Element ;

  GetDP_Begin("Geo_CreateEdgesOfElement");

  D_Element = Geo_GetIM_Den(Geo_Element->Type, &Nbr_Entities2) ;
  Geo_CreateEntitiesOfElement(Nbr_Entities2, D_Element,
			      Geo_Element->NbrNodes, Geo_Element->NumNodes,
			      &Geo_Element->NbrEdges, &Geo_Element->NumEdges,
			      &CurrentGeoData->NbrElementsWithEdges,
			      &CurrentGeoData->NumCurrentEdge,
			      CurrentGeoData->EdgesXNodes) ;
  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N o d e s O f E d g e I n E l e m e n t                   */
/* ------------------------------------------------------------------------ */

int  * Geo_GetNodesOfEdgeInElement(struct Geo_Element * Geo_Element, int Num_Edge) {
  int  Nbr_Entities2 ;

  GetDP_Begin("Geo_GetBodesOfEdgeInElement");

  GetDP_Return( Geo_GetIM_Den(Geo_Element->Type, &Nbr_Entities2) + 
		Num_Edge * NBR_MAX_SUBENTITIES_IN_ELEMENT ) ;
}

/* ------------------------------------------------------------------------ */
/*  G e o _ G e t N o d e s O f F a c e t I n E l e m e n t                 */
/* ------------------------------------------------------------------------ */

int  * Geo_GetNodesOfFacetInElement(struct Geo_Element * Geo_Element, int Num_Facet) {
  int  Nbr_Entities2 ;

  GetDP_Begin("Geo_GetNodesOfFacetInElement");
  
  GetDP_Return( Geo_GetIM_Dfn(Geo_Element->Type, &Nbr_Entities2) + 
		Num_Facet * NBR_MAX_SUBENTITIES_IN_ELEMENT ) ;
}

/* ------------------------------------------------------------------------ */
/*  G e o _ C r e a t e F a c e t s O f E l e m e n t                       */
/* ------------------------------------------------------------------------ */

void  Geo_CreateFacetsOfElement(struct Geo_Element * Geo_Element) {
  int  Nbr_Entities2, * D_Element ;

  GetDP_Begin("Geo_CreatefacetsOfElement");

  D_Element = Geo_GetIM_Dfe(Geo_Element->Type, &Nbr_Entities2) ;
  Geo_CreateEntitiesOfElement(Nbr_Entities2, D_Element,
			      Geo_Element->NbrEdges, Geo_Element->NumEdges,
			      &Geo_Element->NbrFacets, &Geo_Element->NumFacets,
			      &CurrentGeoData->NbrElementsWithFacets,
			      &CurrentGeoData->NumCurrentFacet,
			      CurrentGeoData->FacetsXEdges) ;
  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t I M _ D e n                                               */
/* ------------------------------------------------------------------------ */

int  * Geo_GetIM_Den(int Type_Element, int * Nbe) {

  GetDP_Begin("Geo_GetIM_Den");

  switch (Type_Element) {
  case LINE :        *Nbe = NbrEdges_Line        ; GetDP_Return(*Den_Line) ;
  case TRIANGLE :    *Nbe = NbrEdges_Triangle    ; GetDP_Return(*Den_Triangle) ;
  case QUADRANGLE :  *Nbe = NbrEdges_Quadrangle  ; GetDP_Return(*Den_Quadrangle) ;
  case TETRAHEDRON : *Nbe = NbrEdges_Tetrahedron ; GetDP_Return(*Den_Tetrahedron) ;
  case HEXAHEDRON :  *Nbe = NbrEdges_Hexahedron  ; GetDP_Return(*Den_Hexahedron) ;
  case PRISM :       *Nbe = NbrEdges_Prism       ; GetDP_Return(*Den_Prism) ;
  case PYRAMID :     *Nbe = NbrEdges_Pyramid     ; GetDP_Return(*Den_Pyramid) ;
  default : 
    Msg(GERROR, "Unknown incidence matrix for element type %d", Type_Element);
    GetDP_Return(NULL) ;
  }
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t I M _ D f e                                               */
/* ------------------------------------------------------------------------ */

int  * Geo_GetIM_Dfe(int Type_Element, int * Nbf) {

  GetDP_Begin("Geo_GetIM_Dfe");

  switch (Type_Element) {
  case LINE :        *Nbf = NbrFacets_Line        ; GetDP_Return(NULL) ;
  case TRIANGLE :    *Nbf = NbrFacets_Triangle    ; GetDP_Return(*Dfe_Triangle) ;
  case QUADRANGLE :  *Nbf = NbrFacets_Quadrangle  ; GetDP_Return(*Dfe_Quadrangle) ;
  case TETRAHEDRON : *Nbf = NbrFacets_Tetrahedron ; GetDP_Return(*Dfe_Tetrahedron) ;
  case HEXAHEDRON :  *Nbf = NbrFacets_Hexahedron  ; GetDP_Return(*Dfe_Hexahedron) ;
  case PRISM :       *Nbf = NbrFacets_Prism       ; GetDP_Return(*Dfe_Prism) ;
  case PYRAMID :     *Nbf = NbrFacets_Pyramid     ; GetDP_Return(*Dfe_Pyramid) ;
  default :
    Msg(GERROR, "Unknown incidence matrix for element type %d", Type_Element);
    GetDP_Return(NULL) ;
  }
}

/* ------------------------------------------------------------------------ */
/*  G e o _ G e t I M _ D f n                                               */
/* ------------------------------------------------------------------------ */

int  * Geo_GetIM_Dfn(int Type_Element, int * Nbf) {

  GetDP_Begin("Geo_GetIM_Dfn");

  switch (Type_Element) {
  case LINE :       *Nbf = NbrFacets_Line        ; GetDP_Return(NULL) ;
  case TRIANGLE :   *Nbf = NbrFacets_Triangle    ; GetDP_Return(*Dfn_Triangle) ;
  case QUADRANGLE : *Nbf = NbrFacets_Quadrangle  ; GetDP_Return(*Dfn_Quadrangle) ;
  case TETRAHEDRON :*Nbf = NbrFacets_Tetrahedron ; GetDP_Return(*Dfn_Tetrahedron) ;
  case HEXAHEDRON : *Nbf = NbrFacets_Hexahedron  ; GetDP_Return(*Dfn_Hexahedron) ;
  case PRISM :	    *Nbf = NbrFacets_Prism       ; GetDP_Return(*Dfn_Prism) ;
  case PYRAMID :    *Nbf = NbrFacets_Pyramid     ; GetDP_Return(*Dfn_Pyramid) ;
  default :
    Msg(GERROR, "Unknown incidence matrix for element type %d", Type_Element);
    GetDP_Return(NULL) ;
  }
}

/* ------------------------------------------------------------------------ */
/*  G e o _ G e t I M _ D e n _ X p                                         */
/* ------------------------------------------------------------------------ */

int * Geo_GetIM_Den_Xp(int Type_Element, int * Nbe, int * Nbn) {

  GetDP_Begin("Geo_GetIM_Den_Xp");

  switch (Type_Element) {
  case LINE :
    *Nbe = NbrEdges_Line ; *Nbn = NbrNodes_Line ; 
    GetDP_Return(Den_Line_Xp) ;
  case TRIANGLE :
    *Nbe = NbrEdges_Triangle ; *Nbn = NbrNodes_Triangle ; 
    GetDP_Return(Den_Triangle_Xp) ;
  case QUADRANGLE :
    *Nbe = NbrEdges_Quadrangle ; *Nbn = NbrNodes_Quadrangle ; 
    GetDP_Return(Den_Quadrangle_Xp) ;
  case TETRAHEDRON :
    *Nbe = NbrEdges_Tetrahedron ; *Nbn = NbrNodes_Tetrahedron ; 
    GetDP_Return(Den_Tetrahedron_Xp) ;
  case HEXAHEDRON :
    *Nbe = NbrEdges_Hexahedron ; *Nbn = NbrNodes_Hexahedron ; 
    GetDP_Return(Den_Hexahedron_Xp) ;
  case PRISM :
    *Nbe = NbrEdges_Prism ; *Nbn = NbrNodes_Prism ; 
    GetDP_Return(Den_Prism_Xp) ;
  case PYRAMID :
    *Nbe = NbrEdges_Pyramid ; *Nbn = NbrNodes_Pyramid ; 
    GetDP_Return(Den_Pyramid_Xp) ;
  default :
    Msg(GERROR, "Unknown incidence matrix for element type %d", Type_Element);
    GetDP_Return(NULL) ;
  }
}


/* ------------------------------------------------------------------------ */
/*  G e o _ G e t I M _ D f e _ X p                                         */
/* ------------------------------------------------------------------------ */

int * Geo_GetIM_Dfe_Xp(int Type_Element, int * nbf, int * nbe) {

  GetDP_Begin("Geo_GetIM_Dfe_Xp");

  switch (Type_Element) {
  case LINE :
    *nbf = NbrFacets_Line ; *nbe = NbrEdges_Line ; 
    GetDP_Return(NULL) ;
  case TRIANGLE :
    *nbf = NbrFacets_Triangle ; *nbe = NbrEdges_Triangle ; 
    GetDP_Return(Dfe_Triangle_Xp) ;
  case QUADRANGLE :
    *nbf = NbrFacets_Quadrangle ; *nbe = NbrEdges_Quadrangle ; 
    GetDP_Return(Dfe_Quadrangle_Xp) ;
  case TETRAHEDRON :
    *nbf = NbrFacets_Tetrahedron ; *nbe = NbrEdges_Tetrahedron ; 
    GetDP_Return(Dfe_Tetrahedron_Xp) ;
  case HEXAHEDRON :
    *nbf = NbrFacets_Hexahedron ; *nbe = NbrEdges_Hexahedron ; 
    GetDP_Return(Dfe_Hexahedron_Xp) ;
  case PRISM :
    *nbf = NbrFacets_Prism ; *nbe = NbrEdges_Prism ; 
    GetDP_Return(Dfe_Prism_Xp) ;
  case PYRAMID :
    *nbf = NbrFacets_Pyramid ; *nbe = NbrEdges_Pyramid ; 
    GetDP_Return(Dfe_Pyramid_Xp) ;
  default :
    Msg(GERROR, "Unknown incidence matrix for element type %d", Type_Element);
    GetDP_Return(NULL) ;
  }
}

/* ------------------------------------------------------------------------ */
/*  G e o _ C r e a t e E n t i t i e s O f E l e m e n t                   */
/* ------------------------------------------------------------------------ */

void  Geo_CreateEntitiesOfElement
  (int Nbr_Entities2, int * D_Element,
   int   Geo_Element_NbrEntities1, int  * Geo_Element_NumEntities1,
   int * Geo_Element_NbrEntities2, int ** Geo_Element_NumEntities2,
   int * Geo_NbrElementsWithEntities2,
   int * Geo_NumCurrentEntity2, Tree_T * Geo_Entities2XEntities1) {

  int  i, j, Nbr_Entities1, Num_Entities1[NBR_MAX_SUBENTITIES_IN_ELEMENT], Sign_Entity2 ;
  int  * Entity_P, Entity ;
  struct Entity2XEntity1  Entity2XEntities1, * Entity2XEntities1_P ;

  GetDP_Begin("Geo_CreateEntitiesOfElement");

  *Geo_Element_NbrEntities2 = Nbr_Entities2 ;
  *Geo_Element_NumEntities2 = (int *)Malloc(Nbr_Entities2 * sizeof(int)) ;
  (*Geo_NbrElementsWithEntities2)++ ;

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

    Entity_P = D_Element + i * NBR_MAX_SUBENTITIES_IN_ELEMENT ;