get_dofofelement.c

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

#define RCSID "$Id: Get_DofOfElement.c,v 1.27 2006/02/26 00:42:54 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 "Get_DofOfElement.h"
#include "ExtendedGroup.h"
#include "Cal_Quantity.h"
#include "GeoData.h"
#include "CurrentData.h"
#include "Tools.h"

struct BasisFunction    * BasisFunction_P ;
int                     Nbr_ElementaryBF, Flag_SubSpace ;
struct Group            * GroupSupport_P, * GroupEntity_P ;

/* ------------------------------------------------------------------------ */
/*  G e t _ I n i t D o f O f E l e m e n t                                 */
/* ------------------------------------------------------------------------ */

void  Get_InitDofOfElement(struct Element * Element) {

  GetDP_Begin("Get_InitDofOfElement");

  Element->ElementTrace = NULL ;
  Element->NumLastElementForNodesCoordinates      = -1 ;
  Element->NumLastElementForGroupsOfEntities      = -1 ;
  Element->NumLastElementForSolidAngle            = -1 ;
  Element->NumLastElementForSortedNodesByFacet    = -1 ;

  GetDP_End ;
}


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

void  Get_DofOfElement(struct Element          * Element,
		       struct FunctionSpace    * FunctionSpace_P,
		       struct QuantityStorage  * QuantityStorage_P,
		       List_T                  * BasisFunctionIndex_L) {

  struct BasisFunction  * BasisFunction_P0 ;
  int  Nbr_BasisFunction, Nbr_BasisFunctionAll, i_BFunction, StartingIndex, i ;
  int  * BasisFunctionIndex_P0 = NULL ;

  GetDP_Begin("Get_DofOfElement");

  Current.Element = Element ;
  Nbr_ElementaryBF = 0 ;

  /* Get the SubSpace */

  Nbr_BasisFunctionAll = List_Nbr(FunctionSpace_P->BasisFunction) ;
  BasisFunction_P0 = (Nbr_BasisFunctionAll) ? 
    (struct BasisFunction*)List_Pointer(FunctionSpace_P->BasisFunction, 0) : NULL ;

  if (!BasisFunctionIndex_L) {
    Flag_SubSpace = 0 ;
    Nbr_BasisFunction = Nbr_BasisFunctionAll ;
  }
  else {
    Flag_SubSpace = 1 ;
    Nbr_BasisFunction = List_Nbr(BasisFunctionIndex_L) ;
    BasisFunctionIndex_P0 = (Nbr_BasisFunction) ? 
      (int*)List_Pointer(BasisFunctionIndex_L, 0) : NULL ;
  }

  /* Set the DofData if explicitely specified */
  
  switch (TreatmentStatus) {
  case _CAL :  case _POS :
    if(QuantityStorage_P->DefineQuantity->DofData)
      FunctionSpace_P->DofData = QuantityStorage_P->DefineQuantity->DofData ;
    else
      FunctionSpace_P->DofData = FunctionSpace_P->MainDofData ;
    break;
  }

  /*  For each subset of Basis Functions */

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

    i_BFunction = (!Flag_SubSpace)? i : BasisFunctionIndex_P0[i] ;

    BasisFunction_P = BasisFunction_P0 + i_BFunction ;
    GroupSupport_P = (struct Group*)
      List_Pointer(Problem_S.Group, BasisFunction_P->SupportIndex) ;

    /*  If the BasisFunction exists for this kind of element
	   the interpolation order is lower or equal to the maximum order allowed
	   the element is in the support of the BasisFunction */

    if ( ( BasisFunction_P->ElementType & Current.Element->Type ) 
	 &&
	 ( Flag_ORDER < 0. || BasisFunction_P->Order <= Flag_ORDER )
	 &&
	 ( (GroupSupport_P->Type == REGIONLIST  &&
	    List_Search(GroupSupport_P->InitialList, &Element->Region, fcmp_int))
	   ||
	   (GroupSupport_P->Type == ELEMENTLIST  &&
	    Check_IsEntityInExtendedGroup(GroupSupport_P, Element->Num, 0)) ) ) {

      GroupEntity_P = (struct Group*)
	List_Pointer(Problem_S.Group, BasisFunction_P->EntityIndex) ;

      switch (GroupEntity_P->FunctionType) {

      case NODESOF :
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->GeoElement->NbrNodes, Element->GeoElement->NumNodes,
	   0,  i_BFunction, NODESOF, NULL) ;
	break ;

      case EDGESOF :  case EDGESOFTREEIN :
	if (Element->GeoElement->NbrEdges == 0)
	  Geo_CreateEdgesOfElement(Element->GeoElement) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->GeoElement->NbrEdges, Element->GeoElement->NumEdges,
	   0,  i_BFunction, EDGESOF, NULL) ;
	break ;

      case FACETSOF :  case FACETSOFTREEIN :
	if (Element->GeoElement->NbrEdges == 0)
	  Geo_CreateEdgesOfElement(Element->GeoElement) ;
	if (Element->GeoElement->NbrFacets == 0)
	  Geo_CreateFacetsOfElement(Element->GeoElement) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->GeoElement->NbrFacets, Element->GeoElement->NumFacets,
	   0,  i_BFunction, FACETSOF, NULL) ;
	break ;

      case VOLUMESOF :
	Get_CodesOfElement(FunctionSpace_P, QuantityStorage_P,
			   1, &Element->GeoElement->Num,
			   0,  i_BFunction, VOLUMESOF, NULL) ;
	break ;

      case GROUPSOFNODESOF :
	Get_GroupsOfElementaryEntitiesOfElement
	  (Element, &StartingIndex,
	   Element->GeoElement->NbrNodes, Element->GeoElement->NumNodes,
	   BasisFunction_P) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->NbrGroupsOfEntities, Element->NumGroupsOfEntities,
	   StartingIndex,  i_BFunction, GROUPSOFNODESOF, Element->NumSubFunction[1]) ;
	break ;

      case GROUPSOFEDGESONNODESOF :
	if (Element->GeoElement->NbrEdges == 0)
	  Geo_CreateEdgesOfElement(Element->GeoElement) ;
	Get_GroupsOfEdgesOnNodesOfElement(Element, &StartingIndex) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->NbrGroupsOfEntities, Element->NumGroupsOfEntities,
	   StartingIndex,  i_BFunction, GROUPSOFEDGESONNODESOF, NULL) ;
	break ;

      case GROUPSOFEDGESOF :
	if (Element->GeoElement->NbrEdges == 0)
	  Geo_CreateEdgesOfElement(Element->GeoElement) ;
	Get_GroupsOfElementaryEntitiesOfElement
	  (Element, &StartingIndex,
	   Element->GeoElement->NbrEdges, Element->GeoElement->NumEdges,
	   BasisFunction_P) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->NbrGroupsOfEntities, Element->NumGroupsOfEntities,
	   StartingIndex,  i_BFunction, GROUPSOFEDGESOF, NULL) ;
	break ;

      case REGION :
	Get_RegionForElement(Element, &StartingIndex, BasisFunction_P) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->NbrGroupsOfEntities, Element->NumGroupsOfEntities,
	   StartingIndex,  i_BFunction, REGION, Element->NumSubFunction[1]) ;
	break ;

      case GLOBAL :
	Get_GlobalForElement(Element, &StartingIndex, BasisFunction_P) ;
	Get_CodesOfElement
	  (FunctionSpace_P, QuantityStorage_P,
	   Element->NbrGroupsOfEntities, Element->NumGroupsOfEntities,
	   StartingIndex,  i_BFunction, GLOBAL, NULL) ;
	break ;
      }

    }  /* if Region ... */

  }  /* for i ... */

  QuantityStorage_P->NbrElementaryBasisFunction = Nbr_ElementaryBF ;

  GetDP_End ;
}


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

