wz2dmaps.cxx

Delaunay三角形的网格剖分程序

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

void wz2Dcoordinates::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  X[0]=U[0];  X[1]=U[1];
  wzMapPostprocessing(x);
}

void wz2Dcoordinates::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  U[0]=X[0];  U[1]=X[1];
  wzMapInversePostprocessing(u);
}

void wz2Dpolar::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  wzFloat y = (U[0]>Epsilon)?U[0]:Epsilon;
  X[0] = y*cos(U[1]);
  X[1] = y*sin(U[1]);
  wzMapPostprocessing(x);
}
void wz2Dpolar::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  U[0] = sqrt(X[0]*X[0]+X[1]*X[1]);
  U[1] = atan2(X[1],X[0]);
  wzMapInversePostprocessing(u);
}

void wz2Dexp::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  X[0] = exp(U[0])*cos(U[1]);
  X[1] = exp(U[0])*sin(U[1]);
  wzMapPostprocessing(x);
}
void wz2Dexp::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  U[0] = log(X[0]*X[0]+X[1]*X[1])/2;
  U[1] = atan2(X[1],X[0]);
  wzMapInversePostprocessing(u);
}
void wz2Dlog::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  X[0] = log(U[0]*U[0]+U[1]*U[1])/2;
  X[1] = atan2(U[1],U[0]);
  wzMapPostprocessing(x);
}
void wz2Dlog::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  U[0] = exp(X[0])*cos(X[1]);
  U[1] = exp(X[0])*sin(X[1]);
  wzMapInversePostprocessing(u);
}

void wz2Dcosh::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  wzFloat r= exp(U[0]),ri=1/r,c=cos(U[1]),s=sin(U[1]);
  X[0] = (r+ri)*c/2;
  X[1] = (r-ri)*s/2;
  wzMapPostprocessing(x);
}

// maps on r>0 (outside of unit sphere in exp-coordinates)
void wz2Dcosh::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  double u0,u1,a0,a1,sq,s0,s1;
  double v0,v1;
  a1 = (u0=X[0])*(u1=X[1]);
  a0 = (u0*u0-u1*u1-1)/2;
  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){  where (0,lambda) is the centre of the circle
  if(u0*s0<(-u1)*s1){
    v0 = u0 - s0;
    v1 = u1 - s1;
  }else{
    v0 = u0 + s0;
    v1 = u1 + s1;
  }
  //  if(a1<0)  v0 = u0 - s0;
  //  else      v0 = u0 + s0;
  //  v1 = u1 + s1;
  U[0] = log(v0*v0+v1*v1)/2;
  U[1] = atan2(v1,v0);
  wzMapInversePostprocessing(u);
}

void wz2Dshukovski::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  wzFloat u0=U[0],u1=U[1],r2 = u0*u0+u1*u1;
  //  wzFloat r= exp(U[0]),ri=1/r,c=cos(U[1]),s=sin(U[1]);
  //  X[0] = (r+ri)*c/2;
  //  X[1] = (r-ri)*s/2;
  X[0] = u0*(0.5 + 0.5/r2);
  X[1] = u1*(0.5 - 0.5/r2);
  wzMapPostprocessing(x);
}

// maps on r>0 (outside of unit sphere in exp-coordinates)
void wz2Dshukovski::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  double u0,u1,a0,a1,sq,s0,s1;
  a1 = (u0=X[0])*(u1=X[1]);
  a0 = (u0*u0-u1*u1-1)/2;
  sq = sqrt(a0*a0+a1*a1);
  s0 = sqrt(sq+a0);
  s1 = sqrt(sq-a0);
  if(a1<0) s0 = -s0;
  if(u0*s0<(yc-u1)*s1){  //where (0,yc) is the centre of the critical circle
    U[0] = u0 - s0;
    U[1] = u1 - s1;
  }else{
    U[0] = u0 + s0;
    U[1] = u1 + s1;
  }
  wzMapInversePostprocessing(u);
}

void wz2Dacosh::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  wzFloat r= exp(X[0]),ri=1/r,c=cos(X[1]),s=sin(X[1]);
  U[0] = (r+ri)*c/2;
  U[1] = (r-ri)*s/2;
  wzMapInversePostprocessing(u);
}

// maps on r>0 (outside of unit sphere in exp-coordinates)
void wz2Dacosh::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  double u0,u1,a0,a1,sq,s0,s1;
  double v0,v1;
  a1 = (u0=U[0])*(u1=U[1]);
  a0 = (u0*u0-u1*u1-1)/2;
  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){  where (0,lambda) is the centre of the circle
  if(u0*s0<(-u1)*s1){
    v0 = u0 - s0;
    v1 = u1 - s1;
  }else{
    v0 = u0 + s0;
    v1 = u1 + s1;
  }
  X[0] = log(v0*v0+v1*v1)/2;
  X[1] = atan2(v1,v0);
  wzMapPostprocessing(x);
}
/*
void wz2Dasin::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  wzFloat e = exp(X[1]), i=1/e, c = cos(X[0]), s = sin(X[0]);
  U[0] = (e+i)*s/2;
  U[1] = (e-i)*c/2;
  wzMapInversePostprocessing(u);
}

void wz2Dasin::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  double u0,u1,a0,a1,sq,s0,s1;
  double r0,r1,v0,v1;
  u0 = U[0]/2;
  u1 = U[1]/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)  v0 = u0 - s0;
  else      v0 = u0 + s0;
  v1 = u1 + s1;
  X[0] = atan2(v1,v0);
  X[1] = log(v0*v0+v1*v1)/2;
  wzMapPostprocessing(x);
}

*/

