Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/bsplines/conmat.tpl
1a965ce17d3bfeature/refactoring

Last Change
Raw File

conmat.tpl
View Options

!>
!> @file conmat.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/>.
!>
!> @author
!> (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, 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
LOGICAL, OPTIONAL :: nat_order ! Natural ordering for 2d-1d mapping
!
INTEGER :: n1, nidbas1, ndim1, n1e
INTEGER :: n2, nidbas2, ndim2, n2e
INTEGER :: ng1, ng2
INTEGER :: i1, i2, ig1, ig2
INTEGER :: igt1, igt2, igw1, igw2, irow, jcol
INTEGER, ALLOCATABLE :: left1(:), left2(:)
!
LOGICAL :: nlper1, nlper2, nlnat
!
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, POINTER :: xg1(:,:), xg2(:,:), wg1(:,:), wg2(:,:)
DOUBLE PRECISION, ALLOCATABLE :: fun1(:,:,:,:), fun2(:,:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: mata(:,:,:,:), matc(:,:)
DOUBLE PRECISION, ALLOCATABLE :: matg(:,:,:), matf(:,:,:), matcg(:,:,:)
!===========================================================================
! 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
!
n1e = n1+nidbas1 ! Number of elements in 1st coordinate
n2e = n2+nidbas2 ! Number of elements in 2nd coordinate
iF(nlper2) n2e = n2
!
! 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)
!
! Splines on all intervals
!
kmaxder = 1
IF(PRESENT(maxder)) kmaxder = maxder(1)
ALLOCATE(fun1(0:nidbas1,0:kmaxder,ng1,n1))
ALLOCATE(left1(ng1))
DO i1=1,n1
left1 = i1
CALL basfun(xg1(:,i1), spl%sp1, fun1(:,:,:,i1), left1)
END DO
DEALLOCATE(left1)
!
kmaxder = 1
IF(PRESENT(maxder)) kmaxder = maxder(2)
ALLOCATE(fun2(0:nidbas2,0:kmaxder,ng2,n2))
ALLOCATE(left2(ng2))
DO i2=1,n2
left2 = i2
CALL basfun(xg2(:,i2), spl%sp2, fun2(:,:,:,i2), left2)
END DO
DEALLOCATE(left2)
!
! Ordering in local to global matrix mapping
!
nlnat = .FALSE.
IF(PRESENT(nat_order)) nlnat = nat_order
!===========================================================================
! 2.0 Assembly loop
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms,2))
ALLOCATE(iderw(kterms,2))
ALLOCATE(coefs(kterms,ng1,ng2))
!
! Allocate local matrices
!
ALLOCATE(mata(0:nidbas1,0:nidbas1,0:nidbas2,0:nidbas2))
ALLOCATE(matc(ng1,ng2))
ALLOCATE(matg(0:nidbas2,0:nidbas2,ng2))
ALLOCATE(matf(0:nidbas1,0:nidbas1,ng1))
ALLOCATE(matcg(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)
!
DO it1=0,nidbas1
DO iw1=0,nidbas1
DO ig1=1,ng1
matf(it1,iw1,ig1) = wg1(ig1,i1) * &
& fun1(it1,idert(k,1),ig1,i1) * &
& fun1(iw1,iderw(k,1),ig1,i1)
END DO
END DO
END DO
!
DO it2=0,nidbas2
DO iw2=0,nidbas2
DO ig2=1,ng2
matg(it2,iw2,ig2) = wg2(ig2,i2) * &
& fun2(it2,idert(k,2),ig2,i2) * &
& fun2(iw2,iderw(k,2),ig2,i2)
END DO
END DO
END DO
!
CALL dgemm('N', 'T', ng1, (nidbas2+1)*(nidbas2+1), ng2, 1.0d0, &
& matc, ng1, matg, (nidbas2+1)*(nidbas2+1), 0.0d0, &
& matcg, ng1)
CALL dgemm('N', 'N', (nidbas1+1)*(nidbas1+1), (nidbas2+1)*(nidbas2+1), &
& ng1, 1.0d0, matf, (nidbas1+1)*(nidbas1+1), matcg, 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 = glmap(igt1, igt2, n1e, n2e)
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 = glmap(igw1, igw2, n1e, n2e)
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(fun1)
DEALLOCATE(fun2)
DEALLOCATE(idert, iderw, coefs)
DEALLOCATE(mata)
DEALLOCATE(matc)
DEALLOCATE(matg)
DEALLOCATE(matcg)
DEALLOCATE(matf)
!
CONTAINS
INTEGER FUNCTION glmap(i,j,n1,n2)
INTEGER, INTENT(in) :: i,j,n1,n2
IF(nlnat) THEN
glmap = (j-1)*n1 + i
ELSE
glmap = (i-1)*n2 + j
END IF
END FUNCTION glmap
c4science · Help