
!! This routine is executed exactly once.  This occurs before the domain
!!   is partitioned into chunks.  
SUBROUTINE init

  IMPLICIT none
  !! Without femf.h, all the framework calls confuse the Fortran compiler !!  
  INCLUDE 'femf.h'

  !! Variable declarations !!
  INTEGER :: i, j, nelems, nnodes, ctype, esize
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: conn

  CALL FEM_Print('init called')
  OPEN(20, file='fmesh.dat')
  READ(20,*) nelems, nnodes, ctype

  !! Read in # of nodes/element
  IF (ctype .EQ. FEM_TRIANGULAR) THEN
    esize = 3
  ELSE
    IF(ctype .EQ. FEM_HEXAHEDRAL) THEN
      esize = 8
    ELSE
      esize = 4
    ENDIF
  ENDIF

  !! Read in connectivity info
  ALLOCATE(conn(nelems, esize))
  DO i=1,nelems
    READ(20,*) (conn(i,j),j=1,esize)
  ENDDO
  CLOSE(20)

  
  !! Give framework the connectivity info
  CALL FEM_Set_Mesh(nelems, nnodes, ctype, conn)


END SUBROUTINE init


!! This routine is executed for every chunk of the initial domain !!
!! The parameters contain the info for just those 
!!   elements/nodes in this chunk.
SUBROUTINE driver(numNodes, arrGlobNodeNums, numElems, &
                  arrGlobElemNums, numNodesPerElem, arrConnectivity)


