Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
I don't see where any serious roundoff error could be generated here. I can't find an error in the code, either. It's noticeable that some of the circumcentres come out exactly right and others don't. Are you sure you've initialized the $z$ components to zero? The in-class version of the code assumes that (in forming the cross products) whereas the other one doesn't. –  joriki Nov 18 '11 at 8:24
2  
Your 'in-class' formula is ambiguous because the cross product isn't associative, that is $a\times (b\times c) \neq (a\times b)\times c$. Which one is intended? –  Chris Taylor Nov 18 '11 at 8:26
    
Yes, @joriki, I set the z-components to zero by passing [P,0] to the cross(...) function (where P==[Px,Py]). –  tyblu Nov 18 '11 at 8:53
    
@Chris Taylor, I'm not really sure -- that is a serious omission. It gives me an idea as to how to answer my own question... –  tyblu Nov 18 '11 at 8:56
    
@tyblu: I see, sorry, I don't know the language you're using; should have guessed that that's what that's for... –  joriki Nov 18 '11 at 9:03

1 Answer 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!

share|improve this answer
    
Would be great if someone could double-check this, so I or another can correct the Wikipedia page (perhaps with this or similar derivation...). –  tyblu Nov 18 '11 at 16:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.