Stress calculation

The stress calculation is an uneventful subroutine getstress() . It performs the double inner product of the stiffness matrix with the strain tensor and yields the stress tensor. Note again the nearly mathematical formulation as in (5). Since the stiffness matrix is generally sparse we only use nonzero elements in the calculation, the scale applied to certain components is due to the $\nabla$ calculation according to formula (2).

# getstress:    Hooke's Law
# get stresses from strain
#%
subroutine getstress(strain,stiff,stress,bnd)
implicit none
#include "commons"
TENSOR(strain)
TENSORLAY(strain)
TENSOR(stress)
TENSORLAY(stress)
STIFF(stiff)
STIFFLAY(stiff)
#TENSOR(stiff)
#TENSORLAY(stiff)
BOUND(bnd)
BOUNDLAY(bnd)
integer i,j,k,l

stress=0.

do i=1,ntcomp {
do j=1,ntcomp {               k = ccin(tcomp(i),tcomp(j))
   if (k!=0)  {  
      stress(i,:,:,:) = stiff(k,:,:,:) * stiffmask(tcomp(i),tcomp(j))* 
                        strain(j,:,:,:) + stress(i,:,:,:)
      }
}}

return
end