get_constraintofelement.c

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

  }  /* for i_Constraint ... */

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G e t _ P r e R e s o l u t i o n F o r C o n s t r a i n t             */
/* ------------------------------------------------------------------------ */

void  Get_PreResolutionForConstraint(struct ConstraintInFS * Constraint_P,
				     int * Index_TimeFunction) {

  struct PreResolutionInfo  PreResolutionInfo_S ;

  int  fcmp_Resolution_Name(const void * a, const void * b) ;

  GetDP_Begin("Get_PreResolutionForConstraint");

  *Index_TimeFunction = Constraint_P->ConstraintPerRegion->TimeFunctionIndex ;

  if (Constraint_P->Active.ResolutionIndex < 0)
    if ((Constraint_P->Active.ResolutionIndex = 
	 List_ISearchSeq(Problem_S.Resolution,
			 Constraint_P->ConstraintPerRegion->
			 Case.Solve.ResolutionName, fcmp_Resolution_Name)) < 0) {
      Msg(GERROR, "Unknown ResolutionName '%s' in Constraint",
	  Constraint_P->ConstraintPerRegion->Case.Solve.ResolutionName) ;
    }
  if(List_ISearchSeq(PreResolutionIndex_L, &Constraint_P->Active.ResolutionIndex, 
		     fcmp_int) < 0) {
    PreResolutionInfo_S.Index = Constraint_P->Active.ResolutionIndex ;
    PreResolutionInfo_S.Type  = PR_CONSTRAINT ;
    List_Add(PreResolutionIndex_L, &PreResolutionInfo_S) ;
    Msg(INFO, "  Adding Resolution '%s' for Pre-Resolution (Constraint)",
	Constraint_P->ConstraintPerRegion->Case.Solve.ResolutionName) ;
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */

/* ------------------------------------------------------------------------ */
/*  G e t _ L i n k F o r C o n s t r a i n t   &   C o                     */
/* ------------------------------------------------------------------------ */

struct TwoIntOneDouble { int  Int1, Int2 ; double Double, Double2 ; } ;


void  Get_LinkForConstraint(struct ConstraintInFS * Constraint_P,
			    int Num_Entity,
			    int * CodeEntity_Link, double Value[]) {

  struct TwoIntOneDouble  * TwoIntOneDouble_P ;
  List_T * Couples_L ;

  void  Generate_Link(struct ConstraintInFS * Constraint_P, int Flag_New) ;

  GetDP_Begin("Get_LinkForConstraint");

  if (!Constraint_P->Active.Active)
    Generate_Link(Constraint_P, 1) ;
  else if (Constraint_P->Active.Active->TimeStep != (int)Current.TimeStep)
    Generate_Link(Constraint_P, 0) ;
  else if (Constraint_P->Active.Active->SubTimeStep != Current.SubTimeStep)
    Generate_Link(Constraint_P, 0) ; /* +++ */


  TwoIntOneDouble_P = (struct TwoIntOneDouble *)
    ((Couples_L = Constraint_P->Active.Active->Case.Link.Couples)?
     List_PQuery(Couples_L, &Num_Entity, fcmp_absint) : NULL) ;
  /* old
    List_PQuery(Constraint_P->Active.Active->Case.Link.Couples,
		&Num_Entity, fcmp_absint) ;
  if (!TwoIntOneDouble_P)  Msg(GERROR, "Constraint Link: bad definition") ;
  */

  if (TwoIntOneDouble_P) {
    *CodeEntity_Link = abs(TwoIntOneDouble_P->Int2) ;
    Value[0] = TwoIntOneDouble_P->Double ;
    if (TwoIntOneDouble_P->Int1 < 0)  Value[0] *= -1. ;
    Value[1] = TwoIntOneDouble_P->Double2 ;  /* LinkCplx */
    if (TwoIntOneDouble_P->Int1 < 0)  Value[1] *= -1. ;
  }

  GetDP_End ;
}


/* Data... */

struct NodeXYZ { int NumNode ; double x, y, z ; } ;
struct EdgeNN { int NumEdge ; int Node1, Node2 ; double Coef, Coef2 ; } ;
void  Generate_LinkNodes(struct ConstraintInFS * Constraint_P,
			 List_T * ExtendedList_L, List_T * ExtendedSuppList_L,
			 struct Group * RegionRef_P, struct Group * SubRegionRef_P,
			 int Index_Filter, int Index_Function, int Index_Coef,
			 List_T * Couples_L) ;
void  Generate_LinkEdges(struct ConstraintInFS * Constraint_P,
			 struct Group * Group_P,
			 struct Group * RegionRef_P, struct Group * SubRegionRef_P,
			 List_T * Couples_L) ;
int fcmp_XYZ(const void * a, const void * b) ;
int fcmp_NN(const void * a, const void * b) ;
void  Generate_LinkRegions(struct ConstraintInFS * Constraint_P,
			   List_T * Region_L, List_T * RegionRef_L,
			   int Index_Coef, List_T * Couples_L) ;
void  Generate_ElementaryEntities_EdgeNN
  (List_T * InitialList, List_T ** ExtendedList, int Type_Entity) ;

/* ----- */


void  Generate_Link(struct ConstraintInFS * Constraint_P, int Flag_New) {

  struct ConstraintActive * Active ;
  struct Group  * Group_P, * RegionRef_P, * SubRegionRef_P ;
  int  Nbr_Entity ;

  GetDP_Begin("Generate_Link");

  Msg(DEBUG2, "C o n s t r a i n t   ( L i n k )\n") ;

  if (Flag_New)
    Constraint_P->Active.Active =
      (struct ConstraintActive *)Malloc(sizeof(struct ConstraintActive)) ;

  Active = Constraint_P->Active.Active ;
  Active->TimeStep = (int)Current.TimeStep ;
  Active->SubTimeStep = Current.SubTimeStep ;

  Group_P = (struct Group*)
    List_Pointer(Problem_S.Group, Constraint_P->EntityIndex) ;
  RegionRef_P = (struct Group*)
    List_Pointer(Problem_S.Group,
		 Constraint_P->ConstraintPerRegion->Case.Link.RegionRefIndex) ;
  SubRegionRef_P =
    (Constraint_P->ConstraintPerRegion->Case.Link.SubRegionRefIndex >= 0)?
    (struct Group*)
    List_Pointer(Problem_S.Group,
		 Constraint_P->ConstraintPerRegion->Case.Link.SubRegionRefIndex) :
    NULL ;

  if (Group_P->FunctionType == REGION)
    Nbr_Entity = List_Nbr(Group_P->InitialList) ;
  else
    Nbr_Entity = List_Nbr(Group_P->ExtendedList) ;

  if (Nbr_Entity) {
    if (Flag_New)
      Active->Case.Link.Couples =
	List_Create(Nbr_Entity, 1, sizeof(struct TwoIntOneDouble)) ;
    else
      List_Reset(Active->Case.Link.Couples) ;
  }
  else {
    Active->Case.Link.Couples = NULL ;
    GetDP_End ;
  }

  switch (Group_P->FunctionType) {
  case NODESOF :
    Generate_LinkNodes(Constraint_P,
		       Group_P->ExtendedList, Group_P->ExtendedSuppList,
		       RegionRef_P, SubRegionRef_P,
		       Constraint_P->ConstraintPerRegion->Case.Link.FilterIndex,
		       Constraint_P->ConstraintPerRegion->Case.Link.FunctionIndex,
		       Constraint_P->ConstraintPerRegion->Case.Link.CoefIndex,
		       Active->Case.Link.Couples) ;
    break ;
  case EDGESOF :
    Generate_LinkEdges(Constraint_P, Group_P,
		       RegionRef_P, SubRegionRef_P,
		       Active->Case.Link.Couples) ;
    break ;
  case FACETSOF :
    Msg(GERROR, "Link not yet implemented for FACETSOF") ;
    break ;
  case REGION :
    Generate_LinkRegions(Constraint_P,
			 Group_P->InitialList, RegionRef_P->InitialList,
			 Constraint_P->ConstraintPerRegion->Case.Link.CoefIndex,
			 Active->Case.Link.Couples) ;
    break ;
  default :
    Msg(GERROR, "Bad function type for Constraint Link: %d", Group_P->FunctionType) ;
    break ;
  }

  GetDP_End ;
}


/*  G e n e r a t e _ L i n k N o d e s  */

void  Generate_LinkNodes(struct ConstraintInFS * Constraint_P,
			 List_T * ExtendedList_L, List_T * ExtendedSuppList_L,
			 struct Group * RegionRef_P, struct Group * SubRegionRef_P,
			 int Index_Filter, int Index_Function, int Index_Coef,
			 List_T * Couples_L) {

  int  Nbr_Entity, i, Nbr_EntityRef, Flag_Filter ;
  double TOL=Current.GeoData->CharacteristicLength * 1.e-8;
  struct TwoIntOneDouble  TwoIntOneDouble ;
  struct NodeXYZ  NodeXYZ, NodeXYZRef ;
  List_T  * NodeXYZ_L, * NodeXYZRef_L ;
  List_T  * ExtendedListRef_L, * ExtendedSuppListRef_L ;
  struct Value  Value ;

  GetDP_Begin("Generate_LinkNodes");

  /* Nodes with Constraint */

  Nbr_Entity = List_Nbr(ExtendedList_L) ;
  NodeXYZ_L = List_Create(Nbr_Entity, 1, sizeof(struct NodeXYZ)) ;
  for (i = 0 ; i < Nbr_Entity ; i++) {
    List_Read(ExtendedList_L, i, &NodeXYZ.NumNode) ;
    if (!(ExtendedSuppList_L &&
	  List_Search(ExtendedSuppList_L, &NodeXYZ.NumNode, fcmp_int))) {
      Geo_GetNodesCoordinates( 1, &NodeXYZ.NumNode,
			       &Current.x, &Current.y, &Current.z) ;
      Get_ValueOfExpressionByIndex(Index_Function, NULL, 0., 0., 0., &Value) ;

      Current.x = Value.Val[0] ; Current.y = Value.Val[1] ;
      Current.z = Value.Val[2] ;
      if (Index_Filter < 0)  Flag_Filter = 1 ;
      else {
	Get_ValueOfExpressionByIndex(Index_Filter, NULL, 0., 0., 0., &Value) ;
	Flag_Filter = (int)Value.Val[0] ;
      }

      if (Flag_Filter) {
	NodeXYZ.x = Current.x ; NodeXYZ.y = Current.y ; NodeXYZ.z = Current.z ;
	List_Add(NodeXYZ_L, &NodeXYZ) ;
      }
    }
  }
  Nbr_Entity = List_Nbr(NodeXYZ_L) ;

  /* Nodes of reference (Link) */

  Generate_ElementaryEntities
    (RegionRef_P->InitialList, &ExtendedListRef_L, NODESOF) ;
  if (SubRegionRef_P)
    Generate_ElementaryEntities
      (SubRegionRef_P->InitialList, &ExtendedSuppListRef_L, NODESOF) ;
  else
    ExtendedSuppListRef_L = NULL ;

  Nbr_EntityRef = List_Nbr(ExtendedListRef_L) ;
  NodeXYZRef_L = List_Create(Nbr_EntityRef, 1, sizeof(struct NodeXYZ)) ;
  for (i = 0 ; i < Nbr_EntityRef ; i++) {
    List_Read(ExtendedListRef_L, i, &NodeXYZRef.NumNode) ;
    if (!(ExtendedSuppListRef_L &&
	  List_Search(ExtendedSuppListRef_L, &NodeXYZRef.NumNode, fcmp_int))) {
      Geo_GetNodesCoordinates( 1, &NodeXYZRef.NumNode,
			       &Current.x, &Current.y, &Current.z) ;
      /*
      Get_ValueOfExpressionByIndex(Index_Function, NULL, 0., 0., 0., &Value) ;

      Current.x = Value.Val[0] ; Current.y = Value.Val[1] ;
      Current.z = Value.Val[2] ;
      */
      if (Index_Filter < 0)  Flag_Filter = 1 ;
      else {
	Get_ValueOfExpressionByIndex(Index_Filter, NULL, 0., 0., 0., &Value) ;
	Flag_Filter = (int)Value.Val[0] ;
      }

      if (Flag_Filter) {
	NodeXYZRef.x = Current.x ; NodeXYZRef.y = Current.y ;
	NodeXYZRef.z = Current.z ;
	List_Add(NodeXYZRef_L, &NodeXYZRef) ;
      }
    }
  }
  Nbr_EntityRef = List_Nbr(NodeXYZRef_L) ;

  /*
  Msg(DEBUG2, "Before sorting\n") ;
  Msg(DEBUG2, "- Other\n") ;
  for (i = 0 ; i < Nbr_Entity ; i++) {
    List_Read(NodeXYZ_L, i, &NodeXYZ) ;
    Msg(DEBUG2, "%d -> %d: %e %e :: %e\n",
        i, NodeXYZ.NumNode, NodeXYZ.x, NodeXYZ.y,
	atan2(NodeXYZ.y,NodeXYZ.x)/3.1415926535897*180.) ;
  }
  Msg(DEBUG2, "- Reference (after rotation)\n") ;
  for (i = 0 ; i < Nbr_EntityRef ; i++) {
    List_Read(NodeXYZRef_L, i, &NodeXYZ) ;
    Msg(DEBUG2, "%d -> %d: %e %e :: %e\n",
	i, NodeXYZ.NumNode, NodeXYZ.x, NodeXYZ.y,
	atan2(NodeXYZ.y,NodeXYZ.x)/3.1415926535897*180.) ;
  }
  Msg(DEBUG2, "\n") ;
  */

  List_Sort(NodeXYZ_L   , fcmp_XYZ) ;
  List_Sort(NodeXYZRef_L, fcmp_XYZ) ;

  Msg(DEBUG2, "After sorting\n") ;
  Msg(DEBUG2, "- Other\n") ;
  for (i = 0 ; i < Nbr_Entity ; i++) {
    List_Read(NodeXYZ_L, i, &NodeXYZ) ;
    Msg(DEBUG2, "%d -> %d: %e %e %e :: %e\n",
	    i, NodeXYZ.NumNode, NodeXYZ.x, NodeXYZ.y, NodeXYZ.z,
	    atan2(NodeXYZ.y,NodeXYZ.x)/3.1415926535897*180.) ;
  }
  Msg(DEBUG2, "- Reference (after rotation)\n") ;
  for (i = 0 ; i < Nbr_EntityRef ; i++) {
    List_Read(NodeXYZRef_L, i, &NodeXYZ) ;
    Msg(DEBUG2, "%d -> %d: %e %e %e :: %e\n",
	    i, NodeXYZ.NumNode, NodeXYZ.x, NodeXYZ.y, NodeXYZ.z,
	    atan2(NodeXYZ.y,NodeXYZ.x)/3.1415926535897*180.) ;
  }

  if (Nbr_EntityRef != Nbr_Entity)
    Msg(GERROR, "Constraint Link: bad correspondance of number of Nodes (%d, %d)",
	Nbr_Entity, Nbr_EntityRef) ;

  Msg(DEBUG2, "==> List of link for nodes\n") ;

  for (i = 0 ; i < Nbr_Entity ; i++) {
    List_Read(NodeXYZ_L, i, &NodeXYZ) ;
    List_Read(NodeXYZRef_L, i, &NodeXYZRef) ;

    /* Attention: tolerance !!! */
    if ((fabs(NodeXYZ.x-NodeXYZRef.x) > TOL) ||
	(fabs(NodeXYZ.y-NodeXYZRef.y) > TOL) ||
	(fabs(NodeXYZ.z-NodeXYZRef.z) > TOL))
      Msg(GERROR, "Constraint Link: bad correspondance of Nodes (%d, %d)"
	  " (%e %e %e)",
	  NodeXYZ.NumNode, NodeXYZRef.NumNode,
	  fabs(NodeXYZ.x-NodeXYZRef.x), fabs(NodeXYZ.y-NodeXYZRef.y),
	  fabs(NodeXYZ.z-NodeXYZRef.z)) ;

    TwoIntOneDouble.Int1 = NodeXYZ.NumNode ;