From 8bf1df3423b209948333f1491c1e8d7d11f72574 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Markus=20M=C3=BCtzel?= <markus.muetzel@gmx.de>
Date: Sun, 4 Aug 2024 20:22:50 +0200
Subject: [PATCH] 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 | 14 +++++++++-----
 fem/src/SolverUtils.F90  | 22 +++++++---------------
 2 files changed, 16 insertions(+), 20 deletions(-)

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