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

@@ -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.

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,:), &