v0.9.1
tetcircumcenter.c
Go to the documentation of this file.
1 // Jonathan R Shewchuk <jrs+@cs.cmu.edu>
2 /*****************************************************************************/
3 /* */
4 /* tetcircumcenter() Find the circumcenter of a tetrahedron. */
5 /* */
6 /* The result is returned both in terms of xyz coordinates and xi-eta-zeta */
7 /* coordinates, relative to the tetrahedron's point a' (that is, a' is */
8 /* the origin of both coordinate systems). Hence, the xyz coordinates */
9 /* returned are NOT absolute; one must add the coordinates of a' to */
10 /* find the absolute coordinates of the circumcircle. However, this means */
11 /* that the result is frequently more accurate than would be possible if */
12 /* absolute coordinates were returned, due to limited floating-point */
13 /* precision. In general, the circumradius can be computed much more */
14 /* accurately. */
15 /* */
16 /* The xi-eta-zeta coordinate system is defined in terms of the */
17 /* tetrahedron. Point a' is the origin of the coordinate system. */
18 /* The edge ab' extends one unit along the xi axis. The edge ac' */
19 /* extends one unit along the eta axis. The edge ad' extends one unit */
20 /* along the zeta axis. These coordinate values are useful for linear */
21 /* interpolation. */
22 /* */
23 /* If xi' is NULL on input, the xi-eta-zeta coordinates will not be */
24 /* computed. */
25 /* */
26 /*****************************************************************************/
27
28 #include <stdlib.h>
29
30 void tetcircumcenter_tp(a, b, c, d, circumcenter, xi, eta, zeta)
31 double a[3];
32 double b[3];
33 double c[3];
34 double d[3];
35 double circumcenter[3];
36 double *xi;
37 double *eta;
38 double *zeta;
39 {
40  double xba, yba, zba, xca, yca, zca, xda, yda, zda;
41  double balength, calength, dalength;
42  double xcrosscd, ycrosscd, zcrosscd;
43  double xcrossdb, ycrossdb, zcrossdb;
44  double xcrossbc, ycrossbc, zcrossbc;
45  double denominator;
46  double xcirca, ycirca, zcirca;
47
48  /* Use coordinates relative to point a' of the tetrahedron. */
49  xba = b[0] - a[0];
50  yba = b[1] - a[1];
51  zba = b[2] - a[2];
52  xca = c[0] - a[0];
53  yca = c[1] - a[1];
54  zca = c[2] - a[2];
55  xda = d[0] - a[0];
56  yda = d[1] - a[1];
57  zda = d[2] - a[2];
58  /* Squares of lengths of the edges incident to a'. */
59  balength = xba * xba + yba * yba + zba * zba;
60  calength = xca * xca + yca * yca + zca * zca;
61  dalength = xda * xda + yda * yda + zda * zda;
62  /* Cross products of these edges. */
63  xcrosscd = yca * zda - yda * zca;
64  ycrosscd = zca * xda - zda * xca;
65  zcrosscd = xca * yda - xda * yca;
66  xcrossdb = yda * zba - yba * zda;
67  ycrossdb = zda * xba - zba * xda;
68  zcrossdb = xda * yba - xba * yda;
69  xcrossbc = yba * zca - yca * zba;
70  ycrossbc = zba * xca - zca * xba;
71  zcrossbc = xba * yca - xca * yba;
72
73  /* Calculate the denominator of the formulae. */
74 #ifdef EXACT
75  /* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html */
76  /* to ensure a correctly signed (and reasonably accurate) result, */
77  /* avoiding any possibility of division by zero. */
78  denominator = 0.5 / orient3d(b, c, d, a);
79 #else
80  /* Take your chances with floating-point roundoff. */
81  denominator = 0.5 / (xba * xcrosscd + yba * ycrosscd + zba * zcrosscd);
82 #endif
83
84  /* Calculate offset (from a') of circumcenter. */
85  xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
86  denominator;
87  ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
88  denominator;
89  zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
90  denominator;
91  circumcenter[0] = xcirca;
92  circumcenter[1] = ycirca;
93  circumcenter[2] = zcirca;
94
95  if (xi != (double *) NULL) {
96  /* To interpolate a linear function at the circumcenter, define a */
97  /* coordinate system with a xi-axis directed from a' to b', */
98  /* an eta-axis directed from a' to c', and a zeta-axis directed */
99  /* from a' to `d'. The values for xi, eta, and zeta are computed */
100  /* by Cramer's Rule for solving systems of linear equations. */
101  *xi = (xcirca * xcrosscd + ycirca * ycrosscd + zcirca * zcrosscd) *
102  (2.0 * denominator);
103  *eta = (xcirca * xcrossdb + ycirca * ycrossdb + zcirca * zcrossdb) *
104  (2.0 * denominator);
105  *zeta = (xcirca * xcrossbc + ycirca * ycrossbc + zcirca * zcrossbc) *
106  (2.0 * denominator);
107  }
108 }
109
void tetcircumcenter_tp(a, b, c, d, circumcenter, double *xi, double *eta, double *zeta)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27