constructor Function

  function constructor(template, boundary_locs, boundary_types) result(this)
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! Build a tridiagonal matrix block for finite differences. See the
    ! end of the documentation of the [[fin_diff_block(type)]] type
    ! for a description of how boundary conditions are treated.
    !
    class(abstract_field), intent(in)           :: template
      !! A scalar field with the same grid as any fields passed as
      !! arguments to the [[fin_diff_block(type):solve_for]] method.
    integer, dimension(:), optional, intent(in) :: boundary_locs
      !! The locations in the raw representation of `rhs` for which
      !! boundary conditions are specified. Defaults to there being
      !! none.
    integer, dimension(:), optional, intent(in) :: boundary_types
      !! Integers specifying the type of boundary condition. The type
      !! of boundary condition corresponding to a given integer is
      !! specified in [[boundary_types_mod]]. Only Dirichlet and
      !! Neumann conditions are supported. Defaults to Dirichlet. The
      !! order in which they are stored must match that of
      !! `boundary_locs`.
    type(fin_diff_block)                        :: this
      !! A new finite difference operator

    integer                                   :: i, n, pos
    logical                                   :: use_cached
    real(r8), dimension(:), allocatable, save :: cached_dx_c
    real(r8), save                            :: cached_upper_bound, &
                                                 cached_lower_bound
    class(vector_field), pointer              :: grid
    class(abstract_field), allocatable, save  :: cached_field_type
    real(r8), dimension(:,:), allocatable     :: domain

    if (present(boundary_locs)) then
      this%boundary_locs = boundary_locs
    else
      allocate(this%boundary_locs(0))
    end if
    if (present(boundary_types)) then
      this%boundary_types = boundary_types
    else
      allocate(this%boundary_types(size(this%boundary_locs)))
      this%boundary_types = dirichlet
    end if

    n = template%raw_size()
    domain = template%domain()
    ! Try to use cached copy of inverse grid spacing, if available and suitable 
    if (allocated(cached_field_type)) then
      use_cached = same_type_as(template, cached_field_type) .and. &
                   abs(domain(1,1) - cached_lower_bound) < 1e-10 .and. &
                   abs(domain(1,2) - cached_upper_bound) < 1e-10 .and. &
                   n==size(cached_dx_c)
    else
      use_cached = .false.
    end if
    ! Construct an array containing the distance between grid
    ! points, if necessary
    if (.not. use_cached) then
      if (allocated(cached_field_type)) then
        deallocate(cached_field_type)
        deallocate(cached_dx_c)
      end if
#ifdef DEBUG
      call logger%debug('fin_diff_block','Calculating and caching '// &
                       'grid spacings.')
#endif
      cached_lower_bound = domain(1,1)
      cached_upper_bound = domain(1,2)
      allocate(cached_field_type, mold=template)
      allocate(cached_dx_c(n))
      call template%allocate_vector_field(grid)
      grid = template%grid_spacing()
      ! Use this array to temporarily hold the grid spacings
      cached_dx_c = grid%raw()
      cached_dx_c(2:n-1) = 0.5_r8/cached_dx_c(2:n-1)
      cached_dx_c(1) = 1._r8/cached_dx_c(1)
      cached_dx_c(n) = 1._r8/cached_dx_c(n)
    end if
    ! Create tridiagonal matrix
    allocate(this%diagonal(n))
    this%diagonal = 0._r8
    this%diagonal(1) = -cached_dx_c(1)
    this%diagonal(n) = cached_dx_c(n)
    this%super_diagonal = cached_dx_c(1:n-1)
    this%sub_diagonal = -cached_dx_c(2:n)
    ! Apply boundary conditions
    this%diagonal(1) = this%diagonal(1) + 1e-7
    do i = 1, size(this%boundary_locs)
      pos = this%boundary_locs(i)
      select case(this%boundary_types(i))
      case(dirichlet)
        this%diagonal(pos) = 1._r8
        if (pos < n) this%super_diagonal(pos) = 0._r8
        if (pos > 1) this%sub_diagonal(pos-1) = 1.e-7_r8
      case(neumann)
        if (pos == n) then
          this%diagonal(n) = cached_dx_c(n)
          this%sub_diagonal(n-1) = -cached_dx_c(n)
        else if (pos == 1) then
          this%diagonal(1) = -cached_dx_c(1)
          this%super_diagonal(1) = cached_dx_c(1)
        else
          this%super_diagonal(pos) = cached_dx_c(pos)
          this%sub_diagonal(pos-1) = cached_dx_c(pos)
          this%diagonal(pos) = 0._r8
        end if
      case(free_boundary)
        continue
      case default
        call logger%fatal('fin_diff_block','Boundary condition of '// &
                          'type other than Dirichlet or Neumann encountered.')
        error stop
      end select
    end do
    do i = 1, size(this%boundary_locs)
      pos = this%boundary_locs(i)
      select case(this%boundary_types(i))
      case(dirichlet)
        this%diagonal(pos) = 1._r8
        if (pos < n) this%super_diagonal(pos) = 0._r8
        if (pos > 1) this%sub_diagonal(pos-1) = 1.e-7_r8
      case(neumann)
        if (pos == n) then
          this%diagonal(n) = cached_dx_c(n)
          this%sub_diagonal(n-1) = -cached_dx_c(n)
        else if (pos == 1) then
          this%diagonal(1) = -cached_dx_c(1)
          this%super_diagonal(1) = cached_dx_c(1)
        else
          this%super_diagonal(pos) = cached_dx_c(pos)
          this%sub_diagonal(pos-1) = cached_dx_c(pos)
          this%diagonal(pos) = 0._r8
        end if
      case(free_boundary)
        continue
      case default
        call logger%fatal('fin_diff_block','Boundary condition of '// &
                          'type other than Dirichlet or Neumann encountered.')
        error stop
      end select
    end do
    !this%magnitude = dnrm2(n-1, this%super_diagonal, 1)
#ifdef DEBUG
    call logger%debug('fin_diff_block','Instantiated a finite difference '// &
                      'block object.')
#endif
  end function constructor