diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90
index 6aa6db48d..d379a6255 100644
--- a/src/FEM_mech.f90
+++ b/src/FEM_mech.f90
@@ -327,6 +327,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
    materialpoint_F, &
    materialpoint_P
  use math, only: &
+   math_I3, &
    math_det33, &
    math_inv33  
  use FEsolving, only: &
@@ -362,7 +363,8 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
  CHKERRQ(ierr)
  call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
  call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
- call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
+ !call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
+ call VecCopy(xx_local,x_local,ierr); CHKERRQ(ierr)
  do field = 1, dimPlex; do face = 1, mesh_Nboundaries
    if (params%fieldBC%componentBC(field)%Mask(face)) then
      call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
@@ -395,6 +397,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
        enddo
      enddo
      materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
+       math_I3 + &
        reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
    enddo    
    if (BBarStabilisation) then
diff --git a/src/FEM_utilities.f90 b/src/FEM_utilities.f90
index 1db950e63..6e2eed2fb 100644
--- a/src/FEM_utilities.f90
+++ b/src/FEM_utilities.f90
@@ -411,7 +411,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
    call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
    CHKERRQ(ierr)
    do dof = offset+comp+1, offset+numDof, numComp
-     localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
+     localArray(dof) = BCValue + BCDotValue*timeinc
    enddo
  enddo    
  call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)