wzcoordinates.cxx

有限元学习研究用源代码(老外的),供科研人员参考

#include "wzcoordinates.hxx"
#include <math.h>

void wzCoordinates::dx(wzFloat *d, const wzFloat *yy, wzIndex dir) const
{
  wzFloat y0 = yy[dir], h = 1.e-3*y0+1.0638e-7;
  wzFloat d0[wzPointDim];
  wzIndex i;
  x(d0,yy); 
  ((wzFloat *)yy)[dir] += h;
  x(d,yy);
  for(i=0;i<dim;i++){
    d[i] -= d0[i]; d[i] /= h;
  }
  ((wzFloat *)yy)[dir] = y0;
}

void wzCoordinates::dy(wzFloat *d, const wzFloat *xx, wzIndex dir) const
{
  wzFloat y0 = xx[dir], h = 1.e-3*y0+1.0638e-7;
  wzFloat d0[wzPointDim];
  wzIndex i;
  y(d0,xx); 
  ((wzFloat *)xx)[dir] += h;
  y(d,xx);
  for(i=0;i<dim;i++){
    d[i] -= d0[i]; d[i] /= h;
  }
  ((wzFloat *)xx)[dir] = y0;
}

const wzFloat	wzPi = 4*atan(1);

void wzCoordinates::test()
{
  int i,j,k,l;
  wzFloat xx[3],yy[3],zz[3],beg[3],end[3],d[3],dd,d0,d0max;
  d0max = 0;
  for(i=0;i<dim;i++){
    if(ymin[i] == -wzInfty) beg[i] = -2; else beg[i] = ymin[i]+0.1;
    if(ymax[i] ==  wzInfty) end[i] = +2; else end[i] = ymax[i]-0.1;
    d[i] = (end[i]-beg[i])/11;
    d0max += 0.01*d[i]*d[i];
  }
  xx[0] = beg[0];
  for(i=0;i<11;i++){
    if(dim>1){
      xx[1] = beg[1];
      for(j=0;j<11;j++){
	if(dim>2){
	  xx[2] = beg[2];
	  for(k=0;k<11;k++){
	    x(yy,xx); y(zz,yy);d0=0;
	    for(l=0;l<dim;l++){
	      dd = zz[l]-xx[l]; d0 += dd*dd;
	    }
	    wzAssert(d0<d0max);
	  }
	  xx[2] += d[2];
	}else{
	  x(yy,xx); y(zz,yy);d0=0;
	  for(l=0;l<dim;l++){
	    dd = zz[l]-xx[l]; d0 += dd*dd;
	  }
	  wzAssert(d0<d0max);
	}
	xx[1] += d[1];
      }
    }else{
      x(yy,xx); y(zz,yy);d0=0;
      for(l=0;l<dim;l++){
	dd = zz[l]-xx[l]; d0 += dd*dd;
      }
      wzAssert(d0<d0max);
    }
    xx[0] += d[0];
  }
}

wzCoordinates::wzCoordinates(wzIndex d)
{
  dim = d;
  ymin[0]=ymin[1]=ymin[2]= -wzInfty;
  ymax[0]=ymax[1]=ymax[2]=  wzInfty;
}

void wzCoordinates::x1(wzFloat *x) const {}

void wzCoordinates::y1(wzFloat *x) const {}

void wzCoordinates::x(wzFloat *x, const wzFloat *y) const
{
  for(int i=0;i<dim;i++) x[i]=y[i];
}
void wzCoordinates::y(wzFloat *y, const wzFloat *x) const
{
  for(int i=0;i<dim;i++) y[i]=x[i];
}

void wz3DCoordinates::x1(wzFloat *xx) const
{
  wzFloat X[3]={xx[0],xx[1],xx[2]}; x(xx,X);
}
void wz3DCoordinates::y1(wzFloat *xx) const
{
  wzFloat X[3]={xx[0],xx[1],xx[2]}; y(xx,X);
}
void wz2DCoordinates::x1(wzFloat *xx) const
{
  wzFloat X[2]={xx[0],xx[1]}; this->x(xx,X);
}
void wz2DCoordinates::y1(wzFloat *xx) const
{
  wzFloat X[2]={xx[0],xx[1]}; y(xx,X);
}
void wz1DCoordinates::x1(wzFloat *xx) const
{
  wzFloat X[1]={xx[0]}; x(xx,X);
}
void wz1DCoordinates::y1(wzFloat *xx) const
{
  wzFloat X[1]={xx[0]}; y(xx,X);
}

wzPolarCoordinates::wzPolarCoordinates()
{
  Epsilon = 1.e-7;
  ymin[0] = Epsilon;
  ymin[1] = - wzPi;
  ymax[1] =   wzPi;
  test();
}
void wzPolarCoordinates::setMinimalRadius(wzFloat eps)
{
  Epsilon = eps;
  ymin[0] = Epsilon;
}

