pos_element.c

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

#define RCSID "$Id: Pos_Element.c,v 1.28 2006/02/26 00:42:59 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):
 *   Christophe Trophime
 */

#include "GetDP.h"
#include "GeoData.h"
#include "Get_Geometry.h"
#include "Get_DofOfElement.h"
#include "Data_Active.h"
#include "Pos_Element.h"
#include "Cal_Value.h"
#include "CurrentData.h"

/* ------------------------------------------------------------------------ */
/*  Create/Destroy/Compare                                                  */
/* ------------------------------------------------------------------------ */

void Alloc_PostElement(struct PostElement * PostElement){
  PostElement->NumNodes = (int *) Malloc(PostElement->NbrNodes * sizeof(int)) ;
  /* allocate as much as possible in one step */
  PostElement->u = (double *) Malloc(6 * PostElement->NbrNodes * sizeof(double)) ;
  PostElement->v = &PostElement->u[PostElement->NbrNodes] ;
  PostElement->w = &PostElement->u[2*PostElement->NbrNodes] ;
  PostElement->x = &PostElement->u[3*PostElement->NbrNodes] ;
  PostElement->y = &PostElement->u[4*PostElement->NbrNodes] ;
  PostElement->z = &PostElement->u[5*PostElement->NbrNodes] ;
  PostElement->Value = (struct Value *) 
    Malloc(PostElement->NbrNodes * sizeof(struct Value)) ;
}

struct PostElement * Create_PostElement(int Index, int Type, 
					int NbrNodes, int Depth){
  struct PostElement * PostElement ;

  GetDP_Begin("Create_PostElement");

  PostElement = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
  PostElement->Index = Index ; 
  PostElement->Type = Type ;
  PostElement->Depth = Depth ;
  PostElement->NbrNodes = NbrNodes ;
  if(NbrNodes > 0) Alloc_PostElement(PostElement);

  GetDP_Return(PostElement) ;
}

void Destroy_PostElement(struct PostElement * PostElement){

  GetDP_Begin("Destroy_PostElement");

  if(PostElement->NbrNodes > 0){
    Free(PostElement->NumNodes) ;
    if(PostElement->u) Free(PostElement->u); /* normal case */
    else if(PostElement->x) Free(PostElement->x); /* partial copy */
    Free(PostElement->Value) ;
  }
  Free(PostElement) ;

  GetDP_End ;
}

struct PostElement * NodeCopy_PostElement(struct PostElement *PostElement){
  struct PostElement * Copy ;
  int i ;

  GetDP_Begin("Copy_PostElement");

  Copy = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
  Copy->Index = PostElement->Index ; 
  Copy->Type = PostElement->Type ;
  Copy->Depth = PostElement->Depth ;
  Copy->NbrNodes = PostElement->NbrNodes ;
  if(Copy->NbrNodes > 0){
    Alloc_PostElement(Copy);
    for(i = 0 ; i < Copy->NbrNodes ; i++){
      Copy->NumNodes[i] = PostElement->NumNodes[i];
      Copy->u[i] = PostElement->u[i] ;
      Copy->v[i] = PostElement->v[i] ;
      Copy->w[i] = PostElement->w[i] ;
      Copy->x[i] = PostElement->x[i] ;
      Copy->y[i] = PostElement->y[i] ;
      Copy->z[i] = PostElement->z[i] ;
    }
  }

  GetDP_Return(Copy) ;
}

struct PostElement * PartialCopy_PostElement(struct PostElement *PostElement){
  struct PostElement * Copy ;
  int i ;

  GetDP_Begin("PartialCopy_PostElement");

  Copy = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
  Copy->Index = PostElement->Index ; 
  Copy->Type = PostElement->Type ;
  Copy->Depth = PostElement->Depth ;
  Copy->NbrNodes = PostElement->NbrNodes ;
  if(Copy->NbrNodes > 0){
    Copy->NumNodes = NULL ;
    Copy->u = Copy->v = Copy->w = NULL ;
    /* allocate as much as possible in one step */
    Copy->x = (double *) Malloc(3* Copy->NbrNodes * sizeof(double)) ;
    Copy->y = &Copy->x[Copy->NbrNodes];
    Copy->z = &Copy->x[2 * Copy->NbrNodes];
    Copy->Value = (struct Value *) Malloc(Copy->NbrNodes * sizeof(struct Value)) ;
    for(i = 0 ; i < Copy->NbrNodes ; i++){
      Copy->x[i] = PostElement->x[i] ;
      Copy->y[i] = PostElement->y[i] ;
      Copy->z[i] = PostElement->z[i] ;
      Cal_CopyValue(&PostElement->Value[i], &Copy->Value[i]);
    }
  }

  GetDP_Return(Copy) ;
}

/* 2 PostElements never have the same barycenter unless they are identical */

