dofdata.c

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

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ D e f i n e I n i t S o l v e D o f                             */
/* ------------------------------------------------------------------------ */

void  Dof_DefineInitSolveDof(int D1, int D2, int NbrHar) {
  struct Dof  Dof ;
  int         k ;

  GetDP_Begin("Dof_DefineInitSolveDof");

  Dof.NumType = D1 ;  Dof.Entity = D2 ; 

  for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
    Dof.Harmonic = k ;
    if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
      Dof.Type = DOF_UNKNOWN_INIT ;
      Dof.Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ;
      Tree_Add(CurrentDofData->DofTree, &Dof) ;
    }
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ D e f i n e L i n k D o f                                       */
/* ------------------------------------------------------------------------ */

void  Dof_DefineLinkDof(int D1, int D2, int NbrHar, double Value[], int D2_Link) {
  struct Dof  Dof ;
  int         k ;

  GetDP_Begin("Dof_DefineLinkDof");

  Dof.NumType = D1 ;  Dof.Entity = D2 ;  
  
  for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
    Dof.Harmonic = k ;
    if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
      Dof.Type = DOF_LINK ;
      Dof.Case.Link.Coef = Value[0] ;
      Dof.Case.Link.EntityRef = D2_Link ;
      Dof.Case.Link.Dof = NULL ;
      Tree_Add(CurrentDofData->DofTree, &Dof) ;
    }
  }

  GetDP_End ;
}

void  Dof_DefineLinkCplxDof(int D1, int D2, int NbrHar, double Value[], int D2_Link) {
  struct Dof  Dof ;
  int         k ;

  GetDP_Begin("Dof_DefineLinkCplxDof");

  Dof.NumType = D1 ;  Dof.Entity = D2 ;  
  
  for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
    Dof.Harmonic = k ;
    if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
      Dof.Type = DOF_LINKCPLX ;
      Dof.Case.Link.Coef = Value[0] ;
      Dof.Case.Link.Coef2 = Value[1] ;
      Dof.Case.Link.EntityRef = D2_Link ;
      Dof.Case.Link.Dof = NULL ;
      Tree_Add(CurrentDofData->DofTree, &Dof) ;
    }
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  D o f _ D e f i n e U n k n o w n D o f                                 */
/* ------------------------------------------------------------------------ */

void  Dof_DefineUnknownDof(int D1, int D2, int NbrHar) {
  struct Dof  Dof ;
  int         k ;

  GetDP_Begin("Dof_DefineUnknownDof");

  Dof.NumType = D1 ;  Dof.Entity = D2 ;

  for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
    Dof.Harmonic = k ;
    if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
      Dof.Type = DOF_UNKNOWN ;
      Dof.Case.Unknown.NumDof = -1 ;
      /* Dof.Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ; */
      Tree_Add(CurrentDofData->DofTree, &Dof) ;
    }
  }

  GetDP_End ;
}

void NumberUnknownDof (void *a, void *b) {
  struct Dof * Dof_P ;
  
  GetDP_Begin("NumberSymmetricalDof");

  Dof_P = (struct Dof *)a ;

  if(Dof_P->Type == DOF_UNKNOWN && Dof_P->Case.Unknown.NumDof == -1)
    Dof_P->Case.Unknown.NumDof = ++(CurrentDofData->NbrDof) ;

  GetDP_End ;
}

void  Dof_NumberUnknownDof(void) {

  GetDP_Begin("Dof_NumberUnknownDof");

  if(CurrentDofData->DofTree)
    Tree_Action(CurrentDofData->DofTree, NumberUnknownDof) ;
  else
    List_Action(CurrentDofData->DofList, NumberUnknownDof) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ D e f i n e A s s o c i a t e D o f                             */
/* ------------------------------------------------------------------------ */

void  Dof_DefineAssociateDof(int E1, int E2, int D1, int D2, int NbrHar) {
  struct Dof  Dof, Equ, * Equ_P ;
  int         k ;

  GetDP_Begin("Dof_DefineAssociateDof");

  Equ.NumType = E1 ;  Equ.Entity = E2 ;

  for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
    Equ.Harmonic = k ;
    if ((Equ_P = (struct Dof*)Tree_PQuery(CurrentDofData->DofTree, &Equ))) {
      switch (Equ_P->Type) {
      case DOF_FIXED :
	Equ_P->Type = DOF_FIXEDWITHASSOCIATE ;
	Equ_P->Case.FixedAssociate.NumDof = ++(CurrentDofData->NbrDof) ;
	Dof.NumType = D1 ; Dof.Entity = D2 ; Dof.Harmonic = k ;
	if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
	  Dof.Type = DOF_UNKNOWN ;
	  Dof.Case.Unknown.NumDof = CurrentDofData->NbrDof ;
	  Tree_Add(CurrentDofData->DofTree, &Dof) ;
	}
	break ;
      case DOF_FIXED_SOLVE :
	Equ_P->Type = DOF_FIXEDWITHASSOCIATE_SOLVE ;
	Equ_P->Case.FixedAssociate.NumDof = ++(CurrentDofData->NbrDof) ;
	Dof.NumType = D1 ; Dof.Entity = D2 ; Dof.Harmonic = k ;
	if (!Tree_PQuery(CurrentDofData->DofTree, &Dof)) {
	  Dof.Type = DOF_UNKNOWN ;
	  Dof.Case.Unknown.NumDof = CurrentDofData->NbrDof ;
	  Tree_Add(CurrentDofData->DofTree, &Dof) ;
	}
	break ;
      case DOF_UNKNOWN :  case DOF_UNKNOWN_INIT :
	Dof_DefineUnknownDof(D1, D2, NbrHar) ;
	break ;
      }
    }
  }

  GetDP_End ;
}


