Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/src/conmat2.tpl
f1da25c60953feature/refactoring

Last Change
Raw File

conmat2.tpl
View Options

!>
!> @file conmat2.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 local matrices E and D are precalculted
! (on all n1/n2 intervals and nterms weak-form terms
!
INTERFACE
SUBROUTINE coefeq(x, y, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:)
END SUBROUTINE coefeq
END INTERFACE
INTEGER, OPTIONAL :: maxder(2) ! maximum oder of derivatives
!
INTEGER :: n1, nidbas1, ndim1
INTEGER :: n2, nidbas2, ndim2
INTEGER :: ng1, ng2
INTEGER :: i1, i2, ig1, ig2
INTEGER :: igt1, igt2, igw1, igw2, irow, jcol
INTEGER, ALLOCATABLE :: left1(:), left2(:)
!
LOGICAL :: nlper1, nlper2
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER :: k, kmaxder, it1, iw1, it2, iw2
INTEGER, ALLOCATABLE :: idert(:,:), iderw(:,:) ! Derivative order
DOUBLE PRECISION, ALLOCATABLE :: coefs(:,:,:) ! Terms in weak form
DOUBLE PRECISION :: dummy(mat%nterms)
!
DOUBLE PRECISION, POINTER :: xg1(:,:), xg2(:,:), wg1(:,:), wg2(:,:)
DOUBLE PRECISION, ALLOCATABLE :: fun1(:,:,:), fun2(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: mata(:,:,:,:), matc(:,:), matcd(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: matd(:,:,:,:,:), mate(:,:,:,:,:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
CALL get_dim(spl%sp2, ndim2, n2, nidbas2)
nlper1 = spl%sp1%period
nlper2 = spl%sp2%period
!
! Gauss points and weights on all intervals
!
xg1 => spl%sp1%gausx ! xg1(ng1,n1)
wg1 => spl%sp1%gausw ! wg1(ng1,n1)
ng1 = SIZE(xg1,1)
xg2 => spl%sp2%gausx
wg2 => spl%sp2%gausw
ng2 = SIZE(xg2,1)
!
! Derivative orders in the weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms,2))
ALLOCATE(iderw(kterms,2))
CALL coefeq(xg1(1,1), xg2(1,1), idert, iderw, dummy)
!
! Precalc matrix E in dimension 1
!
kmaxder = 1
IF(PRESENT(maxder)) kmaxder = maxder(1)
ALLOCATE(fun1(0:nidbas1,0:kmaxder,ng1))
ALLOCATE(left1(ng1))
ALLOCATE(mate(0:nidbas1,0:nidbas1,ng1,kterms,n1))
DO i1=1,n1
left1 = i1
CALL basfun(xg1(:,i1), spl%sp1, fun1, left1)
DO k=1,kterms
DO ig1=1,ng1
DO iw1=0,nidbas1
DO it1=0,nidbas1
mate(it1,iw1,ig1,k,i1) = wg1(ig1,i1) * &
& fun1(it1,idert(k,1),ig1) * &
& fun1(iw1,iderw(k,1),ig1)
END DO
END DO
END DO
END DO
END DO
DEALLOCATE(fun1)
DEALLOCATE(left1)
!
! Precalc matrix D in dimension 2
!
kmaxder = 1
IF(PRESENT(maxder)) kmaxder = maxder(2)
ALLOCATE(fun2(0:nidbas2,0:kmaxder,ng2))
ALLOCATE(left2(ng2))
ALLOCATE(matd(0:nidbas2,0:nidbas2,ng2,kterms,n2))
DO i2=1,n2
left2 = i2
CALL basfun(xg2(:,i2), spl%sp2, fun2, left2)
DO k=1,kterms
DO ig2=1,ng2
DO iw2=0,nidbas2
DO it2=0,nidbas2
matd(it2,iw2,ig2,k,i2) = wg2(ig2,i2) * &
& fun2(it2,idert(k,2),ig2) * &
& fun2(iw2,iderw(k,2),ig2)
END DO
END DO
END DO
END DO
END DO
DEALLOCATE(fun2)
DEALLOCATE(left2)
!===========================================================================
! 2.0 Assembly loop
!
! Physical coefficients in Weak form
!
ALLOCATE(coefs(kterms,ng1,ng2))
ALLOCATE(matc(ng1,ng2))
!
! Allocate local matrix A
!
ALLOCATE(mata(0:nidbas1,0:nidbas1,0:nidbas2,0:nidbas2))
ALLOCATE(matcd(ng1,0:nidbas2,0:nidbas2))
!
DO i1=1,n1
DO i2=1,n2
!
! Coefficients of the weak form
!
DO ig1=1,ng1
DO ig2=1,ng2
CALL coefeq(xg1(ig1,i1), xg2(ig2,i2), &
& idert, iderw, coefs(:,ig1,ig2))
END DO
END DO
!
! Compute local matrix: A <- E*(C*D^T) + A
!
mata = 0.0d0
DO k=1,kterms
!
matc(1:ng1,1:ng2) = coefs(k,1:ng1,1:ng2)
!
CALL dgemm('N', 'T', ng1, (nidbas2+1)*(nidbas2+1), ng2, 1.0d0, &
& matc, ng1, matd(0,0,1,k,i2), (nidbas2+1)*(nidbas2+1), 0.0d0, &
& matcd, ng1)
CALL dgemm('N', 'N', (nidbas1+1)*(nidbas1+1), (nidbas2+1)*(nidbas2+1), &
& ng1, 1.0d0, mate(0,0,1,k,i1), (nidbas1+1)*(nidbas1+1), matcd, ng1, 1.0d0, &
& mata, (nidbas1+1)*(nidbas1+1))
!
END DO
!
! Map local matrix A to global matrix
!
DO it1=0,nidbas1
igt1 = i1+it1; IF(nlper1) igt1 = MODULO(igt1-1,n1) + 1
DO it2=0,nidbas2
igt2 = i2+it2; IF(nlper2) igt2 = MODULO(igt2-1, n2) + 1
irow = igt2 + (igt1-1)*n2
DO iw1=0,nidbas1
igw1 = i1+iw1; IF(nlper1) igw1 = MODULO(igw1-1,n1) + 1
DO iw2=0,nidbas2
igw2 = i2+iw2; IF(nlper2) igw2 = MODULO(igw2-1, n2) + 1
jcol = igw2 + (igw1-1)*n2
CALL updtmat(mat, irow, jcol, mata(it1,iw1,it2,iw2))
END DO
END DO
END DO
END DO
!
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(idert, iderw, coefs)
DEALLOCATE(mata)
DEALLOCATE(matc)
DEALLOCATE(matd)
DEALLOCATE(matcd)
DEALLOCATE(mate)
c4science · Help