pos_feminterpolation.c

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

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

#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Get_DofOfElement.h"
#include "Get_Geometry.h"
#include "GeoData.h"
#include "DofData.h"
#include "Pos_Search.h"
#include "CurrentData.h"
#include "Tools.h"

/* ------------------------------------------------------------------------ */
/*  P o s _ F e m I n t e r p o l a t i o n                                 */
/* ------------------------------------------------------------------------ */

int SearchType(int type){

  GetDP_Begin("SearchType");
  
  switch(type){
  case POINT       : GetDP_Return(_ALL) ;
  case LINE        : GetDP_Return(_1D) ;
  case TRIANGLE    :
  case QUADRANGLE  : GetDP_Return(_2D) ;
  case TETRAHEDRON :
  case HEXAHEDRON  :
  case PRISM       :
  case PYRAMID     : GetDP_Return(_3D) ;
  default          : Msg(GERROR, "Unknown Type in 'SearchType'"); GetDP_Return(-1) ;
  }
}

void  Pos_FemInterpolation(struct Element * Element, 
			   struct QuantityStorage * QuantityStorage_P0,
			   struct QuantityStorage * QuantityStorage_P,
			   int Type_Quantity, int Type_Operator, 
			   int Type_Dimension, int UseXYZ, 
			   double u, double v, double w, 
			   double x, double y, double z, 
			   double Val[], int * Type_Value,
			   int Flag_ChangeOfCoordinates) {

  void (*xFunctionBF[NBR_MAX_BASISFUNCTIONS])
    (struct Element *, int, double, double, double, double []) ;
  void (*xChangeOfCoordinates) () = 0;

  struct IntegralQuantityActive    IQA ;
  struct Value                     vBFxDof[NBR_MAX_BASISFUNCTIONS] ;
  struct GeoData                 * GeoData_P ;
  struct Element                   TheElement, * TheElement_P ;
  struct QuantityStorage         * QS_P ;

  double  vBFu[NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
  double  Val_Dof, Val_Dof_r, Val_Dof_i ;

  int  Type_DefineQuantity, SubType_DefineQuantity, Type_Form  ;
  int  i, j, k, Nbr_Dof = 0 ;
  int  GeoDataNum = 0, UseNewGeo = 0 ;

  GetDP_Begin("Pos_FemInterpolation");

  /* -------------
     Quantity Type 
     ------------- */

  Type_DefineQuantity = QuantityStorage_P->DefineQuantity->Type ;

  if(Type_DefineQuantity == INTEGRALQUANTITY){
    
    if(QuantityStorage_P->DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){
      SubType_DefineQuantity = NODOF ; 
    }
    else{
      SubType_DefineQuantity = INTEGRALQUANTITY ;
    }
  }
  else{
    SubType_DefineQuantity = Type_DefineQuantity ;
  }


  /* ---------------
     Get The Element 
     --------------- */

  if(SubType_DefineQuantity != NODOF) {

    if(!QuantityStorage_P->FunctionSpace)
      Msg(GERROR, "No available function space for quantity");

    if(!QuantityStorage_P->FunctionSpace->DofData)
      Msg(GERROR, "No available data to interpolate quantity");

    GeoDataNum = QuantityStorage_P->FunctionSpace->DofData->GeoDataIndex;
    UseNewGeo = (GeoDataNum != Current.GeoData->Num) ;

    if(UseXYZ || UseNewGeo){
      if(UseNewGeo){
	GeoData_P = (struct GeoData *)List_Pointer(GeoData_L, GeoDataNum);
	GeoDataNum = Current.GeoData->Num ;
	Geo_SetCurrentGeoData(Current.GeoData = GeoData_P) ;
      }
      if(!UseXYZ){
	x = y = z = 0. ;
	for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
	  x += Element->x[i] * Element->n[i] ;
	  y += Element->y[i] * Element->n[i] ;
	  z += Element->z[i] * Element->n[i] ;
	}
      }
      Init_SearchGrid(&Current.GeoData->Grid);
      InWhichElement(Current.GeoData->Grid, NULL, &TheElement, 
		     (Type_Dimension >= 0) ? Type_Dimension : SearchType(Element->Type),
		     x, y, z, &u, &v, &w) ;
      TheElement_P = &TheElement ;
      Get_InitDofOfElement(&TheElement) ;
      Get_DofOfElement
	(&TheElement, QuantityStorage_P->FunctionSpace, QuantityStorage_P,
	 QuantityStorage_P->DefineQuantity->IndexInFunctionSpace) ;
    }
    else{
      TheElement_P = Element;
    }
  }
  else{
    TheElement_P = Element ;
  }


  /* ------------------
     Init LocalQuantity
     ------------------ */

  if (Type_DefineQuantity == LOCALQUANTITY) {
    if (TheElement_P->Num != NO_ELEMENT) {
      Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction ;      
      Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBF, Type_Operator, 
			QuantityStorage_P, &Type_Form) ;
      xChangeOfCoordinates =
	(void (*)())Get_ChangeOfCoordinates
	((Flag_ChangeOfCoordinates && TheElement_P->Num != NO_ELEMENT),
	 Type_Form) ;
    }
    else {
      Msg(WARNING, "No element found in mesh for LocalQuantity interpolation");
      Nbr_Dof = 0 ;
      Type_Form = VECTOR ;
    }
  }

  /* ---------------------
     Init IntegralQuantity
     --------------------- */

  else if (Type_DefineQuantity == INTEGRALQUANTITY) {

    if(Type_Operator != NOOP){
      Msg(GERROR, "Operator acting on Integral Quantity");
    }
    Type_Form = VECTOR ;
    Get_InitElementSource(TheElement_P,
			  QuantityStorage_P->DefineQuantity->IntegralQuantity.InIndex) ;

    IQA.IntegrationCase_L =
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    QuantityStorage_P->DefineQuantity->
		    IntegralQuantity.IntegrationMethodIndex)) ->IntegrationCase ;
    IQA.CriterionIndex =
      ((struct IntegrationMethod *)
       List_Pointer(Problem_S.IntegrationMethod,
		    QuantityStorage_P->DefineQuantity->
		    IntegralQuantity.IntegrationMethodIndex)) ->CriterionIndex ;
    IQA.JacobianCase_L =
      ((struct JacobianMethod *)
       List_Pointer(Problem_S.JacobianMethod,
		    QuantityStorage_P->DefineQuantity->
		    IntegralQuantity.JacobianMethodIndex)) ->JacobianCase ;

    xChangeOfCoordinates = (void (*)())Get_ChangeOfCoordinates(0, Type_Form) ;
  }
  

  /* ----------------------
     Compute GlobalQuantity
     ---------------------- */
  
  if (Type_DefineQuantity == GLOBALQUANTITY) {
    if(Current.NbrHar==1){
      if (Type_Quantity == QUANTITY_BF)
	Val[0] = (QuantityStorage_P->BasisFunction[0].Dof->Entity ==
		  Current.SubRegion)? 1. : 0. ;
      else
	Dof_GetRealDofValue
	  (QuantityStorage_P->FunctionSpace->DofData,
	   QuantityStorage_P->BasisFunction[0].Dof,
	   &Val[0]) ;
    }
    else{	      
      for (k = 0 ; k < Current.NbrHar ; k+=2) {
	if (Type_Quantity == QUANTITY_BF) {
	  Val[MAX_DIM*k] = (QuantityStorage_P->BasisFunction[0].Dof->Entity ==
			    Current.SubRegion)? 1. : 0. ;
	  Val[MAX_DIM*(k+1)] = 0. ;
	}
	else {
	  Dof_GetComplexDofValue
	    (QuantityStorage_P->FunctionSpace->DofData,
	     QuantityStorage_P->BasisFunction[0].Dof + k/2*gCOMPLEX_INCREMENT,
	     &Val[MAX_DIM*k], &Val[MAX_DIM*(k+1)]) ;
	}
      }
    }
    *Type_Value = SCALAR ;
    GetDP_End ;
  }
  

  /* -----------------------------------
     Compute Local / Integral Quantities
     ----------------------------------- */

  i = Current.NbrHar * MAX_DIM ;
  for (k = 0 ; k < i ; k++)  Val[k] =  0. ;
  
  while (1) {
    
    if (Type_DefineQuantity == INTEGRALQUANTITY) {
      
      if (Get_NextElementSource(TheElement_P->ElementSource)) {
	
	Get_NodesCoordinatesOfElement(TheElement_P->ElementSource) ;
	
	if(SubType_DefineQuantity != NODOF){
	  Get_DofOfElement(TheElement_P->ElementSource, 
			   QuantityStorage_P->FunctionSpace,
			   QuantityStorage_P,
			   QuantityStorage_P->DefineQuantity->IndexInFunctionSpace) ;	  
	  Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction ;
	  Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBF,
			    QuantityStorage_P->DefineQuantity->IntegralQuantity.TypeOperatorDof,
			    QuantityStorage_P, &IQA.Type_FormDof) ;
	  Type_Form = IQA.Type_FormDof ; /* good form */
	}
	else{
	  Nbr_Dof = 1 ; 
	  xFunctionBF[0] = NULL ; /* for analytic integration tests */
	  Type_Form = IQA.Type_FormDof = VECTOR ; /* form type unknown */

	  for (j = 0 ; j < QuantityStorage_P->DefineQuantity->IntegralQuantity.NbrQuantityIndex ; j++) {
	    QS_P = QuantityStorage_P0 + 
	      QuantityStorage_P->DefineQuantity->IntegralQuantity.QuantityIndexTable[j] ;
	    Get_DofOfElement(TheElement_P->ElementSource, 
			     QS_P->FunctionSpace, 
			     QS_P,
			     QS_P->DefineQuantity->IndexInFunctionSpace) ;
	  }
	}