wzmap.cxx

Delaunay三角形的网格剖分程序

  TmpR[0][0] = TmpR[1][1] = c;
  TmpF[0][1] = TmpR[1][0] = -s;
  TmpR[0][1] = TmpF[1][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 wz2Dmap::rotate(wzFloat phi)
{
  int i,j,k;
  wzFloat TmpR[3][3],TmpF[3][3],OldF[3][3],OldR[3][3];
  Properties &= ~wzMapPropertyRectangular;
  Mode = wzAffineGeneral;
  SingleMode = 0;
  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;
  }
  wzFloat c = cos(phi), s = sin(phi);
  TmpF[0][0] = TmpF[1][1] = c;
  TmpR[0][0] = TmpR[1][1] = c;
  TmpF[0][1] = TmpR[1][0] = -s;
  TmpR[0][1] = TmpF[1][0] = s;
  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 wz3Dmap::switchYZ()
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  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;
  }
  TmpF[1][1] = TmpF[2][2] = 0;
  TmpR[1][1] = TmpR[2][2] = 0;
  TmpF[1][2] = TmpR[2][1] = 1;
  TmpR[1][2] = TmpF[2][1] = 1;
  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::switchXZ()
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  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;
  }
  TmpF[0][0] = TmpF[2][2] = 0;
  TmpR[0][0] = TmpR[2][2] = 0;
  TmpF[0][2] = TmpR[2][0] = 1;
  TmpR[0][2] = TmpF[2][0] = 1;
  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::switchXY()
{
  int i,j,k;
  wzFloat TmpR[4][4],TmpF[4][4],OldF[4][4],OldR[4][4];
  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;
  }
  TmpF[0][0] = TmpF[1][1] = 0;
  TmpR[0][0] = TmpR[1][1] = 0;
  TmpF[0][1] = TmpR[1][0] = 1;
  TmpR[0][1] = TmpF[1][0] = 1;
  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::switchXY()
{
  int i,j,k;
  wzFloat TmpR[3][3],TmpF[3][3],OldF[3][3],OldR[3][3];
  Mode = wzAffineGeneral;
  SingleMode = 0;
  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;
  }
  TmpF[0][0] = TmpF[1][1] = 0;
  TmpR[0][0] = TmpR[1][1] = 0;
  TmpF[0][1] = TmpR[1][0] = 1;
  TmpR[0][1] = TmpF[1][0] = 1;
  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 wzMap::print() const
{
  wzOutput& p=wzOutputDefault;
  p("<>D <>\n invertible in "),dimension(),name;
  chart.print();
  p(" (<>rectangular, <>conformal, <>orthogonal)\n")
    ,isRectangular()?"":"not "
    ,isConformal()?"":"not "
    ,isOrthogonal()?"":"not ";
}

void wz1Dmap::print() const
{
  wzMap::print();
  wzOutput& p=wzOutputDefault;
  switch(Mode){
  case wzAffineStretchShift:
    p("<>\n")," Affine transformation:";
    p(" x = <> u + <>\n" ),Factor,Shift;
    break;
  case wzAffineShift:
    p("<>\n")," Affine transformation:";
    p(" x = u + <>\n"),Shift;
    break;
  case wzAffineNone:
    p("<>\n")," No affine transformation";
  }
}

