Merge commit 'af9ce3271a6aa309a840f6eaba6fb310fbddb788' of https://gi…

…thub.com/elmercsc/elmerfem into devel

Conflicts:
	fem/tests/Zirka/ZirkaTest.F90

fem/src/DefUtils.F90

@@ -3581,9 +3581,13 @@ RECURSIVE FUNCTION DefaultSolve( USolver, BackRotNT ) RESULT(Norm)
 
     CALL Info('DefaultSolve','Calling SolveSystem for linear solution',Level=20)
 
+print*,'a: ', associated(solver), associated(Solver % matrix); flush(6)
     A => Solver % Matrix
+print*,'a: ', associated(solver % variable); flush(6)
     x => Solver % Variable    
+print*,'a: ', associated(A), associated(A % RHS); flush(6)
     b => A % RHS
+print*,'a: ', associated(x % values), SIZE(x % values); flush(6)
     SOL => x % Values
 
     ! Debugging stuff activated only when "Max Output Level" >= 20
@@ -3594,6 +3598,7 @@ RECURSIVE FUNCTION DefaultSolve( USolver, BackRotNT ) RESULT(Norm)
 
 10  CONTINUE
 
+print*,'going sol'; flush(6)
     CALL SolveSystem(A,ParMatrix,b,SOL,x % Norm,x % DOFs,Solver)
 
     IF( InfoActive( 20 ) ) THEN

fem/src/ElementUtils.F90

@@ -1904,6 +1904,7 @@ FUNCTION CreateMatrix( Model, Solver, Mesh, Perm, DOFs, MatrixFormat, &
      IF ( PRESENT(Equation) ) THEN
        CALL Info(Caller,'Creating initial permutation',Level=14)
        k = InitialPermutation( Perm,Model,Solver,Mesh,Eq,DG,GB )
+print*,'aft k', k
        IF ( k <= 0 ) THEN
          IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder)
          RETURN

fem/src/ListMatrix.F90

@@ -482,6 +482,113 @@ SUBROUTINE List_AddMatrixIndexes(List,k1,nk2,Ind)
    END SUBROUTINE List_AddMatrixIndexes
 !-------------------------------------------------------------------------------
 
