The "Discrete Shells" paper1 defines a bending energy for triangle meshes by summing up an energy contribution for each interior edge $i$: $$ E = \sum_{i} \frac{k_i 3 \|\bar{\mathbf{e}}_i\|^2}{\bar{A}_i} \left( \theta_i - \bar{\theta}_i \right)^2 $$ where $k_i$ is the bending stiffness, $\mathbf{e}_i$ is the edge vector, $\bar{A}_i$ is the sum of areas of the two incident triangles, $\bar{\theta}_i$ is the dihedral angle formed by the two incident triangles at rest, and $\theta_i$ is the dihedral angle in the deformed state. I found the notation in "Discrete Bending Forces and Their Jacobians" to be clearest.2
The local information required for computing the contribution at each edge $i$ are the four vertices of the two incident triangles, forming a diamond: $\mathbf{d}_i = \{d_{i1} , d_{i2} , d_{i3} , d_{i4}\}$. From their rest positions, we can easily compute $\bar{A}_i$ and $\|\bar{\mathbf{e}}_i\|^2$: $$ \|\bar{\mathbf{e}}_i\|^2 = \|\bar{\mathbf{v}}_4 - \bar{\mathbf{v}}_3\|^2, \quad\quad \bar{A}_i = \frac{1}{2}\|\bar{\mathbf{e}} \times (\bar{\mathbf{v}}_1 - \bar{\mathbf{v}}_3)\| + \frac{1}{2}\|(\bar{\mathbf{v}}_2 - \bar{\mathbf{v}}_3) \times \mathbf{e} \| $$ using the convention that $\mathbf{v}_3$ and $\mathbf{v}_4$ are the vertices of shared edge $i$ between triangles formed with vertices $\mathbf{v}_1$ and $\mathbf{v}_2$. These unsigend quantities don't require deciding any notion of orientation for the triangles.
The dihedral angle on the other hand could either be the "interior" or "exterior" angle between the two triangle planes. Meanwhile, different formulas or conventions for computing the dihedral angle could lead to different ranges of values: $(0,2\pi)$ or $(-\pi,\pi)$. Since the "Discrete shells" energy is just a squared difference in rest and deformed angles, it would be convenient if our formula and convention for the angle does not depend on the order of the vertices in the diamond.
That is, we would like the squared differences of rest and deformed dihedral angles to be the same regardless of diamond order: $$ (\theta_{1234} - \bar{\theta}_{1234})^2 = (\theta_{1243} - \bar{\theta}_{1243})^2 = (\theta_{2134} - \bar{\theta}_{2134})^2 = (\theta_{2143} - \bar{\theta}_{2143})^2. $$ The particular $\theta,\bar{\theta}$ values might be different (e.g., in sign) but if we use the same order and formula for rest and deformed we only need to worry about the squared difference.
To achieve this, we'll use a range of values $\theta = (-\pi,\pi)$ and ensure that $\theta = 0$ corresponds to the case where the two triangles are coplanar and non-overlapping: flat. And $\theta \rightarrow \pm \pi$ corresponds to the triangles making an extremely sharp wedge.
The formula we'll use is: $$ \theta_{1234} = \text{atan2}\left(\ \hat{\mathbf{e}} \cdot \left(\mathbf{n}_1 \times \mathbf{n}_2\right),\ \mathbf{n}_1 \cdot \mathbf{n}_2\right) $$ where we specifically make use the sign conventions of atan2 and observe that only $\hat{\mathbf{e}}$ needs to be normalized to unit length. The ``normals'' (which are specifically oriented for this diamond formula but need not match the normal convention of the triangle mesh) are computed matching the subterms of $\bar{A}_i$ above: $$ \mathbf{n}_1 = \mathbf{e} \times \left(\mathbf{v}_1 - \mathbf{v}_3\right), \quad\quad \mathbf{n}_2 = \left(\mathbf{v}_2 - \mathbf{v}_3\right) \times \mathbf{e}. $$
The point of these formulas is not to change the bending energy, but rather make sure that it doesn't depend on the order of indices defining the same diamond. As a result it makes it convenient for non-manifold or non-orientable triangle meshes.
Finally, here's some minimal matlab code to demonstrate this in action:
for phi0 = pi*[-0.9 -0.5 -0.1 0 0.1 0.5 0.9]
% rest positions
V0 = [ 0 0 0; 1+cos(phi0) 0 sin(phi0); 1 -1 0; 1 1 0; ];
% All possible diamond orders.
D = [3 4 1 2;
4 3 1 2;
3 4 2 1;
4 3 2 1];
clf;
hold on;
tsurf([3 4 1;3 4 2],V0,'FaceColor','k',falpha(1.0,1),'LineWidth',3);
for phi = linspace(-pi+eps,pi-eps,40)
% deformed positions
V = [ 0 0 0; 1+cos(phi) 0 sin(phi); 1 -1 0; 1 1 0; ];
[E,th0,th] = discrete_shell(V0,V,D);
assert(var(E)<1e-10);
tsurf([3 4 2],V,'CData',mean(E),falpha(0.25,0.5));
end
hold off;
axis equal;
view(-6,4);camproj('persp');
axis([0 2 -1 1 -1 1]);
drawnow;
pause
end
function [E,th0,th] = discrete_shell(V0,V,D)
[th0,e,n1,n2] = dihedral_angle(V0,D);
A = vecnorm(n1,2,2) + vecnorm(n2,2,2);
th = dihedral_angle(V,D);
E = 3*sum(e.^2,2)./A.*(th-th0).^2;
end
function [th,e,n1,n2] = dihedral_angle(V,D)
e = V(D(:,4),:)-V(D(:,3),:);
u = e./vecnorm(e,2,2);
n1 = cross(u,V(D(:,1),:)-V(D(:,3),:),2);
n2 = cross(V(D(:,2),:)-V(D(:,3),:),u,2);
th = atan2(dot(u,cross(n1,n2,2),2), dot(n1,n2,2));
end
1 There's a slightly different version of the same "Discrete Shells" paper online with a diffferent formula: $$ E = \sum_{i} \|\bar{\mathbf{e}}_i\| \left( \theta_i - \bar{\theta}_i \right)^2. $$ This version does not appear to be considered canon. <
2 As a side note, Tamstorf & Grinspun have a formula for determining the bending stiffness $k$ from a Young's modulus $Y$, Poisson's ratio $\nu$ and thickness $h$: $$ k = \frac{Yh^3}{24(1-\nu^2)}, $$ which can be converted to Lamé parameters $\lambda,\mu$: $$ k = \frac{h^3 \mu (\lambda + \mu)}{6\lambda + 12\mu}. $$