if(rheology.planeStress())
  {
    //- add higher order terms
    volScalarField higherTerms = -0.5*volTensorField(gradDU&gradDU.T()).component(tensor::ZZ);

    forAll(gradDU.internalField(), celli)
      {
	gradDU.internalField()[celli][tensor::ZZ] =
	  ((-C.internalField()[celli][symmTensor4thOrder::XXZZ]*DEpsilon.internalField()[celli][symmTensor::XX]
	    - C.internalField()[celli][symmTensor4thOrder::YYZZ]*DEpsilon.internalField()[celli][symmTensor::YY])
                    /
	   C.internalField()[celli][symmTensor4thOrder::ZZZZ])
	  -higherTerms.internalField()[celli];
      }
    forAll(gradDU.boundaryField(), patchi)
      {
	forAll(gradDU.boundaryField()[patchi], facei)
	  {
	    gradDU.boundaryField()[patchi][facei][tensor::ZZ] =
	      ((-C.boundaryField()[patchi][facei][symmTensor4thOrder::XXZZ]*DEpsilon.boundaryField()[patchi][facei][symmTensor::XX]
		- C.boundaryField()[patchi][facei][symmTensor4thOrder::YYZZ]*DEpsilon.boundaryField()[patchi][facei][symmTensor::YY])
                        /
	       C.boundaryField()[patchi][facei][symmTensor4thOrder::ZZZZ])
	      - higherTerms.boundaryField()[patchi][facei];
	  }
      }
  }