int fcmp_PostElement(const void *a, const void *b){
  struct PostElement *PE1, *PE2 ;
  double s1, s2, TOL=Current.GeoData->CharacteristicLength * 1.e-12 ;
  int i;

  PE1 = *(struct PostElement**)a; PE2 = *(struct PostElement**)b;

  if(PE1->NbrNodes != PE2->NbrNodes) 
    return PE1->NbrNodes - PE2->NbrNodes;

  s1 = s2 = 0.0 ;  
  for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->x[i]; s2 += PE2->x[i]; }
  if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;

  s1 = s2 = 0.0 ;
  for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->y[i]; s2 += PE2->y[i]; }
  if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;

  s1 = s2 = 0.0 ;
  for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->z[i]; s2 += PE2->z[i]; }
  if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;

  return 0;
}  

int fcmp_PostElement_v0(const void *a, const void *b){
  return (int)( (*(struct PostElement**)a)->v[0] - 
		(*(struct PostElement**)b)->v[0] ) ;
}  

int fcmp_PostElement_absu0(const void *a, const void *b){
  return (int)( fabs((*(struct PostElement**)b)->u[0]) - 
		fabs((*(struct PostElement**)a)->u[0]) ) ;
}  

/* ------------------------------------------------------------------------ */
/*  C u t _ P o s t E l e m e n t                                           */
/* ------------------------------------------------------------------------ */