void wzPolarCoordinates::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat yy = (y[0]>Epsilon)?y[0]:Epsilon;
  x[0] = yy*cos(y[1]);
  x[1] = yy*sin(y[1]);
}
void wzPolarCoordinates::y(wzFloat *y, const wzFloat *x) const
{
  try{
    y[1] = atan2(x[1],x[0]);
  }catch(...){
    y[1] = 0;  // throw wzCoordinateFailure;
  }
  y[0] = sqrt(x[0]*x[0]+x[1]*x[1]);
  wzAssert(y[1]<=ymax[1]);
  wzAssert(y[1]>=ymin[1]);
}

wzEllipticCoordinates::wzEllipticCoordinates(wzFloat a0)
{
  polar = new wzPolarCoordinates;
  shukowski = new wzComplexShukowski;
  ymin[0] = 1;
  ymin[1] = - wzPi;
  ymax[1] =   wzPi;
  test();
}
void wzEllipticCoordinates::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat z[2];
  //  x[0] = (y[0]+(radius/y[0]))*cos(y[1]);
  //  x[1] = (y[0]-(radius/y[0]))*sin(y[1]);
  polar->x(z,y);
  shukowski->x(x,z);
}
void wzEllipticCoordinates::y(wzFloat *y, const wzFloat *x) const
{
  wzFloat z[2];
  shukowski->y(z,x);
  polar->y(y,z);
}


wzComplexExp::wzComplexExp()
{
  ymin[1] = - wzPi;
  ymax[1] =   wzPi;
  test();
}
void wzComplexExp::x(wzFloat *x, const wzFloat *y) const
{
  x[0] = exp(y[0])*cos(y[1]);
  x[1] = exp(y[0])*sin(y[1]);
}
void wzComplexExp::y(wzFloat *y, const wzFloat *x) const
{
  try{
    y[1] = atan2(x[1],x[0]);
  }catch(...){
    y[1] = 0;  // throw wzCoordinateFailure;
  }
  y[0] = log(x[0]*x[0]+x[1]*x[1])/2;
  wzAssert(y[1]<=ymax[1]);
  wzAssert(y[1]>=ymin[1]);
}

wzComplexShukowski::wzComplexShukowski(wzcoordinates b)
:base(b)
{
  setShiftPoints(0,-1,0,0);
  setLambda(0);
}

void wzComplexShukowski::setLambda(wzFloat l)
{
  lambda = l;
}

void wzComplexShukowski::setShiftPoints(wzFloat x0,wzFloat y0,wzFloat x1,wzFloat y1)
{
  ymax[1] =  0;
  m0 = x1;  m1 = y1;
  d0 = x0-x1;  
  d1 = y0-y1;
  wzFloat r2 = d0*d0+d1*d1;
  d0 /= r2; d1 /= r2;
  r2 = d0*d0+d1*d1;
  rr=sqrt(r2);
  e0 = -2*d0/r2;
  e1 = -2*d1/r2;
  c0 = -d1/rr;
  c1 =  d0/rr;
  test();
}

void wzComplexShukowski::x(wzFloat *x, const wzFloat *y) const
{
  double rs,r0,r1,r,xr[2],q0,q1;
  r0 = y[0] - m0; 
  r1 = y[1] - m1;
  rs = r0*d0+r1*d1;
  q0 = r0+rs*e0;
  q1 = r1+rs*e1;
  r = (r0*r0+r1*r1)*rr*rr;
  if(base!=0){
    xr[0] = y[0]+q0/r;
    xr[1] = y[1]+q1/r;
    base->x(x,xr);
  }else{
    x[0] = y[0]+q0/r;
    x[1] = y[1]+q1/r;
  }
}
void wzComplexShukowski::y(wzFloat *y, const wzFloat *x) const
{
  double u0,u1,a0,a1,sq,s0,s1;
  double r0,r1,xr[2],v0,v1;
  if(base!=0){
    base->y(xr,x);
    r0 = xr[0] - m0; 
    r1 = xr[1] - m1;
  }else{
    r0 = x[0] - m0; 
    r1 = x[1] - m1;
  }
  //  u0 = x[0]/2; u1=x[1]/2;
  u0 = (r0*c0+r1*c1)*rr/2;
  u1 = (r1*c0-r0*c1)*rr/2;
  a0 = (u0*u0-u1*u1-1)/2;
  a1 = u0*u1;
  sq = sqrt(a0*a0+a1*a1);
  s0 = sqrt(sq+a0);
  s1 = sqrt(sq-a0);
  if(a1<0) s0 = -s0;
  if(u0*s0<(lambda-u1)*s1){
    v0 = u0 - s0;
    v1 = u1 - s1;
  }else{
    v0 = u0 + s0;
    v1 = u1 + s1;
  }
  y[0] = m0 + (v0*c0 - v1*c1)/rr;
  y[1] = m1 + (v0*c1 + v1*c0)/rr;
  //  wzAssert(y[1]<=ymax[1]);
}

