Dropped the rotation of the shape function gradients

For a finite element, the discretised weak form terms can be calculated in any local Cartesian coordinate system defined on the element. For example, the Laplace term of an inclined element is computed as

    \int (\nabla \text{N})^{T} \mathbf{K}_l (\nabla \text{N}) \text{d}\Omega

where \nabla \text{N} denotes the gradients of the shape functions with respect to the local coordinates defined on the element (not the coordinates on the reference element for numerical quadrature), and \mathbf{K}_l is the coefficient matrix in the local coordinate system. If the the coefficient matrix is given in the global coordinates system, say \mathbf{K}_g, the local coefficient matrix \mathbf{K}_l is computed as

 \mathbf{K}_g = \mathbf{R} \mathbf{K}_l  \mathbf{R}^{T}

with \mathbf{R} the rotation matrix from the global to the local coordinate system. That means that the Laplace term can be written as

    \int (\nabla \text{N})^{T} \mathbf{R} \mathbf{K}_g  \mathbf{R}^{T} (\nabla \text{N}) \text{d}\Omega

or

    \int ( \mathbf{R} \nabla\text{N})^{T}  \mathbf{K}_g  (\mathbf{R} \nabla \text{N}) \text{d}\Omega

For the later, the coordinate rotation is applied to the matrix of the gradients of the shape functions.

In the existing implementation of LIE, the Laplace term of computation is computed as

    \int ( \mathbf{R}\nabla \text{N} )^{T} \mathbf{R}  \mathbf{K}_g \mathbf{R}^{T} (\mathbf{R} \nabla \text{N}) \text{d}\Omega

which is wrong. Luckily, since \mathbf{R} is orthogonal and normalized, and the local coefficient is always a scalar number of k, the final computation result is correct as:

  \int ( \mathbf{R}\nabla \text{N} )^{T} \mathbf{R}  \mathbf{K}_g \mathbf{R}^{T} (\mathbf{R} \nabla \text{N}) \text{d}\Omega =  \int k ( \nabla \text{N} )^{T}   \mathbf{K}_g  \nabla \text{N} \text{d}\Omega

This MR presents the correction of the coordinate rotation as:

1. dropped the rotation of the matrix of the gradients of the shape functions, which is computed in the element based local coordinate system. This influences all local assemblers.
2. removed the rotation of the coefficient matrix in LIE. The Laplace coefficient of the inclined element, or the related parameter of fractures, is always isotropic.
3. corrected the size of the velocity vector in LIE, i.e from 3 to GlobalDim.
4. moved the application of rotation from dNdx to the velocity in the advection term of the BHE process.

The changes in LIE also correct the orientation of the velocity in LIE.

With the changes, the Laplace term of an inclined element is computed as the regular element in the local assembler as

   \int k ( \nabla \text{N} )^{T}  \nabla \text{N} \text{d}\Omega

except that it is computed in the local Cartesian coordinate system defined on the element.

This MR is a follow-up of MR !3672 (merged), and it closes !3655 (closed).

1. Feature description was added to the changelog
2. Tests covering your feature were added? Yes, with the exiting benchmarks. The reference files of some HM/LIE are updated due to the correction of the flow orientation (inversed), and the change of the velocity vector size.
3. Any new feature or behavior change was documented?