From f6c727dc3f680a16273560da8b0904358dae1625 Mon Sep 17 00:00:00 2001
From: Peter R <peter.raback@csc.fi>
Date: Thu, 14 Nov 2024 14:38:09 +0200
Subject: [PATCH] Fix for retriving data from a given node when there are
 p-degrees of freedom available. This does not affect nodal quantities and
 hence local coordinates are not needed. This affected SaveLine when using
 p-bubbles.

---
 fem/src/ElementUtils.F90              |  9 +++++++--
 fem/src/modules/SaveData/SaveLine.F90 | 18 ++++++++++++++++++
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/fem/src/ElementUtils.F90 b/fem/src/ElementUtils.F90
index 4bd46e9140..d1b4474b61 100644
--- a/fem/src/ElementUtils.F90
+++ b/fem/src/ElementUtils.F90
@@ -3519,9 +3519,14 @@ END SUBROUTINE Ip2DgFieldInElement
     DGVar = ( Var % TYPE == variable_on_nodes_on_elements ) 
     IpVar = ( Var % TYPE == variable_on_gauss_points )
     ElemVar = ( Var % TYPE == Variable_on_elements )       
+
+    ! We do not need P-elements if the value is to be found in node since
+    ! the higher order p-basis does not have any effect there.
     pElem = .FALSE.
-    IF(ASSOCIATED(Var % Solver)) THEN
-      pElem = isActivePElement(Element, Var % Solver)                              
+    IF(.NOT. PRESENT(LocalNode)) THEN
+      IF(ASSOCIATED(Var % Solver)) THEN
+        pElem = isActivePElement(Element, Var % Solver) 
+      END IF
     END IF
       
     PiolaVersion = .FALSE.
diff --git a/fem/src/modules/SaveData/SaveLine.F90 b/fem/src/modules/SaveData/SaveLine.F90
index 37a00804de..7ab0256a59 100644
--- a/fem/src/modules/SaveData/SaveLine.F90
+++ b/fem/src/modules/SaveData/SaveLine.F90
@@ -1336,6 +1336,7 @@ SUBROUTINE SaveExistingLines()
     TYPE(Element_t), POINTER :: Parent
     LOGICAL :: BreakLoop, ParallelComm
     REAL(KIND=dp) :: linepos
+    INTEGER, ALLOCATABLE :: NodeToElement(:)
     
     MaskName = ListGetString(Params,'Save Mask',GotIt) 
     IF(.NOT. GotIt) MaskName = 'Save Line'
@@ -1401,6 +1402,19 @@ SUBROUTINE SaveExistingLines()
         END IF
       END DO
       
+      ! Create a table where from each node we have something pointing to an element.
+      ALLOCATE(NodeToElement(Mesh % NumberOfNodes))
+      NodeToElement = 0
+      DO t = 1,  Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements                
+        CurrentElement => Mesh % Elements(t)
+        IF( ParEnv % PEs > 1 ) THEN
+          IF( CurrentElement % PartIndex /= ParEnv % MyPe ) CYCLE
+        END IF
+        NodeIndexes => CurrentElement % NodeIndexes        
+        NodeToElement(NodeIndexes) = CurrentElement % ElementIndex
+      END DO
+
+      
       IF(CalculateFlux) THEN
         CALL Info(Caller,'Calculating nodal fluxes',Level=8)
         ALLOCATE(PointFluxes(SaveNodes(1),3),PointWeight(SaveNodes(1)), STAT=istat)    
@@ -1565,6 +1579,10 @@ SUBROUTINE SaveExistingLines()
         linepos = -1.0_dp
         DO t = 1, SaveNodes(1)    
           node = InvPerm(t)
+          
+          ! Get some element which may be usefull in evaluating the field.
+          CurrentElement => Mesh % Elements(NodeToElement(node))
+          
           IF( CalculateFlux ) THEN
             CALL WriteFieldsAtElement( CurrentElement, BoundaryIndex(t), node, &
                 dgnode, UseNode = .TRUE., NodalFlux = PointFluxes(t,:), &