void Cut_PostElement(struct PostElement * PE, struct Geo_Element * GE,
		     List_T * PE_L, int Index, int Depth, int Skin,
		     int DecomposeInSimplex){

  struct  Element      E ;
  struct  PostElement  * C[8] ;  
  
  double  u01, u02, u03, u12, u13, u23 ;
  double  v01, v02, v03, v12, v13, v23 ;
  double  w01, w02, w03, w12, w13, w23 ;
  int     i, j, NbCut = 0 ;

  GetDP_Begin("Cut_PostElement");

  /* Recursive division */

  if(PE->Depth < Depth){

    switch(PE->Type){

    case POINT :
      Msg(GERROR, "Impossible to divide a Point recursively");
      break;
      
    case LINE :
      u01 = .5 * (PE->u[0] + PE->u[1]);
      v01 = .5 * (PE->v[0] + PE->v[1]);
      w01 = .5 * (PE->w[0] + PE->w[1]);
      
      C[0] = Create_PostElement(Index, LINE, 2, PE->Depth) ;
      C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
      C[0]->u[1] = u01      ; C[0]->v[1] = v01      ; C[0]->w[1] = w01      ;
      
      C[1] = PE ;
      C[1]->u[0] = u01      ; C[1]->v[0] = v01      ; C[1]->w[0] = w01      ;
      
      NbCut = 2 ;
      break;
      
    case TRIANGLE :
      u01 = .5 * (PE->u[0] + PE->u[1]); u02 = .5 * (PE->u[0] + PE->u[2]); 
      v01 = .5 * (PE->v[0] + PE->v[1]);	v02 = .5 * (PE->v[0] + PE->v[2]); 
      w01 = .5 * (PE->w[0] + PE->w[1]);	w02 = .5 * (PE->w[0] + PE->w[2]); 

      u12 = .5 * (PE->u[1] + PE->u[2]);
      v12 = .5 * (PE->v[1] + PE->v[2]);
      w12 = .5 * (PE->w[1] + PE->w[2]);
           
      C[0] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
      C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
      C[0]->u[1] = u01      ; C[0]->v[1] = v01      ; C[0]->w[1] = w01      ;
      C[0]->u[2] = u02      ; C[0]->v[2] = v02      ; C[0]->w[2] = w02      ;

      C[1] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
      C[1]->u[0] = u01      ; C[1]->v[0] = v01      ; C[1]->w[0] = w01      ;
      C[1]->u[1] = PE->u[1] ; C[1]->v[1] = PE->v[1] ; C[1]->w[1] = PE->w[1] ;
      C[1]->u[2] = u12      ; C[1]->v[2] = v12      ; C[1]->w[2] = w12      ;
      
      C[2] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
      C[2]->u[0] = u02      ; C[2]->v[0] = v02      ; C[2]->w[0] = w02      ;
      C[2]->u[1] = u12      ; C[2]->v[1] = v12      ; C[2]->w[1] = w12      ;
      C[2]->u[2] = PE->u[2] ; C[2]->v[2] = PE->v[2] ; C[2]->w[2] = PE->w[2] ;

      C[3] = PE ;
      C[3]->u[0] = u01      ; C[3]->v[0] = v01      ; C[3]->w[0] = w01      ;
      C[3]->u[1] = u12      ; C[3]->v[1] = v12      ; C[3]->w[1] = w12      ;
      C[3]->u[2] = u02      ; C[3]->v[2] = v02      ; C[3]->w[2] = w02      ;
      
      NbCut = 4 ;
      break;

    case TETRAHEDRON :
      u01 = .5 * (PE->u[0] + PE->u[1]); u02 = .5 * (PE->u[0] + PE->u[2]); 
      v01 = .5 * (PE->v[0] + PE->v[1]);	v02 = .5 * (PE->v[0] + PE->v[2]); 
      w01 = .5 * (PE->w[0] + PE->w[1]);	w02 = .5 * (PE->w[0] + PE->w[2]); 

      u03 = .5 * (PE->u[0] + PE->u[3]); u12 = .5 * (PE->u[1] + PE->u[2]); 
      v03 = .5 * (PE->v[0] + PE->v[3]); v12 = .5 * (PE->v[1] + PE->v[2]); 
      w03 = .5 * (PE->w[0] + PE->w[3]); w12 = .5 * (PE->w[1] + PE->w[2]); 
    
      u13 = .5 * (PE->u[1] + PE->u[3]); u23 = .5 * (PE->u[2] + PE->u[3]);
      v13 = .5 * (PE->v[1] + PE->v[3]);	v23 = .5 * (PE->v[2] + PE->v[3]);
      w13 = .5 * (PE->w[1] + PE->w[3]);	w23 = .5 * (PE->w[2] + PE->w[3]);

      C[0] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
      C[0]->u[1] = u01     ; C[0]->v[1] = v01     ; C[0]->w[1] = w01     ;
      C[0]->u[2] = u02     ; C[0]->v[2] = v02     ; C[0]->w[2] = w02     ;
      C[0]->u[3] = u03     ; C[0]->v[3] = v03     ; C[0]->w[3] = w03     ;

      C[1] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[1]->u[0] = PE->u[1] ; C[1]->v[0] = PE->v[1] ; C[1]->w[0] = PE->w[1] ;
      C[1]->u[1] = u01     ; C[1]->v[1] = v01     ; C[1]->w[1] = w01     ;
      C[1]->u[2] = u12     ; C[1]->v[2] = v12     ; C[1]->w[2] = w12     ;
      C[1]->u[3] = u13     ; C[1]->v[3] = v13     ; C[1]->w[3] = w13     ;

      C[2] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[2]->u[0] = PE->u[2] ; C[2]->v[0] = PE->v[2] ; C[2]->w[0] = PE->w[2] ;
      C[2]->u[1] = u02     ; C[2]->v[1] = v02     ; C[2]->w[1] = w02     ;
      C[2]->u[2] = u12     ; C[2]->v[2] = v12     ; C[2]->w[2] = w12     ;
      C[2]->u[3] = u23     ; C[2]->v[3] = v23     ; C[2]->w[3] = w23     ;

      C[3] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[3]->u[0] = PE->u[3] ; C[3]->v[0] = PE->v[3] ; C[3]->w[0] = PE->w[3] ;
      C[3]->u[1] = u03     ; C[3]->v[1] = v03     ; C[3]->w[1] = w03     ;
      C[3]->u[2] = u13     ; C[3]->v[2] = v13     ; C[3]->w[2] = w13     ;
      C[3]->u[3] = u23     ; C[3]->v[3] = v23     ; C[3]->w[3] = w23     ;

      C[4] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[4]->u[0] = u01     ; C[4]->v[0] = v01     ; C[4]->w[0] = w01     ;
      C[4]->u[1] = u02     ; C[4]->v[1] = v02     ; C[4]->w[1] = w02     ;
      C[4]->u[2] = u03     ; C[4]->v[2] = v03     ; C[4]->w[2] = w03     ;
      C[4]->u[3] = u23     ; C[4]->v[3] = v23     ; C[4]->w[3] = w23     ;

      C[5] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[5]->u[0] = u01     ; C[5]->v[0] = v01     ; C[5]->w[0] = w01     ;
      C[5]->u[1] = u02     ; C[5]->v[1] = v02     ; C[5]->w[1] = w02     ;
      C[5]->u[2] = u12     ; C[5]->v[2] = v12     ; C[5]->w[2] = w12     ;
      C[5]->u[3] = u23     ; C[5]->v[3] = v23     ; C[5]->w[3] = w23     ;

      C[6] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
      C[6]->u[0] = u01     ; C[6]->v[0] = v01     ; C[6]->w[0] = w01     ;
      C[6]->u[1] = u12     ; C[6]->v[1] = v12     ; C[6]->w[1] = w12     ;
      C[6]->u[2] = u13     ; C[6]->v[2] = v13     ; C[6]->w[2] = w13     ;
      C[6]->u[3] = u23     ; C[6]->v[3] = v23     ; C[6]->w[3] = w23     ;

      C[7] = PE ;
      C[7]->u[0] = u01     ; C[7]->v[0] = v01     ; C[7]->w[0] = w01     ;
      C[7]->u[1] = u03     ; C[7]->v[1] = v03     ; C[7]->w[1] = w03     ;
      C[7]->u[2] = u13     ; C[7]->v[2] = v13     ; C[7]->w[2] = w13     ;
      C[7]->u[3] = u23     ; C[7]->v[3] = v23     ; C[7]->w[3] = w23     ;

      NbCut = 8 ;
      break ;

    default :
      Msg(GERROR, "Recursive division not implemented for Quadrangles, Hexahedra, "
	  "Prisms and Pyramids") ;
    }
    
    for(i = 0 ; i < NbCut ; i++){
      C[i]->Depth ++ ;