v0.9.1
tetcircumcenter.c File Reference
#include <stdlib.h>

Go to the source code of this file.

## Functions

void tetcircumcenter_tp (a, b, c, d, circumcenter, double *xi, double *eta, double *zeta)

## ◆ tetcircumcenter_tp()

 void tetcircumcenter_tp ( a , b , c , d , circumcenter , double * xi, double * eta, double * zeta )

Definition at line 30 of file tetcircumcenter.c.

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 }
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