Avoid out-of-bound index.

The test `rotflow` occasionally fails and valgrind flags an out-of-range
error when running that test.

That error occurs because bubble degrees of freedom with indices larger
than the number of nodes are used. Avoid that error by optionally not
including the indices of bubble degrees of freedom in `mGetElementDOFs`.
Only use non-bubble DOFs in `UpdateGlobalEquations`.

fem/src/ElementUtils.F90

@@ -3894,20 +3894,20 @@ END SUBROUTINE EvaluateVariableAtGivenPoint
 
      !> Return number of degrees of freedom and their indexes.
    !------------------------------------------------------------------------------
-   FUNCTION mGetElementDOFs( Indexes, UElement, USolver, NotDG, UMesh )  RESULT(nd)
+   FUNCTION mGetElementDOFs( Indexes, UElement, USolver, NotDG, UMesh, NotBubble )  RESULT(nd)
    !------------------------------------------------------------------------------
      INTEGER :: Indexes(:)
      TYPE(Element_t), OPTIONAL, TARGET :: UElement
      TYPE(Solver_t),  OPTIONAL, TARGET :: USolver
-     LOGICAL, OPTIONAL :: NotDG
+     LOGICAL, OPTIONAL :: NotDG, NotBubble
      TYPE(Mesh_t), OPTIONAL, TARGET :: UMesh
      INTEGER :: nd
 
      TYPE(Solver_t),  POINTER :: Solver
      TYPE(Element_t), POINTER :: Element, Parent, Face
      TYPE(Mesh_t), POINTER :: Mesh
 
-     LOGICAL :: Found, GB, DGDisable, NeedEdges
+     LOGICAL :: Found, GB, DGDisable, BubbleDisable, NeedEdges
      INTEGER :: i,j,k,id, nb, p, NDOFs, MaxNDOFs, EDOFs, MaxEDOFs, FDOFs, MaxFDOFs, BDOFs
      INTEGER :: Ind, ElemFamily, ParentFamily, face_type, face_id
      INTEGER :: NodalIndexOffset, EdgeIndexOffset, FaceIndexOffset
@@ -3992,7 +3992,7 @@ FUNCTION mGetElementDOFs( Indexes, UElement, USolver, NotDG, UMesh )  RESULT(nd)
 
      NDOFs = Solver % Def_Dofs(ElemFamily,id,1)
      IF (NDOFs > 0) THEN
-       DO i=1,Element % TYPE % NumberOfNodes
+       DO i=1,MIN(SIZE(Element % NodeIndexes),Element % TYPE % NumberOfNodes)
          DO j=1,NDOFs
            nd = nd + 1
            Indexes(nd) = MaxNDOFs * (Element % NodeIndexes(i)-1) + j
@@ -4259,7 +4259,11 @@ FUNCTION mGetElementDOFs( Indexes, UElement, USolver, NotDG, UMesh )  RESULT(nd)
          END IF
        END SELECT
      ELSE
-       IF (ASSOCIATED(Element % BubbleIndexes) .AND. Solver % GlobalBubbles) THEN
+       BubbleDisable=.FALSE.
+       IF (PRESENT(NotBubble)) BubbleDisable=NotBubble
+
+       IF (.NOT. BubbleDisable .AND. ASSOCIATED(Element % BubbleIndexes) &
+           .AND. Solver % GlobalBubbles) THEN
          BDOFs = 0
          nb = Solver % Def_Dofs(ElemFamily,id,5)
          p = Solver % Def_Dofs(ElemFamily,id,6)

fem/src/SolverUtils.F90

@@ -760,10 +760,10 @@ SUBROUTINE UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
      IF ( Rotate ) THEN
        NormalIndexes = 0
 
-       np = mGetElementDOFs(pIndexes,Element)
-       np = MIN(np, n)
+       np = mGetElementDOFs(pIndexes,Element,NotBubble=.TRUE.)
+       np = MIN(np,n)
        NormalIndexes(1:np) = NT % BoundaryReorder(pIndexes(1:np))
-       
+
        CALL RotateMatrix( LocalStiffMatrix, LocalForce, n, dim, NDOFs, &
           NormalIndexes, NT % BoundaryNormals, NT % BoundaryTangent1, NT % BoundaryTangent2 )
      END IF
@@ -850,18 +850,10 @@ SUBROUTINE UpdateGlobalEquationsVec( Gmtr, Lmtr, Gvec, Lvec, n, &
      IF ( Rotate ) THEN
        Ind = 0
 
-       np = mGetElementDOFs(pIndexes,Element)
+       np = mGetElementDOFs(pIndexes,Element,NotBubble=.TRUE.)
        np = MIN(np,n)
+       Ind(1:np) = NT % BoundaryReorder(pIndexes(1:np))
 
-       DO i=1,np
-         j = pIndexes(i)
-         IF(j > 0 .AND. j <= SIZE(NT % BoundaryReorder) ) THEN
-           Ind(i) = NT % BoundaryReorder(j)
-         ELSE
-           Ind(i) = j
-         END IF
-       END DO
-
        ! TODO: See that RotateMatrix is vectorized
        CALL RotateMatrix( Lmtr, Lvec, n, dim, NDOFs, Ind, NT % BoundaryNormals, &
                     NT % BoundaryTangent1, NT % BoundaryTangent2 )
@@ -14491,7 +14483,7 @@ END SUBROUTINE BlockSolveExt
           END IF
         END DO
       END IF
-      
+
       ! This is temporal (?) fix for a glitch where the complex eigen vector
       ! is expanded to one where real and complex parts follow each other. 
       IF( ListGetLogical( Solver % Values,'Expand Eigen Vectors', Stat ) ) THEN
@@ -19713,7 +19705,7 @@ SUBROUTINE ControlLinearSystem(Solver,PreSolve)
                 MAXVAL(f(i::dofs,iControl)),SUM(f(i::dofs,iControl))
           END DO
         END  IF
-       
+
         IF( ABS(cAmp(iControl)) > 1.0e-20 ) THEN
           b(1:nsize) = b(1:nsize) + cAmp(iControl) * f(1:nsize,iControl)
         END IF