diff --git a/fem/src/ElementUtils.F90 b/fem/src/ElementUtils.F90
index dcc8e5374f..fc44cfcd2f 100644
--- a/fem/src/ElementUtils.F90
+++ b/fem/src/ElementUtils.F90
@@ -3894,12 +3894,12 @@ 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
 
@@ -3907,7 +3907,7 @@ FUNCTION mGetElementDOFs( Indexes, UElement, USolver, NotDG, UMesh )  RESULT(nd)
      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)
diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90
index 6828d09339..8d6d0b1de1 100644
--- a/fem/src/SolverUtils.F90
+++ b/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