Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/dismat.f90
03540012fedefeature/refactoring

Last Change
Raw File

dismat.f90
View Options

!>
!> @file dismat.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
SUBROUTINE dismat(spl, mat)
!
! Assembly of FE matrix mat using spline spl
!
USE bsplines
USE matrix
IMPLICIT NONE
TYPE(spline2d), INTENT(in) :: spl
TYPE(gbmat), INTENT(inout) :: mat
!
INTEGER :: n1, nidbas1, ndim1, ng1
INTEGER :: n2, nidbas2, ndim2, ng2
INTEGER :: i, j, ig1, ig2
INTEGER :: iterm, iw1, iw2, igw1, igw2, it1, it2, igt1, igt2, irow, jcol
DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:,:)
DOUBLE PRECISION:: contrib
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER, ALLOCATABLE :: idert(:,:,:,:), iderw(:,:,:,:) ! Derivative order
DOUBLE PRECISION, ALLOCATABLE :: coefs(:,:,:) ! Terms in weak form
INTEGER, ALLOCATABLE :: left1(:), left2(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
CALL get_dim(spl%sp2, ndim2, n2, nidbas2)
WRITE(*,'(/a, 5i6)') 'n1, nidbas1, ndim1 =', n1, nidbas1, ndim1
WRITE(*,'(a, 5i6)') 'n2, nidbas2, ndim2 =', n2, nidbas2, ndim2
!
! Gauss quadature
!
CALL get_gauss(spl%sp1, ng1)
CALL get_gauss(spl%sp2, ng2)
ALLOCATE(xg1(ng1), wg1(ng1))
ALLOCATE(xg2(ng1), wg2(ng1))
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms,2,ng1,ng2), iderw(kterms,2,ng1,ng2))
ALLOCATE(coefs(kterms,ng1,ng2))
!
ALLOCATE(fun1(0:nidbas1,0:1,ng1)) ! Spline and 1st derivative
ALLOCATE(fun2(0:nidbas2,0:1,ng2)) !
!===========================================================================
! 2.0 Assembly loop
!
ALLOCATE(left1(ng1))
ALLOCATE(left2(ng2))
DO i=1,n1
CALL get_gauss(spl%sp1, ng1, i, xg1, wg1)
left1 = i
CALL basfun(xg1, spl%sp1, fun1, left1)
DO j=1,n2
CALL get_gauss(spl%sp2, ng2, j, xg2, wg2)
left2 = j
CALL basfun(xg2, spl%sp2, fun2, left2)
!
DO ig1=1,ng1
DO ig2=1,ng2
CALL coefeq(xg1(ig1), xg2(ig2), &
& idert(:,:,ig1,ig2), &
& iderw(:,:,ig1,ig2), &
& coefs(:,ig1,ig2))
END DO
END DO
!
DO iw1=0,nidbas1 ! Weight function in dir 1
igw1 = i+iw1
DO iw2=0,nidbas2 ! Weight function in dir 2
igw2 = MODULO(j+iw2-1, n2) + 1
irow = igw2 + (igw1-1)*n2
DO it1=0,nidbas1 ! Test function in dir 1
igt1 = i+it1
DO it2=0,nidbas2 ! Test function in dir 2
igt2 = MODULO(j+it2-1, n2) + 1
jcol = igt2 + (igt1-1)*n2
!-------------
contrib = 0.0d0
DO ig1=1,ng1
DO ig2=1,ng2
DO iterm=1,kterms
contrib = contrib + &
& fun1(iw1,iderw(iterm,1,ig1,ig2),ig1) * &
& fun2(iw2,iderw(iterm,2,ig1,ig2),ig2) * &
& coefs(iterm,ig1,ig2) * &
& fun2(it2,idert(iterm,2,ig1,ig2),ig2) * &
& fun1(it1,idert(iterm,1,ig1,ig2),ig1) * &
& wg1(ig1) * wg2(ig2)
END DO
END DO
END DO
CALL updtmat(mat, irow, jcol, contrib)
!-------------
END DO
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(xg1, wg1, fun1)
DEALLOCATE(xg2, wg2, fun2)
DEALLOCATE(idert, iderw, coefs)
DEALLOCATE(left1,left2)
!
CONTAINS
SUBROUTINE coefeq(x, y, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt,1))
!
! Weak form = Int(x*dw/dx*dt/dx + 1/x*dw/dy*dt/dy)dxdy
!
c(1) = x !
idt(1,1) = 1
idt(1,2) = 0
idw(1,1) = 1
idw(1,2) = 0
!
c(2) = 1.d0/x
idt(2,1) = 0
idt(2,2) = 1
idw(2,1) = 0
idw(2,2) = 1
END SUBROUTINE coefeq
END SUBROUTINE dismat
c4science · Help