wz2droundedbox.cxx

Delaunay三角形的网格剖分程序

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

wz2DroundedBox::wz2DroundedBox(wzFloat xl, wzFloat yl, wzFloat rr)
  :lx(xl),ly(yl)
{
  chart.setBorder(0,wzInfty,- wzPi*phi-2*lx-2*ly, wzPi*phi+2*lx+2*ly);
  phi = rr; // lx+ly;
  wzAssert(phi>0);
  //  ymin[1] = - wzPi*phi-2*lx-2*ly;
  //  ymax[1] =   wzPi*phi+2*lx+2*ly;
  st1 = ly;
  st2 = ly + phi*wzPi/2;
  up  = ly + phi*wzPi/2 + lx;
  st3 = ly + phi*wzPi/2 + 2*lx;
  st4 = ly + phi*wzPi   + 2*lx;
  st5 = wzPi*phi+2*lx+2*ly;
  test();
}

void wz2DroundedBox::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  wzFloat sig,f;
  if(U[1]>0) {sig=1; f=U[1];} else {sig=-1; f=-U[1];}
  if(f<st2){
    if(f<st1){
      X[0] = lx+U[0];
      X[1] = f;
    }else{
      X[0] = lx+U[0]*cos((f-st1)/phi);
      X[1] = ly+U[0]*sin((f-st1)/phi);
    }
  }else if(f<st4){
    if(f<st3){
      X[0] = -lx-f+st3;
      X[1] = ly+U[0];
    }else{
      X[0] = -lx-U[0]*sin((f-st3)/phi);
      X[1] =  ly+U[0]*cos((f-st3)/phi);
    }
  }else{
    X[0] = -lx-U[0];
    X[1] = st5-f;
  }
  X[1] *= sig;
  wzMapPostprocessing(x);
}

void wz2DroundedBox::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  wzFloat xx = X[0], yy=X[1],f,sig;
  if(yy>0){sig = 1; f = yy;}else{sig = -1;f = -yy;}
  if(f<=ly){
    if(xx>=0){
      U[0] = xx-lx;
      U[1] = f;
    }else{
      U[0] = -lx-xx;
      U[1] = st5-f;
    }
  }else if(xx >  lx){
    f -= ly; xx -= lx;
    try{U[1] = atan2(f,xx);}catch(...){U[1] = 0;}
    U[1] *= phi;
    U[1] += st1;
    U[0] = sqrt(f*f+xx*xx);
  }else if(xx < -lx){
    f -= ly; xx = -lx-xx;
    try{U[1] = atan2(xx,f);}catch(...){U[1] = 0;}
    U[1] *= phi;
    U[1] += st3;
    U[0] = sqrt(f*f+xx*xx);
  }else{
    U[0] =  f - ly;
    U[1] = up - xx;
  }
  U[1] *= sig;
  wzMapInversePostprocessing(u);
}