void  Get_GroupsOfElementaryEntitiesOfElement
  (struct Element * Element,
   int * StartingIndex, int Nbr_ElementaryEntities, int Num_ElementaryEntities[],
   struct BasisFunction * BasisFunction_P) {

  /* external input/output :  GroupEntity_P     : In  */

  int            i, j, Num_Entity, Nbr_SubFunction, i_SF ;
  struct TwoInt  * Key_P, Key ;
  int iFirst;

  GetDP_Begin("Get_GroupsOfElementaryEntitiesOfElement");

  if (Element->NumLastElementForGroupsOfEntities != Element->Num) {
    Element->NumLastElementForGroupsOfEntities = Element->Num ;
    Element->NbrGroupsOfEntities = 0 ;
  }

  *StartingIndex = Element->NbrGroupsOfEntities ;

  if (GroupEntity_P->ExtendedList == NULL) Generate_ExtendedGroup(GroupEntity_P) ;

  for (i = 0 ; i < Nbr_ElementaryEntities ; i++) {
    Num_Entity = abs(Num_ElementaryEntities[i]) ;

    /*
    if ((Key_P = (struct TwoInt*)List_PQuery(GroupEntity_P->ExtendedList, 
					     &Num_Entity, fcmp_absint ))) {
    */
    iFirst = 1;
    Key.Int1 = Num_Entity;
    Key_P = &Key;
    while (List_LQuery(GroupEntity_P->ExtendedList, &Key, fcmp_absint, iFirst)) {
      iFirst = 0;

      j = *StartingIndex ;
      while ((j < Element->NbrGroupsOfEntities) &&
	     (Element->NumGroupsOfEntities[j] != Key_P->Int2))  j++ ;

      if (!BasisFunction_P->SubFunction) {
	if (j == Element->NbrGroupsOfEntities) {
	  Element->NumSubFunction[1][j] = 0 ;
	  Element->NumSubFunction[0][j] = -1 ;
	  Element->NumGroupsOfEntities[j] = Key_P->Int2 ;
	  Element->NbrEntitiesInGroups[Element->NbrGroupsOfEntities++] = 0 ;
	}
	Element->NumEntitiesInGroups[j][Element->NbrEntitiesInGroups[j]++] =
	  (Key_P->Int1 > 0)?  (i+1) : -(i+1) ;
      }
      else { /* For SubFunctions (basis functions for a global function) */
	Nbr_SubFunction = List_Nbr(BasisFunction_P->SubFunction) ;

	if (j == Element->NbrGroupsOfEntities) {
	  for (i_SF = 0 ; i_SF < Nbr_SubFunction ; i_SF++) {
	    Element->NumSubFunction[1][j+i_SF] = i_SF ;
	    Element->NumSubFunction[0][j+i_SF] =
	      *((int *)List_Pointer(BasisFunction_P->SubFunction, i_SF)) ;
	    if (BasisFunction_P->SubdFunction)
	      Element->NumSubFunction[2][j+i_SF] =
		*((int *)List_Pointer(BasisFunction_P->SubdFunction, i_SF)) ;
	    Element->NumGroupsOfEntities[j+i_SF] = Key_P->Int2 ;
	    Element->NbrEntitiesInGroups[Element->NbrGroupsOfEntities++] = 0 ;
	  }
	}
	for (i_SF = 0 ; i_SF < Nbr_SubFunction ; i_SF++)
	  Element->NumEntitiesInGroups[j+i_SF][Element->NbrEntitiesInGroups[j+i_SF]++] =
	    (Key_P->Int1 > 0)?  (i+1) : -(i+1) ;
      }

    } /* if Key_P */
  }

  GetDP_End ;
}

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

void  Get_GroupsOfEdgesOnNodesOfElement(struct Element * Element,
					int * StartingIndex) {
  /* external input/output :  GroupEntity_P     : In  */

  int  i, j, Num_Edge, * Num_Nodes, Num_Node ;

  GetDP_Begin("Get_GroupsOfEdgesOnNodesOfElement");
  
  if (Element->NumLastElementForGroupsOfEntities != Element->Num) {
    Element->NumLastElementForGroupsOfEntities = Element->Num ;