void Dof_GetDof_Four(struct DofData * DofData_P, struct DofData * DofData2_P) {
  int NbrHar2, NbrDof2, i ;
  struct Dof * Dof_P ;

  NbrHar2 = DofData2_P->NbrHar ;
  NbrDof2 = List_Nbr(DofData2_P->DofList) ;

  for (i=0 ; i<NbrDof2 ; i+=NbrHar2) {

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

    printf("i %d, dof %d \n", i/NbrHar2, List_ISearch(DofData_P->DofList, Dof_P, fcmp_Dof)) ; 
	  
  }
}



/* ------------------------------------------------------------------------ */
/*  P R O C E S S I N G   ( C o d e s   i n   L i s t )                     */
/* ------------------------------------------------------------------------ */

/* ------------------------------------------------------------------------ */
/*  D o f _ G e t D o f S t r u c t                                         */
/* ------------------------------------------------------------------------ */

struct Dof  * Dof_GetDofStruct(struct DofData * DofData_P, int D1, int D2, int D3) {
  struct Dof  Dof ;

  GetDP_Begin("Dof_getDofStruct");

  Dof.NumType = D1 ;  Dof.Entity = D2 ;  Dof.Harmonic = D3 ;

  GetDP_Return((struct Dof *)List_PQuery(DofData_P->DofList, &Dof, fcmp_Dof)) ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ U p d a t e L i n k D o f                                       */
/* ------------------------------------------------------------------------ */

void  Dof_UpdateLinkDof(struct Dof *Dof_P, int NbrHar, double Value[], int D2_Link) {
  int         k ;

  GetDP_Begin("Dof_UpdateLinkDof");

  if (Dof_P->Type == DOF_LINK || Dof_P->Type == DOF_LINKCPLX) {
    /*
    fprintf(stderr,"===> %d %d %.16g\n", Dof_P->NumType, Dof_P->Entity, Value[0]) ;
    */
    for(k=0 ; k<NbrHar ; k+=gSCALAR_SIZE){
      Dof_P->Case.Link.Coef = Value[0] ;

      if (Dof_P->Type == DOF_LINKCPLX)
	Dof_P->Case.Link.Coef2 = Value[1] ;

      Dof_P->Case.Link.EntityRef = D2_Link ;
      Dof_P->Case.Link.Dof = NULL ;
    }

  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ A s s e m b l e I n M a t                                       */
/* ------------------------------------------------------------------------ */

void  Dof_AssembleInMat(struct Dof * Equ_P, struct Dof * Dof_P, int NbrHar, 
			double * Val, gMatrix * Mat, gVector * Vec) {
  gScalar tmp, tmp2 ;
  double  valtmp[2], d1, d2 ;

  GetDP_Begin("Dof_AssembleInMat");

  switch (Equ_P->Type) {
  case DOF_UNKNOWN :
  case DOF_FIXEDWITHASSOCIATE :

    switch (Dof_P->Type) {
      
    case DOF_UNKNOWN :
      if(NbrHar==1){
	if(Val[0])
	  LinAlg_AddDoubleInMatrix
	    (Val[0], Mat, 
	     Equ_P->Case.Unknown.NumDof-1, Dof_P->Case.Unknown.NumDof-1) ;
      }
      else
	LinAlg_AddComplexInMatrix
	  (Val[0], Val[1], Mat,
	   Equ_P->Case.Unknown.NumDof-1, Dof_P->Case.Unknown.NumDof-1,
	   (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1,
	   (gSCALAR_SIZE==1)?((Dof_P+1)->Case.Unknown.NumDof-1):-1) ;
      break ;

    case DOF_FIXED :
    case DOF_FIXEDWITHASSOCIATE :
      if(Vec){
	if(NbrHar==1){
	  if(Val[0]){
	    LinAlg_ProdScalarDouble
	      (&Dof_P->Val, 
	       CurrentDofData->CurrentSolution->
	       TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
	       &tmp);
	    LinAlg_ProdScalarDouble(&tmp, -Val[0], &tmp) ;
	    LinAlg_AddScalarInVector(&tmp, Vec, Equ_P->Case.Unknown.NumDof-1) ;
	  }
	}
	else{
	  LinAlg_ProdScalarDouble
	    (&Dof_P->Val, 
	     CurrentDofData->CurrentSolution->
	     TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
	     &tmp);
	  if(gSCALAR_SIZE == 2){
	    LinAlg_ProdScalarComplex(&tmp, -Val[0], -Val[1], &valtmp[0], &valtmp[1]) ;
	  }
	  else{
	    LinAlg_GetDoubleInScalar(&d1, &tmp);
	    LinAlg_ProdScalarDouble
	      (&(Dof_P+1)->Val, 
	       CurrentDofData->CurrentSolution->
	       TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex],
	       &tmp2);
	    LinAlg_GetDoubleInScalar(&d2, &tmp2);
	    valtmp[0] = -d1*Val[0] + d2*Val[1] ;
	    valtmp[1] = -d1*Val[1] - d2*Val[0] ;
	  }
	  LinAlg_AddComplexInVector
	    (valtmp[0], valtmp[1], Vec,
	     Equ_P->Case.Unknown.NumDof-1,
	     (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
	}
      }
      break ;

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

    case DOF_LINKCPLX :
      if(NbrHar==1)
	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_AssembleInMat(Equ_P, Dof_P->Case.Link.Dof, NbrHar, valtmp, Mat, 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_AssembleInMat(Equ_P->Case.Link.Dof, Dof_P, NbrHar, valtmp, Mat, 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_AssembleInMat(Equ_P->Case.Link.Dof, Dof_P, NbrHar, valtmp, Mat, Vec) ;
    break ;

  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  D o f _ A s s e m b l e I n V e c                                       */
/* ------------------------------------------------------------------------ */

void  Dof_AssembleInVec(struct Dof * Equ_P, struct Dof * Dof_P, int NbrHar, 
			double * Val, struct Solution * OtherSolution,
			gVector * Vec0, gVector * Vec) {
  gScalar tmp ;
  double valtmp[2] ;
  double a, b, c, d ;

  GetDP_Begin("Dof_AssembleInVec");

  switch (Equ_P->Type) {
  case DOF_UNKNOWN :
  case DOF_FIXEDWITHASSOCIATE :

    switch (Dof_P->Type) {

    case DOF_UNKNOWN :
      if(NbrHar==1){	
	if(Val[0]){
	  LinAlg_GetDoubleInVector(&a, Vec0, Dof_P->Case.Unknown.NumDof-1) ;
	  a *= Val[0] ;
	  LinAlg_AddDoubleInVector(a, Vec, Equ_P->Case.Unknown.NumDof-1) ;
	}
      }
      else{
	LinAlg_GetComplexInVector(&a, &b, Vec0,
				  Dof_P->Case.Unknown.NumDof-1, 
				  (gSCALAR_SIZE==1)?((Dof_P+1)->Case.Unknown.NumDof-1):-1) ;
	c = a * Val[0] - b * Val[1] ; 
	d = a * Val[1] + b * Val[0] ;
	LinAlg_AddComplexInVector(c, d, Vec,
				  Equ_P->Case.Unknown.NumDof-1, 
				  (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
      }
      break ;

    case DOF_FIXED :
    case DOF_FIXEDWITHASSOCIATE :
      if(NbrHar==1){	
	if(Val[0]){
	  LinAlg_ProdScalarDouble
	    (&Dof_P->Val, 
	     Val[0] * OtherSolution->
	     TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex], 
	     &tmp) ;
	  LinAlg_AddScalarInVector(&tmp, Vec, Equ_P->Case.Unknown.NumDof-1) ;
	}
      }
      else{
	if(gSCALAR_SIZE == 2){
	  LinAlg_ProdScalarComplex
	    (&Dof_P->Val, 
	     Val[0] * OtherSolution->
	     TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex], 
	     Val[1] * OtherSolution->
	     TimeFunctionValues[Dof_P->Case.FixedAssociate.TimeFunctionIndex], 
	     &a, &b) ;
	  LinAlg_AddComplexInVector(a, b, Vec,
				    Equ_P->Case.Unknown.NumDof-1, 
				    (gSCALAR_SIZE==1)?((Equ_P+1)->Case.Unknown.NumDof-1):-1) ;
	}
	else{
	  Msg(GERROR, "Assemby in vectors with more than one harmonic not yet implemented") ;
	}
      }
      break ;

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

    case DOF_LINKCPLX :
      if(NbrHar==1)