rbox.c

关于网格剖分的

	}
      }
      for (j= 0; j<numpoints; j++) {
        if (iswidth) 
          apex= qh_RANDOMint % (dim+1);
        else
          apex= -1;
        for (k= 0; k<dim; k++)
	  coord[k]= 0.0;
	norm= 0.0;
        for (i= 0; i<dim+1; i++) {
    	  randr= qh_RANDOMint;
          factor= randr/randmax;
          if (i == apex)
            factor *= width;
          norm += factor;
          for (k= 0; k<dim; k++) {
	    simplexp= simplex + i*dim + k;
            coord[k] += factor * (*simplexp);
          }
        }
        for (k= 0; k<dim; k++)
          coord[k] /= norm;
        if (iscdd)
          out1( 1.0);
        for (k=0; k < dim; k++) 
	  out1( coord[k] * box);
        fprintf (fp, "\n");
      }
      isregular= 0; /* continue with isbox */
      numpoints= 0;
    }
    /* ============= mesh distribution =============== */
    if (ismesh) {
      nthroot= pow (numpoints, 1.0/dim) + 0.99999;
      for (k= dim; k--; )
	mult[k]= 0;
      for (i= 0; i < numpoints; i++) {
	for (k= 0; k < dim; k++) {
	  if (k == 0)
	    out1( mult[0] * meshn + mult[1] * (-meshm));
	  else if (k == 1)
	    out1( mult[0] * meshm + mult[1] * meshn);
	  else
	    out1( mult[k] * meshr );
	}
        fprintf (fp, "\n");
	for (k= 0; k < dim; k++) {
	  if (++mult[k] < nthroot)
	    break;
	  mult[k]= 0;
	}
      }
      if (iscdd)
	fprintf (fp, "end\nhull\n");
      return 0;
    }

    /* ============= regular points for 's' =============== */
    if (isregular && !islens) {
      if (dim != 2 && dim != 3) {
	fprintf(stderr, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
	exit(1);
      }
      if (!isaxis || radius == 0.0) {
	isaxis= 1;
	radius= 1.0;
      }
      if (dim == 3) {
        if (iscdd)
          out1( 1.0);
	out3n( 0.0, 0.0, -box);
	if (!isgap) {
          if (iscdd)
            out1( 1.0);
  	  out3n( 0.0, 0.0, box);
  	}
      }
      angle= 0.0;
      anglediff= 2.0 * PI/numpoints;
      for (i=0; i < numpoints; i++) {
	angle += anglediff;
	x= radius * cos (angle);
	y= radius * sin (angle);
	if (dim == 2) {
          if (iscdd)
            out1( 1.0);
	  out2n( x*box, y*box);
	}else {
	  norm= sqrt (1.0 + x*x + y*y);
          if (iscdd)
            out1( 1.0);
	  out3n( box*x/norm, box*y/norm, box/norm);
	  if (isgap) {
	    x *= 1-gap;
	    y *= 1-gap;
	    norm= sqrt (1.0 + x*x + y*y);
            if (iscdd)
              out1( 1.0);
	    out3n( box*x/norm, box*y/norm, box/norm);
	  }
	}
      }
      if (iscdd)
	fprintf (fp, "end\nhull\n");
      return 0;
    }
    /* ============= regular points for 'r Ln D2' =============== */
    if (isregular && islens && dim == 2) {
      double cos_0;
      
      angle= lensangle;
      anglediff= 2 * lensangle/(numpoints - 1);
      cos_0= cos (lensangle);
      for (i=0; i < numpoints; i++, angle -= anglediff) {
	x= radius * sin (angle);
  	y= radius * (cos (angle) - cos_0);
        if (iscdd)
          out1( 1.0);
	out2n( x*box, y*box);
	if (i != 0 && i != numpoints - 1) {
          if (iscdd)
            out1( 1.0);
	  out2n( x*box, -y*box);
	}
      }
      if (iscdd)
	fprintf (fp, "end\nhull\n");
      return 0;
    }
    /* ============= regular points for 'r Ln D3' =============== */
    if (isregular && islens && dim != 2) {
      if (dim != 3) {
	fprintf(stderr, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
	exit(1);
      }
      angle= 0.0;
      anglediff= 2* PI/numpoints;
      if (!isgap) {
	isgap= 1;
        gap= 0.5;
      }
      offset= sqrt (radius * radius - (1-gap)*(1-gap)) - lensbase;
      for (i=0; i < numpoints; i++, angle += anglediff) {
  	x= cos (angle);
	y= sin (angle);
        if (iscdd)
          out1( 1.0);
	out3n( box*x, box*y, 0);
	x *= 1-gap;
	y *= 1-gap;
        if (iscdd)
          out1( 1.0);
	out3n( box*x, box*y, box * offset);
        if (iscdd)
          out1( 1.0);
	out3n( box*x, box*y, -box * offset);
      }
      if (iscdd)
	fprintf (fp, "end\nhull\n");
      return 0;
    }
    /* ============= apex of 'Zn' distribution + gendim =============== */
    if (isaxis) {
      gendim= dim-1;
      if (iscdd)
        out1( 1.0);
      for (j=0; j < gendim; j++)
	out1( 0.0);
      out1( -box);
      fprintf (fp, "\n");
    }else if (islens) 
      gendim= dim-1;
    else
      gendim= dim;
    /* ============= generate random point in unit cube =============== */
    for (i=0; i < numpoints; i++) {
      norm= 0.0;
      for (j=0; j < gendim; j++) {
	randr= qh_RANDOMint;
	coord[j]= 2.0 * randr/randmax - 1.0;
	norm += coord[j] * coord[j];
      }
      norm= sqrt (norm);
      /* ============= dim-1 point of 'Zn' distribution ========== */
      if (isaxis) {
	if (!isgap) {
	  isgap= 1;
	  gap= 1.0;
	}
	randr= qh_RANDOMint;
	rangap= 1.0 - gap * randr/randmax;
	factor= radius * rangap / norm;
	for (j=0; j<gendim; j++)
	  coord[j]= factor * coord[j];
      /* ============= dim-1 point of 'Ln s' distribution =========== */
      }else if (islens && issphere) {
        if (!isgap) {
	  isgap= 1;
	  gap= 1.0;
	}
	randr= qh_RANDOMint;
	rangap= 1.0 - gap * randr/randmax;
	factor= rangap / norm;
	for (j=0; j<gendim; j++)
	  coord[j]= factor * coord[j];
      /* ============= dim-1 point of 'Ln' distribution ========== */
      }else if (islens && !issphere) {
        if (!isgap) {
	  isgap= 1;
	  gap= 1.0;
	}
	j= qh_RANDOMint % gendim;
	if (coord[j] < 0)
	  coord[j]= -1.0 - coord[j] * gap;
	else
	  coord[j]= 1.0 - coord[j] * gap;
      /* ============= point of 'l' distribution =============== */
      }else if (isspiral) {
	if (dim != 3) {
	  fprintf(stderr, "rbox error: spiral distribution is available only in 3d\n\n");
	  exit(1);
	}
	coord[0]= cos(2*PI*i/(numpoints - 1));
	coord[1]= sin(2*PI*i/(numpoints - 1));
	coord[2]= 2.0*(double)i/(double)(numpoints-1) - 1.0;
      /* ============= point of 's' distribution =============== */
      }else if (issphere) {
	factor= 1.0/norm;
	if (iswidth) {
  	  randr= qh_RANDOMint;
	  factor *= 1.0 - width * randr/randmax;
	}
	for (j=0; j<dim; j++)
	  coord[j]= factor * coord[j];
      }
      /* ============= project 'Zn s' point in to sphere =============== */
      if (isaxis && issphere) {
	coord[dim-1]= 1.0;
	norm= 1.0;
	for (j=0; j<gendim; j++)
	  norm += coord[j] * coord[j];
	norm= sqrt (norm);
	for (j=0; j<dim; j++)
	  coord[j]= coord[j] / norm;
	if (iswidth) {
  	  randr= qh_RANDOMint;
	  coord[dim-1] *= 1 - width * randr/randmax;
	}
      /* ============= project 'Zn' point onto cube =============== */
      }else if (isaxis && !issphere) {  /* not very interesting */
        randr= qh_RANDOMint;
	coord[dim-1]= 2.0 * randr/randmax - 1.0;
      /* ============= project 'Ln' point out to sphere =============== */
      }else if (islens) {
	coord[dim-1]= lensbase;
	for (j=0, norm= 0; j<dim; j++)
	  norm += coord[j] * coord[j];
	norm= sqrt (norm);
	for (j=0; j<dim; j++)
	  coord[j]= coord[j] * radius/ norm;
	coord[dim-1] -= lensbase;
	if (iswidth) {
  	  randr= qh_RANDOMint;
	  coord[dim-1] *= 1 - width * randr/randmax;
	}
	if (qh_RANDOMint > randmax/2)
	  coord[dim-1]= -coord[dim-1];
      /* ============= project 'Wn' point toward boundary =============== */
      }else if (iswidth && !issphere) {
	j= qh_RANDOMint % gendim;
	if (coord[j] < 0)
	  coord[j]= -1.0 - coord[j] * width;
	else
	  coord[j]= 1.0 - coord[j] * width;
      }
      /* ============= write point =============== */
      if (iscdd)
        out1( 1.0);
      for (k=0; k < dim; k++) 
	out1( coord[k] * box);
      fprintf (fp, "\n");
    }
    /* ============= write cube vertices =============== */
    if (addcube) {
      for (j=0; j<cubesize; j++) {
        if (iscdd)
          out1( 1.0);
	for (k=dim-1; k>=0; k--) {
	  if (j & ( 1 << k))
	    out1( cube);
	  else
	    out1( -cube);
	}
	fprintf (fp, "\n");
      }
    }
    /* ============= write diamond vertices =============== */
    if (adddiamond) {
      for (j=0; j<diamondsize; j++) {
        if (iscdd)
          out1( 1.0);
	for (k=dim-1; k>=0; k--) {
	  if (j/2 != k)
	    out1( 0.0);
	  else if (j & 0x1)
	    out1( diamond);
	  else
	    out1( -diamond);
	}
	fprintf (fp, "\n");
      }
    }
    if (iscdd)
      fprintf (fp, "end\nhull\n");
    return 0;
  } /* rbox */

/*------------------------------------------------
-outxxx - output functions
*/
int roundi( double a) {
  if (a < 0.0) {
    if (a - 0.5 < INT_MIN) {
      fprintf(stderr, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
      exit (1);
    }
    return a - 0.5;
  }else {
    if (a + 0.5 > INT_MAX) {
      fprintf(stderr, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
      exit (1);
    }
    return a + 0.5;
  }
} /* roundi */

void out1(double a) {

  if (isinteger) 
    fprintf(fp, "%d ", roundi( a+out_offset));
  else
    fprintf(fp, qh_REAL_1, a+out_offset);
} /* out1 */

void out2n( double a, double b) {

  if (isinteger)
    fprintf(fp, "%d %d\n", roundi(a+out_offset), roundi(b+out_offset));
  else
    fprintf(fp, qh_REAL_2n, a+out_offset, b+out_offset);
} /* out2n */

void out3n( double a, double b, double c) { 

  if (isinteger)
    fprintf(fp, "%d %d %d\n", roundi(a+out_offset), roundi(b+out_offset), roundi(c+out_offset));
  else
    fprintf(fp, qh_REAL_3n, a+out_offset, b+out_offset, c+out_offset);
} /* out3n */

/*-------------------------------------------------
-rand & srand- generate pseudo-random number between 1 and 2^31 -2
  from Park & Miller's minimimal standard random number generator
  Communications of the ACM, 31:1192-1201, 1988.
notes:
  does not use 0 or 2^31 -1
  this is silently enforced by qh_srand()
  copied from geom2.c
*/
static int seed = 1;  /* global static */

int qh_rand( void) {
#define qh_rand_a 16807
#define qh_rand_m 2147483647
#define qh_rand_q 127773  /* m div a */
#define qh_rand_r 2836    /* m mod a */
  int lo, hi, test;

  hi = seed / qh_rand_q;  /* seed div q */
  lo = seed % qh_rand_q;  /* seed mod q */
  test = qh_rand_a * lo - qh_rand_r * hi;
  if (test > 0)
    seed= test;
  else
    seed= test + qh_rand_m;
  return seed;
} /* rand */

void qh_srand( int newseed) {
  if (newseed < 1)
    seed= 1;
  else if (newseed >= qh_rand_m)
    seed= qh_rand_m - 1;
  else
    seed= newseed;
} /* qh_srand */