Solution (source code)

= Solution

The surface element relation is <Nanson's formula>. The <force> on a corresponding pair of current and reference surface elements must agree:
$$
t\,da=(S\circ y)n\,da=(S\circ y)\operatorname{cof}F\,N\,dA=T N\,dA.
$$
Therefore the first <Piola-Kirchhoff stress tensor> is
$$
\boxed{T=(S\circ y)\operatorname{cof}F=J(S\circ y)F^{-T}.}
$$
It measures <force> per reference <area>, with a current <force> direction and a reference normal direction.

In components, differentiate $T_{iA}=S_{ij}(y(x))(\operatorname{cof}F)_{jA}$. The <Piola identity> says $\partial_A(\operatorname{cof}Dy)_{jA}=0$. One proof writes the <cofactor matrix> as $\frac12\varepsilon_{jkl}\varepsilon_{ABC}y_{k,B}y_{l,C}$; differentiating in $A$ produces symmetric second <derivatives> contracted against antisymmetric epsilon symbols, hence zero. The <chain rule> consequently gives
$$
\partial_AT_{iA}=\frac{\partial S_{ij}}{\partial y_k}(y)F_{kA}(\operatorname{cof}F)_{jA}
=J\frac{\partial S_{ij}}{\partial y_j}(y),
$$
since $F_{kA}(\operatorname{cof}F)_{jA}=J\delta_{kj}$. Thus, with the coordinate composition made explicit,
$$
\boxed{\operatorname{Div}_xT=J(\operatorname{div}_yS)\circ y.}
$$
The distinction between reference <divergence> and current <divergence> is essential to this identity.