gf_laplacexform.c

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

    
    switch (Element->ElementSource->Type) {
      
    case POINT :
      Val->Type = VECTOR ;
      
      if (Element->Num == Element->ElementSource->Num) {
	Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
	GetDP_End ;
      }
      
      xxs = x - Element->ElementSource->x[0] ;
      yys = y - Element->ElementSource->y[0] ;
      r2 = SQU(xxs)+SQU(yys) ;
      
      if (r2 > EPSILON2) {
	Val->Val[0] = - ONE_OVER_TWO_PI * xxs / r2 ;
	Val->Val[1] = - ONE_OVER_TWO_PI * yys / r2 ;
	Val->Val[2] = 0. ;
      }
      else {
	Val->Val[0] = Val->Val[1] = Val->Val[2] = 0. ;
      }
      break ;
      
    case TRIANGLE :
      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;

      if(xFunctionBF == NULL) {
	cglt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Type = VECTOR ;
	Val->Val[0] = valr ;
	Val->Val[1] = vali ;
	Val->Val[2] = 0. ;
      }
      if(xFunctionBF == (CAST)BF_GradNode) {
	DetJ = (xs[2]-xs[0]) * (ys[1]-ys[0]) - (xs[1]-xs[0]) * (ys[2]-ys[0]) ;
	switch(NumEntity){
	case 1 : mx = (ys[2]-ys[1])/DetJ ; my = (xs[1]-xs[2])/DetJ ; break;
	case 2 : mx = (ys[0]-ys[2])/DetJ ; my = (xs[2]-xs[0])/DetJ ; break;
	case 3 : mx = (ys[1]-ys[0])/DetJ ; my = (xs[0]-xs[1])/DetJ ; break;
	}
	cglt2d_(&x,&y,&xs[0],&xs[1],&xs[2],&ys[0],&ys[1],&ys[2],&valr,&vali);
	Val->Type = SCALAR ;
	Val->Val[0] = my*valr - mx*vali ; 
      }
      else{
	Msg(GERROR, "Unknown Basis Function Type for 'GF_GradLaplacexForm'");
      }
      break ;

    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    break ;
    
    
  case _3D :
    
    switch (Element->ElementSource->Type) {
      
    case LINE :
      Val->Type = VECTOR ;

      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      zs[0] = Element->ElementSource->z[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      zs[1] = Element->ElementSource->z[1] ;
      
      a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) + SQU(zs[0]-zs[1]) ;
      b = 2. * ((x-xs[0])*(xs[0]-xs[1]) + 
		(y-ys[0])*(ys[0]-ys[1]) + 
		(z-zs[0])*(zs[0]-zs[1])) ;
      c = SQU(x-xs[0]) + SQU(y-ys[0]) + SQU(z-zs[0]) + SQU(RADIUS) ;
      
      I1 = 2./(4.*a*c-b*b) * ( (2.*a+b)/sqrt(a+b+c) - b/sqrt(c) ) ;
      I2 = 2./(-4.*a*c+b*b) * ( (2.*c+b)/sqrt(a+b+c) - 2.*sqrt(c) ) ;
      a2 = sqrt(a) ;
      
      Val->Val[0] = ONE_OVER_FOUR_PI * ( (xs[0]-x) * I1 + (xs[1]-xs[0]) * I2 ) * a2 ;
      Val->Val[1] = ONE_OVER_FOUR_PI * ( (ys[0]-y) * I1 + (ys[1]-ys[0]) * I2 ) * a2 ;
      Val->Val[2] = ONE_OVER_FOUR_PI * ( (zs[0]-z) * I1 + (zs[1]-zs[0]) * I2 ) * a2 ;
      break ;

    case TRIANGLE : 
      Val->Type = VECTOR ;

      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      zs[0] = Element->ElementSource->z[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      zs[1] = Element->ElementSource->z[1] ;
      xs[2] = Element->ElementSource->x[2] ; ys[2] = Element->ElementSource->y[2] ;
      zs[2] = Element->ElementSource->z[2] ;
      

      if(xFunctionBF == (CAST)BF_Volume) Type_Int = 1 ;
      if(xFunctionBF == (CAST)BF_Node)   Type_Int = 2 ;

      /* triangle side lengths */
      a = sqrt(SQU(xs[1]-xs[0]) + SQU(ys[1]-ys[0]) +  SQU(zs[1]-zs[0]));
      b = sqrt(SQU(xs[2]-xs[1]) + SQU(ys[2]-ys[1]) +  SQU(zs[2]-zs[1]));
      c = sqrt(SQU(xs[2]-xs[0]) + SQU(ys[2]-ys[0]) +  SQU(zs[2]-zs[0]));
            
      /* local system (u,v,w) centered at (xs[0],ys[0],zs[0]) */
      u[0] = (xs[1]-xs[0])/a;
      u[1] = (ys[1]-ys[0])/a;
      u[2] = (zs[1]-zs[0])/a;
      
      /* triangle normal */
      Geo_CreateNormal(Element->ElementSource->Type,xs,ys,zs,n);

      v[0] = n[1]*u[2]-n[2]*u[1];
      v[1] = n[2]*u[0]-n[0]*u[2];
      v[2] = n[0]*u[1]-n[1]*u[0];
	
      u2 = (xs[2]-xs[0])*u[0] + (ys[2]-ys[0])*u[1] + (zs[2]-zs[0])*u[2]; /* u2 coordinate */
      v2 = (xs[2]-xs[0])*v[0] + (ys[2]-ys[0])*v[1] + (zs[2]-zs[0])*v[2];  /* triangle height, v2 coordinate*/

      /* local coordinates of the observation point (xl, yl, zl)*/ 
      xl = u[0] * (x-xs[0]) +  u[1] * (y-ys[0]) + u[2] * (z-zs[0]);
      yl = v[0] * (x-xs[0]) +  v[1] * (y-ys[0]) + v[2] * (z-zs[0]);
      zl = n[0] * (x-xs[0]) +  n[1] * (y-ys[0]) + n[2] * (z-zs[0]);

      area = a * v2/2 ;/* Triangle area */
      if (!zl) zl = sqrt(area) * 1e-15 ;

      s0m = -( (a-xl) * (a-u2) + yl*v2 ) / b;
      s0p = s0m + b;
      s1p = ( xl * u2 + yl * v2 ) / c;
      s1m = s1p - c;
      s2m = - xl;
      s2p = a - xl;
  
      /* distance observation point projection on triangle plane to triangle local vertices*/
      t00 = (yl * (u2-a) + v2 * (a-xl)) / b;
      t10 = (xl * v2 - yl * u2) / c;
      t20 = yl;

      t0m_2 = ((a-xl)*(a-xl) + yl*yl); 
      t0p_2 = ((u2-xl)*(u2-xl) + (v2-yl)*(v2-yl)); 
      t1p_2 = (xl*xl + yl*yl);
      
      /* minimum distances^2 from the observation point to each triangle side*/
      zl_2 = SQU(zl) ;
      r00_2 = SQU(t00) + zl_2 ;
      r10_2 = SQU(t10) + zl_2 ;
      r20_2 = SQU(t20) + zl_2 ;
     
      r00 = sqrt(r00_2);
      r10 = sqrt(r10_2);
      r20 = sqrt(r20_2);

      /* distances from observation point to the vertices*/
      r0p = sqrt(t0p_2 + zl_2);      
      r0m = sqrt(t0m_2 + zl_2);
      r1p = sqrt(t1p_2 + zl_2);
      
      EPS = EPSILON*(fabs(s0m)+fabs(s0p));

      B0 = (r00 <= EPS) ? 0. : atan(t00*s0p/(r00_2+fabs(zl)*r0p))-atan(t00*s0m/(r00_2+fabs(zl)*r0m));
      f20 = ((r0m + s0m) <= EPS) ? log(s0m/s0p) : log((r0p + s0p) / (r0m + s0m)) ;

      EPS = EPSILON*(fabs(s1m)+fabs(s1p)) ; 
      B1 = (r10 <=EPS) ? 0. : atan(t10*s1p/(r10_2+fabs(zl)*r1p))-atan(t10*s1m/(r10_2+fabs(zl)*r0p));
      f21 = ((r0p + s1m)<=EPS) ? log(s1m/s1p) : log((r1p + s1p) / (r0p + s1m));

      EPS = EPSILON*(fabs(s2m)+fabs(s2p)) ; 
      B2 =  (r20 <= EPS) ? 0. : atan(t20*s2p/(r20_2+fabs(zl)*r0m))-atan(t20*s2m/(r20_2+fabs(zl)*r1p));
      f22 = ((r1p + s2m)< EPS) ? log(s2m/s2p): log((r0m + s2p) / (r1p + s2m));
      
      B = B0 + B1 + B2 ; 

      s0[0] = (xs[2] - xs[1])/b ;
      s0[1] = (ys[2] - ys[1])/b ;
      s0[2] = (zs[2] - zs[1])/b ;

      s1[0] = (xs[0] - xs[2])/c ;
      s1[1] = (ys[0] - ys[2])/c ;
      s1[2] = (zs[0] - zs[2])/c ;

      m0[0] = s0[1] * n[2] - s0[2]* n[1] ;
      m0[1] = s0[2] * n[0] - s0[0]* n[2] ;
      m0[2] = s0[0] * n[1] - s0[1]* n[0] ;

      m1[0] = s1[1]* n[2] - s1[2] * n[1] ;
      m1[1] = s1[2]* n[0] - s1[0] * n[2] ;
      m1[2] = s1[0]* n[1] - s1[1] * n[0] ;
      
      m2[0] = u[1] * n[2] - u[2]* n[1] ;
      m2[1] = u[2] * n[0] - u[0]* n[2] ;
      m2[2] = u[0] * n[1] - u[1]* n[0] ;

      /* Grad(1/r) integral solution*/
      I[0] = -n[0] * THESIGN(zl) * B - (m0[0]*f20 + m1[0]*f21 + m2[0]*f22) ; 
      I[1] = -n[1] * THESIGN(zl) * B - (m0[1]*f20 + m1[1]*f21 + m2[1]*f22) ;
      I[2] = -n[2] * THESIGN(zl) * B - (m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

      switch ( Type_Int ){
      case 1 : /* BF_Volume */	
	Val->Val[0] = I[0]/area ; 	
	Val->Val[1] = I[1]/area ; 	
	Val->Val[2] = I[2]/area ; 	
	break;
    
      case 2 : /* BF_Node */
	if (!v2 ) 	
	  Msg(GERROR,"1/0 in GF_LaplacexForm (case _3D TRIANGLE) v2 %e", v2);

	f0[0] = s0[0] * t00 * f20 - m0[0]*(r0p-r0m) ; /* fi */
	f0[1] = s0[1] * t00 * f20 - m0[1]*(r0p-r0m) ;
	f0[2] = s0[2] * t00 * f20 - m0[2]*(r0p-r0m) ;

	f1[0] = s1[0] * t10 * f21 - m1[0]*(r1p-r0p) ;
	f1[1] = s1[1] * t10 * f21 - m1[1]*(r1p-r0p) ;
	f1[2] = s1[2] * t10 * f21 - m1[2]*(r1p-r0p) ;

	f2[0] = u[0] * t20 * f22 - m2[0]*(r0m-r1p) ;
	f2[1] = u[1] * t20 * f22 - m2[1]*(r0m-r1p) ;
	f2[2] = u[2] * t20 * f22 - m2[2]*(r0m-r1p) ;  

	umf2i = u[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) +
	        u[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) + 
	        u[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

	us0 = u[0] * s0[0] + u[1] * s0[1] + u[2] * s0[2] ; 
	us1 = u[0] * s1[0] + u[1] * s1[1] + u[2] * s1[2] ; 
	us2 = u[0] *  u[0] + u[1] *  u[1] + u[2] *  u[2] ; 

	vmf2i = v[0]*(m0[0]*f20 + m1[0]*f21 + m2[0]*f22) + 
	        v[1]*(m0[1]*f20 + m1[1]*f21 + m2[1]*f22) + 
	        v[2]*(m0[2]*f20 + m1[2]*f21 + m2[2]*f22) ;

	vs0 = v[0] * s0[0] + v[1] * s0[1] + v[2] * s0[2] ; 
	vs1 = v[0] * s1[0] + v[1] * s1[1] + v[2] * s1[2] ; 
	vs2 = v[0] *  u[0] + v[1] *  u[1] + v[2] *  u[2] ; 

	B *= fabs(zl);
	umf2i *= zl ;
	vmf2i *= zl ;

	Iua[0] = n[0] * umf2i - B * u[0] + f0[0] * us0 + f1[0] * us1 + f2[0] * us2 ;
	Iua[1] = n[1] * umf2i - B * u[1] + f0[1] * us0 + f1[1] * us1 + f2[1] * us2 ;
	Iua[2] = n[2] * umf2i - B * u[2] + f0[2] * us0 + f1[2] * us1 + f2[2] * us2 ;
	
	Iva[0] = n[0] * vmf2i - B * v[0] + f0[0] * vs0 + f1[0] * vs1 + f2[0] * vs2 ;
	Iva[1] = n[1] * vmf2i - B * v[1] + f0[1] * vs0 + f1[1] * vs1 + f2[1] * vs2 ;
	Iva[2] = n[2] * vmf2i - B * v[2] + f0[2] * vs0 + f1[2] * vs1 + f2[2] * vs2 ;
	
	switch(NumEntity){
	case 1 :  
	  N10 = 1 - xl/a + (u2/a -1) * yl/v2 ;
	  Val->Val[0] = N10 * I[0]  - Iua[0]/a + (u2/a-1) * Iva[0]/v2 ;
	  Val->Val[1] = N10 * I[1]  - Iua[1]/a + (u2/a-1) * Iva[1]/v2 ;
	  Val->Val[2] = N10 * I[2]  - Iua[2]/a + (u2/a-1) * Iva[2]/v2 ;
	  break;
	case 2 : 
	  N20 = xl/a - u2/a * yl/v2 ;
	  Val->Val[0] = N20 * I[0] + Iua[0]/a - u2/a * Iva[0]/v2 ;
	  Val->Val[1] = N20 * I[1] + Iua[1]/a - u2/a * Iva[1]/v2 ;
	  Val->Val[2] = N20 * I[2] + Iua[2]/a - u2/a * Iva[2]/v2 ;
	  break;
	case 3 :  
	  N30 = yl/v2 ;
	  Val->Val[0] = N30 * I[0] + Iva[0]/v2 ;
	  Val->Val[1] = N30 * I[1] + Iva[1]/v2 ;
	  Val->Val[2] = N30 * I[2] + Iva[2]/v2 ;
	  break;
	}
	break;
      }
      
      Val->Val[0] *= ONE_OVER_FOUR_PI ;
      Val->Val[1] *= ONE_OVER_FOUR_PI ;
      Val->Val[2] *= ONE_OVER_FOUR_PI ;      
      break ;
      
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_GradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));
    }
    break ;
    
  default :
    Msg(GERROR, "Unknown Dimension (%d) for 'GF_GradLaplacexForm'",
	(int)Fct->Para[0]);
    
  }

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G F _ N P x G r a d L a p l a c e x F o r m                             */
/* ------------------------------------------------------------------------ */

void GF_NPxGradLaplacexForm (F_ARG2) {

  double  xs[MAX_NODES], ys[MAX_NODES] ;
  double  xp[MAX_NODES], yp[MAX_NODES], N[3] ;
  int     Type_Int;
  double  a, b, c, d, m, n, Jp, i1, Is, I1=0 ;
  struct Value ValGrad ;

  GetDP_Begin("GF_NPxGradLaplacexForm");

  Val->Type = SCALAR ;
  Val->Val[MAX_DIM] = 0.0 ;
  
  if (Element->Num == Element->ElementSource->Num) {
    Val->Val[0] = 0.0 ;
    GetDP_End ;
  }
  
  switch ((int)Fct->Para[0]) {
    
  case _2D : 

    switch (Element->ElementSource->Type) {

    case LINE :	  
      if (Element->Type != LINE)
	Msg(GERROR, "GF_NPxGradLaplacexForm not ready for mixed geometrical elements");

      xs[0] = Element->ElementSource->x[0] ; ys[0] = Element->ElementSource->y[0] ;
      xs[1] = Element->ElementSource->x[1] ; ys[1] = Element->ElementSource->y[1] ;
      
      if(xFunctionBF == (CAST)BF_Volume) {
	if ((x == xs[0]) && (y == ys[0]))
	  Type_Int = 1 ; 
	else if ((x == xs[1]) && (y == ys[1]))
	  Type_Int = 2 ; 
	else                                    
	  Type_Int = 3 ; 
	
	xp[0] = Element->x[0] ; yp[0] = Element->y[0] ;
	xp[1] = Element->x[1] ; yp[1] = Element->y[1] ;
	
	a = SQU(xs[0]-xs[1]) + SQU(ys[0]-ys[1]) ;
	b = 2. * ((x-xs[0]) * (xs[0]-xs[1]) + (y-ys[0]) * (ys[0]-ys[1])) ;
	c = SQU(x-xs[0]) + SQU(y-ys[0]) ;     	
	d = 4.*a*c - b*b ;
	
	switch (Type_Int) {	  
	case 1 :  
	case 2 :  
	  Msg(GERROR, "Degenerate case not done in 'GF_NPxGradLaplacexForm'");
	  break ;	
	case 3 :  
	  if (fabs(d) < EPSILON2) {
	    I1 = 0.0 ;
	  }
	  else {
	    if(d<0) Msg(GERROR, "Unexpected value in 'GF_NPxGradLaplacexForm'");
	    i1 = sqrt(d) ;
	    Is = 2. / i1 * (atan((2.*a+b)/i1) -  atan(b/i1)) ;
	    Jp = sqrt(SQU(xp[0]-xp[1])+SQU(yp[0]-yp[1])) ;
	    m = ((ys[0]-ys[1]) * (xp[0]-xp[1]) + (xs[0]-xs[1]) * (yp[1]-yp[0])) / Jp ;
	    n = ((yp[1]-yp[0]) * (x-xs[0]) + (xp[0]-xp[1]) * (y-ys[0])) / Jp ; 	
	    I1 = m /(2.*a) * log((a+b)/c+1.) + (n - m*b/(2.*a)) * Is ;	   
	  }
	  break ;
	}
	Val->Val[0] = - ONE_OVER_TWO_PI * I1 ; 
      }
      else {	
	Msg(GERROR, "Unknown Basis Function Type for 'GF_NPxGradLaplacexForm'");
      }
      break ;

    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));     
    }
    break ;

  case _3D:
    switch (Element->ElementSource->Type) {
    case TRIANGLE :
      Geo_CreateNormal(Element->Type,
		       Element->x,Element->y,Element->z, N);

      GF_GradLaplacexForm(Element, Fct, xFunctionBF, NumEntity, x, y, z, &ValGrad) ;

      Val->Val[0] = N[0]*ValGrad.Val[0] + N[1]*ValGrad.Val[1] + N[2]*ValGrad.Val[2] ;
      break ;
    default :
      Msg(GERROR, "Unknown Element Type (%s) for 'GF_NPxGradLaplacexForm'",
	  Get_StringForDefine(Element_Type, Element->ElementSource->Type));  
    }
    break ;

  default :
    Msg(GERROR, "Unknown Dimension (%d) for 'GF_NPxGradLaplacexForm'",
	(int)Fct->Para[0]);
  }

  GetDP_End ;
}



/* ------------------------------------------------------------------------ */
/*  G F _ N S x G r a d L a p l a c e x F o r m                             */
/* ------------------------------------------------------------------------ */

void GF_NSxGradLaplacexForm (F_ARG2) {
    
  GetDP_Begin("GF_NSxGradLaplacexForm");

  Msg(GERROR, "Not done: 'GF_NSxGradLaplacexForm'");
  
  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  G F _ A p p r o x i m a t e L a p l a c e x F o r m                     */
/* ------------------------------------------------------------------------ */

void GF_ApproximateLaplacexForm (F_ARG2) {

  GetDP_Begin("GF_ApproxilateLaplacexForm");

  switch ((int)Fct->Para[1]) {

  case 0 :
    GF_LaplacexForm(Element, Fct, (CAST)BF_Volume, 1, x, y, z, Val);
    break ;

  default :
    Msg(GERROR, "Bad Parameter Value in 'GF_ApproximateLaplacexForm'");
    break;

  }

  GetDP_End ;
}


#undef F_ARG2
#undef CAST