I'm trying to understand why integrating \begin{align} dw =c_{ijkl}\epsilon_{kl}\,d\epsilon_{ij}, \end{align} gives \begin{align} w =\frac{1}{2}c_{ijkl}\epsilon_{kl}\epsilon_{ij}, \tag{1} \end{align} given that the constants have the properties $$c_{ijkl}=c_{klij}=c_{lkij}=c_{klji}=c_{jikl}=c_{jilk}.$$
If the indices are $i,j,k,l=1,2$ then there are terms like $$ c_{1111}\epsilon_{11}\,d\epsilon_{11}, \ \ \ c_{2222}\epsilon_{22}\,d\epsilon_{22},\ \ \ \ 4\times c_{1212}\epsilon_{12}\,d\epsilon_{12}, $$ which seem to agree with (1). However, the terms (multiples due to terms being the same through the symmetries of $c_{ijkl}$) $$ 4\times c_{1222}\epsilon_{22}\,d\epsilon_{12},\ \ \ 4\times c_{1112}\epsilon_{12}\,d\epsilon_{11}, \ \ \ 2\times c_{1122}\epsilon_{22}\,d\epsilon_{11}, $$ would not seem to give the factor of 1/2 in (1) when integrated (integrating constants).
Where am I going wrong?
By the product rule $$d(\epsilon_{ij}\epsilon_{kl})=d\epsilon_{ij}\epsilon_{kl}+\epsilon_{ij}d\epsilon_{kl}$$ contracting with $c_{ijkl}$ $$c_{ijkl}d(\epsilon_{ij}\epsilon_{kl})=c_{ijkl}d\epsilon_{ij}\epsilon_{kl}+c_{ijkl}\epsilon_{ij}d\epsilon_{kl}$$ Then because of the symmetry $c_{ijkl}=c_{klij}$ the two terms in the RHS are equal, so that $$c_{ijkl}d(\epsilon_{ij}\epsilon_{kl})=2c_{ijkl}d\epsilon_{ij}\epsilon_{kl}$$ $$\implies c_{ijkl}\epsilon_{ij}d\epsilon_{kl}=\frac{1}{2}d(c_{ijkl}\epsilon_{ij}\epsilon_{kl})$$ Where $c_{ijkl}$ can get inside $d$ because it is constant. Integrate both sides to get the desired result.