0
$\begingroup$

I'm writing a finite volume program over a 2D triangular mesh, and at one point I need to calculate the circumcenters of the triangles. The equation given in class and that on Wikipedia give different results.

In class:

$ \vec{r}_c = \vec{r}_i + \frac{( {\lVert \vec{r}_{ik} \rVert}^2 \vec{r}_{ij} - {\lVert \vec{r}_{ij} \rVert}^2 \vec{r}_{ik}) \times \vec{r}_{ik} \times \vec{r}_{ij}} { 2 {\lVert \vec{r}_{ik} \times \vec{r}_{ij} \rVert}^2} $

Wikipedia:

$ x_c = \frac{{\lVert A \rVert}^2(B_y-C_y)+{\lVert B \rVert}^2(C_y-A_y)+{\lVert C \rVert}^2(A_y-B_y)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}$

$y_c = \frac{{\lVert A \rVert}^2(C_x-B_x)+{\lVert B \rVert}^2(A_x-C_x)+{\lVert C \rVert}^2(B_x-A_x)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))} $

Where $\vec{r}_{ab} = \vec{r}_b - \vec{r}_a$, and $(\vec{r}_i,\vec{r}_j,\vec{r}_k)$ and $(A,B,C)$ are the triangle vertices.

I implented both after running into problems with the final solution (cross section of fluid velocity in a strangely-shaped tube), and found they gave different solutions.

Plotting the triangles and calculated circumcenters gives the following.

Using in-class formula:
MATLAB plot 1

Using Wikipedia formula:
MATLAB plot 2

Circumcenters were calculated with the following lines of code. In-class formula:

C = [I,0]+(cross(cross(((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0]),[K-I,0]),[J-I,0]))/(2*(norm(cross([K-I,0],[J-I,0])))^2); C = C(1,1:2); 

Wikipedia:

D = 2*(I(1)*(J(2)-K(2))+J(1)*(K(2)-I(2))+K(1)*(I(2)-J(2))); xc = ((I(2)^2+I(1)^2)*(J(2)-K(2))+(J(2)^2+J(1)^2)*(K(2)-I(2))+(K(2)^2+K(1)^2)*(I(2)-J(2)))/D; yc = ((I(1)^2+I(2)^2)*(J(1)-K(1))+(J(1)^2+J(2)^2)*(K(1)-I(1))+(K(1)^2+K(2)^2)*(I(1)-J(1)))/D; C = [xc -yc]; 

I figure that there must be something wrong with my implementation of the in-class formula, as many acute triangles don't have a calculated circumcenter inside of them, but I can't find the error. Can you?

PS: Yeah, I see the more elegant methods of implementing the Wikipedia formula. Meh.

  • 0
    @tyblu: I see, sorry, I don't know the language you're using; should have guessed that that's what that's for...2019-02-26

1 Answers 1

1

Substituting variables $(\vec{r}_i,\vec{r}_j,\vec{r}_k) = (A,B,C)$, with all z-coordinates equal to zero, into the 'in-class' formula gives the following:

$\vec{r}_c = A+\frac{({\lVert C-A \rVert}^2(B-A)-{\lVert B-A \rVert}^2(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$

$\vec{r}_c = A+\frac{(({\lVert A \rVert}^2+{\lVert C \rVert}^2-2A \cdot C)(B-A)-({\lVert A \rVert}^2+{\lVert B \rVert}^2-2A \cdot B)(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2} $

$\vec{r}_c = A+\frac{\left({\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]\right) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2} $

Taking a hint from the denominator, I assume the cross product follows the following order. (The terms to the right combine to form a unit vector.)

$\vec{r}_c = A+\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]}{2{\lVert (C-A)\times(B-A) \rVert}} \times \frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}$

$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(A\cdot C)(B-A)-(A\cdot B)(C-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$

Let $Y=B-A$, $Z=C-A$, $A\cdot B=X\cdot Y$, and $A\cdot C=X\cdot Z$, giving the following:

$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(X\cdot Z)Y-(X\cdot Y)Z}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$

$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{X\times(Y\times Z)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$

$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}+X\times\frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$

$\vec{r}_c = A+(X\times\hat z)\times\hat z +\frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$

Now, if I knew how to solve for $-X$, I'd show the following:

$\vec{r}_c = \frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$

In other words, the two formulas are off by a factor of -1, and it appears to be the Wikipedia one that maps improperly, as I had been 'fixing' it by setting yc=-yc -- the x-axis is symmetric, so I didn't notice that it was also flipped.

The 'in-class' formula needs to be implemented like this: (@Chris Taylor hit it on the head!)

C = [I,0]+(cross((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0],cross([K-I,0],[J-I,0])))/(2*(norm(cross([K-I,0],[J-I,0])))^2); C = C(1,1:2); 

I'll have to try it out later, as it's time to head to class!

  • 0
    Would be great if someone could double-check this, so I or another can correct [the Wikipedia page](http://en.wikipedia.org/wiki/Circumcenter#Cartesian_coordinates) (perhaps with this or similar derivation...).2011-11-18