+
+!-------------------------------------------------------------------------------
+   SUBROUTINE List_AddMatrixRow(List,k1,nk2,Ind,Vals)
+   ! Add an array of sorted indeces to a row in ListMatrix_t. "ind" may
+   ! contain duplicate entries.
+!-------------------------------------------------------------------------------
+     IMPLICIT NONE
+
+     TYPE(ListMatrix_t) :: List(:)
+     INTEGER, INTENT(IN) :: k1, nk2
+     INTEGER, INTENT(INOUT) :: Ind(nk2)
+     REAL(KIND=dp), INTENT(INOUT) :: Vals(nk2)
+
+     TYPE(ListMatrixEntry_t), POINTER :: RowPtr, PrevPtr, Entry, Dummy
+!-------------------------------------------------------------------------------
+     INTEGER :: i,k2,k2i,j, k,prevind
+
+     IF (k1>SIZE(List)) THEN
+       CALL Fatal('List_AddMatrixIndexes','Row index out of bounds: '//TRIM(I2S(k1)))
+     END IF
+
+     CALL SortF(nk2, Ind, Vals)
+
+     ! Add each element in Ind to the row list
+     RowPtr => List(k1) % Head
+
+     ! First element needs special treatment as it may modify 
+     ! the list starting point
+     IF (.NOT. ASSOCIATED(RowPtr)) THEN
+       Dummy => NULL() 
+       Entry => List_GetMatrixEntry(Ind(1),Dummy)
+       Entry % Val = Vals(1)
+       List(k1) % Degree = 1
+       List(k1) % Head => Entry
+       k2i = 2
+       prevind = ind(1)
+     ELSE IF (RowPtr % Index > Ind(1)) THEN
+       Entry => List_GetMatrixEntry(Ind(1),RowPtr)
+       Entry % Val = Vals(1)
+       List(k1) % Degree = List(k1) % Degree + 1
+       List(k1) % Head => Entry
+       k2i = 2
+       prevind = ind(1)
+     ELSE IF (RowPtr % Index == Ind(1)) THEN
+        k2i = 2
+        prevind = ind(1)
+     ELSE
+       k2i = 1
+       prevind = -1
+     END IF
+
+     PrevPtr => List(k1) % Head
+     RowPtr  => List(k1) % Head % Next
+
+     DO i=k2i,nk2
+       k2 = Ind(i)
+       if (k2 == prevind) cycle
+
+       ! Find a correct place place to add index to
+       DO WHILE( ASSOCIATED(RowPtr) )
+         IF (RowPtr % Index >= k2) EXIT
+         PrevPtr => RowPtr
+         RowPtr  => RowPtr % Next
+       END DO
+
+       IF (ASSOCIATED(RowPtr)) THEN
+         ! Do not add duplicates
+         IF (RowPtr % Index /= k2) THEN
+           ! Create new element between PrevPtr and RowPtr
+           Entry => List_GetMatrixEntry(k2,RowPtr)
+           Entry % Val = Vals(i)
+           PrevPtr % Next => Entry
+           List(k1) % Degree = List(k1) % Degree + 1
+
+           ! Advance to next element in list
+           PrevPtr => Entry
+         ELSE
+           ! Advance to next element in list
+           RowPtr % Val = RowPtr % Val + Vals(i)
+           PrevPtr => RowPtr
+           RowPtr  => RowPtr % Next
+         END IF
+       ELSE
+         EXIT
+       END IF
+
+       prevind = k2
+     END DO
+
+     DO j=i,nk2
+       k2 = Ind(j)
+       if (k2 == prevind) cycle
+       prevind = k2
+
+       Dummy => NULL()
+       Entry => List_GetMatrixEntry(k2,Dummy)
+       Entry % Val = Vals(j)
+       PrevPtr % Next => Entry
+       PrevPtr => PrevPtr % Next
+       List(k1) % Degree = List(k1) % Degree + 1
+     END DO
+!-------------------------------------------------------------------------------
+   END SUBROUTINE List_AddMatrixRow
+!-------------------------------------------------------------------------------
+
+
+
 !-------------------------------------------------------------------------------
    FUNCTION List_GetMatrixEntry(ind, next) RESULT(ListEntry)
 !-------------------------------------------------------------------------------

fem/src/Lists.F90

@@ -117,13 +117,25 @@ SUBROUTINE SetGetMatcParams(nparams,params,resul)
      REAL(KIND=dp) :: params(:)
      CHARACTER(*), OPTIONAL :: resul
 
-     INTEGER :: i,l
+     INTEGER :: i,j,l
      CHARACTER(LEN=1024) :: pcmd,res
 
      IF(nparams==0) THEN
        pcmd = "tx=0"
      ELSE
+#if 0
        WRITE(pcmd,*)  [(params(i),i=1,nparams)]
+#else
+       ! cray ftn output from above can be somewhat convoluted, do this instead
+       j = 1
+       DO i=1,nparams
+         WRITE(pcmd(j:), *) params(i)
+         DO WHILE(pcmd(j:j) == ' '); j=j+1; END DO
+         DO WHILE(pcmd(j:j) /= ' '); j=j+1; END DO
+         IF(pcmd(j-1:j-1)=='.') pcmd(j-1:j-1) = ' '
+         j = j + 1
+       END DO
+#endif
        IF(PRESENT(resul)) THEN
          pcmd = TRIM(resul)//'='//TRIM(pcmd)
        ELSE
@@ -243,6 +255,7 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, &
          INTEGER :: body_id, MaxGroup, group0, group
          INTEGER, POINTER :: DgMap(:), DgMaster(:), DgSlave(:)
          LOGICAL :: GotDgMap, GotMaster, GotSlave
+!------------------------------------------------------------------------------
 
          DgMap => ListGetIntegerArray( Solver % Values,'DG Reduced Basis Mapping',GotDgMap )
          DgMaster => ListGetIntegerArray( Solver % Values,'DG Reduced Basis Master Bodies',GotMaster )
@@ -326,7 +339,8 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, &
              ' db nodes from bulk hits',Level=15)
 
          IF ( FoundDG ) THEN
-           RETURN ! Discontinuous bodies !!!
+           GOTO 10
+!          RETURN ! Discontinuous bodies !!!
          END IF
        END BLOCK
      END IF
@@ -408,7 +422,8 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, &
            ' nodes from bulk hits',Level=15)
 
        IF ( FoundDG ) THEN
-          RETURN ! Discontinuous galerkin !!!
+          GOTO 10
+!         RETURN ! Discontinuous galerkin !!!
        END IF
      END IF
 
@@ -678,6 +693,9 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, &
 
      IF ( ALLOCATED(EdgeDOFs) ) DEALLOCATE(EdgeDOFs)
      IF ( ALLOCATED(FaceDOFs) ) DEALLOCATE(FaceDOFs)
+
+10   CONTINUE
+
 !------------------------------------------------------------------------------
    END FUNCTION InitialPermutation
 !------------------------------------------------------------------------------

fem/src/LoadMod.F90

@@ -297,7 +297,7 @@ FUNCTION loadfunction(quiet, abort_not_found, library, fname) RESULT(ptr)
             CHARACTER :: library(*), fname(*)
             ! TYPE(C_FUNPTR) :: ptr
             INTEGER(KIND=AddrInt) :: ptr
