Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/zconmat_1d.tpl
f1da25c60953feature/refactoring

Last Change
Raw File

zconmat_1d.tpl
View Options

!>
!> @file zconmat_1d.tpl
!>
!> @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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
!
! In this version s[lines are precalculted
! (on all n1/n2 intervals
!
INTERFACE
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE COMPLEX, INTENT(out) :: c(:)
END SUBROUTINE coefeq
END INTERFACE
INTEGER, OPTIONAL :: maxder ! maximum oder of derivatives
!
INTEGER :: n1, nidbas1, ndim1, ng1
INTEGER :: i1, ig1
INTEGER :: irow, jcol
INTEGER, ALLOCATABLE :: left1(:)
!
LOGICAL :: nlper1
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER :: k, kmaxder, it1, iw1
INTEGER, ALLOCATABLE :: idert(:), iderw(:) ! Derivative order
DOUBLE COMPLEX :: one=(1.0d0,0.0d0), zero=(0.0d0,0.0d0)
DOUBLE COMPLEX, ALLOCATABLE :: coefs(:,:) ! Terms in weak form
!
DOUBLE PRECISION, POINTER :: xg1(:,:), wg1(:,:)
DOUBLE PRECISION, ALLOCATABLE :: fun1(:,:,:,:), fun2(:,:,:,:)
DOUBLE COMPLEX, ALLOCATABLE :: mata(:,:), matc(:)
DOUBLE COMPLEX, ALLOCATABLE :: matf(:,:,:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, ndim1, n1, nidbas1)
nlper1 = spl%period
!
! Gauss points and weights on all intervals
!
xg1 => spl%gausx ! xg1(ng1,n1)
wg1 => spl%gausw ! wg1(ng1,n1)
ng1 = SIZE(xg1,1)
!
! Splines on all intervals
!
kmaxder = 1
IF(PRESENT(maxder)) kmaxder = maxder
ALLOCATE(fun1(0:nidbas1,0:kmaxder,ng1,n1))
ALLOCATE(left1(ng1))
DO i1=1,n1
left1 = i1
CALL basfun(xg1(:,i1), spl, fun1(:,:,:,i1), left1)
END DO
DEALLOCATE(left1)
!===========================================================================
! 2.0 Assembly loop
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms))
ALLOCATE(iderw(kterms))
ALLOCATE(coefs(kterms,ng1))
!
! Allocate local matrices
!
ALLOCATE(mata(0:nidbas1,0:nidbas1))
ALLOCATE(matc(ng1))
ALLOCATE(matf(0:nidbas1,0:nidbas1,ng1))
!
DO i1=1,n1
!
! Coefficients of the weak form
!
DO ig1=1,ng1
CALL coefeq(xg1(ig1,i1), idert, iderw, coefs(:,ig1))
END DO
!
! Compute local matrix: A <- F*c + A
!
mata = 0.0d0
DO k=1,kterms
!
matc(1:ng1) = coefs(k,1:ng1)
!
DO it1=0,nidbas1
DO iw1=0,nidbas1
DO ig1=1,ng1
matf(it1,iw1,ig1) = wg1(ig1,i1) * &
& fun1(it1,idert(k),ig1,i1) * &
& fun1(iw1,iderw(k),ig1,i1)
END DO
END DO
END DO
!
CALL zgemv('N', (nidbas1+1)*(nidbas1+1), ng1, one, matf, &
& (nidbas1+1)*(nidbas1+1), matc, 1, one, mata, 1)
END DO
!
! Map local matrix A to global matrix
!
DO it1=0,nidbas1
irow = i1+it1; IF(nlper1) irow = MODULO(irow-1,n1) + 1
DO iw1=0,nidbas1
jcol = i1+iw1; IF(nlper1) jcol = MODULO(jcol-1,n1) + 1
CALL updtmat(mat, irow, jcol, mata(it1,iw1))
END DO
END DO
!
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun1)
DEALLOCATE(idert, iderw, coefs)
DEALLOCATE(mata)
DEALLOCATE(matc)
DEALLOCATE(matf)
c4science · Help