!>
!> @file pde1dp_mod.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE pde1dp_mod
!
  USE bsplines
  USE matrix
  IMPLICIT NONE
  DOUBLE PRECISION, ALLOCATABLE :: bcoef(:)
  TYPE(spline1d), SAVE :: splx
!  
CONTAINS 
  SUBROUTINE meshdist(c, x, nx)
!
!   Construct an 1d  non-equidistant mesh given a
!   mesh distribution function defined in FDIST
!
    DOUBLE PRECISION, INTENT(in) :: c(5)
    INTEGER, INTENT(iN) :: nx
    DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
    INTEGER :: nintg
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
    DOUBLE PRECISION :: a, b, dx, f0, f1, scal
    INTEGER :: i, k
!
    a=x(0)
    b=x(nx)
    nintg = 10*nx
    ALLOCATE(xint(0:nintg), fint(0:nintg))
!
!  Mesh distribution
!
    dx = (b-a)/REAL(nintg)
    xint(0) = a
    fint(0) = 0.0d0
    f1 = fdist(xint(0))
    DO i=1,nintg
       f0 = f1
       xint(i) = xint(i-1) + dx
       f1 = fdist(xint(i))
       fint(i) = fint(i-1) + 0.5*(f0+f1)
    END DO
!
!  Normalization
!
    scal = REAL(nx) / fint(nintg)
    fint(0:nintg) = fint(0:nintg) * scal
!
!  Obtain mesh point by (inverse) interpolation
!
    k = 1
    DO i=1,nintg-1
       IF( fint(i) .GE. REAL(k) ) THEN
          x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
               &   (k-fint(i))
          k = k+1
       END IF
    END DO
!
    DEALLOCATE(xint, fint)
  CONTAINS
    DOUBLE PRECISION FUNCTION fdist(x)
      DOUBLE PRECISION, INTENT(in) :: x
      fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
    END FUNCTION fdist
  END SUBROUTINE meshdist
!+++
  SUBROUTINE dismat(spl, mat)
!
!   Assembly FE matrix (with periodic BC) mat using spline spl
!
    TYPE(spline1d), INTENT(in) :: spl
    TYPE(periodic_mat), INTENT(inout) :: mat
    INTEGER :: dim, nx, nidbas, ngauss
    INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
    DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
    DOUBLE PRECISION :: contrib
!
    INTEGER :: kterms                ! Number of terms in weak form
    INTEGER, ALLOCATABLE :: idert(:), iderw(:)  ! Derivative order
    DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
    CALL get_dim(spl, dim, nx, nidbas)
    ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
!   Weak form
!
    kterms = mat%mat%nterms
    ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
!   Gauss quadature
!
    CALL get_gauss(spl, ngauss)
    ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
    DO i=1,nx
       CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
       DO igauss=1,ngauss
          CALL basfun(xgauss(igauss), spl, fun, i)
          CALL coefeq(xgauss(igauss), idert, iderw, coefs)
          DO iterm=1,kterms
             DO jt=0,nidbas
                DO iw=0,nidbas
                   contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
                        &    fun(iw,iderw(iterm)) * wgauss(igauss)
                   irow=MODULO(i+iw-1,nx) + 1   ! Periodic BC
                   jcol=MODULO(i+jt-1,nx) + 1
                   CALL updtmat(mat, irow, jcol, contrib)
                END DO
             END DO
          END DO
       END DO
    END DO
!===========================================================================
!              9.0  Epilogue
!
    DEALLOCATE(fun)
    DEALLOCATE(xgauss, wgauss)
    DEALLOCATE(iderw, idert, coefs)
  END SUBROUTINE dismat
!
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Mass matrix
!
    c(1) = 1.0d0
    idt(1) = 0
    idw(1) = 0
  END SUBROUTINE coefeq
!+++
  SUBROUTINE disrhs(spl, rhs)
!
!   Assenbly the RHS using spline spl
!
    TYPE(spline1d), INTENT(in) :: spl
    DOUBLE PRECISION, INTENT(out) :: rhs(:)
    INTEGER :: dim, nrank, nx, nidbas, ngauss
    INTEGER :: i, igauss, it, irow
    DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
    DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
    DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
    nrank = SIZE(rhs)
    CALL get_dim(spl, dim, nx, nidbas)
!
    ALLOCATE(fun(0:nidbas,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
    CALL get_gauss(spl, ngauss)
    ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
    rhs(:) = 0.0d0
!
    DO i=1,nx
       CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
       DO igauss=1,ngauss
          CALL basfun(xgauss(igauss), spl, fun, i)
          contrib = wgauss(igauss) * rhseq(xgauss(igauss))
          DO it=0,nidbas
             irow=MODULO(i+it-1,nx) + 1   ! Periodic BC
             rhs(irow) = rhs(irow) + contrib*fun(it,1)
          END DO
       END DO
    END DO
!===========================================================================
!              9.0  Epilogue
!
    DEALLOCATE(fun)
    DEALLOCATE(xgauss, wgauss)
  END SUBROUTINE disrhs
!
DOUBLE PRECISION FUNCTION rhseq(x)
  DOUBLE PRECISION, INTENT(in) :: x
  DOUBLE PRECISION :: xarr(1), farr(1)
  INTEGER, SAVE :: icall =0
  xarr(1) = x
  IF( icall.EQ.0 ) THEN
     icall = icall+1
     CALL gridval(splx, xarr, farr, 0, bcoef)
  ELSE
     CALL gridval(splx, xarr, farr, 0)
  END IF
  rhseq = farr(1)
END FUNCTION rhseq

END MODULE pde1dp_mod