wzmap.cxx

Delaunay三角形的网格剖分程序

#define wzDebug
#include "wzmap.hxx"
#include "wzoutput.hxx"
#include <math.h>

wzMap::wzMap(wzIndex d, int dl)
  :Dim(d),DimLocal(dl<0?d:dl),Properties(0),base(0),Base(0)
  ,full(DimLocal),chart(DimLocal),image(Dim)
{;}

void wzMap::inverse(wzFloat *u, const wzFloat *x) const
{int i; for(i=0;i<DimLocal;i++) u[i]=0;}
void wz1Dmap::inverse(wzFloat *u, const wzFloat *x) const
{u[0]=0;}
void wz2Dmap::inverse(wzFloat *u, const wzFloat *x) const
{u[0]=u[1]=0;}
void wz3Dmap::inverse(wzFloat *u, const wzFloat *x) const
{u[0]=u[1]=u[2]=0;}

void wzMap::compose(wzmap b){base=b; Base=b;}
void wz1Dmap::compose(wz1dmap b){Base=b; base=(wz1Dmap*)b;}
void wz2Dmap::compose(wz2dmap b){Base=b; base=(wz2Dmap*)b;}
void wz3Dmap::compose(wz3dmap b){Base=b; base=(wz3Dmap*)b;}

wzBoolean wzMap::isRectangular() const
{
  if(!(Properties&wzMapPropertyRectangular)) return wzFalse;
  if(Base) return Base->isRectangular();
  return wzTrue;
}
wzBoolean wzMap::isConformal()   const
{
  if(!(Properties&wzMapPropertyConformal)) return wzFalse;
  if(Base) return Base->isConformal();
  return wzTrue;
}
wzBoolean wzMap::isInvertible()      const
{
  if(!(Properties&wzMapPropertyInvertible)) return wzFalse;
  if(Base) return Base->isInvertible();
  return wzTrue;
}
wzBoolean wzMap::isInvertibleLeft()  const
{
  if(!(Properties&wzMapPropertyInvertibleLeft)) return wzFalse;
  if(Base) return Base->isInvertibleLeft();
  return wzTrue;
}
wzBoolean wzMap::isInvertibleRight() const
{
  if(!(Properties&wzMapPropertyInvertibleRight)) return wzFalse;
  if(Base) return Base->isInvertibleRight();
  return wzTrue;
}

wzBoolean wzMap::isOrthogonal()  const
{
  if(!(Properties&wzMapPropertyOrthogonal)) return wzFalse;
  if(!Base) return wzTrue;
  if(Properties&wzMapPropertyConformal)
    return Base->isOrthogonal();
  return Base->isRectangular();
}

wzMap::~wzMap(){base=0;}

wz3Dcoordinates::wz3Dcoordinates(wz3dmap b)
  :wz3Dmap(b,3)
{Properties = wzIsCoordinateCarthesian; name="Carthesian coordinates";}

wz2Dcoordinates::wz2Dcoordinates(wz2dmap b)
  :wz2Dmap(b,2)
{Properties = wzIsCoordinateCarthesian; name="Carthesian coordinates";}

wz1Dcoordinates::wz1Dcoordinates(wz1dmap b)
  :wz1Dmap(b,1)
{Properties = wzIsCoordinateCarthesian; name="Carthesian coordinates";}

void wz3Dmap::callDirect(wzFloat *xx) const
{
  wzFloat X[3]={xx[0],xx[1],xx[2]}; (*(wzMap*)this)(xx,X);
}
void wz3Dmap::callInverse(wzFloat *xx) const
{
  const wzFloat X[3]={xx[0],xx[1],xx[2]}; wzMap::inverse(xx,X);
}

void wz2Dmap::callDirect(wzFloat *xx) const
{
  wzFloat X[2]={xx[0],xx[1]}; (*(wzMap*)this)(xx,X);
}
void wz2Dmap::callInverse(wzFloat *xx) const
{
  const wzFloat X[2]={xx[0],xx[1]}; wzMap::inverse(xx,X);
}

