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} $
$ 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:
Using Wikipedia formula:
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.