!------------------------------------------------------------------------------
!> Simplified advection-diffusion solver with constant coefficients.
!------------------------------------------------------------------------------
SUBROUTINE AdvDiffSolver( Model,Solver,dt,TransientSimulation )
!------------------------------------------------------------------------------
  USE DefUtils

  IMPLICIT NONE
!------------------------------------------------------------------------------
  TYPE(Solver_t) :: Solver
  TYPE(Model_t) :: Model

  REAL(KIND=dp) :: dt
  LOGICAL :: TransientSimulation
!------------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------------
  LOGICAL :: AllocationsDone = .FALSE., Found
  TYPE(Element_t),POINTER :: Element

  REAL(KIND=dp) :: Norm
  INTEGER :: n, nb, nd, t, istat, NoActive
  TYPE(Matrix_t), POINTER :: A
  TYPE(Mesh_t), POINTER :: Mesh
  REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), MASS(:,:), LOAD(:), FORCE(:)
  SAVE STIFF, MASS, LOAD, FORCE, AllocationsDone
!------------------------------------------------------------------------------

  ! Allocate some permanent storage, this is done first time only:
  !--------------------------------------------------------------
  CALL Info( 'AdvDiffSolver','')
  CALL Info( 'AdvDiffSolver','----------------------------------------------------------')
  CALL Info( 'AdvDiffSolver','Solving advection-diffusion equation with material handles')

  Mesh => GetMesh()

  IF ( .NOT. AllocationsDone ) THEN
    N = Solver % Mesh % MaxElementDOFs  ! just big enough for elemental arrays
    ALLOCATE( FORCE(N), LOAD(N), STIFF(N,N), MASS(n,n), STAT=istat )
    IF ( istat /= 0 ) THEN
      CALL Fatal( 'AdvDiffSolver', 'Memory allocation error.' )
    END IF
    AllocationsDone = .TRUE.
   END IF
   
   ! System assembly:
   !-----------------
   CALL ResetTimer('AdvDiffSolver::Assembly')
   CALL DefaultInitialize()

   NoActive = GetNOFActive()
   DO t=1,NoActive

      Element => GetActiveElement(t)
      n  = GetElementNOFNodes()
      nd = GetElementNOFDOFs()
      nb = GetElementNOFBDOFs()

      ! Get element local matrix and rhs vector:
      ! Usually Elmer assembles stiffness and mass matrices and 
      ! adds them directly into global matrix. Here they are
      ! kept separately for later use.
      !------------------------------------------------------------
      CALL LocalMatrix( t==1,STIFF, MASS, FORCE, LOAD, Element, n, nd+nb )

      CALL DefaultUpdateEquations( STIFF, FORCE )
      IF( TransientSimulation ) THEN
        CALL DefaultUpdateMass( MASS )
      END IF
   END DO

   CALL DefaultFinishAssembly()
   CALL DefaultDirichletBCs()

   A => GetMatrix()

   ! Here the damping values are just added to the mass matrix
   ! to obtain the normal Galerkin method.    
   ! If FCT is implemented then something smarter should be done.
   !-------------------------------------------------------------
   IF (TransientSimulation) CALL Default1stOrderTimeGlobal()
   CALL CheckTimer('AdvDiffSolver::Assembly',Delete=.TRUE.)

   ! And finally, solve:
   !--------------------
   Norm = DefaultSolve()
   
   CALL Info( 'AdvDiffSolver','Add done')
   CALL Info( 'AdvDiffSolver','----------------------------------------------------------')


CONTAINS

!------------------------------------------------------------------------------
  SUBROUTINE LocalMatrix(  Initialize, STIFF, MASS, FORCE, LOAD, Element, n, nd )
!------------------------------------------------------------------------------
    LOGICAL :: Initialize 
    REAL(KIND=dp) :: STIFF(:,:), MASS(:,:), FORCE(:), LOAD(:)
    INTEGER :: n, nd
    TYPE(Element_t), POINTER :: Element
