+ if (u > 0)
+ result = partial (node1, node2, v / (double) u);
+ else
+ result = node1;
+
+ gp->nodes[node].position = normalise (result);
+ gp->nodes[node].initial_position = gp->nodes[node].position;
+ gp->nodes[node].normal = zero_vector;
+ node++;
+ }
+ }
+
+ /*
+ * Determine which nodes make up each face. The complexity is caused
+ * by having to determine the correct nodes for the edges of the
+ * tetrahedron since the common nodes on the edges are only calculated
+ * once (see above).
+ */
+ for (u = 0; u < (nodes_on_edge - 1); u++)
+ {
+ for (v = 0; v <= u; v++)
+ {
+ {
+ if (side < 2)
+ {
+ gp->faces[face].node1 = base + ((u * (u + 1)) / 2) + v;
+ gp->faces[face].node2 =
+ base + ((u + 1) * (u + 2)) / 2 + v + 1;
+ gp->faces[face].node3 =
+ base + ((u + 1) * (u + 2)) / 2 + v;
+
+ if ((side == 1) && (u == (nodes_on_edge - 2)))
+ {
+ gp->faces[face].node3 =
+ ((u + 1) * (u + 2)) / 2 +
+ nodes_on_edge - v - 1;
+ gp->faces[face].node2 =
+ ((u + 1) * (u + 2)) / 2 +
+ nodes_on_edge - v - 2;
+ }
+ }
+ else if (side < 3)
+ {
+ gp->faces[face].node1 =
+ base + (((u - 1) * u) / 2) + v - 1;
+ gp->faces[face].node2 = base + ((u) * (u + 1)) / 2 + v;
+ gp->faces[face].node3 =
+ base + ((u) * (u + 1)) / 2 + v - 1;
+
+ if (u == (nodes_on_edge - 2))
+ {
+ int n = nodes_on_edge - v - 1;
+ gp->faces[face].node2 =
+ ((nodes_on_edge *
+ (nodes_on_edge + 1)) / 2) +
+ ((n - 1) * (n + 0)) / 2;
+ gp->faces[face].node3 =
+ ((nodes_on_edge *
+ (nodes_on_edge + 1)) / 2) +
+ ((n + 0) * (n + 1)) / 2;
+ }
+ if (v == 0)
+ {
+ gp->faces[face].node1 = (((u + 1) * (u + 2)) / 2) - 1;
+ gp->faces[face].node3 = (((u + 2) * (u + 3)) / 2) - 1;
+ }
+ }
+ else
+ {
+ gp->faces[face].node1 =
+ base + (((u - 2) * (u - 1)) / 2) + v - 1;
+ gp->faces[face].node2 = base + ((u - 1) * u) / 2 + v;
+ gp->faces[face].node3 = base + ((u - 1) * u) / 2 + v - 1;
+
+ if (v == 0)
+ {
+ gp->faces[face].node1 =
+ base2 + ((u * (u + 1)) / 2) - 1;
+ gp->faces[face].node3 =
+ base2 + ((u + 1) * (u + 2)) / 2 - 1;
+ }
+ if (u == (nodes_on_edge - 2))
+ {
+ gp->faces[face].node3 =
+ ((nodes_on_edge *
+ (nodes_on_edge + 1)) / 2) +
+ ((v + 1) * (v + 2)) / 2 - 1;
+ gp->faces[face].node2 =
+ ((nodes_on_edge *
+ (nodes_on_edge + 1)) / 2) +
+ ((v + 2) * (v + 3)) / 2 - 1;
+ }
+ if (v == u)
+ {
+ gp->faces[face].node1 = (u * (u + 1)) / 2;
+ gp->faces[face].node2 = ((u + 1) * (u + 2)) / 2;
+ }
+ }
+ face++;
+ }
+
+ if (v < u)
+ {
+ if (side < 2)
+ {
+ gp->faces[face].node1 = base + ((u * (u + 1)) / 2) + v;
+ gp->faces[face].node2 =
+ base + ((u * (u + 1)) / 2) + v + 1;
+ gp->faces[face].node3 =
+ base + (((u + 1) * (u + 2)) / 2) + v + 1;
+
+ if ((side == 1) && (u == (nodes_on_edge - 2)))
+ {
+ gp->faces[face].node3 =
+ ((u + 1) * (u + 2)) / 2 +
+ nodes_on_edge - v - 2;
+ }
+ }
+ else if (side < 3)
+ {
+ gp->faces[face].node1 =
+ base + ((u * (u - 1)) / 2) + v - 1;
+ gp->faces[face].node2 = base + ((u * (u - 1)) / 2) + v;
+ gp->faces[face].node3 = base + ((u * (u + 1)) / 2) + v;
+
+ if (u == (nodes_on_edge - 2))
+ {
+ int n = nodes_on_edge - v - 1;
+ gp->faces[face].node3 =
+ ((nodes_on_edge *
+ (nodes_on_edge + 1)) / 2) +
+ ((n + 0) * (n - 1)) / 2;
+ }
+ if (v == 0)
+ {
+ gp->faces[face].node1 = (((u + 1) * (u + 2)) / 2) - 1;
+ }
+ }
+ else
+ {
+ gp->faces[face].node1 =
+ base + (((u - 2) * (u - 1)) / 2) + v - 1;
+ gp->faces[face].node2 =
+ base + (((u - 2) * (u - 1)) / 2) + v;
+ gp->faces[face].node3 = base + (((u - 1) * u) / 2) + v;
+
+ if (v == 0)
+ {
+ gp->faces[face].node1 = base2 + (u * (u + 1)) / 2 - 1;
+ }
+ if (u == (nodes_on_edge - 2))
+ {
+ gp->faces[face].node3 =
+ ((nodes_on_edge * (nodes_on_edge + 1)) / 2) +
+ ((v + 2) * (v + 3)) / 2 - 1;
+ }
+ if (v == (u - 1))
+ {
+ gp->faces[face].node2 = (u * (u + 1)) / 2;
+ }
+ }
+ face++;
+ }