-            TYPE(C_FUNPTR) :: cptr
+            TYPE(C_PTR) :: cptr
 
             INTERFACE
                 FUNCTION loadfunction_c(quiet, abort_not_found, library, fname ) RESULT(cptr) &
@@ -307,7 +307,7 @@ FUNCTION loadfunction_c(quiet, abort_not_found, library, fname ) RESULT(cptr) &
                     INTEGER(C_INT) :: abort_not_found
                     CHARACTER(C_CHAR) :: library(*), fname(*)
                     ! INTEGER(CAddrInt) :: fptr
-                    TYPE(C_FUNPTR) :: cptr
+                    TYPE(C_PTR) :: cptr
                 END FUNCTION loadfunction_c
             END INTERFACE
 

fem/src/MeshUtils.F90

@@ -27468,7 +27468,7 @@ SUBROUTINE ElmerGraphColour(Graph, Colouring, ConsistentColours)
 
     INTEGER :: nc, dualmaxdeg, i, v, w, uci, wci, vli, vti, vcol, wcol, &
             nrc, nunc, nthr, TID, allocstat, gn
-    INTEGER, ALLOCATABLE :: colours(:)
+    INTEGER, POINTER :: colours(:)
     INTEGER, PARAMETER :: VERTEX_PER_THREAD = 100
     LOGICAL :: consistent
 
@@ -27649,7 +27649,8 @@ SUBROUTINE ElmerGraphColour(Graph, Colouring, ConsistentColours)
 
     ! Set up colouring data structure
     Colouring % nc = nc
-    CALL MOVE_ALLOC(colours, Colouring % colours)
+!   CALL MOVE_ALLOC(colours, Colouring % colours)
+    Colouring % colours => colours
   END SUBROUTINE ElmerGraphColour
 
   SUBROUTINE Colouring_Deallocate(Colours)
@@ -27709,7 +27710,7 @@ SUBROUTINE ElmerBoundaryGraphColour(Mesh, Colours, BoundaryColours)
 
     TYPE(Element_t), POINTER :: Element
     INTEGER :: elem, nelem, nbelem, astat, lcolour, rcolour, nbc
-    INTEGER, ALLOCATABLE :: bcolours(:)
+    INTEGER, POINTER :: bcolours(:)
 
     nelem = Mesh % NumberOfBulkElements
     nbelem = Mesh % NumberOfBoundaryElements
@@ -27757,7 +27758,8 @@ SUBROUTINE ElmerBoundaryGraphColour(Mesh, Colours, BoundaryColours)
 
     ! Set up colouring data structure
     BoundaryColours % nc = nbc
-    CALL MOVE_ALLOC(bcolours, BoundaryColours % colours)
+!   CALL MOVE_ALLOC(bcolours, BoundaryColours % colours)
+    BoundaryColours % colours => bcolours
   END SUBROUTINE ElmerBoundaryGraphColour
 
   ! Given CRS indices, referenced indirectly from graph, 

fem/src/SolverUtils.F90

@@ -7365,16 +7365,16 @@ END SUBROUTINE SetDirichletPoint
 !------------------------------------------------------------------------------
   SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec )
 !------------------------------------------------------------------------------
-    TYPE(Model_t), POINTER :: Model  !< The current model structure
-    TYPE(Mesh_t), POINTER :: Mesh    !< The current mesh structure
+    TYPE(Model_t) :: Model  !< The current model structure
+    TYPE(Mesh_t)  :: Mesh    !< The current mesh structure
     CHARACTER(LEN=*) :: SourceName   !< Name of the keyword setting the source term
     INTEGER :: DOFs                  !< The total number of DOFs for this equation
     INTEGER :: Perm(:)               !< The node reordering info
     LOGICAL :: GotSrc                !< Did we get something?
     REAL(KIND=dp) :: SrcVec(:)       !< The assemblied source vector
 !------------------------------------------------------------------------------
     TYPE(Element_t), POINTER :: Element
-    INTEGER :: i,t,n,bc,bf,FirstElem,LastElem,nlen
+    INTEGER :: i,j,k,t,n,bc,bf,FirstElem,LastElem,nlen
     LOGICAL :: Found,AnyBC,AnyBF,Axisymmetric
     REAL(KIND=dp) :: Coeff
     REAL(KIND=dp), ALLOCATABLE :: FORCE(:,:)
