wz3dmaps.cxx

Delaunay三角形的网格剖分程序

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

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

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

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

void wz3Dcylinder::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  wzFloat yy = (U[0]>Epsilon)?U[0]:Epsilon;
  X[0] = yy*cos(U[1]);
  X[1] = yy*sin(U[1]);
  X[2] = U[2];
  wzMapPostprocessing(x);
}

void wz3Dcylinder::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]);
  U[2] = X[2];
  wzMapInversePostprocessing(u);
}

void wz3Dcylindrical::inverse(wzFloat *u, const wzFloat *x) const
{
  wzMapInversePreprocessing(u,x);
  if(XY)
    XY->inverse(&u[0],&X[0]);
  else{
    u[0] = X[0];
    u[1] = X[1];
  }
  if(Z)
    Z->inverse(&u[2],&X[2]);
  else
    u[2] = X[2];
}

void wz3Dcylindrical::operator()(wzFloat *x, const wzFloat *u) const
{
  wzMapPreprocessing(x,u);
  if(XY)
    (*XY)(&X[0],&u[0]);
  else{
    X[0] = u[0];
    X[1] = u[1];
  }
  if(Z)
    (*Z)(&X[2],&u[2]);
  else
    X[2] = u[2];
  wzMapPostprocessing(x);
}


wz3Dspherical::wz3Dspherical()
{
  Epsilon  = 1.e-9;
  Epsilon1 = 1.e-13;
  name = "spherical coordinates";
  Properties = wzIsCoordinateOrthogonal;
  full.setBorder(0);
  chart.setBorder(0,wzInfty,-wzPi/2,wzPi/2,-wzPi,wzPi);
}

void wz3Dspherical::setMinimalRadius(wzFloat eps, wzFloat eps1)
{
  Epsilon  = eps;
  Epsilon1 = eps*eps1;
}

wz3Dcylinder::wz3Dcylinder()
{
  Epsilon = 4.e-9;
  name = "cylinder coordinates";
  Properties = wzIsCoordinateOrthogonal;
  chart.setBorder(0,wzInfty,-wzPi,wzPi);
}

void wz3Dcylinder::setMinimalRadius(wzFloat eps)
{
  Epsilon = eps;
}

wz3Dcylindrical::wz3Dcylindrical(wz2dmap x2, wz1dmap zz)
  :wz3Dcoordinates(),xy(x2),XY(x2),z(zz),Z(zz)
{
  base=0;Base=0;
  name = "cylindrical coordinates";
  if(!x2 || x2->isOrthogonal())
    Properties = wzIsCoordinateOrthogonal;
  else
    Properties = wzIsCoordinate;
}

void wz3Dcylindrical::composeXY(wz2dmap b){xy=b; XY=b;}
void wz3Dcylindrical::composeZ(wz1dmap b){z=b; Z=b;}
void wz3Dcylindrical::compose(wz3dmap b){;}

void wz3Drectangular::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];
  if(zz)
    zz->inverse(&u[2],&X[2]);
  else
    u[2] = X[2];
}

void wz3Drectangular::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];
  if(zz)
    (*zz)(&X[2],&u[2]);
  else
    X[2] = u[2];
  wzMapPostprocessing(x);
}

wz3Drectangular::wz3Drectangular(wz1dmap x0, wz1dmap y0, wz1dmap z0)
  :wz3Dcoordinates(),x(x0),xx(x0),y(y0),yy(y0),z(z0),zz(z0)
{
  base=0;Base=0;
  name = "rectangular coordinates";
  Properties = wzIsCoordinateRectangular;
}

void wz3Drectangular::compose(wz3dmap x){;}
void wz3Drectangular::composeX(wz1dmap b){x=b; xx=b;}
void wz3Drectangular::composeY(wz1dmap b){y=b; yy=b;}
void wz3Drectangular::composeZ(wz1dmap b){z=b; zz=b;}