void wz2Dmap::print() const
{
  wzMap::print();
  wzOutput& p=wzOutputDefault;
  switch(Mode){
  case wzAffineGeneral:
    p("<>\n")," Affine transformation:";
    p(" x = <8> u + <8> v + <8>\n"),Full[0][0],Full[0][1],Full[0][2];
    p(" y = <8> u + <8> v + <8>\n"),Full[1][0],Full[1][1],Full[1][2];
    p(" u = <8> x + <8> y + <8>\n")
      ,Reverse[0][0],Reverse[0][1],Reverse[0][2];
    p(" v = <8> x + <8> y + <8>\n")
      ,Reverse[1][0],Reverse[1][1],Reverse[1][2];
    break;
  case wzAffineStretchShift:
    p("<>\n")," Affine transformation:";
    p(" x = <8> u + <8>," ),Factor[0],Shift[0];
    p(" y = <8> v + <8>\n"),Factor[1],Shift[1];
    break;
  case wzAffineShift:
    p("<>\n")," Affine transformation:";
    p(" x = u + <8>, y = v + <8>\n"),Shift[0],Shift[1];
    break;
  case wzAffineNone:
    p("<>\n")," No affine transformation";
  }
}
void wz3Dmap::print() const
{
  wzMap::print();
  wzOutput& p=wzOutputDefault;
  switch(Mode){
  case wzAffineGeneral:
    p("<>\n")," Affine transformation:";
    p(" x = <8> u + <8> v + <8> w + <8>\n")
      ,Full[0][0],Full[0][1],Full[0][2],Full[0][3];
    p(" y = <8> u + <8> v + <8> w + <8>\n")
      ,Full[1][0],Full[1][1],Full[1][2],Full[1][3];
    p(" z = <8> u + <8> v + <8> w + <8>\n")
      ,Full[2][0],Full[2][1],Full[2][2],Full[2][3];
    p(" u = <8> x + <8> y + <8> z + <8>\n")
      ,Reverse[0][0],Reverse[0][1],Reverse[0][2],Reverse[0][3];
    p(" v = <8> x + <8> y + <8> z + <8>\n")
      ,Reverse[1][0],Reverse[1][1],Reverse[1][2],Reverse[1][3];
    p(" w = <8> x + <8> y + <8> z + <8>\n")
      ,Reverse[2][0],Reverse[2][1],Reverse[2][2],Reverse[2][3];
    break;
  case wzAffineStretchShift:
    p("<>\n")," Affine transformation:";
    p(" x = <8> u + <8>, y = <8> v + <8>, z = <8> w + <8>\n")
      ,Factor[0],Shift[0],Factor[1],Shift[1],Factor[2],Shift[2];
    break;
  case wzAffineShift:
    p("<>\n")," Affine transformation:";
    p(" x = u + <8>, y = v + <8>, z = w + <8>\n")
      ,Shift[0],Shift[1],Shift[2];
    break;
  case wzAffineNone:
    p("<>\n")," No affine transformation";
  }
}

void wz1Dcoordinates::print() const
{
  wz1Dmap::print();
  wzOutput& p=wzOutputDefault;
  p("<>\n")," Examples:";
  wzFloat x[1],r[1];
  wzFloat lm[]={-1};  call(x,lm); inverse(r,x);
  p("[<>] => [<>]\n"),r[0],x[0];
  wzFloat l0[]={0};  call(x,l0); inverse(r,x);
  p("[<>] => [<>]\n"),r[0],x[0];
  wzFloat l1[]={0.5};  call(x,l1); inverse(r,x);
  p("[<>] => [<>]\n"),r[0],x[0];
  wzFloat l2[]={1};  call(x,l2); inverse(r,x);
  p("[<>] => [<>]\n"),r[0],x[0];
  wzFloat l3[]={2.0};  call(x,l3); inverse(r,x);
  p("[<>] => [<>]\n"),r[0],x[0];
}

void wz2Dcoordinates::print() const
{
  wz2Dmap::print();
  wzOutput& p=wzOutputDefault;
  p("<>\n")," Examples:";
  wzFloat x[2],r[2];
  wzFloat l0[]={0,0};  call(x,l0); inverse(r,x);
  p("[<>,<>] => [<>,<>]\n"),r[0],r[1],x[0],x[1];
  wzFloat l1[]={1,0};  call(x,l1); inverse(r,x);
  p("[<>,<>] => [<>,<>]\n"),r[0],r[1],x[0],x[1];
  wzFloat l2[]={0,1};  call(x,l2); inverse(r,x);
  p("[<>,<>] => [<>,<>]\n"),r[0],r[1],x[0],x[1];
}

void wz3Dcoordinates::print() const
{
  wz3Dmap::print();
  wzOutput& p=wzOutputDefault;
  p("<>\n")," Examples:";
  wzFloat x[3],r[3];
  wzFloat l0[]={0,0,0};  call(x,l0); inverse(r,x);
  p("[<>,<>,<>] => [<>,<>,<>]\n"),r[0],r[1],r[2],x[0],x[1],x[2];
  wzFloat l1[]={1,0,0};  call(x,l1); inverse(r,x);
  p("[<>,<>,<>] => [<>,<>,<>]\n"),r[0],r[1],r[2],x[0],x[1],x[2];
  wzFloat l2[]={0,1,0};  call(x,l2); inverse(r,x);
  p("[<>,<>,<>] => [<>,<>,<>]\n"),r[0],r[1],r[2],x[0],x[1],x[2];
  wzFloat l3[]={0,0,1};  call(x,l3); inverse(r,x);
  p("[<>,<>,<>] => [<>,<>,<>]\n"),r[0],r[1],r[2],x[0],x[1],x[2];
}