!------------------------------------------------------------------------------
    REAL(KIND=dp) :: Basis(nd),dBasisdx(nd,3),DetJ
    REAL(KIND=dp) :: Velo(3,nd), Velo_AtIp(3)
    REAL(KIND=dp) :: D_AtIp, C_AtIp, M_AtIp, S_AtIp
    LOGICAL :: Stat
    INTEGER :: i,t,p,q,dim
    TYPE(GaussIntegrationPoints_t) :: IP
    TYPE(ValueList_t), POINTER :: BodyForce, Material
    TYPE(Nodes_t) :: Nodes
    TYPE(ValueHandle_t) :: S_handle, D_handle, C_handle, M_handle

    SAVE Nodes, S_handle, D_handle, C_handle, M_handle

!------------------------------------------------------------------------------

    ! Initialize the handles for speedier execution.
    ! The names of the keywords are derived from the variable name.
    !--------------------------------------------------------------
    IF( Initialize ) THEN
      CALL ListInitRealAtIp(S_Handle,'Body Force',&
          TRIM(ComponentName(Solver % Variable))//' Source')      
      CALL ListInitRealAtIp(D_Handle,'Material',&
          TRIM(ComponentName(Solver % Variable))//' Diffusivity')
!      CALL ListInitRealAtIp(C_Handle,'Material',&
!          TRIM(ComponentName(Solver % Variable))//' Capacity')
!      CALL ListInitRealAtIp(M_Handle,'Material',&
!          TRIM(ComponentName(Solver % Variable))//' Density')      
    END IF


    dim = CoordinateSystemDimension()

    CALL GetElementNodes( Nodes )

    STIFF = 0._dp
    MASS = 0._dp
    FORCE = 0._dp
    LOAD = 0._dp

    Velo=0._dp; Velo_atip=0._dp
    CALL GetScalarLocalSolution(Velo(1,:),'Velocity 1')
    CALL GetScalarLocalSolution(Velo(2,:),'Velocity 2')
    CALL GetScalarLocalSolution(Velo(3,:),'Velocity 3')

    BodyForce => GetBodyForce()
    Material => GetMaterial()

    ! Numerical integration:
    !-----------------------
    IP = GaussPoints( Element )

    DO t=1,IP % n
      ! Basis function values & derivatives at the integration point:
      !--------------------------------------------------------------
      stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
              IP % W(t), detJ, Basis, dBasisdx )

      ! The source term at the integration point:
      !------------------------------------------
      S_AtIp = ListGetRealAtIp( S_Handle,BodyForce,Basis)

      ! Material parameters at integration point
      !--------------------------------------------------
      D_AtIp = ListGetRealAtIp( D_Handle,Material,Basis )
!      C_AtIp = ListGetRealAtIp( C_Handle,Material,Basis )
!      M_AtIp = ListGetRealAtIp( M_Handle,Material,Basis )
      C_AtIp = 1.0_dp
      M_ATIp = 1.0_dp

      ! Note that if the Handles would not be initialized
      ! an extra parameter for their name should be given as below
      !-----------------------------------------------------------
      ! M_AtIp = ListGetRealAtIp( M_Handle,Material,Basis, &
      ! TRIM(ComponentName(Solver % Variable))//' Density')

      ! Velocity at integration point
      !------------------------------
      DO i=1,dim
        Velo_AtIp(i) = SUM(Velo(i,1:n) * Basis(1:n))
      END DO

      ! Finally, the elemental matrix & vector:
      ! First the diffusive part
      !----------------------------------------
      STIFF(1:nd,1:nd) = STIFF(1:nd,1:nd) + IP % s(t) * DetJ * &
             D_AtIp * MATMUL( dBasisdx, TRANSPOSE( dBasisdx ) )

      ! Convective part
      DO p=1,nd
        DO q=1,nd
          STIFF(p,q) = STIFF(p,q) + IP % s(t) * DetJ * &
             C_AtIp * SUM(Velo_AtIp * dBasisdx(q,:)) * Basis(p)
        END DO
      END DO

      ! Mass matrix
      IF( TransientSimulation ) THEN
        DO p=1,nd
          DO q=1,nd
            MASS(p,q) = MASS(p,q) + IP % s(t) * DetJ * &
                M_AtIp * Basis(q) * Basis(p)
          END DO
        END DO
      END IF

      ! Source term
      FORCE(1:nd) = FORCE(1:nd) + IP % s(t) * DetJ * S_AtIP * Basis(1:nd)
    END DO
!------------------------------------------------------------------------------
  END SUBROUTINE LocalMatrix
!------------------------------------------------------------------------------

END SUBROUTINE AdvDiffSolver
