sll_m_ellipt_2d_fem_sps_stencil_new_projector.F90

Go to the documentation of this file.

module sll_m_ellipt_2d_fem_sps_stencil_new_projector
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
 
   use sll_m_working_precision, only: f64
 
   use sll_m_vector_space_c1_block, only: sll_t_vector_space_c1_block
 
   use sll_m_linear_operator_matrix_stencil_to_stencil, only: sll_t_linear_operator_matrix_stencil_to_stencil
 
   use sll_m_linear_operator_matrix_c1_block_new, only: sll_t_linear_operator_matrix_c1_block_new
 
   implicit none
 
   public :: sll_t_ellipt_2d_fem_sps_stencil_new_projector
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   ! Working precision
   integer, parameter :: wp = f64
 
   type :: sll_t_ellipt_2d_fem_sps_stencil_new_projector
 
      integer :: n1
      integer :: n2
      integer :: p1
      integer :: p2
 
      ! Temporary storage needed for projection
      real(wp), allocatable :: qp_temp(:, :, :)
      real(wp), allocatable :: l(:, :, :)
 
   contains
 
      procedure :: init => s_ellipt_2d_fem_sps_stencil_new_projector__init
      procedure :: change_basis_matrix => s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_matrix
      procedure :: change_basis_vector => s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vector
      procedure :: change_basis_vector_inv => s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vecinv
      procedure :: free => s_ellipt_2d_fem_sps_stencil_new_projector__free
 
   end type sll_t_ellipt_2d_fem_sps_stencil_new_projector
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   ! Initializer
   subroutine s_ellipt_2d_fem_sps_stencil_new_projector__init(self, n1, n2, p1, p2, L)
      class(sll_t_ellipt_2d_fem_sps_stencil_new_projector), intent(inout) :: self
      integer, intent(in) :: n1
      integer, intent(in) :: n2
      integer, intent(in) :: p1
      integer, intent(in) :: p2
      real(wp), intent(in) :: L(:, :, :) ! matrix of barycentric coordinates
 
      self%n1 = n1
      self%n2 = n2
      self%p1 = p1
      self%p2 = p2
 
      sll_assert(size(l, 1) == 2)
      sll_assert(size(l, 2) == n2)
      sll_assert(size(l, 3) == 3)
 
      ! Allocate temporary storage
      allocate (self%L(size(l, 1), size(l, 2), size(l, 3)), source=l)
      allocate (self%Qp_temp(size(l, 1), size(l, 2), size(l, 3)))
 
   end subroutine s_ellipt_2d_fem_sps_stencil_new_projector__init
 
   ! Change basis: C1 projection of stiffness and mass matrices
   ! NOTE: 'self' has intent(inout) because temporary storage has to be assigned
   subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_matrix(self, Ql, Qp)
      class(sll_t_ellipt_2d_fem_sps_stencil_new_projector), intent(inout) :: self
      type(sll_t_linear_operator_matrix_stencil_to_stencil), intent(inout) :: Ql
      type(sll_t_linear_operator_matrix_c1_block_new), intent(inout) :: Qp
 
      integer :: i1, i2, j1, j2, k1, k2, ip, jp, ll
 
      associate(n1 => self%n1, &
                 n2 => self%n2, &
                 p1 => self%p1, &
                 p2 => self%p2)
 
         ! Fill blocks
         self%Qp_temp(:, :, :) = 0.0_wp
         qp%block1%A(:, :) = 0.0_wp
         qp%block2%A(:, :, :) = 0.0_wp
         qp%block3%A(:, :, :) = 0.0_wp
         do ll = 1, 3
            do i2 = 1, n2
               do i1 = 1, n1
                  do k2 = -p2, p2
                     do k1 = -p1, p1
 
                        j1 = i1 + k1
                        j2 = i2 + k2
 
                        ! Hand-made modulo operation
                        if (j1 < 1) then
                           j1 = j1 + n1
                        else if (j1 > n1) then
                           j1 = j1 - n1
                        end if
 
                        ! Hand-made modulo operation
                        if (j2 < 1) then
                           j2 = j2 + n2
                        else if (j2 > n2) then
                           j2 = j2 - n2
                        end if
 
                        ! block 1: 3 x 3
                        if (i1 <= 2 .and. j1 <= 2) then
                           self%Qp_temp(i1, i2, ll) = self%Qp_temp(i1, i2, ll) + ql%A(k1, k2, i1, i2)*self%L(j1, j2, ll)
                           ! block 2: 3 x (n1-2)*n2
                        else if (i1 <= 2 .and. j1 <= 2 + p1) then
                           qp%block2%A(ll, j1 - 2, j2) = qp%block2%A(ll, j1 - 2, j2) + self%L(i1, i2, ll)*ql%A(k1, k2, i1, i2)
                           ! block 3: (n1-2)*n2 x 3
                        else if (i1 <= 2 + p1 .and. j1 <= 2) then
                           qp%block3%A(i1 - 2, i2, ll) = qp%block3%A(i1 - 2, i2, ll) + ql%A(k1, k2, i1, i2)*self%L(j1, j2, ll)
                        end if
 
                     end do
                  end do
               end do
            end do
         end do
         ! complete calculation of block 1
         do ip = 1, 3
            do jp = 1, 3
               do i2 = 1, n2
                  do i1 = 1, 2
                     qp%block1%A(ip, jp) = qp%block1%A(ip, jp) + self%L(i1, i2, ip)*self%Qp_temp(i1, i2, jp)
                  end do
               end do
            end do
         end do
 
         ! block 4: (n1-2)*n2 x (n1-2)*n2
         qp%block4%A(:, :, :, :) = 0.0_wp
         do i2 = 1, n2
            do i1 = 1, n1 - 2
               do k2 = -p2, p2
                  do k1 = -p1, p1
                     j1 = modulo(i1 - 1 + k1, n1 - 2)
                     ll = modulo(i1 + 1 + k1, n1)
                     if (ll == 2 + j1) qp%block4%A(k1, k2, i1, i2) = ql%A(k1, k2, i1 + 2, i2)
                  end do
               end do
            end do
         end do
 
      end associate
 
   end subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_matrix
 
   ! Change basis: C1 projection of vectors
   subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vector(self, V, Vp)
      class(sll_t_ellipt_2d_fem_sps_stencil_new_projector), intent(in) :: self
      real(wp), intent(in) :: V(:, :)
      type(sll_t_vector_space_c1_block), intent(inout) :: Vp
 
      integer :: i1, i2, ll
 
      associate(n1 => self%n1, &
                 n2 => self%n2, &
                 p1 => self%p1, &
                 p2 => self%p2)
 
         ! Checks
         sll_assert(size(v, 1) == n1)
         sll_assert(size(v, 2) == n2)
         sll_assert(size(vp%vd%array) == 3)
         sll_assert(size(vp%vs%array, 1) == n1 - 2 + 2*p1)
         sll_assert(size(vp%vs%array, 2) == n2 + 2*p2)
 
         vp%vd%array(:) = 0.0_wp
         do ll = 1, 3
            do i2 = 1, n2
               do i1 = 1, 2
                  vp%vd%array(ll) = vp%vd%array(ll) + self%L(i1, i2, ll)*v(i1, i2)
               end do
            end do
         end do
 
         do i2 = 1, n2
            do i1 = 3, n1
               vp%vs%array(i1 - 2, i2) = v(i1, i2)
            end do
         end do
 
      end associate
 
   end subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vector
 
   ! Change basis: C1 projection of vectors
   subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vecinv(self, Vp, V)
      class(sll_t_ellipt_2d_fem_sps_stencil_new_projector), intent(in) :: self
      type(sll_t_vector_space_c1_block), intent(inout) :: Vp
      real(wp), intent(inout) :: V(:)
 
      integer :: i, i1, i2, j, j1, j2, ll
 
      associate(n1 => self%n1, &
                 n2 => self%n2, &
                 p1 => self%p1, &
                 p2 => self%p2)
 
         ! Checks
         sll_assert(size(vp%vd%array) == 3)
         sll_assert(size(vp%vs%array, 1) == n1 - 2 + 2*p1)
         sll_assert(size(vp%vs%array, 2) == n2 + 2*p2)
         sll_assert(size(v) == n1*n2)
 
         v(:) = 0.0_wp
         do ll = 1, 3
            do i2 = 1, n2
               do i1 = 1, 2
                  i = (i1 - 1)*n2 + i2
                  v(i) = v(i) + self%L(i1, i2, ll)*vp%vd%array(ll)
               end do
            end do
         end do
 
         do j2 = 1, n2
            do j1 = 1, n1 - 2
               j = (j1 - 1)*n2 + j2
               v(2*n2 + j) = vp%vs%array(j1, j2)
            end do
         end do
 
      end associate
 
   end subroutine s_ellipt_2d_fem_sps_stencil_new_projector__change_basis_vecinv
 
   ! Deallocate allocatables
   subroutine s_ellipt_2d_fem_sps_stencil_new_projector__free(self)
      class(sll_t_ellipt_2d_fem_sps_stencil_new_projector), intent(inout) :: self
 
      deallocate (self%Qp_temp)
      deallocate (self%L)
 
   end subroutine s_ellipt_2d_fem_sps_stencil_new_projector__free
 
end module sll_m_ellipt_2d_fem_sps_stencil_new_projector