dofdata.c

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

	Msg(GERROR,"LinkCplx only valid for Complex systems") ;
      else{
	valtmp[0] = Val[0] * Dof_P->Case.Link.Coef - Val[1] * Dof_P->Case.Link.Coef2 ;
	valtmp[1] = Val[1] * Dof_P->Case.Link.Coef + Val[0] * Dof_P->Case.Link.Coef2 ;
      }
      Dof_AssembleInVec(Equ_P, Dof_P->Case.Link.Dof, NbrHar, 
			valtmp, OtherSolution, Vec0, Vec) ;
      break ;

    case DOF_FIXED_SOLVE :  case DOF_FIXEDWITHASSOCIATE_SOLVE :
      Msg(GERROR,"Wrong Constraints: "
	  "remaining Dof(s) waiting to be fixed by a Resolution");
      break;

    case DOF_UNKNOWN_INIT :
      Msg(GERROR,"Wrong Initial Constraints: "
	  "remaining Dof(s) with non-fixed initial conditions");
      break;
    }
    break ;

  case DOF_LINK :
    if(NbrHar==1)
      valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
    else{
      valtmp[0] = Val[0] * Equ_P->Case.Link.Coef ;
      valtmp[1] = Val[1] * Equ_P->Case.Link.Coef ;
    }
    Dof_AssembleInVec(Equ_P->Case.Link.Dof, Dof_P, NbrHar, 
		      valtmp, OtherSolution, Vec0, Vec) ;
    break ;

  case DOF_LINKCPLX :
    if(NbrHar==1)
      Msg(GERROR,"LinkCplx only valid for Complex systems") ;
    else{ /* Warning: conjugate! */
      valtmp[0] = Val[0] * Equ_P->Case.Link.Coef + Val[1] * Equ_P->Case.Link.Coef2 ;
      valtmp[1] = Val[1] * Equ_P->Case.Link.Coef - Val[0] * Equ_P->Case.Link.Coef2 ;
    }
    Dof_AssembleInVec(Equ_P->Case.Link.Dof, Dof_P, NbrHar, 
		      valtmp, OtherSolution, Vec0, Vec) ;
    break ;

  }

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  D o f _ T r a n s f e r S o l u t i o n T o C o n s t r a i n t         */
/* ------------------------------------------------------------------------ */