wz2Dpolar::wz2Dpolar(wz2dmap m)
  :wz2Dcoordinates(m)
{
  Epsilon = 4.e-9;
  name = "polar coordinates";
  Properties = wzIsCoordinateOrthogonal;
  full.setBorder(0);
  chart.setBorder(0,wzInfty,-wzPi,wzPi);
}
void wz2Dpolar::setMinimalRadius(wzFloat eps)
{
  Epsilon = eps;
}

wz2Dexp::wz2Dexp(wz2dmap m)
  :wz2Dcoordinates(m)
{
  name = "complex exponential";
  Properties = wzIsCoordinateConformal;
  chart.setBorder(-wzInfty,wzInfty,-wzPi,wzPi);
}

void wz2Drectangular::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  if(xx)
    xx->inverse(&u[0],&X[0]);
  else
    u[0] = X[0];
  if(yy)
    yy->inverse(&u[1],&X[1]);
  else
    u[1] = X[1];
}

void wz2Drectangular::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  if(xx)
    (*xx)(&X[0],&u[0]);
  else
    X[0] = u[0];
  if(yy)
    (*yy)(&X[1],&u[1]);
  else
    X[1] = u[1];
  wzMapPostprocessing(x);
}

wz2Drectangular::wz2Drectangular(wz1dmap x0, wz1dmap y0)
  :wz2Dcoordinates(),x(x0),xx(x0),y(y0),yy(y0)
{
  base=0;Base=0;
  name = "rectangular coordinates";
  Properties = wzIsCoordinateRectangular;
}

void wz2Drectangular::compose(wz2dmap x){;}
void wz2Drectangular::composeX(wz1dmap b){x=b; xx=b;}
void wz2Drectangular::composeY(wz1dmap b){y=b; yy=b;}

wz2Dcosh::wz2Dcosh(wz2dmap x)
  :wz2Dcoordinates(x)
{
  name = "complex cosinus hyperbolicus";
  Properties = wzIsCoordinateConformal;
  //  chart.setBorder(-3,3,0,wzPi);
  chart.setBorder(0,3,-wzPi,wzPi);
}

wz2Dacosh::wz2Dacosh(wz2dmap x)
  :wz2Dcoordinates(x)
{
  name = "complex arcus cosinus hyperbolicus";
  Properties = wzIsCoordinateConformal;
  image.setBorder(0,wzInfty,-wzPi,wzPi);
  // a nice choice for a picture:
  //  chart.setBorder(0,8,0,4);
}

wz2Dlog::wz2Dlog(wz2dmap x)
  :wz2Dcoordinates(x)
{
  name = "complex logarithm";
  Properties = wzIsCoordinateConformal;
  image.setBorder(0,wzInfty,-wzPi,wzPi);
  // a nice choice for a picture:
  //  chart.setBorder(0,8,0,4);
}

wz2Dshukovski::wz2Dshukovski(wz2dmap x, wzFloat y)
  :wz2Dcoordinates(x)
{
  name = "Shukovski function";
  Properties = wzIsCoordinateConformal;
  image.setBorder(-wzInfty,wzInfty,-wzInfty,wzInfty);
  full.setBorder(-wzInfty,wzInfty,-wzInfty,wzInfty);
  setYofCircleCenter(y);
}

void wz2Dshukovski::setYofCircleCenter(wzFloat y)
{
  yc=y;
  if(yc>=wzInfty/2){
    chart.setBorder(-wzInfty,wzInfty,-wzInfty,0);
  }else if(yc>=0){
    chart.setBorder(-wzInfty,wzInfty,-wzInfty,yc-sqrt(1+yc*yc));
  }else if(yc>-wzInfty/2){
    chart.setBorder(-wzInfty,wzInfty,yc+sqrt(1+yc*yc),wzInfty);
  }else{
    chart.setBorder(-wzInfty,wzInfty,0,wzInfty);
  }
}

wz2DshiftShukovski::wz2DshiftShukovski(wz2dmap b,wzFloat y)
  :wz2Dcoordinates(b)
{
  setShiftPoints(0,-1,0,0);
  name = "Shukovski function";
  Properties = wzIsCoordinateConformal;
  setYofCircleCenter(y);
}

void wz2DshiftShukovski::setYofCircleCenter(wzFloat l)
{
  lambda = l;
}

void wz2DshiftShukovski::setShiftPoints(wzFloat x0,wzFloat y0,
					wzFloat x1,wzFloat y1)
{
  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;
}

void wz2DshiftShukovski::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  double rs,r0,r1,r,q0,q1;
  r0 = U[0] - m0;
  r1 = U[1] - m1;
  rs = r0*d0+r1*d1;
  q0 = r0+rs*e0;
  q1 = r1+rs*e1;
  r = (r0*r0+r1*r1)*rr*rr;
  X[0] = U[0]+q0/r;
  X[1] = U[1]+q1/r;
  wzMapPostprocessing(x);
}

void wz2DshiftShukovski::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  double u0,u1,a0,a1,sq,s0,s1;
  double r0,r1,v0,v1;
  r0 = X[0] - m0;
  r1 = X[1] - m1;
  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;
  }
  U[0] = m0 + (v0*c0 - v1*c1)/rr;
  U[1] = m1 + (v0*c1 + v1*c0)/rr;
  wzMapInversePostprocessing(u);
}