@@ -7389,8 +7389,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec
     CALL Info(Caller,'Checking for generalized source terms: '&
         //SourceName(1:nlen),Level=15)
 
-    ALLOCATE( ActiveBC(Model % NumberOfBCs ), &
-        ActiveBF(Model % NumberOfBodyForces) )
+    ALLOCATE( ActiveBC(Model % NumberOfBCs), ActiveBF(Model % NumberOfBodyForces) )
 
     ! First make a quick test going through the short boundary condition and
     ! body force lists.
@@ -7439,6 +7438,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec
     DO t=FirstElem, LastElem 
       Element => Mesh % Elements(t)
       Indexes => Element % NodeIndexes
+      n = Element % Type % NumberOfNodes
 
       IF( t > Mesh % NumberOfBulkElements ) THEN
         Found = .FALSE.
@@ -7465,8 +7465,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec
       CALL LocalSourceAssembly(Element, dofs, FORCE )
 
       DO i=1,dofs
-        SrcVec(dofs*(Perm(Indexes)-1)+i) = SrcVec(dofs*(Perm(Indexes)-1)+i) + &
-            Coeff * FORCE(i,1:n)
+        SrcVec(dofs*(Perm(Indexes)-1)+i) = SrcVec(dofs*(Perm(Indexes)-1)+i) + FORCE(i,1:n)
       END DO
     END DO
 
@@ -16187,12 +16186,11 @@ SUBROUTINE FinalizeLumpedMatrix( Solver )
           END DO
           CLOSE(11)
         END IF
-
+
+        WRITE(Message,*) 'Normalization checksum: ',CheckSum
+        CALL Info(Caller,Message) 
       END IF
       CALL Info( Caller,'Constraint modes fluxes was saved to file '//TRIM(MatrixFile),Level=5)
-
-      WRITE(Message,*) 'Normalization checksum: ',CheckSum
-      CALL Info(Caller,Message) 
     END IF
   END IF
 
@@ -19700,21 +19698,20 @@ SUBROUTINE ControlLinearSystem(Solver,PreSolve)
         ! We need to add the control source here in order to be able to use
         ! standard means for convergence monitoring. 
         CALL Info(Caller,'Computing source term for: '//TRIM(str),Level=7)
-        CALL SetNodalSources( CurrentModel,Mesh,str, &
-            dofs, Perm, GotF, f(:,iControl) )
+        CALL SetNodalSources( CurrentModel,Mesh,str,dofs, Perm, GotF, f(:,iControl) )
 
-       ! The additional source needs to be nullified for Dirichlet conditions
-       IF( ALLOCATED( A % ConstrainedDOF ) ) THEN
-         WHERE( A % ConstrainedDOF ) f(:,iControl) = 0.0_dp
-       END IF
+        ! The additional source needs to be nullified for Dirichlet conditions
+        IF( ALLOCATED( A % ConstrainedDOF ) ) THEN
+          WHERE( A % ConstrainedDOF ) f(:,iControl) = 0.0_dp
+        END IF
 
-       IF(InfoActive(10)) THEN
-         DO i=1,dofs
-           PRINT *,'ranges b:',i,MINVAL(b(i::dofs)),MAXVAL(b(i::dofs)),SUM(b(i::dofs))
-           PRINT *,'ranges f:',i,MINVAL(f(i::dofs,iControl)),&
-               MAXVAL(f(i::dofs,iControl)),SUM(f(i::dofs,iControl))
-         END DO
-       END  IF
+        IF(InfoActive(10)) THEN
+          DO i=1,dofs
+            PRINT *,'ranges b:',i,MINVAL(b(i::dofs)),MAXVAL(b(i::dofs)),SUM(b(i::dofs))
+            PRINT *,'ranges f:',i,MINVAL(f(i::dofs,iControl)),&
+                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)
@@ -25138,6 +25135,17 @@ END SUBROUTINE P2LagrangeSwapper
    !-------------------------------------------------------------------------------------
    FUNCTION EvalFieldAtElem( Mesh, Var, Element, Basis, dofi, eigeni, imVal, GotVal ) RESULT ( Val )
 
+     INTERFACE 
+       SUBROUTINE Ip2DgFieldInElement( Mesh, Parent, nip, fip, np, fdg )
+         USE Types
+         IMPLICIT NONE
+         TYPE(Mesh_t), POINTER :: Mesh
+         TYPE(Element_t), POINTER :: Parent
+         INTEGER :: nip, np
+         REAL(KIND=dp) :: fip(:), fdg(:)
+       END SUBROUTINE Ip2DgFieldInElement
+     END INTERFACE
+
      TYPE(Mesh_t), POINTER :: Mesh
      TYPE(Variable_t), POINTER :: Var
      TYPE(Element_t), POINTER :: Element

fem/src/Types.F90

@@ -854,7 +854,7 @@ MODULE Types
 
    TYPE Graphcolour_t
      INTEGER :: nc
-     INTEGER, ALLOCATABLE :: colours(:)
+     INTEGER, POINTER :: colours(:) => Null()
    END TYPE Graphcolour_t
 
    TYPE MortarBC_t 

fem/src/modules/SteadyPhaseChange.F90

@@ -1,4 +1,4 @@
-/*****************************************************************************/
+!/*****************************************************************************/
 ! *
 ! *  Elmer, A Finite Element Software for Multiphysical Problems
 ! *