!!--------------- Begin Variable Declaration ----------------!!

  IMPLICIT none

  TYPE Node
      DOUBLE PRECISION               :: prev
      DOUBLE PRECISION               :: temp
  END TYPE Node

  TYPE Element
      DOUBLE PRECISION               :: prev
      DOUBLE PRECISION               :: temp
  END TYPE Element

  TYPE BoundaryCondition
      INTEGER			     :: nodeNum
      DOUBLE PRECISION               :: bc
  END TYPE

  !! Parameters !! 
  INTEGER                        :: numNodes, numElems, numNodesPerElem 
  INTEGER, DIMENSION(numNodes)   :: arrGlobNodeNums
  INTEGER, DIMENSION(numElems)   :: arrGlobElemNums
  INTEGER, DIMENSION(numElems,numNodesPerElem)   :: arrConnectivity

  !! Local variables !!
  INTEGER                    :: i,j          !! loop indices
  INTEGER                    :: nodeNum      !! node number
  DOUBLE PRECISION           :: temperature  !! init temp
  INTEGER                    :: numGlobNodes !! global # of
  INTEGER                    :: numGlobBC    !! # of global boundary cond.
  INTEGER                    :: numBC        !! # of local boundary conditions
  INTEGER                    :: ignore       !! field of input to ignore
  DOUBLE PRECISION           :: bc           !! bc value 
  INTEGER                    :: fid          !! field id
  LOGICAL                    :: converged    !! is error small enough to quit?
  REAL                       :: errChunk     !! error value over chunk
  REAL			     :: errDomain    !! error value over entire domain 
  TYPE(BoundaryCondition), DIMENSION(:), ALLOCATABLE       :: arrBCs 
  TYPE(Node), DIMENSION(numNodes)                          :: arrNodes
  TYPE(Element), DIMENSION(numElems)                       :: arrElems

  !!---------------- End Variable Declaration -----------------!!

  !! Without this, all the framework calls confuse the Fortran compiler !!  
  INCLUDE 'femf.h'

  !! Read in intial node positions FOR THIS CHUNK ONLY !!
  OPEN(21, file='grid.dat')
  READ(21,*) 
  READ(21,*) numGlobNodes
  DO i = 1, numGlobNodes
    READ(21,*) nodeNum, temperature
    DO j = 1, numNodes                            !! only if nodenum
      IF (arrGlobNodeNums(j) == nodeNum) THEN     !!  belongs to this chunk
        arrNodes(j)%temp = temperature            !!  do we care about it.
        arrNodes(j)%prev = 0.d0
      ENDIF
    ENDDO
  ENDDO



  !! Read in boundary conditions FOR THIS CHUNK ONLY !!
  !! We put only the local bcs into an array of size globalBcCount
  !!   Efficiency is lacking, but the goal here is clarity.
  READ(21,*)
  READ(21,*) numGlobBC 
  ALLOCATE(arrBCs(numGlobBC))
  numBC = 0                                       !! local bc count
  DO i = 1, numGlobBC
    READ(21,*) ignore, nodeNum, bc
    DO j = 1, numNodes
      IF (arrGlobNodeNums(j) == nodeNum) THEN
        numBC = numBC + 1
        arrBCs(numBC)%nodeNum = nodeNum
        arrBCs(numBC)%bc = bc
      ENDIF
    ENDDO
  ENDDO 
 
  !! Set elements initially to zero !!
  DO i = 1, numElems
    arrElems(i)%temp = 0.d0
    arrElems(i)%prev = 0.d0
  ENDDO

  !! Register the node size and offset of temp field with Framework !!
  !! FEM_Create_Field( FEM_type, # of variables, 
  !!    offset of first var from beginning of Node TYPE, 
  !!    sizeof(Node TYPE) )
  fid = FEM_Create_Field(FEM_DOUBLE, 1, &
           offsetof(arrNodes(1), arrNodes(1)%temp), &
	   offsetof(arrNodes(1), arrNodes(2)))
  
  !! Combine shared nodes
  !! FEM_Update_Field(id returned from FEM_Create_Field, address of array)
  CALL FEM_Update_Field(fid, arrNodes(1))

  converged = .false.
  DO WHILE (.NOT.converged)

    !! Apply the boundary conditions
    DO i = 1, numBC
      arrNodes(arrBCs(i)%nodeNum)%temp = arrBCs(i)%bc
    ENDDO

    !! Do the calculations on the elements: avg of surrounding neighbors !!
    DO i = 1, numElems
      arrElems(i)%prev = arrElems(i)%temp
      arrElems(i)%temp = 0.d0
      DO j = 1, numNodesPerElem
        arrElems(i)%temp = arrElems(i)%temp + &
             arrNodes(arrConnectivity(i,j))%temp
      ENDDO
      arrElems(i)%temp = arrElems(i)%temp / numNodesPerElem
    ENDDO

    !! Prepare the nodes for calculation !!
    DO i = 1, numNodes
      arrNodes(i)%prev = arrNodes(i)%temp
      arrNodes(i)%temp = 0.d0
    ENDDO

    !! Do calculation on nodes !!
    DO i = 1, numElems
      DO j = 1, numNodesPerElem
        arrNodes(arrConnectivity(i,j))%temp =  &
             arrNodes(arrConnectivity(i,j))%temp + &
             arrElems(i)%temp / numNodesPerElem
      ENDDO
    ENDDO

    !! Combine shared nodes !!
    !! FEM_Update_Field(id returned from FEM_Create_Field, address of array)
    CALL FEM_Update_Field(fid, arrNodes(1)) 


    !! Calculate the error on this chunk
    errChunk = 0.d0
    DO i = 1, numElems
      errChunk = errChunk + ABS(arrElems(i)%prev - arrElems(i)%temp)
    ENDDO

    !! Now, sum the errors from all the elements in all the chunks
    !! FEM_Reduce(FEM_type, contrib from this chunk, variable to store result,
    !!            FEM_operation) 
    !! Write(*,'(A,F15.13)') 'err for chunk = ', errChunk
    CALL FEM_Reduce(FEM_DOUBLE, errChunk, errDomain, &
                    FEM_Max)

    !! If we are on the first processor, print error achieved
    !! IF (FEM_My_Partition() == 0) WRITE(*,'(A,F12.10)') 'Error = ',  errDomain
    converged = (errDomain <= 1*10**(-8))
  ENDDO
END SUBROUTINE driver


!! This routine is executed exactly once.  This occurs after the 
!!   calculations on each chunk are complete.
SUBROUTINE finalize()
  INCLUDE 'femf.h'
  CALL FEM_Print('finalize called')
END SUBROUTINE finalize

