anchor.diff:        end forall
       enddo
   
  -    total_mass = 0.0d0
  -    do i = 1, nn*ndof, ndof
  -      do j = 1, nn*ndof, ndof
  -        total_mass = total_mass + mass(j,i)
  -      enddo
  -    enddo
  -
       is_lumped = .true.
  -    if(is_lumped)then
  -      diag_mass = 0.0d0
  -      do i = 1, nn*ndof, ndof
  -        diag_mass = diag_mass + mass(i,i)
  -      enddo
  -
  -      diag_mass = 1.0d0/diag_mass
  -      do i = 1, nn*ndof
  -        lumped(i) = lumped(i) + mass(i,i)*total_mass*diag_mass
  -      enddo
  -
  -      mass = 0.0d0
  -      do i = 1, nn*ndof
  -        mass(i,i) = lumped(i)
  -      enddo
  -    endif
  +    if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
     end subroutine mass_C2
   
     subroutine mass_C3(etype, nn, ecoord, gausses, mass, lumped, temperature)
  
positive.code:       end forall
      enddo
  
      total_mass = 0.0d0
      do i = 1, nn*ndof, ndof
        do j = 1, nn*ndof, ndof
  
negative.code:       enddo
      enddo
  
      diag_mass = 0.0d0
      do i = 1, nn*ndof, ndof
        diag_mass = diag_mass + mass(i,i)
      enddo
  
      diag_mass = 1.0d0/diag_mass
      do i = 1, nn*ndof
        lumped(i) = lumped(i) + mass(i,i)*total_mass*diag_mass
      enddo
  
      mass = 0.0d0
      do i = 1, nn*ndof
        mass(i,i) = lumped(i)
      enddo
    end subroutine get_lumped_mass
  
    function get_length(ecoord)
      use hecmw
  