cal_derhamtermoffemequation.c

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

      Get_JacobianFunction(Element->JacobianCase->TypeJacobian + 
			   Cell_RelativeJacobianType,
			   Cells[0].Type, &Cell_Type_Dimension) ;
  }
  
  /*  ------------------------------------------------------------------------  */
  /*  C o m p u t a t i o n   o f   E l e m e n t a r y   S u b m a t r i x     */
  /*  ------------------------------------------------------------------------  */
  
  /* Loop on source elements (> 1 only if integral quantity) */
  
  while (1) {
    
    if (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) {
      
      if ( Get_NextElementSource(Element->ElementSource) ) {
	
	/* Get DOF of source element */
	
	Get_DofOfElement(Element->ElementSource,
			 QuantityStorageDof_P->FunctionSpace,
			 QuantityStorageDof_P, NULL) ;
	
	/* Get function value for shape function */
	
	Get_NodesCoordinatesOfElement(Element->ElementSource) ;
	Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ;
	Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
			  QuantityStorageDof_P->DefineQuantity->IntegralQuantity.TypeOperatorDof, 
			  QuantityStorageDof_P, &FI->IntegralQuantityActive.Type_FormDof) ;
	
	/* Initialize Integral Quantity (integration & jacobian) */
	
	Cal_InitIntegralQuantity(Element, &FI->IntegralQuantityActive, 
				 QuantityStorageDof_P);
      }
      else { 
        break ;
      }
      
    }

    /* Elementary Submatrix initialization */

    for (i = 0 ; i < Nbr_Equ ; i++)
      for (j = 0 ; j < Nbr_Dof ; j++)
	for (k = 0 ; k < Current.NbrHar ; k++)
	  Ek[i][j][k] = 0. ;

    /* Cell Value initialization */
    
    for(i_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++)
      for (j = 0 ; j < Nbr_Dof ; j++)
	Cal_ZeroValue(&Cell_Value[i_Cell][j]);


    /* 
       printf("Nb Equ= %d  nbdof=%d cells =%d", Nbr_Equ , Nbr_Dof, Nbr_Cells); 
       */

    
    /*  -------------------------------------  */
    /*  P o i n t   c o l l o c a t i o n      */
    /*  -------------------------------------  */
    
    if(Cells[0].Type == POINT){
    
      for(i_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++) {

	Current.Element = Element ;	 
	Current.u = Cells[i_Cell].x[0] ;
	Current.v = Cells[i_Cell].y[0] ;
	Current.w = Cells[i_Cell].z[0] ;

	Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;
	Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;

	if (FI->Flag_InvJac)
	  Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
			    &Element->Jac, &Element->InvJac) ;

	/* Shape Functions (+ surrounding expression) */

	Cal_vBFxDof(EquationTerm_P, FI, 
		    QuantityStorage_P0, QuantityStorageDof_P,
		    Nbr_Dof, xFunctionBFDof, NULL, Cell_Value[i_Cell]);
      }
    }

    /*  -------------------------------------  */
    /*  I n t e g r a t i o n                  */
    /*  -------------------------------------  */
    
    else {

      switch (IntegrationCase_P->Type) {
	
      case GAUSS :  
      case GAUSSLEGENDRE :  
	
	Quadrature_P = (struct Quadrature*)
	  List_PQuery(IntegrationCase_P->Case, &Cells[0].Type, fcmp_int);
	
	if(!Quadrature_P)
	  Msg(GERROR, "Unknown type of Element (%s) for Integration method (%s)",
	      Get_StringForDefine(Element_Type, Cells[0].Type),
	      ((struct IntegrationMethod *)
	       List_Pointer(Problem_S.IntegrationMethod,
			    EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
	
	Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
	Get_IntPoint  = (void(*)(int,int,double*,double*,double*,double*))
	  Quadrature_P->Function ;    
		
	for(i_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++) {
	  
	  for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
	    
	    /* u,v,w in the Cell */
	    
	    Get_IntPoint(Nbr_IntPoints, i_IntPoint, &u, &v, &w, &weight) ;
	    
	    Get_BFGeoElement(&Cells[i_Cell], u, v, w) ;
	    
	    Cells[i_Cell].DetJac = Cell_Get_Jacobian(&Cells[i_Cell], &Cells[i_Cell].Jac) ;
	    
	    /* back to u,v,w in the original element. This is quite
               inefficient (approximatively 3 jacobian inversions
               instead of 1 if the transformation was defined in the
               reverse order... -> to change ! */
	    
	    x2 = y2 = z2 = 0. ;
	    for (i = 0 ; i < Cells[i_Cell].GeoElement->NbrNodes ; i++) {
	      x2 += Cells[i_Cell].x[i] * Cells[i_Cell].n[i];
	      y2 += Cells[i_Cell].y[i] * Cells[i_Cell].n[i];
	      z2 += Cells[i_Cell].z[i] * Cells[i_Cell].n[i];
	    }
	    xyz2uvwInAnElement(Element, x2, y2, z2, &u2, &v2, &w2) ;
	    
	    Current.Element = Element ;
	    Current.u = u2 ;
	    Current.v = v2 ;
	    Current.w = w2 ;

	    /* 
	       Msg(INFO, "Elm=%d  uvw=%g %g %g  xyz=%g %g %g",
	       Current.Element->Num, Current.u, Current.v, Current.w, x2, y2, z2);
	    */

	    /* Shape Functions (+ surrounding expression) */
	    
	    Cal_vBFxDof(EquationTerm_P, FI, 
			QuantityStorage_P0, QuantityStorageDof_P,
			Nbr_Dof, xFunctionBFDof, NULL, vBFxDof);
	    
	    Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Cells[i_Cell].DetJac) : weight ;
	    
	    /* Scalar product (by the tangent or normal to the cell). This is 
	       definitely a hack and should be completely changed */
	    
	    switch(FI->Type_FormDof){
	    case FORM0 :
	    case FORM3 :
	    case SCALAR :
	      break ;
	    case FORM1 :
	    case FORM2 :
	    case VECTOR :
	      for(j = 0 ; j < Nbr_Dof ; j++)
		Cal_ProductValue(&Cell_Vector[i_Cell], &vBFxDof[j], &vBFxDof[j]) ;
	      break;
	    default :
	      Msg(GERROR, "Unknown type of Dof (%s) for scalar product in deRham Map",
		  Get_StringForDefine(Field_Type, FI->Type_FormDof));
	      break;
	    }
	    
	    for (j = 0 ; j < Nbr_Dof ; j++)
	      for (k = 0 ; k < Current.NbrHar ; k++)	  
		Cell_Value[i_Cell][j].Val[MAX_DIM*k] += Factor * vBFxDof[j].Val[MAX_DIM*k];
	    
	    
	  }  /* for i_IntPoint ... */            
	  
	} /* for i_Cell */
	
	break ; /* case GAUSS */
	
      default :
	Msg(GERROR, "Unknown type of Integration method (%s)",
	    ((struct IntegrationMethod *)
	     List_Pointer(Problem_S.IntegrationMethod,
			  EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
	break;
	
      } /* switch integration method */

    } /* if point, else */


    /* Discrete derivative if any */

    switch (Group_P->FunctionType) {
      
    case BOUNDARYOFDUALNODESOF : 
      Product = 1 ;
      Dk = Geo_GetIM_Den_Xp(Element->Type, &Nbr_Ent1, &Nbr_Ent2); break;
      break ;

    case BOUNDARYOFDUALEDGESOF : 
      Product = 1 ;
      Dk = Geo_GetIM_Dfe_Xp(Element->Type, &Nbr_Ent1, &Nbr_Ent2); break;
      break;
      
    default :
      Product = 0 ; 
      break ;
    }

    /*
      if(Nbr_Cells != Nbr_Ent1 || Nbr_Equ != Nbr_Ent2)
      Msg(GERROR, "Error in deRham Map: wrong discrete derivative");
    */
    
    
    if(Flag_VERBOSE == 10) {
      printf("H D      = ") ;
      for (j = 0 ; j < Nbr_Dof ; j++)
	Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof) ; 
      printf("\n") ;
      for (i = 0 ; i < Nbr_Equ ; i++) { 
	Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof) ; 
	printf("[ ") ;
	for (j = 0 ; j < Nbr_Dof ; j++) 
	  printf("% .5e ", Cell_Value[i][j].Val[0]) ;
	printf("]\n") ;
      }
    }

    /* D^T (H D) if product ; H otherwise */

    if(Product){
      for (i = 0 ; i < Nbr_Equ ; i++)
	for (j = 0 ; j < Nbr_Dof ; j++)
	  for (n = 0 ; n < Nbr_Cells ; n++)
	    if((ii = Dk [ n*Nbr_Ent1 + QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement ])){
	      for (k = 0 ; k < Current.NbrHar ; k++)
		Ek[i][j][k] += ii * Cell_Value[n][j].Val[MAX_DIM*k] ;
	    }
      
      if(Flag_VERBOSE == 10) {
	printf("D^T H D  = ") ;
	for (j = 0 ; j < Nbr_Dof ; j++)
	  Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof) ; 
	printf("\n") ;
	for (i = 0 ; i < Nbr_Equ ; i++) { 
	  Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof) ; 
	  printf("[ ") ;
	  for (j = 0 ; j < Nbr_Dof ; j++) 
	    printf("% .5e ", Ek[i][j][0]) ;
	  printf("]\n") ;
	}
      }

    }
    else{ /* ok, this is stupid */
      for (i = 0 ; i < Nbr_Equ ; i++)
	for (j = 0 ; j < Nbr_Dof ; j++)
	  for (k = 0 ; k < Current.NbrHar ; k++)
	    Ek[i][j][k] = Cell_Value[i][j].Val[MAX_DIM*k] ;
    }


    
    
    /* Assembly in global matrix */
    /* ------------------------- */

    for (i = 0 ; i < Nbr_Equ ; i++)  
      for (j = 0 ; j < Nbr_Dof ; j++)
	((void (*)(struct Dof*, struct Dof*, double*)) 
	 FI->Function_AssembleTerm)
	  (QuantityStorageEqu_P->BasisFunction[i].Dof,
	   QuantityStorageDof_P->BasisFunction[j].Dof,  Ek[i][j]) ;


    /* Exit while(1) directly if not integral quantity */

    if (FI->Type_DefineQuantityDof != INTEGRALQUANTITY)  break ;

  }  /* while (1) ... */

  GetDP_End ;
}