void wz1Dmap::callDirect(wzFloat *xx) const
{
  wzFloat X[1]={xx[0]}; (*(wzMap*)this)(xx,X);
}
void wz1Dmap::callInverse(wzFloat *xx) const
{
  const wzFloat X[1]={xx[0]}; wzMap::inverse(xx,X);
}

wz1dmap   wzMap::map1D() const{return Dim==1?(wz1Dmap*)this:0;}
wz2dmap   wzMap::map2D() const{return Dim==2?(wz2Dmap*)this:0;}
wz3dmap   wzMap::map3D() const{return Dim==3?(wz3Dmap*)this:0;}


wz3Dmap::wz3Dmap(wz3dmap b, int dl)
  :wzMap(3,dl)
  ,SingleMode(1),Mode(wzAffineNone)
{
  int i,j;
  compose(b);
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Full[i][j] = 0;
      Reverse[i][j] = 0;
    }
    Full[i][i] = 1;
    Reverse[i][i] = 1;
  }
  for(i=0;i<3;i++){
    Shift[i] = 0;
    Factor[i] = 1;
  }
}

wz2Dmap::wz2Dmap(wz2dmap b, int dl)
  :wzMap(2,dl)
  ,SingleMode(1),Mode(wzAffineNone)
{
  int i,j;
  compose(b);
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      Full[i][j] = 0;
      Reverse[i][j] = 0;
    }
    Full[i][i] = 1;
    Reverse[i][i] = 1;
  }
  for(i=0;i<2;i++){
    Shift[i] = 0;
    Factor[i] = 1;
  }
}

wz1Dmap::wz1Dmap(wz1dmap b, int dl)
  :wzMap(1,dl)
  ,Mode(wzAffineNone)
{
  compose(b);
  Shift = 0;
  Factor = 1;
}