void  Dof_TransferSolutionToConstraint(struct DofData * DofData_P) {

  struct Dof * Dof_P, * Dof_P0 ;
  int  i ;

  GetDP_Begin("Dof_TransferSolutionToConstraint");

  if (!DofData_P->NbrAnyDof){
    GetDP_End;
  }

  Dof_P0 = (struct Dof *)List_Pointer(DofData_P->DofList, 0) ;


  for (i = 0 ; i < DofData_P->NbrAnyDof ; i++) {
    Dof_P = Dof_P0 + i ;

    switch (Dof_P->Type) {

    case DOF_UNKNOWN :
      Dof_P->Type = DOF_FIXED ;
      LinAlg_GetScalarInVector(&Dof_P->Val, 
			       &DofData_P->CurrentSolution->x, 
			       Dof_P->Case.Unknown.NumDof-1) ;
      Dof_P->Case.FixedAssociate.TimeFunctionIndex = 0 ;
      break ;

    case DOF_FIXED :
    case DOF_FIXEDWITHASSOCIATE :
    case DOF_LINK :
    case DOF_LINKCPLX :
      break ;

    default :  break ;
    }
  }

  DofData_P->NbrDof = 0 ;

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  D o f _ G e t D o f V a l u e                                           */
/* ------------------------------------------------------------------------ */

gScalar Dof_GetDofValue(struct DofData * DofData_P, struct Dof * Dof_P) {
  gScalar tmp ;
  
  GetDP_Begin("Dof_GetDofValue");

  switch (Dof_P->Type) {

  case DOF_UNKNOWN :
    if(!DofData_P->CurrentSolution->SolutionExist){
      Msg(GERROR, "Empty solution in DofData %d", DofData_P->Num);
    }
    LinAlg_GetScalarInVector(&tmp, &DofData_P->CurrentSolution->x, 
			     Dof_P->Case.Unknown.NumDof-1) ;
    break ;

  case DOF_FIXED :
  case DOF_FIXEDWITHASSOCIATE :
    LinAlg_ProdScalarDouble(&Dof_P->Val, 
			    ((Dof_P->Case.FixedAssociate.TimeFunctionIndex)?
			     DofData_P->CurrentSolution->TimeFunctionValues
			     [Dof_P->Case.FixedAssociate.TimeFunctionIndex] :
			     1.),
			    &tmp);
    break ;

  case DOF_LINK :
    tmp = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ;
    LinAlg_ProdScalarDouble(&tmp, Dof_P->Case.Link.Coef, &tmp) ;
    break ;

  case DOF_LINKCPLX :
    /* Too soon to treat LinkCplx: we need the real and imaginary parts */
    Msg(GERROR, "Cannot call Dof_GetDofValue for LinkCplx");
    break ;

  default :
    LinAlg_ZeroScalar(&tmp) ;
    break ;
  }

  GetDP_Return(tmp) ;
}

void Dof_GetRealDofValue(struct DofData * DofData_P, struct Dof * Dof_P, double *d) {
  gScalar tmp ;

  GetDP_Begin("Dof_GetRealDofValue");

  if (Dof_P->Type == DOF_LINKCPLX) {
    Msg(GERROR, "Cannot call Dof_GetRealDofValue for LinkCplx");
  }

  tmp = Dof_GetDofValue(DofData_P, Dof_P) ;
  LinAlg_GetDoubleInScalar(d, &tmp) ;

  GetDP_End ;
}

void Dof_GetComplexDofValue(struct DofData * DofData_P, struct Dof * Dof_P, 
			    double *d1, double *d2) {
  gScalar tmp1, tmp2 ;
  double valtmp[2] ;

  GetDP_Begin("Dof_GetComplexDofValue");

  if(gSCALAR_SIZE == 1){
    if(Dof_P->Type == DOF_LINKCPLX) { /* Can only be done here */
      if (Dof_P->Case.Link.Dof->Type == DOF_LINKCPLX) { /* recurse */
	Dof_GetComplexDofValue(DofData_P, Dof_P->Case.Link.Dof, d1, d2);
      }
      else{
	tmp1 = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ;
	tmp2 = Dof_GetDofValue(DofData_P, (Dof_P+1)->Case.Link.Dof) ;
	LinAlg_GetDoubleInScalar(d1, &tmp1) ;
	LinAlg_GetDoubleInScalar(d2, &tmp2) ;
      }
    }
    else{
      tmp1 = Dof_GetDofValue(DofData_P, Dof_P) ;
      tmp2 = Dof_GetDofValue(DofData_P, Dof_P+1) ;
      LinAlg_GetDoubleInScalar(d1, &tmp1) ;
      LinAlg_GetDoubleInScalar(d2, &tmp2) ;
    }
  }
  else{
    if (Dof_P->Type == DOF_LINKCPLX) { /* Can only be done here */
      if (Dof_P->Case.Link.Dof->Type == DOF_LINKCPLX) { /* recurse */
	Dof_GetComplexDofValue(DofData_P, Dof_P->Case.Link.Dof, d1, d2);
      }
      else{
	tmp1 = Dof_GetDofValue(DofData_P, Dof_P->Case.Link.Dof) ; 
	LinAlg_GetComplexInScalar(d1, d2, &tmp1) ;
      }
    }
    else{
      tmp1 = Dof_GetDofValue(DofData_P, Dof_P) ; 
      LinAlg_GetComplexInScalar(d1, d2, &tmp1) ;
    }
  }

  if(Dof_P->Type == DOF_LINKCPLX){
    valtmp[0] = Dof_P->Case.Link.Coef*(*d1) - Dof_P->Case.Link.Coef2*(*d2) ;
    valtmp[1] = Dof_P->Case.Link.Coef*(*d2) + Dof_P->Case.Link.Coef2*(*d1) ;
    *d1 = valtmp[0] ;
    *d2 = valtmp[1] ; 
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------- */
/*  D o f _ D e f i n e Unknown D o f F r o m Solve o r Init D o f           */
/* ------------------------------------------------------------------------- */

void  Dof_DefineUnknownDofFromSolveOrInitDof(struct DofData ** DofData_P) {
  int         i, Nbr_AnyDof ;
  struct Dof  * Dof_P ;

  GetDP_Begin("Dof_DefineUnknownDofFromSolveOrInitDof");

  Nbr_AnyDof = List_Nbr((*DofData_P)->DofList) ;

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

    Dof_P = (struct Dof*)List_Pointer((*DofData_P)->DofList, i) ;

    switch (Dof_P->Type) {
    case DOF_FIXED_SOLVE :
    case DOF_FIXEDWITHASSOCIATE_SOLVE :
      Dof_P->Type = DOF_UNKNOWN ;
      Dof_P->Case.Unknown.NumDof = ++((*DofData_P)->NbrDof) ;
      break ;
    case DOF_UNKNOWN_INIT :
      Dof_P->Type = DOF_UNKNOWN ;
      break ;
    }
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ T r a n s f e r D o f                                           */
/* ------------------------------------------------------------------------ */

void  Dof_TransferDof(struct DofData  * DofData_P1, 
		      struct DofData ** DofData_P2) {

  int              i, Nbr_AnyDof ;
  struct Dof       Dof, * Dof_P ;
  struct Solution  * Solutions_P0 ;

  GetDP_Begin("Dof_TransferDof");

  Nbr_AnyDof   = List_Nbr(DofData_P1->DofList) ;
  Solutions_P0 = (struct Solution*)List_Pointer(DofData_P1->Solutions, 0) ;
  DofData_P1->CurrentSolution = Solutions_P0 ;

  for(i = 0; i < Nbr_AnyDof; i++) {
    Dof = *(struct Dof *)List_Pointer(DofData_P1->DofList, i) ;

    if((Dof_P = (struct Dof*)Tree_PQuery((*DofData_P2)->DofTree, &Dof))){

      switch (Dof_P->Type) {
      case DOF_FIXED_SOLVE :
	Dof_P->Type = DOF_FIXED ;
	Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
	break ;
      case DOF_FIXEDWITHASSOCIATE_SOLVE :
	Dof_P->Type = DOF_FIXEDWITHASSOCIATE ;
	Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
	break ;
      case DOF_UNKNOWN_INIT :
	Dof_P->Val = Dof_GetDofValue(DofData_P1, &Dof) ;
	break ;

	/* Un DOF_UNKNOWN_INIT prendra toujours une valeur obtenue par
	   pre-resolution, meme si on ne demande qu'une simple initialisation ...
	   Pourquoi pas definir un type DOF_UNKNOWN_INIT_SOLVE,
	   propre a Dof_DefineInitSolveDof() ? ... */

      }
    }
  }
  
  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  D o f _ I n i t D o f F o r N o D o f                                   */
/* ------------------------------------------------------------------------ */

void Dof_InitDofForNoDof(struct Dof * DofForNoDof, int NbrHar) {
  int k ;
  double Val[2] = {1.,0.} ;

  GetDP_Begin("Dof_InitDofForNoDof");

  for (k=0 ; k<NbrHar ; k+=gSCALAR_SIZE) {
    DofForNoDof[k].Type = DOF_FIXED ;
    LinAlg_SetScalar(&DofForNoDof[k].Val, &Val[k%2]) ;
    DofForNoDof[k].Case.FixedAssociate.TimeFunctionIndex = 0 ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------- */
/*  P r i n t _  D o f N u m b e r                         */
/* ------------------------------------------------------- */

void Print_DofNumber(struct Dof *Dof_P){

  GetDP_Begin("Print_DofNumber");

  switch(Dof_P->Type){
  case DOF_UNKNOWN : 
    printf("%d ", Dof_P->Case.Unknown.NumDof) ; 
    break ;
  case DOF_FIXED   : 
    printf("Fixed ") ; 
    break ;
  case DOF_FIXEDWITHASSOCIATE : 
    printf("Assoc-%d ", Dof_P->Case.FixedAssociate.NumDof) ; 
    break ;
  case DOF_LINK : 
    printf("Link-");
    Print_DofNumber(Dof_P->Case.Link.Dof);
    break ;
  case DOF_LINKCPLX : 
    printf("LinkCplx-");
    Print_DofNumber(Dof_P->Case.Link.Dof);
    break ;
  default : 
    printf(" ? ") ; 
    break ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------- */
/*  D u m m y  D o f s                                     */
/* ------------------------------------------------------- */

void Dof_GetDummies(struct DefineSystem * DefineSystem_P, struct DofData * DofData_P){

  struct Formulation      * Formulation_P ; 
  struct DefineQuantity   * DefineQuantity_P ;
  struct FunctionSpace    * FunctionSpace_P ;
  struct BasisFunction    * BasisFunction_P ;
  struct GlobalQuantity   * GlobalQuantity_P ;
  struct Dof              * Dof_P ;

  int i, j, k, l, iDof, ii,iit, iNum, iHar;
  int Nbr_Formulation, Index_Formulation ;
  int * DummyDof;
  double DummyFrequency, * Val_Pulsation;

  GetDP_Begin("Dof_GetDummies");

  if (!(Val_Pulsation = Current.DofData->Val_Pulsation))
    Msg(GERROR, "Dof_GetDummies can only be used for harmonic problems");

  DummyDof = DofData_P->DummyDof = (int *)Malloc(DofData_P->NbrDof*sizeof(int)); 
  for (iDof = 0 ; iDof < DofData_P->NbrDof ; iDof++) DummyDof[iDof]=0;

  Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ;

  for (i = 0 ; i < Nbr_Formulation ; i++) {
    List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ;
    Formulation_P = (struct Formulation*)
      List_Pointer(Problem_S.Formulation, Index_Formulation) ;
    for (j = 0 ; j < List_Nbr(Formulation_P->DefineQuantity) ; j++) {
      DefineQuantity_P = (struct DefineQuantity*)
	List_Pointer(Formulation_P->DefineQuantity, j) ;
      for (l = 0 ; l < List_Nbr(DefineQuantity_P->DummyFrequency) ; l++) {
	DummyFrequency = *(double *)List_Pointer(DefineQuantity_P->DummyFrequency, l) ;

	iHar=-1;
	for (k = 0 ; k < Current.NbrHar/2 ; k++)      
	  if (fabs (Val_Pulsation[k]-TWO_PI*DummyFrequency) <= 1e-10 * Val_Pulsation[k]) {
	    iHar = 2*k; break;
	  }
	if(iHar>=0) {
	  FunctionSpace_P = (struct FunctionSpace*)
	    List_Pointer(Problem_S.FunctionSpace, DefineQuantity_P->FunctionSpaceIndex) ;
	  
	  for (k = 0 ; k < List_Nbr(FunctionSpace_P->BasisFunction) ; k++) {
	    BasisFunction_P = (struct BasisFunction *)
	      List_Pointer(FunctionSpace_P->BasisFunction, k) ;
	    iNum = ((struct BasisFunction *)BasisFunction_P)->Num;
	    ii=iit=0;
	    for (iDof = 0 ; iDof < List_Nbr(DofData_P->DofList) ; iDof++) { 
	      Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
	      if (Dof_P->Type == DOF_UNKNOWN && Dof_P->NumType == iNum) {
		iit++;
		if (Dof_P->Harmonic == iHar || Dof_P->Harmonic == iHar+1) { 
		  DummyDof[Dof_P->Case.Unknown.NumDof-1]=1; ii++;
		}
	      }
	    }
	    Msg(INFO, "Freq %e (%d/%d) Form %d  Quant %d  Basis %d  #dummies %d/%d", 
		Val_Pulsation[iHar/2]/TWO_PI, iHar/2, Current.NbrHar/2, 
		i, j, ((struct BasisFunction *)BasisFunction_P)->Num, ii, iit) ;
	  }
	  
	  for (k = 0 ; k < List_Nbr(FunctionSpace_P->GlobalQuantity) ; k++) {
	    GlobalQuantity_P = (struct GlobalQuantity *)
	      List_Pointer(FunctionSpace_P->GlobalQuantity, k) ;
	    iNum = ((struct GlobalQuantity *)GlobalQuantity_P)->Num;
	    ii=iit=0;    
	    for (iDof = 0 ; iDof < List_Nbr(DofData_P->DofList) ; iDof++) { 
	      Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
	      if (Dof_P->Type == DOF_UNKNOWN && Dof_P->NumType == iNum) {
		iit++;
		if (Dof_P->Harmonic == iHar || Dof_P->Harmonic == iHar+1) {  
		  DummyDof[Dof_P->Case.Unknown.NumDof-1]=1; ii++;
		}
	      }
	    }
	    Msg(INFO, "Freq %e (%d/%d) Form %d  Quant %d  Global %d  #dummies %d/%d", 
		Val_Pulsation[iHar/2]/TWO_PI, iHar/2, Current.NbrHar/2, 
		i, j, ((struct GlobalQuantity *)GlobalQuantity_P)->Num, ii, iit) ;
	  }
	    
	}   /*  end DummyFrequency in DofData */ 
      }   /* end DummyFrequency in Quantity */
    }   /* end Quantity */
  }   /* end Formulation */

  i=0;
  for (iDof = 0 ; iDof < DofData_P->NbrDof ; iDof++) {
    if(DummyDof[iDof]) i++;
    /*    
    Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, iDof) ;
    Msg(INFO, "Dof Num iHar, Entity %d %d %d",
	iDof, Dof_P->NumType, Dof_P->Harmonic, Dof_P->Entity);
    */
  }

  Msg(INFO, "Total %d Dummies %d", DofData_P->NbrDof,i) ;
  
  /*  
  Dof_OpenFile(DOF_PRE, "hallopp", "w+") ; 
  Dof_WriteFilePRE(DofData_P) ;
  Dof_CloseFile(DOF_PRE) ; 
  */

  GetDP_End ;
}