void wz3DContinuation::y(wzFloat *y, const wzFloat *x) const
{
  original->y(y,x);
  if(scale==0){
    y[2] = x[2];
  }else{
    scale->y(&y[2],&x[2]);
  }
}
void wz3DContinuation::x(wzFloat *x, const wzFloat *y) const
{
  original->x(x,y);
  if(scale==0){
    x[2] = y[2];
  }else{
    scale->x(&x[2],&y[2]);
  }
}
wz3DContinuation::wz3DContinuation(wzcoordinates orig,wzcoordinates s)
  :original(orig),scale(s)
{}

void wzXYSwitch::y(wzFloat *y, const wzFloat *x) const
{
  wzFloat t[3];
  original->y(t,x); 
  y[0] = t[1];  y[1] = t[0];  y[2] = t[2];
}
void wzXYSwitch::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat t[3];
  t[0] = y[1];  t[1] = y[0];  t[2] = y[2];
  original->x(x,t); 
}
wzXYSwitch::wzXYSwitch(wzcoordinates orig):original(orig){}
void wzYZSwitch::y(wzFloat *y, const wzFloat *x) const
{
  wzFloat t[3];
  t[0] = x[0];  t[1] = x[2];  t[2] = x[1];
  original->y(y,t); 
}
void wzYZSwitch::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat t[3];
  original->x(t,y); 
  x[0] = t[0];  x[1] = t[2];  x[2] = t[1];
}
wzYZSwitch::wzYZSwitch(wzcoordinates orig):original(orig){}
void wzZXSwitch::y(wzFloat *y, const wzFloat *x) const
{
  wzFloat t[3];
  original->y(t,x); 
  y[0] = t[2];  y[1] = t[1];  y[2] = t[0];
}
void wzZXSwitch::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat t[3];
  t[0] = y[2];  t[1] = y[1];  t[2] = y[0];
  original->x(x,t); 
}
wzZXSwitch::wzZXSwitch(wzcoordinates orig):original(orig){}

void wzScaleLog::y(wzFloat *y, const wzFloat *x) const
{
  y[0] = exp(x[0]);
}
void wzScaleLog::x(wzFloat *x, const wzFloat *y) const
{
  if(y[0]<=0)	x[0] = -wzInfty;
  else 		x[0] = log(y[0]);
}
wzScaleExp::wzScaleExp(){}

void wzScaleExp::x(wzFloat *x, const wzFloat *y) const
{
  x[0] = exp(y[0]);
}
void wzScaleExp::y(wzFloat *y, const wzFloat *x) const
{
  if(x[0]<=0) 	y[0] = -wzInfty;
  else		y[0] = log(x[0]);
}
wzScaleLog::wzScaleLog(){}


wzSphericalCoordinates::wzSphericalCoordinates()
{
  Epsilon = 1.e-7;
  ymin[0] = Epsilon;
  ymin[1] = - wzPi/2;
  ymax[1] =   wzPi/2;
  ymin[2] = - wzPi;
  ymax[2] =   wzPi;
}
void wzSphericalCoordinates::setMinimalRadius(wzFloat eps)
{
  Epsilon = eps;
  ymin[0] = Epsilon;
}
void wzSphericalCoordinates::x(wzFloat *x, const wzFloat *y) const
{
  wzFloat yy = (y[0]>Epsilon)?y[0]:Epsilon;
  wzFloat r = yy*cos(y[1]);  if(r<1.e-8) r = 1.e-8;
  x[0] = r*cos(y[2]);
  x[1] = r*sin(y[2]);
  x[2] = yy*sin(y[1]);
}
void wzSphericalCoordinates::y(wzFloat *y, const wzFloat *x) const
{
  try{
    y[2] = atan2(x[1],x[0]);
  }catch(...){
    y[2] = 0;  // throw wzCoordinateFailure;
  }
  wzFloat r = sqrt(x[0]*x[0]+x[1]*x[1]);
  try{
    y[1] = atan2(x[2],r);
  }catch(...){
    y[1] = 0;  // throw wzCoordinateFailure;
  }
  y[0] = sqrt(r*r+x[2]*x[2]);
}
 
/*
wzCylinderCoordinates::wzCylinderCoordinates(wzcoordinates c)
  :coordinates(c)
{}
void wzCylinderCoordinates::x(wzFloat *x, const wzFloat *y) const
{
  coordinates->x(x,y);
  x[2] = y[2];
}
void wzCylinderCoordinates::y(wzFloat *y, const wzFloat *x) const
{
  y[2] = x[2];
  coordinates->y(y,x);
}

*/