void wz3Dmap::shift(wzFloat sx, wzFloat sy, wzFloat sz)
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4],S[3];
  switch(Mode){
  case wzAffineNone:
    Mode = wzAffineShift;
  case wzAffineGeneral:
  case wzAffineStretchShift:
  case wzAffineShift:
    ;
  }
  Shift[0] += S[0] = sx;
  Shift[1] += S[1] = sy;
  Shift[2] += S[2] = sz;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = TmpR[i][i] = 1;
  }
  for(i=0;i<3;i++){
    TmpF[i][3] = S[i];
    TmpR[i][3] = -S[i];
  }
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<4;k++){
	Full[i][j]	+= TmpF[i][k]*OldF[k][j];
	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz2Dmap::shift(wzFloat sx, wzFloat sy)
{
  int i,j,k;
  wzFloat TmpR[3][3],TmpF[3][3],OldF[3][3],OldR[3][3],S[2];
  switch(Mode){
  case wzAffineNone:
    Mode = wzAffineShift;
  case wzAffineGeneral:
  case wzAffineStretchShift:
  case wzAffineShift:
    ;
  }
  Shift[0] += S[0] = sx;
  Shift[1] += S[1] = sy;
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = TmpR[i][i] = 1;
  }
  for(i=0;i<2;i++){
    TmpF[i][2] = S[i];
    TmpR[i][2] = -S[i];
  }
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<3;k++){
	Full[i][j]   	+= TmpF[i][k]*OldF[k][j];
	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz1Dmap::shift(wzFloat sx)
{
  switch(Mode){
  case wzAffineNone:
    Mode = wzAffineShift;
  case wzAffineGeneral:
  case wzAffineStretchShift:
  case wzAffineShift:
    ;
  }
  Shift += sx;
}

void wz3Dmap::stretch(wzFloat fx)
{stretchAnisotropic(fx,fx,fx);}
void wz3Dmap::stretchAnisotropic(wzFloat fx, wzFloat fy, wzFloat fz)
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4],F[4];
  if((fx != fy) || (fx != fz)){
    Properties &= ~wzMapPropertyConformal;
    if(!(Properties & wzMapPropertyRectangular))
      Properties &= ~wzMapPropertyOrthogonal;
  }
  switch(Mode){
  case wzAffineNone:
  case wzAffineShift:
    Mode = wzAffineStretchShift;
  case wzAffineStretchShift:
  case wzAffineGeneral:
    ;
  }
  Factor[0] *= (F[0] = fx);
  Factor[1] *= (F[1] = fy);
  Factor[2] *= (F[2] = fz);
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = 1;
    TmpR[i][i] = 1;
  }
  for(i=0;i<3;i++){
    TmpF[i][i] = F[i];
    TmpR[i][i] = 1/F[i];
    Shift[i]  *= F[i];
  }
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<4;k++){
 	Full[i][j]	+= TmpF[i][k]*OldF[k][j];
 	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz2Dmap::stretch(wzFloat fx)
{stretchAnisotropic(fx,fx);}

void wz2Dmap::stretchAnisotropic(wzFloat fx, wzFloat fy)
{
  int i,j,k;
  wzFloat TmpR[3][3],TmpF[3][3],OldF[3][3],OldR[3][3],F[3];
  if(fx != fy){
    Properties &= ~wzMapPropertyConformal;
    if(!(Properties & wzMapPropertyRectangular))
      Properties &= ~wzMapPropertyOrthogonal;
  }
  switch(Mode){
  case wzAffineNone:
  case wzAffineShift:
    Mode = wzAffineStretchShift;
  case wzAffineStretchShift:
  case wzAffineGeneral:
    ;
  }
  Factor[0] *= (F[0] = fx);
  Factor[1] *= (F[1] = fy);
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = 1;
    TmpR[i][i] = 1;
  }
  for(i=0;i<2;i++){
    TmpF[i][i] = F[i];
    TmpR[i][i] = 1/F[i];
    Shift[i]  *= F[i];
  }
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<3;k++){
 	Full[i][j]   	+= TmpF[i][k]*OldF[k][j];
 	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz1Dmap::stretch(wzFloat fx)
{
  switch(Mode){
  case wzAffineNone:
  case wzAffineShift:
    Mode = wzAffineStretchShift;
  case wzAffineStretchShift:
  case wzAffineGeneral:
    ;
  }
  Factor *= fx;
  Shift  *= fx;
}

void wz3Dmap::rotateX(wzFloat phi)
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  Properties &= ~wzMapPropertyRectangular;
  Mode = wzAffineGeneral;
  SingleMode = 0;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = 1;
    TmpR[i][i] = 1;
  }
  wzFloat c = cos(phi), s = sin(phi);
  TmpF[1][1] = TmpF[2][2] = c;
  TmpR[1][1] = TmpR[2][2] = c;
  TmpF[1][2] = TmpR[2][1] = -s;
  TmpR[1][2] = TmpF[2][1] = s;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<4;k++){
	Full[i][j]	+= TmpF[i][k]*OldF[k][j];
	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz3Dmap::rotateY(wzFloat phi)
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  Properties &= ~wzMapPropertyRectangular;
  Mode = wzAffineGeneral;
  SingleMode = 0;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = 1;
    TmpR[i][i] = 1;
  }
  wzFloat c = cos(phi), s = sin(phi);
  TmpF[0][0] = TmpF[2][2] = c;
  TmpR[0][0] = TmpR[2][2] = c;
  TmpF[0][2] = TmpR[2][0] = s;
  TmpR[0][2] = TmpF[2][0] = -s;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Full[i][j] = Reverse[i][j]=0;
      for(k=0;k<4;k++){
	Full[i][j]	+= TmpF[i][k]*OldF[k][j];
	Reverse[i][j]	+= OldR[i][k]*TmpR[k][j];
      }
    }
  }
}

void wz3Dmap::rotateZ(wzFloat phi)
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  Properties &= ~wzMapPropertyRectangular;
  Mode = wzAffineGeneral;
  SingleMode = 0;
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      TmpF[i][j] = TmpR[i][j] = 0;
      OldF[i][j] = Full[i][j];
      OldR[i][j] = Reverse[i][j];
    }
    TmpF[i][i] = 1;
    TmpR[i][i] = 1;
  }
  wzFloat c = cos(phi), s = sin(phi);
  TmpF[0][0] = TmpF[1][1] = c;