void wzMap::dx(const wzFloat *x, wzFloat *dx0,  wzFloat *dx1,  wzFloat *dx2) const
{
  switch(Dim){
  case 0:
    return;
  case 1:
    ((wz1Dmap *)this)->dx(x,dx0); return;
  case 2:
    ((wz2Dmap *)this)->dx(x,dx0,dx1); return;
  case 3:
    ((wz3Dmap *)this)->dx(x,dx0,dx1,dx2); return;
  }
}
void wz3Dmap::dx(const wzFloat *u, wzFloat *dx0,  wzFloat *dx1,  wzFloat *dx2) const
{
  int i;
  wzFloat x[3],xh[3],uh[wzPointDim];
  wzFloat h=1.e-7;
  (*this)(x,u);
  for(i=0;i<DimLocal;i++) uh[i]=u[i];
  for(i=0;i<DimLocal;i++){
    uh[i] += h;   (*this)(xh,uh); uh[i]=u[i];
    dx0[i]=xh[0]-x[0];
    dx1[i]=xh[1]-x[1];
    dx2[i]=xh[2]-x[2];
  }
}
void wz2Dmap::dx(const wzFloat *u, wzFloat *dx0,  wzFloat *dx1) const
{
  int i;
  wzFloat x[2],xh[2],uh[wzPointDim];
  wzFloat h=1.e-7;
  (*this)(x,u);
  for(i=0;i<DimLocal;i++) uh[i]=u[i];
  for(i=0;i<DimLocal;i++){
    uh[i] += h;   (*this)(xh,uh); uh[i]=u[i];
    dx0[i]=xh[0]-x[0];
    dx1[i]=xh[1]-x[1];
  }
}
void wz1Dmap::dx(const wzFloat *u, wzFloat *dx0) const
{
  int i;
  wzFloat x,xh,uh[wzPointDim];
  wzFloat h=1.e-7;
  (*this)(&x,u);
  for(i=0;i<DimLocal;i++) uh[i]=u[i];
  for(i=0;i<DimLocal;i++){
    uh[i] += h;   (*this)(&xh,uh); uh[i]=u[i];
    dx0[i]=xh-x;
  }
}

/*

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

*/

void wzMap::test() const
{
  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(chart.min[i] == -wzInfty) beg[i] = -2; else beg[i] = chart.min[i]+0.1;
    if(chart.max[i] ==  wzInfty) end[i] = +2; else end[i] = chart.max[i]-0.1;
    d[i] = (end[i]-beg[i])/11;
    d0max += 0.01*d[i]*d[i];
  }
  wzOutput& p=wzOutputDefault;
  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++){
	    call(yy,xx); inverse(zz,yy);d0=0;
	    for(l=0;l<Dim;l++){
	      dd = zz[l]-xx[l]; d0 += dd*dd;
	    }
       	    if(d0>d0max)
	      p("[<>,<>,<>]=>[<>,<>,<>]=>[<>,<>,<>]\n")
		,xx[0],xx[1],xx[2]
       	       	,yy[0],yy[1],yy[2]
		,zz[0],zz[1],zz[2];
	  }
	  xx[2] += d[2];
	}else{
	  call(yy,xx); inverse(zz,yy);d0=0;
	  for(l=0;l<Dim;l++){
	    dd = zz[l]-xx[l]; d0 += dd*dd;
	  }
	  if(d0>d0max)
      	    p("[<>,<>,0]=>[<>,<>,0]=>[<>,<>,0]\n")
      	      ,xx[0],xx[1]
      	      ,yy[0],yy[1]
      	      ,zz[0],zz[1];
      	}
      	xx[1] += d[1];
      }
    }else{
      call(yy,xx); inverse(zz,yy);d0=0;
      for(l=0;l<Dim;l++){
      	dd = zz[l]-xx[l]; d0 += dd*dd;
      }
      if(d0>d0max)
	p("[<>,0,0]=>[<>,0,0]=>[<>,0,0]\n")
	      ,xx[0]
	      ,yy[0]
	      ,zz[0];
    }
    xx[0] += d[0];
  }
}