v0.14.0
Loading...
Searching...
No Matches
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
30void tetcircumcenter_tp(a, b, c, d, circumcenter, xi, eta, zeta)
31double a[3];
32double b[3];
33double c[3];
34double d[3];
35double circumcenter[3];
36double *xi;
37double *eta;
38double *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
constexpr double a
double eta
const double c
speed of light (cm/ns)
double zeta
Viscous hardening.
Definition plastic.cpp:177
void tetcircumcenter_tp(a, b, c, d, circumcenter, double *xi, double *eta, double *zeta)