stencil_stargeneral.F90

Go to the documentation of this file.

!! Copyright (C) 2016 U. De Giovannini, H Huebener; 2019 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program 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 General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module stencil_stargeneral_oct_m
  use affine_coordinates_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use messages_oct_m
  use stencil_oct_m
 
  implicit none
 
  private
  public ::                     &
    stencil_stargeneral_size_lapl, &
    stencil_stargeneral_extent,    &
    stencil_stargeneral_get_lapl,  &
    stencil_stargeneral_pol_lapl,  &
    stencil_stargeneral_get_arms
 
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_get_arms(this, dim, coord_system)
    type(stencil_t),            intent(inout) :: this
    integer,                    intent(in)    :: dim
    class(coordinate_system_t), intent(in)    :: coord_system
 
    real(real64)   :: theta(dim), basis_vectors(dim, dim)
    integer :: arm(dim)
    integer :: ii, jj, kk
    real(real64)   :: norm, min_norm(dim)
 
    real(real64), parameter :: tol = 1e-14_real64
 
    push_sub(stencil_stargeneral_get_arms)
 
    assert(dim <= 3)
 
    this%stargeneral%narms = 0
    this%stargeneral%arms = 0
 
    if (dim == 1) then
      !we are done
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
    select type (coord_system)
    class is (affine_coordinates_t)
      basis_vectors = coord_system%basis%vectors
    class default
      message(1) = "Stencil stargeneral can currently only be used with affine coordinate systems."
      call messages_fatal(1)
    end select
 
    if (dim == 2) then
      !get the angle between the primitive vectors
      theta(1) = acos(clamp(dot_product(basis_vectors(:, 1), basis_vectors(:, 2)), -m_one, m_one))
 
      if (theta(1) < m_pi*m_half) then
        this%stargeneral%narms = this%stargeneral%narms + 1
        arm = (/1,-1/)
        this%stargeneral%arms(this%stargeneral%narms, 1:dim) = arm
      else if (theta(1) > m_pi*m_half) then
        this%stargeneral%narms = this%stargeneral%narms + 1
        arm = (/1,+1/)
        this%stargeneral%arms(this%stargeneral%narms, 1:dim) = arm
      end if
      !if theta == pi/2 we do not need additional arms
 
      ! NB: you have supposed the axis of the 2D system is along z.
      !we are done
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
    ! dim>2
 
    !We first count how many arms we need
    theta(1) = acos(clamp(dot_product(basis_vectors(:, 1), basis_vectors(:, 2)), -m_one, m_one))
    if (abs(theta(1) - m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
    theta(2) = acos(clamp(dot_product(basis_vectors(:, 2), basis_vectors(:, 3)), -m_one, m_one))
    if (abs(theta(2)-m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
    theta(3) = acos(clamp(dot_product(basis_vectors(:, 3), basis_vectors(:, 1)), -m_one, m_one))
    if (abs(theta(3) - m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
 
    !Cubic cell, no extra arms
    if (this%stargeneral%narms == 0) then
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
 
    !We find the three closest neighbours, to define the direction of the arms
    min_norm = m_huge
    do ii= -1,1
      do jj = -1,1
        do kk = -1,1
          !The nearest points could not be along the primitive axis
          if (ii == 0 .and. jj == 0) cycle
          if (ii == 0 .and. kk == 0) cycle
          if (kk == 0 .and. jj == 0) cycle
 
          !Nor in a plane formed by orthogonal vectors
          if (abs(theta(1) - m_pi*m_half) <= 1e-8_real64 .and. kk == 0) cycle
          if (abs(theta(2) - m_pi*m_half) <= 1e-8_real64 .and. ii == 0) cycle
          if (abs(theta(3) - m_pi*m_half) <= 1e-8_real64 .and. jj == 0) cycle
 
          norm = sum((ii*basis_vectors(:, 1) + jj*basis_vectors(:, 2) + kk*basis_vectors(:, 3))**2)
 
          if (min_norm(1) - norm > tol) then
            if (min_norm(2) - min_norm(1) > tol .and. this%stargeneral%narms >= 2) then
              if (min_norm(3) - min_norm(2) > tol .and. this%stargeneral%narms == 3) then
                this%stargeneral%arms(3, 1:3) = this%stargeneral%arms(2, 1:3)
                min_norm(3) =  min_norm(2)
              end if
              this%stargeneral%arms(2, 1:3) = this%stargeneral%arms(1, 1:3)
              min_norm(2) = min_norm(1)
            end if
            min_norm(1) = norm
            this%stargeneral%arms(1, 1:3) = (/ii, jj, kk/)
            cycle
          end if
 
          ! Note: this should not be needed since we now have tolerances
          if (this%stargeneral%arms(1, 1) == -ii &
            .and. this%stargeneral%arms(1, 2) == -jj &
            .and. this%stargeneral%arms(1, 3) == -kk) cycle
 
          if (this%stargeneral%narms == 1) cycle
 
          if (min_norm(2) - norm > tol) then
            if (min_norm(3) - min_norm(2) > tol .and. this%stargeneral%narms == 3) then
              this%stargeneral%arms(3, 1:3) = this%stargeneral%arms(2, 1:3)
              min_norm(3) =  min_norm(2)
            end if
            min_norm(2) = norm
            this%stargeneral%arms(2, 1:3) = (/ii, jj, kk/)
            cycle
          end if
 
          if (this%stargeneral%narms == 2) cycle
 
          ! Note: this should not be needed since we now have tolerances
          if (this%stargeneral%arms(2, 1) == -ii &
            .and. this%stargeneral%arms(2, 2) == -jj &
            .and. this%stargeneral%arms(2, 3) == -kk) cycle
 
          if (min_norm(3) - norm > tol) then
            min_norm(3) = norm
            this%stargeneral%arms(3, 1:3) = (/ii, jj, kk/)
          end if
        end do
      end do
    end do
 
    pop_sub(stencil_stargeneral_get_arms)
  end subroutine stencil_stargeneral_get_arms
 
 
  ! ---------------------------------------------------------
  integer function stencil_stargeneral_size_lapl(this, dim, order) result(n)
    type(stencil_t),     intent(inout) :: this
    integer, intent(in) :: dim
    integer, intent(in) :: order
 
    push_sub(stencil_stargeneral_size_lapl)
 
    !normal star
    n = 2*dim*order + 1
 
    ! star general
    n = n + 2 * order * this%stargeneral%narms
 
 
    pop_sub(stencil_stargeneral_size_lapl)
  end function stencil_stargeneral_size_lapl
 
 
  ! ---------------------------------------------------------
  integer function stencil_stargeneral_extent(dir, order)
    integer, intent(in) :: dir
    integer, intent(in) :: order
 
    integer :: extent
 
    push_sub(stencil_stargeneral_extent)
 
    extent = 0
    if (dir >= 1 .and. dir <= 3) then
      if (order <= 2) then
        extent = 2
      else
        extent = order
      end if
    end if
    stencil_stargeneral_extent = extent
 
    pop_sub(stencil_stargeneral_extent)
  end function stencil_stargeneral_extent
 
 
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_get_lapl(this, dim, order)
    type(stencil_t), intent(inout) :: this
    integer,         intent(in)    :: dim
    integer,         intent(in)    :: order
 
    integer :: i, j, n
    logical :: got_center
 
    push_sub(stencil_stargeneral_get_lapl)
 
    call stencil_allocate(this, dim, stencil_stargeneral_size_lapl(this, dim, order))
 
    n = 1
    select case (dim)
    case (1)
      n = 1
      do i = 1, dim
        do j = -order, order
          if (j == 0) cycle
          n = n + 1
          this%points(i, n) = j
        end do
      end do
    case (2)
      n = 1
      do i = 1, dim
        do j = -order, order
          if (j == 0) cycle
          n = n + 1
          this%points(i, n) = j
        end do
      end do
 
      do j = -order, order
        if (j == 0) cycle
        do i = 1, this%stargeneral%narms
          n = n + 1
          this%points(1:2, n) = this%stargeneral%arms(i, 1:2)*j
        end do
      end do
 
    case (3)
      got_center = .false.
 
      n = 0
      do i = 1, dim
        do j = -order, order
 
          ! count center only once
          if (j == 0) then
            if (got_center) then
              cycle
            else
              got_center = .true.
            end if
 
          end if
          n = n + 1
          this%points(i, n) = j
        end do
      end do
 
      do j = -order, order
        if (j == 0) cycle
        do i = 1, this%stargeneral%narms
          n = n + 1
          this%points(1:3, n) = this%stargeneral%arms(i, 1:3)*j
        end do
      end do
 
    end select
 
    call stencil_init_center(this)
 
    pop_sub(stencil_stargeneral_get_lapl)
  end subroutine stencil_stargeneral_get_lapl
 
 
 
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_pol_lapl(this, dim, order, pol)
    type(stencil_t), intent(in) :: this
    integer, intent(in)          :: dim
    integer, intent(in)          :: order
    integer, intent(out)         :: pol(:,:)
 
    integer :: i, j, n, j1
    logical :: zeros_treated(dim)
 
    push_sub(stencil_stargeneral_pol_lapl)
 
    n = 1
    select case (dim)
    case (1)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
    case (2)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
 
      do j = 1, 2*order
        do i = 1, this%stargeneral%narms
          n = n + 1
          if (sum(this%stargeneral%arms(i,1:dim)) == 0) then
            pol(1:2, n) = (/j,1/)
          else
            pol(1:2, n) = (/1,j/)
          end if
        end do
      end do
 
    case (3)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
 
      ! first, treat arms with zeros: they should have polynomials with zeros in the
      ! corresponding dimension
      zeros_treated = .false.
      do i = 1, this%stargeneral%narms
        if (this%stargeneral%arms(i, 1) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/0, 1, j1/)
          end do
          zeros_treated(1) = .true.
        else if (this%stargeneral%arms(i, 2) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/j1, 0, 1/)
          end do
          zeros_treated(2) = .true.
        else if (this%stargeneral%arms(i, 3) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/1, j1, 0/)
          end do
          zeros_treated(3) = .true.
        end if
      end do ! end loop on number of arms
      ! then treat the remaining arms, use the polynomials that have not been used yet
      do i = 1, this%stargeneral%narms - count(zeros_treated)
        if (.not. zeros_treated(1)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/0, 1, j1/)
          end do
          zeros_treated(1) = .true.
        else if (.not. zeros_treated(2)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/j1, 0, 1/)
          end do
          zeros_treated(2) = .true.
        else if (.not. zeros_treated(3)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/1, j1, 0/)
          end do
          zeros_treated(3) = .true.
        end if
      end do ! end loop on number of arms
    end select
 
    pop_sub(stencil_stargeneral_pol_lapl)
  end subroutine stencil_stargeneral_pol_lapl
 
 
end module stencil_stargeneral_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: