anchor.diff:                    s_start = (n - 1)*3+1
                     s_end   = s_start+2
                                       
  +                     ! note that we multiply by MeshMap%MapMotions(i)%shape_fn(n1) after forming the first two terms of this matrix
                     MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) = OuterProduct( MeshMap%DisplacedPosition(:,i,n1), Src%RotationAcc(:,n) )
                     tmp=dot_product( MeshMap%DisplacedPosition(:,i,n1), Src%RotationAcc(:,n) )
                     do j=0,2
  
positive.code:                   n = MeshMap%MapMotions(i)%OtherMesh_Element
                    d_start = (i-1)*3+1
                    d_end   = d_start+2
                    s_start = (n - 1)*3+1
                    s_end   = s_start+2
                                      
                    tmpVec=cross_product(  Src%RotationVel(:,n), MeshMap%DisplacedPosition(:,i,1) )
                    MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) = MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) + OuterProduct( tmpVec, Src%RotationVel(:,n) )
                    
                    tmp = dot_product( Src%RotationVel(:,n), MeshMap%DisplacedPosition(:,i,1)) 
                    tmpVec = tmp*Src%RotationVel(:,n)
                    MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) = MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) + SkewSymMat( tmpVec )                  
                 end do
              end if
                                 
  
negative.code:                   s_start = (n - 1)*3+1
                    s_end   = s_start+2
                                      
                    tmpVec=cross_product(  Src%RotationVel(:,n), MeshMap%DisplacedPosition(:,i,n1) ) * MeshMap%MapMotions(i)%shape_fn(n1)
                    MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) = MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) + OuterProduct( tmpVec, Src%RotationVel(:,n) )
                    
                    tmp = dot_product( Src%RotationVel(:,n), MeshMap%DisplacedPosition(:,i,n1)) * MeshMap%MapMotions(i)%shape_fn(n1)
                    tmpVec = tmp*Src%RotationVel(:,n) 
                    MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) = MeshMap%dM%ta1( d_start:d_end, s_start:s_end ) - SkewSymMat( tmpVec )                   
                 end do
                    
              end do
  