maxwell_boundary_op.F90

Go to the documentation of this file.

!! Copyright (C) 2019 R. Jestaedt, H. Appel, F. Bonafe, M. Oliveira, 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 maxwell_boundary_op_oct_m
  use accel_oct_m
  use box_sphere_oct_m
  use box_parallelepiped_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use external_waves_oct_m
  use global_oct_m
  use grid_oct_m
  use index_oct_m
  use io_oct_m
  use io_function_oct_m
  use math_oct_m
  use linear_medium_to_em_field_oct_m
  use math_oct_m
  use maxwell_function_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use par_vec_oct_m
  use parser_oct_m
  use profiling_oct_m
  use states_elec_oct_m
  use string_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use space_oct_m
  use states_mxll_oct_m
 
  implicit none
 
  private
  public ::                          &
    bc_mxll_init,                    &
    bc_mxll_end,                     &
    bc_mxll_t,                       &
    bc_mxll_generate_pml_parameters, &
    bc_mxll_initialize_pml_simple,   &
    inner_and_outer_points_mapping,  &
    surface_grid_points_mapping
 
  type pml_t
    real(real64)         :: width
    integer              :: points_number
    integer, allocatable :: points_map(:)
    integer, allocatable :: points_map_inv(:)
    real(real64)         :: power
    real(real64)         :: refl_error
    real(real64),   allocatable :: kappa(:,:)
    real(real64),   allocatable :: sigma_e(:,:)
    real(real64),   allocatable :: sigma_m(:,:)
    real(real64),   allocatable :: a(:,:)
    real(real64),   allocatable :: b(:,:)
    real(real64),   allocatable :: c(:,:)
    real(real64),   allocatable :: mask(:)
    complex(real64),   allocatable :: aux_ep(:,:,:)
    complex(real64),   allocatable :: aux_mu(:,:,:)
    complex(real64),   allocatable :: conv_plus(:,:,:)
    complex(real64),   allocatable :: conv_minus(:,:,:)
    complex(real64),   allocatable :: conv_plus_old(:,:,:)
    complex(real64),   allocatable :: conv_minus_old(:,:,:)
    logical              :: parameters_initialized = .false.
    ! GPU buffers
    type(accel_mem_t)    :: buff_a, buff_b, buff_c, buff_conv_plus, buff_conv_minus, buff_map, buff_conv_plus_old
  end type pml_t
 
  type bc_mxll_t
    integer              :: bc_type(3)
    integer              :: bc_ab_type(3)
    real(real64)         :: bc_bounds(2, 3)
    logical              :: ab_user_def
    real(real64),   allocatable :: ab_ufn(:)
 
    real(real64)         :: ab_width
    real(real64)         :: mask_width
    integer              :: mask_points_number(3)
    integer, allocatable :: mask_points_map(:,:)
    real(real64),   allocatable :: mask(:,:)
 
    integer              :: der_bndry_mask_points_number
    integer, allocatable :: der_bndry_mask_points_map(:)
    real(real64),   allocatable :: der_bndry_mask(:)
 
    type(pml_t)          :: pml
    type(single_medium_box_t)   :: medium(3)
 
    integer              :: constant_points_number
    integer, allocatable :: constant_points_map(:)
    complex(real64),   allocatable :: constant_rs_state(:,:)
    type(accel_mem_t)    :: buff_constant_points_map
 
    integer              :: mirror_points_number(3)
    integer, allocatable :: mirror_points_map(:,:)
 
    logical                :: do_plane_waves = .false.   !! look here afterwards
    type(external_waves_t) :: plane_wave
    logical                :: plane_waves_dims(1:3) = .false.
 
    real(real64)         :: zero_width
    integer              :: zero_points_number(3)
    integer, allocatable :: zero_points_map(:,:)
    real(real64),   allocatable :: zero(:,:)
  end type bc_mxll_t
 
  integer, public, parameter ::   &
    MXLL_BC_ZERO          = 0,    &
    mxll_bc_constant      = 1,    &
    mxll_bc_mirror_pec    = 2,    &
    mxll_bc_mirror_pmc    = 3,    &
    mxll_bc_plane_waves   = 4,    &
    mxll_bc_periodic      = 5,    &
    mxll_bc_medium        = 6
 
  integer, parameter ::                  &
    MXLL_PLANE_WAVES_NEGATIVE_SIDE = -1, &
    mxll_plane_waves_positive_side = 1
 
  integer, public, parameter ::   &
    MXLL_AB_NOT_ABSORBING = 0,    &
    mxll_ab_mask          = 1,    &
    mxll_ab_cpml          = 2,    &
    mxll_ab_mask_zero     = 7
 
contains
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_init(bc, namespace, space, gr, st)
    type(bc_mxll_t),          intent(inout) :: bc
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(states_mxll_t),      intent(inout) :: st
 
    integer             :: idim, nlines, icol, ncols, ab_shape_dim
    real(real64)        :: bounds(2, space%dim), ab_bounds(2, space%dim)
    type(block_t)       :: blk
    character(len=1024) :: string
    character(len=50)   :: ab_type_str
    character(len=1), dimension(3), parameter :: dims = ["x", "y", "z"]
    logical             :: plane_waves_check, ab_mask_check, ab_pml_check
    logical             :: constant_check, zero_check
    real(real64) :: ep_factor, mu_factor, sigma_e_factor, sigma_m_factor
 
    push_sub(bc_mxll_init)
 
    call profiling_in('BC_MXLL_INIT')
 
    plane_waves_check = .false.
    ab_mask_check = .false.
    ab_pml_check = .false.
    constant_check = .false.
    zero_check = .false.
 
    bc%ab_user_def = .false.
    bc%bc_ab_type(:) = mxll_ab_not_absorbing ! default option
 
    !%Variable MaxwellAbsorbingBoundaries
    !%Type block
    !%Section Maxwell
    !%Description
    !% Type of absorbing boundaries used for Maxwell propagation in each direction.
    !%
    !% Example:
    !%
    !% <tt>%MaxwellAbsorbingBoundaries
    !% <br>&nbsp;&nbsp;   cpml | cpml | cpml
    !% <br>%</tt>
    !%
    !%Option not_absorbing 0
    !% No absorbing boundaries.
    !%Option mask 1
    !% A mask equal to the wavefunctions mask is applied to the Maxwell states at the boundaries
    !%Option cpml 2
    !% Perfectly matched layer absorbing boundary
    !%Option mask_zero 7
    !% Absorbing boundary region is set to zero
    !%End
 
    if (parse_block(namespace, 'MaxwellAbsorbingBoundaries', blk) == 0) then
      ! find out how many lines (i.e. states) the block has
      nlines = parse_block_n(blk)
      if (nlines /= 1) then
        message(1) = 'MaxwellAbsorbingBounaries has to consist of one line!'
        call messages_fatal(1, namespace=namespace)
      end if
 
      ncols = parse_block_cols(blk, 0)
      if (ncols /= 3) then
        message(1) = 'MaxwellAbsorbingBoundaries has to consist of three columns!'
        call messages_fatal(1, namespace=namespace)
      end if
 
      do icol = 1, ncols
        call parse_block_integer(blk, 0, icol-1, bc%bc_ab_type(icol))
      end do
 
      call parse_block_end(blk)
    end if
 
    do idim = 1, 3
      if (bc%bc_ab_type(idim) == mxll_ab_mask) ab_mask_check = .true.
      if (bc%bc_ab_type(idim) == mxll_ab_cpml) ab_pml_check = .true.
      if (bc%bc_type(idim) == mxll_bc_constant) constant_check = .true.
      if (bc%bc_type(idim) == mxll_bc_zero) zero_check = .true.
    end do
 
    if (ab_mask_check .or. ab_pml_check) then
 
      call messages_print_with_emphasis(msg='Maxwell Absorbing Boundaries', namespace=namespace)
      write(message(1),'(a)')  "Please keep in mind that"
      write(message(2),'(a)')  "with larger ABWidth, comes great responsibility."
      write(message(3),'(a)')  "AbsorbingBoundaries occupy space in the simulation box,"
      write(message(4),'(a)')  "hence choose your Lsize wisely."
      call messages_info(4, namespace=namespace)
 
    end if
 
    do idim = 1, st%dim
      select case (bc%bc_type(idim))
 
      case (mxll_bc_zero, mxll_bc_mirror_pec, mxll_bc_mirror_pmc)
 
        bounds(1, idim) = (gr%idx%nr(2, idim) - gr%idx%enlarge(idim))*gr%spacing(idim)
        bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
 
      case (mxll_bc_constant, mxll_bc_periodic)
 
        bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
        bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
 
      case (mxll_bc_plane_waves)
 
        bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
        bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
        plane_waves_check = .true.
        bc%do_plane_waves = .true.
 
      case (mxll_bc_medium)
        call bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
        call maxwell_medium_points_mapping(bc, gr, bounds)
        call bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
 
      end select
 
      select type (box => gr%box)
      type is (box_sphere_t)
        ab_shape_dim = 1
        if (space%is_periodic()) then
          message(1) = "Sphere box shape can only work for non-periodic systems"
          call messages_fatal(1, namespace=namespace)
        end if
      type is (box_parallelepiped_t)
        if (bc%bc_type(idim) == mxll_bc_periodic .and. .not. box%axes%orthogonal) then
          message(1) = "Maxwell propagation does not work for non-orthogonal boxes with periodic boundary conditions."
          call messages_fatal(1, namespace=namespace)
        end if
 
        ab_shape_dim = space%dim
        ab_bounds(1, idim) = bounds(1, idim)
        ab_bounds(2, idim) = bounds(1, idim)
      class default
        message(1) = "Box shape for Maxwell propagation not supported yet"
        call messages_fatal(1, namespace=namespace)
      end select
 
      if (bc%bc_ab_type(idim) /= mxll_ab_not_absorbing) then
 
        call messages_print_var_option("MaxwellAbsorbingBoundaries", bc%bc_ab_type(idim), namespace=namespace)
 
        select case (bc%bc_ab_type(idim))
        case (mxll_ab_mask_zero)
          if (bc%bc_type(idim) == mxll_bc_periodic) then
            message(1) = "Zero absorbing boundary conditions do not work in periodic directions"
            call messages_fatal(1, namespace=namespace)
          end if
 
          call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
          bc%zero_width = bc%ab_width
 
        case (mxll_ab_mask)
          call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
          bc%mask_width = bc%ab_width
 
        case (mxll_ab_cpml)
          call bc_mxll_pml_init(bc, gr, namespace, ab_bounds, idim)
 
        case default
          message(1) = "Absorbing boundary type not implemented for Maxwell propagation"
          call messages_fatal(1, namespace=namespace)
        end select
 
      end if
 
      select case (bc%bc_ab_type(idim))
      case (mxll_ab_mask, mxll_ab_cpml, mxll_ab_mask_zero)
        bounds(1, idim) = ab_bounds(1, idim)
        bounds(2, idim) = bounds(2, idim)
        bc%bc_bounds(:, idim) = bounds(:, idim)
      case default
        bc%bc_bounds(:, idim) = bounds(:, idim)
      end select
 
      select type (box => gr%box)
      type is (box_parallelepiped_t)
        select case (bc%bc_ab_type(idim))
        case (mxll_ab_cpml)
          ab_type_str = "PML"
        case (mxll_ab_mask)
          ab_type_str = "Mask"
        case (mxll_ab_mask_zero)
          ab_type_str = "Zero"
        case default
          ab_type_str = ""
        end select
 
        if (bc%bc_ab_type(idim) == mxll_ab_cpml .or. bc%bc_ab_type(idim) == mxll_ab_mask .or. &
          bc%bc_ab_type(idim) == mxll_ab_mask_zero) then
          string = trim(ab_type_str)//" Absorbing Boundary"
          write(string,'(2a, f10.3,3a)') trim(string), " in  "//dims(idim)//' direction spans from: ', &
            units_from_atomic(units_inp%length, ab_bounds(1, idim) ), ' [', &
            trim(units_abbrev(units_inp%length)), ']'
          write(message(1),'(a)') trim(string)
 
          string = "to "
          write(string,'(a,f10.3,3a)') trim(string),&
            units_from_atomic(units_inp%length, ab_bounds(2, idim) ), ' [', &
            trim(units_abbrev(units_inp%length)), ']'
 
          write(message(2),'(a)') trim(string)
          call messages_info(2, namespace=namespace)
        end if
 
      class default
 
        write(message(1),'(a,es10.3,3a)') &
          "  Lower bound = ", units_from_atomic(units_inp%length, ab_bounds(1, idim) ),&
          ' [', trim(units_abbrev(units_inp%length)), ']'
        write(message(2),'(a,es10.3,3a)') &
          "  Upper bound = ", units_from_atomic(units_inp%length, ab_bounds(2, idim) ),&
          ' [', trim(units_abbrev(units_inp%length)), ']'
        call messages_info(2, namespace=namespace)
 
      end select
 
    end do
 
    ! initialization of surfaces
    call maxwell_surfaces_init(gr, st, bounds)
 
    ! mapping of mask boundary points
    if (ab_mask_check) then
      call maxwell_mask_points_mapping(bc, gr, ab_bounds)
    end if
 
    ! mapping of pml boundary points
    if (ab_pml_check) then
      call maxwell_pml_points_mapping(bc, gr, ab_bounds)
    end if
 
    ! mapping of constant boundary points
    if (constant_check) then
      call maxwell_constant_points_mapping(bc, gr, bounds)
    end if
 
    ! mapping of plane waves boundary points
    if (plane_waves_check) then
      call maxwell_plane_waves_points_mapping(bc, gr, bounds, namespace)
      call external_waves_init(bc%plane_wave, namespace)
    end if
 
    ! mapping of zero points
    if (zero_check) then
      call maxwell_zero_points_mapping(bc, gr, bounds)
    end if
 
    if (ab_mask_check) then
      call bc_mxll_generate_mask(bc, gr, ab_bounds)
    end if
 
    if (ab_pml_check) then
      call bc_mxll_generate_pml(bc, space)
    end if
 
    !call bc_generate_zero(bc, gr, ab_bounds)
 
    if (debug%info) call bc_mxll_write_info(bc, gr, namespace, space)
 
    if (ab_mask_check .or. ab_pml_check) then
      call messages_print_with_emphasis(namespace=namespace)
    end if
 
    call profiling_out('BC_MXLL_INIT')
 
    pop_sub(bc_mxll_init)
  end subroutine bc_mxll_init
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_end(bc)
    type(bc_mxll_t),   intent(inout) :: bc
 
    integer :: idim
 
    push_sub(bc_mxll_end)
 
    safe_deallocate_a(bc%ab_ufn)
 
    safe_deallocate_a(bc%mask_points_map)
    safe_deallocate_a(bc%mask)
 
    safe_deallocate_a(bc%der_bndry_mask)
    safe_deallocate_a(bc%der_bndry_mask_points_map)
 
    call pml_end(bc%pml)
    do idim = 1, 3
      call single_medium_box_end(bc%medium(idim))
    end do
 
    safe_deallocate_a(bc%constant_points_map)
    safe_deallocate_a(bc%constant_rs_state)
    if (accel_is_enabled()) then
      call accel_release_buffer(bc%buff_constant_points_map)
    end if
 
    safe_deallocate_a(bc%mirror_points_map)
 
    call external_waves_end(bc%plane_wave)
 
    safe_deallocate_a(bc%zero_points_map)
    safe_deallocate_a(bc%zero)
 
    pop_sub(bc_mxll_end)
  end subroutine bc_mxll_end
 
  ! ---------------------------------------------------------
  subroutine pml_end(pml)
    type(pml_t),   intent(inout) :: pml
 
    push_sub(pml_end)
 
    safe_deallocate_a(pml%points_map)
    safe_deallocate_a(pml%points_map_inv)
    safe_deallocate_a(pml%kappa)
    safe_deallocate_a(pml%sigma_e)
    safe_deallocate_a(pml%sigma_m)
    safe_deallocate_a(pml%a)
    safe_deallocate_a(pml%b)
    safe_deallocate_a(pml%c)
    safe_deallocate_a(pml%mask)
    safe_deallocate_a(pml%aux_ep)
    safe_deallocate_a(pml%aux_mu)
    safe_deallocate_a(pml%conv_plus)
    safe_deallocate_a(pml%conv_minus)
    safe_deallocate_a(pml%conv_plus_old)
    safe_deallocate_a(pml%conv_minus_old)
    if (accel_is_enabled()) then
      call accel_release_buffer(pml%buff_map)
      call accel_release_buffer(pml%buff_a)
      call accel_release_buffer(pml%buff_b)
      call accel_release_buffer(pml%buff_c)
      call accel_release_buffer(pml%buff_conv_plus)
      call accel_release_buffer(pml%buff_conv_minus)
      call accel_release_buffer(pml%buff_conv_plus_old)
    end if
 
 
    pop_sub(pml_end)
  end subroutine pml_end
 
 
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
    type(grid_t),        intent(in)    :: gr
    type(namespace_t),   intent(in)    :: namespace
    real(real64),        intent(inout) :: bounds(:,:)
    integer,             intent(in)    :: idim
    real(real64),        intent(out)   :: ep_factor
    real(real64),        intent(out)   :: mu_factor
    real(real64),        intent(out)   :: sigma_e_factor
    real(real64),        intent(out)   :: sigma_m_factor
 
    real(real64) :: width
 
    push_sub(bc_mxll_medium_init)
 
    call profiling_in('BC_MXLL_MEDIUM_INIT')
 
    !%Variable MediumWidth
    !%Type float
    !%Default 0.0 a.u.
    !%Section Maxwell::Boundaries
    !%Description
    !% Width of the boundary region with medium
    !%End
    call parse_variable(namespace, 'MediumWidth', m_zero, width, units_inp%length)
    bounds(1,idim) = ( gr%idx%nr(2, idim) - gr%idx%enlarge(idim) ) * gr%spacing(idim)
    bounds(1,idim) = bounds(1,idim) - width
    bounds(2,idim) = ( gr%idx%nr(2, idim) ) * gr%spacing(idim)
 
    !%Variable MediumEpsilonFactor
    !%Type float
    !%Default 1.0.
    !%Section Maxwell::Boundaries
    !%Description
    !% Linear medium electric susceptibility.
    !%End
    call parse_variable(namespace, 'MediumpsilonFactor', m_one, ep_factor, unit_one)
 
    !%Variable MediumMuFactor
    !%Type float
    !%Default 1.0
    !%Section Maxwell::Boundaries
    !%Description
    !% Linear medium magnetic susceptibility.
    !%End
    call parse_variable(namespace, 'MediumMuFactor', m_one, mu_factor, unit_one)
 
    !%Variable MediumElectricSigma
    !%Type float
    !%Default 0.
    !%Section Maxwell::Boundaries
    !%Description
    !% Electric conductivity of the linear medium.
    !%End
 
    call parse_variable(namespace, 'MediumElectricSigma', m_zero, sigma_e_factor, unit_one)
    !%Variable MediumMagneticSigma
    !%Type float
    !%Default 0.
    !%Section Maxwell::Boundaries
    !%Description
    !% Magnetic conductivity of the linear medium.
    !%End
    call parse_variable(namespace, 'MediumMagneticSigma', m_zero, sigma_m_factor, unit_one)
 
    call profiling_out('BC_MXLL_MEDIUM_INIT')
 
    pop_sub(bc_mxll_medium_init)
  end subroutine bc_mxll_medium_init
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
    type(bc_mxll_t),     intent(inout) :: bc
    type(grid_t),        intent(in)    :: gr
    type(namespace_t),   intent(in)    :: namespace
    real(real64),        intent(inout) :: ab_bounds(:,:)
    integer,             intent(in)    :: idim
 
    real(real64)        :: default_width
 
    push_sub(bc_mxll_ab_bounds_init)
 
    !%Variable MaxwellABWidth
    !%Type float
    !%Section Maxwell::Boundaries
    !%Description
    !% Width of the region used to apply the absorbing boundaries. The default value is twice
    !% the derivative order.
    !%End
 
    default_width = m_two * gr%der%order * maxval(gr%spacing(1:3))
    call parse_variable(namespace, 'MaxwellABWidth', default_width, bc%ab_width, units_inp%length)
 
    if (bc%ab_width < default_width) then
      bc%ab_width = default_width
      write(message(1),'(a)') 'Absorbing boundary width has to be larger or equal than twice the derivatives order times spacing!'
      write(message(2),'(a,es10.3)') 'Absorbing boundary width is set to: ', default_width
      call messages_info(2, namespace=namespace)
    end if
 
    ab_bounds(1, idim) = ab_bounds(2, idim) - bc%ab_width
 
    pop_sub(bc_mxll_ab_bounds_init)
  end subroutine bc_mxll_ab_bounds_init
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_pml_init(bc, gr, namespace, ab_bounds, idim)
    type(bc_mxll_t),     intent(inout) :: bc
    type(grid_t),        intent(in)    :: gr
    type(namespace_t),   intent(in)    :: namespace
    real(real64),        intent(inout) :: ab_bounds(:,:)
    integer,             intent(in)    :: idim
 
    push_sub(bc_mxll_pml_init)
 
    call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
    bc%pml%width = bc%ab_width
 
    !%Variable MaxwellABPMLPower
    !%Type float
    !%Default 3.5
    !%Section Maxwell::Boundaries
    !%Description
    !% Exponential of the polynomial profile for the non-physical conductivity of the PML.
    !% Should be between 2 and 4
    !%End
    call parse_variable(namespace, 'MaxwellABPMLPower', 3.5_real64, bc%pml%power, unit_one)
 
    !%Variable MaxwellABPMLReflectionError
    !%Type float
    !%Default 1.0e-16
    !%Section Maxwell::Boundaries
    !%Description
    !% Tolerated reflection error for the PML
    !%End
    call parse_variable(namespace, 'MaxwellABPMLReflectionError', 1.0e-16_real64, bc%pml%refl_error, unit_one)
 
    pop_sub(bc_mxll_pml_init)
  end subroutine bc_mxll_pml_init
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_write_info(bc, mesh, namespace, space)
    type(bc_mxll_t),       intent(in) :: bc
    class(mesh_t),         intent(in) :: mesh
    type(namespace_t),     intent(in) :: namespace
    class(space_t),        intent(in) :: space
 
    integer :: err, idim, idim2
    real(real64), allocatable :: tmp(:)
    logical :: mask_check, pml_check, medium_check
    character(1) :: dim_label(3)
 
    push_sub(bc_mxll_write_info)
 
    call profiling_in('BC_MXLL_WRITE_INFO')
 
    mask_check = .false.
    pml_check = .false.
    medium_check = .false.
 
    dim_label = (/'x', 'y', 'z'/)
 
    do idim = 1, 3
      if (bc%bc_ab_type(idim) == mxll_ab_mask) then
        mask_check = .true.
      end if
      if (bc%bc_ab_type(idim) == mxll_ab_cpml) then
        pml_check = .true.
      end if
      if (bc%bc_ab_type(idim) == mxll_bc_medium) then
        medium_check = .true.
      end if
    end do
 
    if (mask_check) then
      safe_allocate(tmp(1:mesh%np))
      do idim = 1, 3
        tmp(:) = m_zero
        call get_mask_io_function(bc%mask, bc, tmp, idim)
        call write_files("maxwell_mask", tmp)
      end do
      safe_deallocate_a(tmp)
    else if (pml_check) then
      safe_allocate(tmp(1:mesh%np))
      do idim = 1, 3
        ! sigma for electric field dim = idim
        tmp(:) = m_one
        call get_pml_io_function(bc%pml%sigma_e(:, idim), bc, tmp)
        call write_files("maxwell_sigma_e-"//dim_label(idim), tmp)
 
        ! sigma for magnetic dim = idim
        tmp(:) = m_zero
        call get_pml_io_function(bc%pml%sigma_m(:, 1), bc, tmp)
        call write_files("maxwell_sigma_m-"//dim_label(idim), tmp)
 
        ! pml_a for electric field dim = idim
        tmp(:) = m_zero
        call get_pml_io_function(real(bc%pml%a(:, idim), real64), bc, tmp)
        call write_files("maxwell_sigma_pml_a_e-"//dim_label(idim), tmp)
      end do
      safe_deallocate_a(tmp)
    end if
 
    if (medium_check) then
      safe_allocate(tmp(1:mesh%np))
      do idim = 1, 3
        ! medium epsilon
        tmp(:) = p_ep
        call get_medium_io_function(bc%medium(idim)%ep, bc, tmp, idim)
        call write_files("maxwell_ep"//dim_label(idim), tmp)
        ! medium mu
        tmp(:) = p_mu
        call get_medium_io_function(bc%medium(idim)%mu, bc, tmp, idim)
        call write_files("maxwell_mu"//dim_label(idim), tmp)
        ! medium epsilon
        tmp(:) = p_c
        call get_medium_io_function(bc%medium(idim)%c, bc, tmp, idim)
        call write_files("maxwell_c"//dim_label(idim), tmp)
        do idim2 = 1, 3
          ! medium epsilon aux field dim = idim
          tmp(:) = m_zero
          call get_medium_io_function(bc%medium(idim)%aux_ep(:, idim2), bc, tmp, idim)
          call write_files("maxwell_aux_ep-"//dim_label(idim)//"-"//dim_label(idim2), tmp)
 
          ! medium mu aux field dim = idim
          tmp(:) = m_zero
          call get_medium_io_function(bc%medium(idim)%aux_mu(:, idim2), bc, tmp, idim)
          call write_files("maxwell_aux_mu-"//dim_label(idim)//"-"//dim_label(idim2), tmp)
        end do
      end do
 
      safe_deallocate_a(tmp)
    end if
 
    call profiling_out('BC_MXLL_WRITE_INFO')
 
    pop_sub(bc_mxll_write_info)
  contains
 
    subroutine get_pml_io_function(pml_func, bc, io_func)
      real(real64),       intent(in)    :: pml_func(:)
      type(bc_mxll_t),    intent(in)    :: bc
      real(real64),       intent(inout) :: io_func(:)
 
      integer :: ip, ip_in
 
      do ip_in = 1, bc%pml%points_number
        ip          = bc%pml%points_map(ip_in)
        io_func(ip) = pml_func(ip_in)
      end do
 
    end subroutine get_pml_io_function
 
    subroutine get_mask_io_function(mask_func, bc, io_func, idim)
      real(real64),       intent(in)    :: mask_func(:,:)
      type(bc_mxll_t),    intent(in)    :: bc
      real(real64),       intent(inout) :: io_func(:)
      integer,            intent(in)    :: idim
 
      integer :: ip, ip_in
 
      do ip_in = 1, bc%mask_points_number(idim)
        ip          = bc%mask_points_map(ip_in, idim)
        io_func(ip) = mask_func(ip_in, idim)
      end do
 
    end subroutine get_mask_io_function
 
    subroutine get_medium_io_function(medium_func, bc, io_func, idim)
      real(real64),       intent(in)    :: medium_func(:)
      type(bc_mxll_t),    intent(in)    :: bc
      real(real64),       intent(inout) :: io_func(:)
      integer,            intent(in)    :: idim
 
      integer :: ip, ip_in
 
      do ip_in = 1, bc%medium(idim)%points_number
        ip          = bc%medium(idim)%points_map(ip_in)
        io_func(ip) = medium_func(ip_in)
      end do
 
    end subroutine get_medium_io_function
 
    subroutine write_files(filename, tmp)
      character(len=*), intent(in) :: filename
      real(real64),     intent(in) :: tmp(:)
 
      call dio_function_output(io_function_fill_how("VTK"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("AxisX"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("AxisY"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("AxisZ"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneX"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneY"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneZ"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
        unit_one, err)
    end subroutine write_files
 
  end subroutine bc_mxll_write_info
 
  ! ---------------------------------------------------------
  subroutine maxwell_mask_points_mapping(bc, mesh, bounds)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, ip_in_max, point_info, idim
 
    push_sub(maxwell_mask_points_mapping)
 
    call profiling_in('MAXWELL_MASK_POINTS_MAPPING')
    ip_in_max = 1
    do idim = 1, 3
      if (bc%bc_ab_type(idim) == mxll_ab_mask) then
        ! allocate mask points map
        ip_in = 0
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
          end if
        end do
        if (ip_in > ip_in_max) ip_in_max = ip_in
        bc%mask_points_number(idim) = ip_in
      end if
    end do
    safe_allocate(bc%mask(1:ip_in_max, 1:idim))
    safe_allocate(bc%mask_points_map(1:ip_in_max, 1:idim))
 
    do idim = 1, 3
      if (bc%bc_ab_type(idim) == mxll_ab_mask) then
        ! mask points mapping
        ip_in = 0
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
            bc%mask_points_map(ip_in, idim) = ip
          end if
        end do
      end if
    end do
 
    call profiling_out('MAXWELL_MASK_POINTS_MAPPING')
    pop_sub(maxwell_mask_points_mapping)
  end subroutine maxwell_mask_points_mapping
 
  ! ---------------------------------------------------------
  subroutine maxwell_pml_points_mapping(bc, mesh, bounds)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, point_info
 
    push_sub(maxwell_pml_points_mapping)
 
    call profiling_in('MXWL_PML_POINTS_MAPPING')
 
    ! allocate pml points map
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
      end if
    end do
    bc%pml%points_number = ip_in
    safe_allocate(bc%pml%points_map(1:ip_in))
    safe_allocate(bc%pml%points_map_inv(1:mesh%np))
    bc%pml%points_map_inv = 0
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
        bc%pml%points_map(ip_in) = ip
        bc%pml%points_map_inv(ip) = ip_in
      end if
    end do
 
    call profiling_out('MXWL_PML_POINTS_MAPPING')
 
    pop_sub(maxwell_pml_points_mapping)
  end subroutine maxwell_pml_points_mapping
 
  ! ---------------------------------------------------------
  subroutine maxwell_constant_points_mapping(bc, mesh, bounds)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, point_info
 
    push_sub(maxwell_constant_points_mapping)
 
    call profiling_in('MAXWELL_CONSTANT_POINTS_MAP')
 
    ! allocate constant points map
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
      end if
    end do
    bc%constant_points_number = ip_in
    safe_allocate(bc%constant_points_map(1:ip_in))
    safe_allocate(bc%constant_rs_state(1:ip_in, 1:3))
 
    ! zero constant mapping
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
        bc%constant_points_map(ip_in) = ip
      end if
    end do
 
    if (accel_is_enabled()) then
      call accel_create_buffer(bc%buff_constant_points_map, accel_mem_read_only, type_integer, bc%constant_points_number)
      call accel_write_buffer(bc%buff_constant_points_map, bc%constant_points_number, bc%constant_points_map)
    end if
 
    call profiling_out('MAXWELL_CONSTANT_POINTS_MAP')
 
    pop_sub(maxwell_constant_points_mapping)
  end subroutine maxwell_constant_points_mapping
 
  ! ---------------------------------------------------------
  subroutine maxwell_plane_waves_points_mapping(bc, mesh, bounds, namespace)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
    type(namespace_t),   intent(in)    :: namespace
 
    integer :: ip, ip_in, point_info
 
    push_sub(maxwell_plane_waves_points_mapping)
 
    call profiling_in('MXLL_PW_POINTS_MAP')
 
    bc%plane_waves_dims = (bc%bc_type(1:3) == mxll_bc_plane_waves)
 
    ! allocate map
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
      end if
    end do
    bc%plane_wave%points_number = ip_in
    safe_allocate(bc%plane_wave%points_map(1:ip_in))
 
    ! mapping
    ip_in = 0
    do ip = 1, mesh%np
      call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
      if (point_info == 1) then
        ip_in = ip_in + 1
        bc%plane_wave%points_map(ip_in) = ip
      end if
    end do
 
    if (accel_is_enabled()) then
      call accel_create_buffer(bc%plane_wave%buff_map, accel_mem_read_only, type_integer, bc%plane_wave%points_number)
      call accel_write_buffer(bc%plane_wave%buff_map, bc%plane_wave%points_number, bc%plane_wave%points_map)
    end if
 
    call profiling_out('MXLL_PW_POINTS_MAP')
 
    pop_sub(maxwell_plane_waves_points_mapping)
  end subroutine maxwell_plane_waves_points_mapping
 
  ! ---------------------------------------------------------
  subroutine maxwell_zero_points_mapping(bc, mesh, bounds)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, ip_in_max, point_info, idim
 
    push_sub(maxwell_zero_points_mapping)
 
    call profiling_in('MXLL_ZERO_POINTS_MAPPING')
 
    ip_in_max = 0
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_zero) then
        ! allocate zero points map
        ip_in = 0
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
          end if
        end do
        if (ip_in > ip_in_max) ip_in_max = ip_in
      end if
    end do
    bc%zero_points_number = ip_in
    safe_allocate(bc%zero(1:ip_in_max,3))
    safe_allocate(bc%zero_points_map(1:ip_in_max, 1:3))
 
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_zero) then
        ! zero points mapping
        ip_in = 0
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
            bc%zero_points_map(ip_in, idim) = ip
          end if
        end do
      end if
    end do
 
    call profiling_out('MXLL_ZERO_POINTS_MAPPING')
 
    pop_sub(maxwell_zero_points_mapping)
  end subroutine maxwell_zero_points_mapping
 
  ! ---------------------------------------------------------
  subroutine maxwell_medium_points_mapping(bc, mesh, bounds)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, ip_in_max, ip_bd, ip_bd_max, point_info, boundary_info, idim
 
    push_sub(maxwell_medium_points_mapping)
 
    call profiling_in('MXLL_MEDIUM_POINTS_MAPPING')
 
    ip_in_max = 0
    ip_bd_max = 0
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_medium) then
        ! allocate pml points map
        ip_in = 0
        ip_bd = 0
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          call maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
          end if
          if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_bd = ip_bd + 1
          end if
        end do
        bc%medium(idim)%points_number = ip_in
        bc%medium(idim)%bdry_number = ip_bd
      end if
    end do
    do idim = 1, 3
      safe_allocate(bc%medium(idim)%points_map(1:ip_in_max))
      safe_allocate(bc%medium(idim)%bdry_map(1:ip_bd_max))
    end do
 
    ip_in = 0
    ip_bd = 0
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_medium) then
        do ip = 1, mesh%np
          call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
          call maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
          if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_in = ip_in + 1
            bc%medium(idim)%points_map(ip_in) = ip
          end if
          if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
            ip_bd = ip_bd + 1
            bc%medium(idim)%bdry_map(ip_bd) = ip
          end if
        end do
      end if
    end do
 
    call profiling_out('MXLL_MEDIUM_POINTS_MAPPING')
 
    pop_sub(maxwell_medium_points_mapping)
  end subroutine maxwell_medium_points_mapping
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_generate_pml(bc, space)
    type(bc_mxll_t),    intent(inout) :: bc
    class(space_t),     intent(in)    :: space
 
 
    push_sub(bc_mxll_generate_pml)
 
    call profiling_in('BC_MXLL_GENERATE_PML')
 
    safe_allocate(bc%pml%kappa(1:bc%pml%points_number, 1:space%dim))
    safe_allocate(bc%pml%sigma_e(1:bc%pml%points_number, 1:space%dim))
    safe_allocate(bc%pml%sigma_m(1:bc%pml%points_number, 1:space%dim))
    safe_allocate(bc%pml%a(1:bc%pml%points_number, 1:space%dim))
    safe_allocate(bc%pml%b(1:bc%pml%points_number, 1:space%dim))
    safe_allocate(bc%pml%c(1:bc%pml%points_number, 1:3))
    safe_allocate(bc%pml%mask(1:bc%pml%points_number))
    safe_allocate(bc%pml%conv_plus(1:bc%pml%points_number, 1:3, 1:3))
    safe_allocate(bc%pml%conv_minus(1:bc%pml%points_number, 1:3, 1:3))
    safe_allocate(bc%pml%conv_plus_old(1:bc%pml%points_number, 1:3, 1:3))
    safe_allocate(bc%pml%conv_minus_old(1:bc%pml%points_number, 1:3, 1:3))
    safe_allocate(bc%pml%aux_ep(1:bc%pml%points_number, 1:3, 1:3))
    safe_allocate(bc%pml%aux_mu(1:bc%pml%points_number, 1:3, 1:3))
 
    bc%pml%kappa                 = m_one
    bc%pml%sigma_e               = m_zero
    bc%pml%sigma_m               = m_zero
    bc%pml%a                     = m_zero
    bc%pml%b                     = m_zero
    bc%pml%c                     = m_zero
    bc%pml%mask                  = m_one
    bc%pml%conv_plus             = m_z0
    bc%pml%conv_minus            = m_z0
    bc%pml%conv_plus_old         = m_z0
    bc%pml%conv_minus_old        = m_z0
 
    if (accel_is_enabled()) then
      call accel_create_buffer(bc%pml%buff_map, accel_mem_read_only, type_integer, bc%pml%points_number)
      call accel_create_buffer(bc%pml%buff_a, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
      call accel_create_buffer(bc%pml%buff_b, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
      call accel_create_buffer(bc%pml%buff_c, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
      call accel_create_buffer(bc%pml%buff_conv_plus, accel_mem_read_write, type_cmplx, &
        int(bc%pml%points_number, int64)*space%dim*space%dim)
      call accel_create_buffer(bc%pml%buff_conv_minus, accel_mem_read_write, type_cmplx, &
        int(bc%pml%points_number, int64)*space%dim*space%dim)
      call accel_create_buffer(bc%pml%buff_conv_plus_old, accel_mem_read_write, type_cmplx, &
        int(bc%pml%points_number, int64)*space%dim*space%dim)
 
      call accel_write_buffer(bc%pml%buff_a, int(bc%pml%points_number, int64)*space%dim, bc%pml%a)
      call accel_write_buffer(bc%pml%buff_b, int(bc%pml%points_number, int64)*space%dim, bc%pml%b)
      call accel_write_buffer(bc%pml%buff_c, int(bc%pml%points_number, int64)*space%dim, bc%pml%c)
      call accel_write_buffer(bc%pml%buff_conv_plus, int(bc%pml%points_number, int64)*space%dim*space%dim, bc%pml%conv_plus)
      call accel_write_buffer(bc%pml%buff_conv_minus, int(bc%pml%points_number, int64)*space%dim*space%dim, bc%pml%conv_minus)
      call accel_write_buffer(bc%pml%buff_conv_plus_old, int(bc%pml%points_number, int64)*space%dim*space%dim, bc%pml%conv_plus_old)
    end if
 
    call profiling_out('BC_MXLL_GENERATE_PML')
 
    pop_sub(bc_mxll_generate_pml)
  end subroutine bc_mxll_generate_pml
 
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_generate_pml_parameters(bc, space, gr, c_factor, dt)
    type(bc_mxll_t),    intent(inout) :: bc
    class(space_t),     intent(in)    :: space
    type(grid_t),       intent(in)    :: gr
    real(real64),       intent(in)    :: c_factor
    real(real64), optional,    intent(in)    :: dt
 
    integer :: ip, ip_in, idim
    real(real64)   :: width(3)
    real(real64)   :: ddv(3), ss_e, ss_m, ss_max, aa_e, aa_m, bb_e, bb_m, gg, hh, kk, ll_e, ll_m
    real(real64), allocatable  :: tmp(:), tmp_grad(:,:)
 
 
    push_sub(bc_mxll_generate_pml_parameters)
    safe_allocate(tmp(gr%np_part))
    safe_allocate(tmp_grad(gr%np, 1:space%dim))
    ! here bounds are ab_bounds, for which ab_bounds(1,:) = bc%bounds(1,:)
    ! width is stored in bc%pml%width
    ! assuming that the width is the same in the 3 dimensions (only case implemented now), we can change the following line:
    ! width(1:3) = bounds(2, 1:3) - bounds(1, 1:3) !! original line
    ! by
    width(1:3) = (/ bc%pml%width, bc%pml%width, bc%pml%width /)
 
    ! PML variables for all boundary points
    do ip_in = 1, bc%pml%points_number
      ip = bc%pml%points_map(ip_in)
      ddv(1:3) = abs(gr%x(ip, 1:3)) - bc%bc_bounds(1, 1:3)
      do idim = 1, space%dim
        if (ddv(idim) >= m_zero) then
          gg     = (ddv(idim)/bc%pml%width)**bc%pml%power
          hh     = (m_one-ddv(idim)/bc%pml%width)**bc%pml%power
          kk     = m_one
          ss_max = -(bc%pml%power + m_one)*p_c*c_factor*p_ep*log(bc%pml%refl_error)/(m_two * bc%pml%width)
          ss_e   = gg * ss_max
          ss_m   = gg * ss_max
          ll_e   = ss_e*kk
          ll_m   = ss_m*kk
          bb_e   = exp(-(ss_e/(p_ep))*dt)
          bb_m   = exp(-(ss_m/(p_ep))*dt)
          aa_e   = (bb_e - 1)
          aa_m   = (bb_m - 1)
          if (abs(ll_e) < m_epsilon) aa_e = m_zero
          if (abs(ll_m) < m_epsilon) aa_m = m_zero
          bc%pml%sigma_e(ip_in, idim) = ss_e
          bc%pml%sigma_m(ip_in, idim) = ss_m
          bc%pml%a(ip_in, idim)       = aa_e
          bc%pml%b(ip_in, idim)       = bb_e
          bc%pml%kappa(ip_in, idim)   = kk
          bc%pml%mask(ip_in)          = bc%pml%mask(ip_in) * (m_one - sin(ddv(idim)*m_pi/(m_two*(width(idim))))**2)
        else
          bc%pml%kappa(ip_in, idim)   = m_one
          bc%pml%sigma_e(ip_in, idim) = m_zero
          bc%pml%sigma_m(ip_in, idim) = m_zero
          bc%pml%a(ip_in, idim)       = m_zero
          bc%pml%b(ip_in, idim)       = m_zero
          bc%pml%mask(ip_in)          = m_one
        end if
      end do
    end do
 
    ! PML auxiliary epsilon for all boundary points
    do idim = 1, space%dim
      tmp = p_ep
      do ip_in = 1, bc%pml%points_number
        ip = bc%pml%points_map(ip_in)
        tmp(ip) = p_ep / bc%pml%kappa(ip_in, idim)
      end do
      call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
      do ip_in = 1, bc%pml%points_number
        ip = bc%pml%points_map(ip_in)
        bc%pml%aux_ep(ip_in, :, idim) = tmp_grad(ip, :)/(m_four*p_ep*bc%pml%kappa(ip_in, idim))
      end do
    end do
 
    ! PML auxiliary mu
    do idim = 1, space%dim
      ! P_mu is proportional to 1/P_c**2, so we need to devide by c_factor
      tmp = p_mu/c_factor**2
      do ip_in = 1, bc%pml%points_number
        ip = bc%pml%points_map(ip_in)
        tmp(ip) = p_mu/c_factor**2 / bc%pml%kappa(ip_in, idim)
      end do
      call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
      do ip_in = 1, bc%pml%points_number
        ip = bc%pml%points_map(ip_in)
        bc%pml%aux_mu(ip_in, :, idim) = tmp_grad(ip, :)/(m_four*p_mu/c_factor**2 * bc%pml%kappa(ip_in, idim))
      end do
    end do
 
    ! PML auxiliary c for all boundary points
    do idim = 1, space%dim
      do ip_in = 1, bc%pml%points_number
        bc%pml%c(ip_in, idim) = p_c*c_factor/bc%pml%kappa(ip_in, idim)
      end do
    end do
    safe_deallocate_a(tmp)
    safe_deallocate_a(tmp_grad)
 
    if (accel_is_enabled()) then
      call accel_write_buffer(bc%pml%buff_map, bc%pml%points_number, bc%pml%points_map)
      call accel_write_buffer(bc%pml%buff_a, int(bc%pml%points_number, int64)*space%dim, bc%pml%a)
      call accel_write_buffer(bc%pml%buff_b, int(bc%pml%points_number, int64)*space%dim, bc%pml%b)
      call accel_write_buffer(bc%pml%buff_c, int(bc%pml%points_number, int64)*space%dim, bc%pml%c)
      call accel_write_buffer(bc%pml%buff_conv_plus, int(bc%pml%points_number, int64)*space%dim*space%dim, bc%pml%conv_plus)
      call accel_write_buffer(bc%pml%buff_conv_minus, int(bc%pml%points_number, int64)*space%dim*space%dim, bc%pml%conv_minus)
    end if
 
    bc%pml%parameters_initialized = .true.
 
    pop_sub(bc_mxll_generate_pml_parameters)
 
  end subroutine bc_mxll_generate_pml_parameters
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_initialize_pml_simple(bc, space, gr, c_factor, dt)
    type(bc_mxll_t),    intent(inout) :: bc
    class(space_t),     intent(in)    :: space
    type(grid_t),       intent(in)    :: gr
    real(real64),       intent(in)    :: c_factor
    real(real64),       intent(in)    :: dt
 
    integer :: ip_in, ip, idir
    real(real64) :: ddv(1:space%dim), kappa, sigma, alpha
    push_sub(bc_mxll_initialize_pml_simple)
 
    ! PML variables for all boundary points
    do ip_in = 1, bc%pml%points_number
      ip = bc%pml%points_map(ip_in)
      ddv(1:3) = abs(gr%x(ip, 1:3)) - bc%bc_bounds(1, 1:3)
      do idir = 1, space%dim
        if (ddv(idir) >= m_zero) then
          sigma = (ddv(idir)/bc%pml%width)**bc%pml%power * &
            -(bc%pml%power + m_one)*p_c*c_factor*p_ep*log(bc%pml%refl_error)/(m_two * bc%pml%width)
          ! kappa and alpha could be set to different values, but these seem to be fine
          kappa = m_one
          ! non-zero values for alpha can be used to modify the low-frequency behavior
          alpha = m_zero
          bc%pml%b(ip_in, idir) = exp(-(alpha + sigma/kappa)/p_ep*dt)
          if (abs(sigma) < m_epsilon) then
            bc%pml%c(ip_in, idir) = m_zero
          else
            bc%pml%c(ip_in, idir) = m_one/(kappa + kappa**2*alpha/sigma) * &
              (bc%pml%b(ip_in, idir) - 1)
          end if
        else
          bc%pml%b(ip_in, idir) = m_zero
          bc%pml%c(ip_in, idir) = m_zero
        end if
      end do
    end do
    if (accel_is_enabled()) then
      call accel_write_buffer(bc%pml%buff_map, &
        bc%pml%points_number, bc%pml%points_map)
      call accel_write_buffer(bc%pml%buff_b, &
        int(bc%pml%points_number, int64)*space%dim, bc%pml%b)
      call accel_write_buffer(bc%pml%buff_c, &
        int(bc%pml%points_number, int64)*space%dim, bc%pml%c)
    end if
    pop_sub(bc_mxll_initialize_pml_simple)
  end subroutine bc_mxll_initialize_pml_simple
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_generate_mask(bc, mesh, bounds)
    type(bc_mxll_t),    intent(inout) :: bc
    class(mesh_t),      intent(in)    :: mesh
    real(real64),       intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, idim, ip_in_max
    real(real64)   :: ddv(3), tmp(3), width(3)
    real(real64), allocatable :: mask(:)
 
    push_sub(bc_mxll_generate_mask)
 
    call profiling_in('BC_MXLL_GENERATE_MASK')
 
    ip_in_max = maxval(bc%mask_points_number(:))
 
    safe_allocate(mask(1:mesh%np))
 
    mask(:) = m_one
 
    width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
    tmp(:)   = m_zero
 
    do ip = 1, mesh%np
      tmp = m_one
      mask(ip) = m_one
      ddv(1:3) = abs(mesh%x(ip, 1:3)) - bounds(1, 1:3)
      do idim = 1, mesh%box%dim
        if (ddv(idim) >= m_zero) then
          if (ddv(idim)  <=  width(idim)) then
            tmp(idim) = m_one - sin(ddv(idim) * m_pi / (m_two * (width(idim))))**2
          else
            tmp(idim) = m_one
          end if
        end if
        mask(ip) = mask(ip) * tmp(idim)
      end do
    end do
 
    do idim = 1, mesh%box%dim
      do ip_in = 1, bc%mask_points_number(idim)
        ip = bc%mask_points_map(ip_in, idim)
        bc%mask(ip_in,idim) = mask(ip)
      end do
    end do
 
    safe_deallocate_a(mask)
 
    call profiling_out('BC_MXLL_GENERATE_MASK')
 
    pop_sub(bc_mxll_generate_mask)
  end subroutine bc_mxll_generate_mask
 
  ! ---------------------------------------------------------
  subroutine bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
    type(bc_mxll_t),         intent(inout) :: bc
    class(space_t),          intent(in)    :: space
    type(grid_t),            intent(in)    :: gr
    real(real64),            intent(in)    :: bounds(:,:)
    real(real64),            intent(in)    :: ep_factor
    real(real64),            intent(in)    :: mu_factor
    real(real64),            intent(in)    :: sigma_e_factor
    real(real64),            intent(in)    :: sigma_m_factor
 
    integer :: ip, ipp, ip_in, ip_in_max, ip_bd, idim, point_info
    real(real64)   :: dd, dd_min, dd_max, xx(3), xxp(3)
    real(real64), allocatable  :: tmp(:), tmp_grad(:,:)
 
    push_sub(bc_mxll_generate_medium)
 
    call profiling_in('BC_MXLL_GENERATE_MEDIUM')
 
    ip_in_max = max(bc%medium(1)%points_number, bc%medium(2)%points_number, bc%medium(3)%points_number)
    dd_max = max(2*gr%spacing(1), 2*gr%spacing(2), 2*gr%spacing(3))
 
    do idim = 1, 3
      call single_medium_box_allocate(bc%medium(idim), ip_in_max)
      safe_allocate(tmp(gr%np_part))
      safe_allocate(tmp_grad(1:gr%np_part,1:space%dim))
      bc%medium(idim)%aux_ep = m_zero
      bc%medium(idim)%aux_mu = m_zero
      bc%medium(idim)%c = p_c
 
      tmp = p_ep
      do  ip = 1, gr%np_part
        call maxwell_box_point_info(bc, gr, ip, bounds, point_info)
        if ((point_info /= 0) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim))) then
          xx(:) = gr%x(ip, :)
          dd_min = m_huge
          do ip_bd = 1, bc%medium(idim)%bdry_number
            ipp = bc%medium(idim)%bdry_map(ip_bd)
            xxp(:) = gr%x(ipp, :)
            dd = norm2(xx(1:3) - xxp(1:3))
            if (dd < dd_min) dd_min = dd
          end do
          tmp(ip) = p_ep * (m_one + ep_factor * m_one/(m_one + exp(-m_five/dd_max * (dd_min-m_two*dd_max))))
        end if
      end do
      call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
      do ip_in = 1, bc%medium(idim)%points_number
        ip = bc%medium(idim)%points_map(ip_in)
        bc%medium(idim)%aux_ep(ip_in, :) = &
          tmp_grad(ip, :)/(m_four*p_ep*ep_factor * m_one/(m_one + exp(-m_five/dd_max-dd)))
      end do
    end do
 
    do idim = 1, 3
      tmp = p_mu
      do ip = 1, gr%np_part
        call maxwell_box_point_info(bc, gr, ip, bounds, point_info)
        if ((point_info == 1) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim))) then
          xx(:) = gr%x(ip, :)
          dd_min = m_huge
          do ip_bd = 1, bc%medium(idim)%bdry_number
            ipp = bc%medium(idim)%bdry_map(ip_bd)
            xxp(:) = gr%x(ipp,:)
            dd = norm2(xx(1:3) - xxp(1:3))
            if (dd < dd_min) dd_min = dd
          end do
          tmp(ip) = p_mu * (m_one + mu_factor * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
        end if
      end do
      call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
      do ip_in = 1, bc%medium(idim)%points_number
        ip = bc%medium(idim)%points_map(ip_in)
        bc%medium(idim)%aux_mu(ip_in, :) = &
          tmp_grad(ip, :)/(m_four * p_mu * mu_factor * m_one/(m_one + exp(-m_five/dd_max-dd)))
      end do
    end do
 
    do idim = 1, 3
      do ip_in = 1, bc%medium(idim)%points_number
        ip = bc%medium(idim)%points_map(ip_in)
        xx(:) = gr%x(ip, :)
        dd_min = m_huge
        do ip_bd = 1, bc%medium(idim)%bdry_number
          ipp = bc%medium(idim)%bdry_map(ip_bd)
          xxp(:) = gr%x(ipp, :)
          dd = norm2(xx(1:3) - xxp(1:3))
          if (dd < dd_min) dd_min = dd
        end do
        bc%medium(idim)%ep(ip_in) = p_ep * (m_one + ep_factor &
          * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
        bc%medium(idim)%mu(ip_in) = p_mu * (m_one + mu_factor &
          * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
        bc%medium(idim)%sigma_e(ip_in) = (m_one + sigma_e_factor &
          * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
        bc%medium(idim)%sigma_m(ip_in) = (m_one + sigma_m_factor &
          * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
        bc%medium(idim)%c(ip_in) = m_one / sqrt(bc%medium(idim)%ep(ip_in) * bc%medium(idim)%mu(ip_in))
      end do
    end do
 
    safe_deallocate_a(tmp)
    safe_deallocate_a(tmp_grad)
 
    call profiling_out('BC_MXLL_GENERATE_MEDIUM')
 
    pop_sub(bc_mxll_generate_medium)
  end subroutine bc_mxll_generate_medium
 
  ! ---------------------------------------------------------
  subroutine maxwell_surfaces_init(mesh, st, bounds)
    class(mesh_t),            intent(in)    :: mesh
    type(states_mxll_t),      intent(inout) :: st
    real(real64),             intent(in)    :: bounds(:,:)
 
 
    push_sub(maxwell_surfaces_init)
 
    call profiling_in('MAXWELL_SURFACES_INIT')
 
    ! y-z surface at -x boundary
    st%surface(1, 1)%spacing   = m_half*(mesh%spacing(2) + mesh%spacing(3))
    st%surface(1, 1)%origin(:) = m_zero
    st%surface(1, 1)%origin(1) = -bounds(1, 1)
    st%surface(1, 1)%n(:) = m_zero
    st%surface(1, 1)%n(1) = -m_one
    st%surface(1, 1)%u(:) = m_zero
    st%surface(1, 1)%u(2) = -m_one
    st%surface(1, 1)%v(:) = m_zero
    st%surface(1, 1)%v(3) = m_one
    st%surface(1, 1)%nu   = -int(bounds(1, 2)/mesh%spacing(2))
    st%surface(1, 1)%mu   =  int(bounds(1, 2)/mesh%spacing(2))
    st%surface(1, 1)%nv   = -int(bounds(1, 3)/mesh%spacing(3))
    st%surface(1, 1)%mv   =  int(bounds(1, 3)/mesh%spacing(3))
 
    ! y-z surface at +x boundary
    st%surface(2, 1)%spacing   = m_half*(mesh%spacing(2) + mesh%spacing(3))
    st%surface(2, 1)%origin(:) = m_zero
    st%surface(2, 1)%origin(1) = bounds(1, 1)
    st%surface(2, 1)%n(:) = m_zero
    st%surface(2, 1)%n(1) = m_one
    st%surface(2, 1)%u(:) = m_zero
    st%surface(2, 1)%u(2) = m_one
    st%surface(2, 1)%v(:) = m_zero
    st%surface(2, 1)%v(3) = m_one
    st%surface(2, 1)%nu   = -int(bounds(1, 2)/mesh%spacing(2))
    st%surface(2, 1)%mu   =  int(bounds(1, 2)/mesh%spacing(2))
    st%surface(2, 1)%nv   = -int(bounds(1, 3)/mesh%spacing(3))
    st%surface(2, 1)%mv   =  int(bounds(1, 3)/mesh%spacing(3))
 
    ! x-z surface at -y boundary
    st%surface(1, 2)%spacing   = m_half*(mesh%spacing(1) + mesh%spacing(3))
    st%surface(1, 2)%origin(:) = m_zero
    st%surface(1, 2)%origin(2) = -bounds(1, 2)
    st%surface(1, 2)%n(:) = m_zero
    st%surface(1, 2)%n(2) = -m_one
    st%surface(1, 2)%u(:) = m_zero
    st%surface(1, 2)%u(1) = m_one
    st%surface(1, 2)%v(:) = m_zero
    st%surface(1, 2)%v(3) = m_one
    st%surface(1, 2)%nu   = -int(bounds(1, 1)/mesh%spacing(1))
    st%surface(1, 2)%mu   =  int(bounds(1, 1)/mesh%spacing(1))
    st%surface(1, 2)%nv   = -int(bounds(1, 3)/mesh%spacing(3))
    st%surface(1, 2)%mv   =  int(bounds(1, 3)/mesh%spacing(3))
 
    ! x-z surface at +y boundary
    st%surface(2, 2)%spacing   = m_half*(mesh%spacing(1) + mesh%spacing(3))
    st%surface(2, 2)%origin(:) = m_zero
    st%surface(2, 2)%origin(2) = bounds(1, 2)
    st%surface(2, 2)%n(:) = m_zero
    st%surface(2, 2)%n(2) = m_one
    st%surface(2, 2)%u(:) = m_zero
    st%surface(2, 2)%u(1) = m_one
    st%surface(2, 2)%v(:) = m_zero
    st%surface(2, 2)%v(3) = -m_one
    st%surface(2, 2)%nu   = -int(bounds(1, 1)/mesh%spacing(1))
    st%surface(2, 2)%mu   =  int(bounds(1, 1)/mesh%spacing(1))
    st%surface(2, 2)%nv   = -int(bounds(1, 3)/mesh%spacing(3))
    st%surface(2, 2)%mv   =  int(bounds(1, 3)/mesh%spacing(3))
 
    ! x-y surface at -z boundary
    st%surface(1, 3)%spacing   = m_half*(mesh%spacing(1) + mesh%spacing(2))
    st%surface(1, 3)%origin(:) = m_zero
    st%surface(1, 3)%origin(3) = -bounds(1, 3)
    st%surface(1, 3)%n(:) = m_zero
    st%surface(1, 3)%n(3) = -m_one
    st%surface(1, 3)%u(:) = m_zero
    st%surface(1, 3)%u(1) = m_one
    st%surface(1, 3)%v(:) = m_zero
    st%surface(1, 3)%v(2) = -m_one
    st%surface(1, 3)%nu   = -int(bounds(1, 1)/mesh%spacing(1))
    st%surface(1, 3)%mu   =  int(bounds(1, 1)/mesh%spacing(1))
    st%surface(1, 3)%nv   = -int(bounds(1, 2)/mesh%spacing(2))
    st%surface(1, 3)%mv   =  int(bounds(1, 2)/mesh%spacing(2))
 
    ! x-y surface at +z boundary
    st%surface(2, 3)%spacing   = m_half*(mesh%spacing(1) + mesh%spacing(2))
    st%surface(2, 3)%origin(:) = m_zero
    st%surface(2, 3)%origin(3) = bounds(1, 3)
    st%surface(2, 3)%n(:) = m_zero
    st%surface(2, 3)%n(3) = m_one
    st%surface(2, 3)%u(:) = m_zero
    st%surface(2, 3)%u(1) = m_one
    st%surface(2, 3)%v(:) = m_zero
    st%surface(2, 3)%v(2) = m_one
    st%surface(2, 3)%nu   = -int(bounds(1, 1)/mesh%spacing(1))
    st%surface(2, 3)%mu   =  int(bounds(1, 1)/mesh%spacing(1))
    st%surface(2, 3)%nv   = -int(bounds(1, 2)/mesh%spacing(2))
    st%surface(2, 3)%mv   =  int(bounds(1, 2)/mesh%spacing(2))
 
    call profiling_out('MAXWELL_SURFACES_INIT')
 
    pop_sub(maxwell_surfaces_init)
  end subroutine maxwell_surfaces_init
 
  ! ---------------------------------------------------------
  subroutine maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
    type(bc_mxll_t),     intent(inout) :: bc
    class(mesh_t),       intent(in)    :: mesh
    integer,             intent(in)    :: ip
    real(real64),        intent(in)    :: bounds(:,:)
    integer,             intent(out)   :: point_info
 
    real(real64)   :: rr, dd, xx(3), width(3)
 
    point_info = 0
 
    width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
    xx = m_zero
    xx(1:mesh%box%dim) = mesh%x(ip, 1:mesh%box%dim)
 
    if (bc%ab_user_def) then
 
      dd = bc%ab_ufn(ip) - bounds(1, 1)
      if (dd > m_zero) then
        if (bc%ab_ufn(ip) < bounds(2, 1)) then
          point_info = 1
        end if
      end if
 
    else ! bc%ab_user_def == .false.
 
      select type (box => mesh%box)
      type is (box_sphere_t)
        rr = norm2(xx -  box%center)
        dd = rr -  bounds(1, 1)
        if (dd > m_zero) then
          if (dd  <  width(1)) then
            point_info = 1
          end if
        end if
 
      type is (box_parallelepiped_t)
        ! Limits of boundary region
        if (all(abs(xx(1:3) - box%center) <= bounds(2, 1:3))) then
          if (any(abs(xx(1:3) - box%center) > bounds(1, 1:3))) then
            point_info = 1
          else
            point_info = 0
          end if
        else
          point_info = -1
        end if
 
      class default
        ! Other box shapes are not supported.
        assert(.false.)
      end select
    end if
 
  end subroutine maxwell_box_point_info
 
  ! ---------------------------------------------------------
  subroutine maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
    class(mesh_t),       intent(in)    :: mesh
    integer,             intent(in)    :: ip
    real(real64),        intent(in)    :: bounds(:,:)
    integer,             intent(out)   :: boundary_info
 
    real(real64)   :: xx(3)
 
    boundary_info = 0
 
    xx = m_zero
    xx(1:mesh%box%dim) = mesh%x(ip, 1:mesh%box%dim)
    if (is_close(abs(xx(1)), bounds(1, 1)) .and. (all(abs(xx(2:3)) <= bounds(1, 2:3))) .or. &
      is_close(abs(xx(2)), bounds(1, 2)) .and. (all(abs(xx(1:3:2)) <= bounds(1, 1:3:2))) .or. &
      is_close(abs(xx(3)), bounds(1, 3)) .and. (all(abs(xx(1:2)) <= bounds(1, 1:2)))) then
      boundary_info = 1
    else
      boundary_info = 0
    end if
 
  end subroutine maxwell_boundary_point_info
 
  ! ---------------------------------------------------------
  subroutine inner_and_outer_points_mapping(mesh, st, bounds)
    class(mesh_t),       intent(in)    :: mesh
    type(states_mxll_t), intent(inout) :: st
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ip, ip_in, ip_bd, point_info
    real(real64)   :: xx(mesh%box%dim)
 
    push_sub(inner_and_outer_points_mapping)
 
    call profiling_in('IN_OUT_POINTS_MAP')
 
    ! allocate inner and boundary points points map
    ip_in = 0
    ip_bd = 0
    do ip = 1, mesh%np
      xx = mesh%x(ip, :)
      if (all(abs(xx) <= bounds(2,1:mesh%box%dim))) then
        if (any(abs(xx) > bounds(1,1:mesh%box%dim))) then
          point_info = 1
        else
          point_info = 0
        end if
      else
        point_info = -1
      end if
      if (point_info == 0) then
        ip_in = ip_in + 1
      else
        ip_bd = ip_bd + 1
      end if
    end do
    st%inner_points_number = ip_in
    safe_allocate(st%inner_points_map(1:ip_in))
    safe_allocate(st%inner_points_mask(1:mesh%np))
    st%boundary_points_number = ip_bd
    safe_allocate(st%boundary_points_map(1:ip_bd))
    safe_allocate(st%boundary_points_mask(1:mesh%np))
    st%inner_points_mask = .false.
    st%boundary_points_mask = .false.
 
    ! inner and boundary points mapping
    ip_in = 0
    ip_bd = 0
    do ip = 1, mesh%np
      xx = mesh%x(ip, :)
      if (all(abs(xx) <= bounds(2,1:mesh%box%dim))) then
        if (any(abs(xx) > bounds(1,1:mesh%box%dim))) then
          point_info = 1
        else
          point_info = 0
        end if
      else
        point_info = -1
      end if
      if (point_info == 0) then
        ip_in = ip_in + 1
        st%inner_points_map(ip_in) = ip
        st%inner_points_mask(ip) = .true.
      else
        ip_bd = ip_bd + 1
        st%boundary_points_map(ip_bd) = ip
        st%boundary_points_mask(ip) = .true.
      end if
    end do
 
    if (accel_is_enabled()) then
      call accel_create_buffer(st%buff_inner_points_map, accel_mem_read_only, type_integer, st%inner_points_number)
      call accel_write_buffer(st%buff_inner_points_map, st%inner_points_number, st%inner_points_map)
      call accel_create_buffer(st%buff_boundary_points_map, accel_mem_read_only, type_integer, st%boundary_points_number)
      call accel_write_buffer(st%buff_boundary_points_map, st%boundary_points_number, st%boundary_points_map)
    end if
 
    call profiling_out('IN_OUT_POINTS_MAP')
 
    pop_sub(inner_and_outer_points_mapping)
  end subroutine inner_and_outer_points_mapping
 
  ! ---------------------------------------------------------
  subroutine surface_grid_points_mapping(mesh, st, bounds)
    class(mesh_t),       intent(in)    :: mesh
    type(states_mxll_t), intent(inout) :: st
    real(real64),        intent(in)    :: bounds(:,:)
 
    integer :: ix, ix_max, iix, iy, iy_max, iiy, iz, iz_max, iiz, idx1, idx2, nn_max
    integer(int64) :: ip_global
    integer, allocatable :: nn(:,:,:,:)
    real(real64)   :: rr(3), delta(3), vec(2), min_1(3), max_1(3), min_2(3), max_2(3)
 
    push_sub(surface_grid_points_mapping)
 
    call profiling_in('SURF_GRID_POINTS_MAPPING')
 
    st%surface_grid_rows_number(1) = 3
    ix_max  = st%surface_grid_rows_number(1)
    st%surface_grid_rows_number(2) = 3
    iy_max  = st%surface_grid_rows_number(2)
    st%surface_grid_rows_number(3) = 3
    iz_max  = st%surface_grid_rows_number(3)
 
    delta(1) = m_two * abs(bounds(1,1)) / real(ix_max, real64)
    delta(2) = m_two * abs(bounds(1,2)) / real(iy_max, real64)
    delta(3) = m_two * abs(bounds(1,3)) / real(iz_max, real64)
 
    st%surface_grid_element(1) = delta(2) * delta(3)
    st%surface_grid_element(2) = delta(1) * delta(3)
    st%surface_grid_element(3) = delta(1) * delta(2)
 
    safe_allocate(nn(1:2, 1:3, 1:3, 1:3))
 
    st%surface_grid_center(1, 1, :, :) = -int(bounds(1,1))
    do iy = 1, iy_max
      do iz = 1, iz_max
        rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
        rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
        st%surface_grid_center(1, 2, iy, iz) = int(rr(2))
        st%surface_grid_center(1, 3, iy, iz) = int(rr(3))
      end do
    end do
    st%surface_grid_center(2, 1, :, :) = int(bounds(1,1))
    do iy = 1, iy_max
      do iz = 1, iz_max
        rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
        rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
        st%surface_grid_center(2, 2, iy, iz) = int(rr(2))
        st%surface_grid_center(2, 3, iy, iz) = int(rr(3))
      end do
    end do
 
    st%surface_grid_center(1, 2, :, :) = -int(bounds(1,2))
    do ix = 1, ix_max
      do iz = 1, iz_max
        rr(1) = -bounds(1,1) + delta(1)/m_two + (ix-1) * delta(1)
        rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
        st%surface_grid_center(1, 1, ix, iz) = int(rr(1))
        st%surface_grid_center(1, 3, ix, iz) = int(rr(3))
      end do
    end do
    st%surface_grid_center(2, 2, :, :) = int(bounds(1,2))
    do ix = 1, ix_max
      do iz = 1, iz_max
        rr(1) = -bounds(1,2) + delta(1)/m_two + (ix-1) * delta(1)
        rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
        st%surface_grid_center(2, 1, ix, iz) = int(rr(1))
        st%surface_grid_center(2, 3, ix, iz) = int(rr(3))
      end do
    end do
 
    st%surface_grid_center(1, 3, :, :) = -int(bounds(1,3))
    do ix = 1, ix_max
      do iy = 1, iy_max
        rr(1) = -bounds(1,1) + delta(1)/m_two + (ix-1) * delta(1)
        rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
        st%surface_grid_center(1, 1, ix, iy) = int(rr(1))
        st%surface_grid_center(1, 2, ix, iy) = int(rr(2))
      end do
    end do
    st%surface_grid_center(2, 3, :, :) = int(bounds(1,3))
    do ix = 1, ix_max
      do iy = 1, iy_max
        rr(1) = -bounds(1,2) + delta(1)/m_two + (ix-1) * delta(1)
        rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
        st%surface_grid_center(2, 1, ix, iy) = int(rr(1))
        st%surface_grid_center(2, 2, ix, iy) = int(rr(2))
      end do
    end do
 
    st%surface_grid_points_number(:,:,:) = 0
 
    nn_max = 0
 
    do iy = 1, iy_max
      do iz = 1, iz_max
        min_1(iy) = -bounds(1,2) + (iy-1) * delta(2)
        max_1(iy) = -bounds(1,2) + iy * delta(2)
        min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
        max_2(iz) = -bounds(1,3) + iz * delta(3)
      end do
    end do
    do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
      do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
        vec(1) = iiy * mesh%spacing(2)
        vec(2) = iiz * mesh%spacing(3)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          st%surface_grid_points_number(1, idx1, idx2) = st%surface_grid_points_number(1, idx1, idx2) + 1
          if (nn_max < st%surface_grid_points_number(1, idx1, idx2)) then
            nn_max = st%surface_grid_points_number(1, idx1, idx2)
          end if
        end if
      end do
    end do
 
    do ix = 1, ix_max
      do iz = 1, iz_max
        min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
        max_1(ix) = -bounds(1,1) + ix * delta(1)
        min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
        max_2(iz) = -bounds(1,3) + iz * delta(3)
      end do
    end do
    do iix = mesh%idx%nr(1, 1), mesh%idx%nr(2, 1)
      do iiz = mesh%idx%nr(1, 3), mesh%idx%nr(2, 3)
        vec(1)     = iix * mesh%spacing(1)
        vec(2)     = iiz * mesh%spacing(3)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          st%surface_grid_points_number(2,idx1,idx2) = st%surface_grid_points_number(2,idx1,idx2)+1
          if (nn_max < st%surface_grid_points_number(2, idx1, idx2)) then
            nn_max = st%surface_grid_points_number(2, idx1, idx2)
          end if
        end if
      end do
    end do
 
    do ix = 1, ix_max
      do iy = 1, iy_max
        min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
        max_1(ix) = -bounds(1,1) + ix * delta(1)
        min_2(iy) = -bounds(1,2) + (iy-1) * delta(2)
        max_2(iy) = -bounds(1,2) + iy * delta(2)
      end do
    end do
    do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
      do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
        vec(1) = iix * mesh%spacing(1)
        vec(2) = iiy * mesh%spacing(2)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          st%surface_grid_points_number(3, idx1, idx2) = st%surface_grid_points_number(3, idx1, idx2) + 1
          if (nn_max < st%surface_grid_points_number(3, idx1, idx2)) then
            nn_max = st%surface_grid_points_number(3, idx1, idx2)
          end if
        end if
      end do
    end do
 
    ! originally there were three allocated of the same pointer here
    safe_allocate(st%surface_grid_points_map(1:2, 1:st%dim, 1:ix_max, 1:iy_max, 1:nn_max))
 
    nn(:,:,:,:) = 0
 
    do iy = 1, iy_max
      do iz = 1, iz_max
        min_1(iy) = -bounds(1,2) + (iy-1) * delta(2)
        max_1(iy) = -bounds(1,2) + iy * delta(2)
        min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
        max_2(iz) = -bounds(1,3) + iz * delta(3)
      end do
    end do
    do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
      do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
        vec(1) = iiy * mesh%spacing(2)
        vec(2) = iiz * mesh%spacing(3)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          nn(1, 1, idx1, idx2) = nn(1, 1, idx1, idx2) + 1
          rr(1) = -bounds(1, 1)
          rr(2) = iiy * mesh%spacing(2)
          rr(3) = iiz * mesh%spacing(3)
          iix = int(-bounds(1,1)/mesh%spacing(1))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(1, 1, idx1, idx2, nn(1, 1, idx1, idx2)) = ip_global
          nn(2, 1, idx1, idx2) = nn(2, 1, idx1, idx2) + 1
          rr(1) = bounds(1,1)
          rr(2) = iiy * mesh%spacing(2)
          rr(3) = iiz * mesh%spacing(3)
          iix = int(bounds(1,1)/mesh%spacing(1))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(2, 1, idx1, idx2, nn(2, 1, idx1, idx2)) = ip_global
        end if
      end do
    end do
 
    do ix = 1, ix_max
      do iz = 1, iz_max
        min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
        max_1(ix) = -bounds(1,1) + ix * delta(1)
        min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
        max_2(iz) = -bounds(1,3) + iz * delta(3)
      end do
    end do
    do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
      do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
        vec(1) = iix * mesh%spacing(1)
        vec(2) = iiz * mesh%spacing(3)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          nn(1, 2, idx1, idx2) = nn(1, 2, idx1, idx2) + 1
          rr(1) = iix * mesh%spacing(1)
          rr(2) = -bounds(1, 2)
          rr(3) = iiz * mesh%spacing(3)
          iiy = int(-bounds(1,2)/mesh%spacing(2))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(1, 2, idx1, idx2, nn(1, 2, idx1, idx2)) = ip_global
          nn(2, 2, idx1, idx2) = nn(2, 2, idx1, idx2) + 1
          rr(1) = iix * mesh%spacing(1)
          rr(2) = bounds(1,2)
          rr(3) = iiz * mesh%spacing(3)
          iiy = int(bounds(1,2)/mesh%spacing(2))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(2, 2, idx1, idx2, nn(2, 2, idx1, idx2)) = ip_global
        end if
      end do
    end do
 
    do ix = 1, ix_max
      do iy = 1, iy_max
        min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
        max_1(ix) = -bounds(1,1) + ix * delta(1)
        min_2(iy) = -bounds(1,2) + (iy-1) * delta(2)
        max_2(iy) = -bounds(1,2) + iy * delta(2)
      end do
    end do
    do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
      do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
        vec(1) = iix * mesh%spacing(1)
        vec(2) = iiy * mesh%spacing(2)
        call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
        if (idx1 /= 0 .and. idx2 /= 0) then
          nn(1, 3, idx1, idx2) = nn(1, 3, idx1, idx2) + 1
          rr(1) = iix * mesh%spacing(1)
          rr(2) = iiy * mesh%spacing(2)
          rr(3) = -bounds(1,3)
          iiz = int(-bounds(1,3)/mesh%spacing(3))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(1, 3, idx1, idx2, nn(1, 3, idx1, idx2)) = ip_global
          nn(2, 3, idx1, idx2) = nn(2, 3, idx1, idx2) + 1
          rr(1) = iix * mesh%spacing(1)
          rr(2) = iiy * mesh%spacing(2)
          rr(3) = bounds(1,3)
          iiz = int(bounds(1,3)/mesh%spacing(3))
          ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
          st%surface_grid_points_map(2, 3, idx1, idx2, nn(2, 3, idx1, idx2)) = ip_global
        end if
      end do
    end do
 
    safe_deallocate_a(nn)
 
    call profiling_out('SURF_GRID_POINTS_MAPPING')
 
    pop_sub(surface_grid_points_mapping)
  contains
 
    subroutine get_surface_indices(vec, min_1, max_1, min_2, max_2, index_1, index_2)
      real(real64),   intent(in)  :: vec(:)
      real(real64),   intent(in)  :: min_1(:)
      real(real64),   intent(in)  :: max_1(:)
      real(real64),   intent(in)  :: min_2(:)
      real(real64),   intent(in)  :: max_2(:)
      integer, intent(out) :: index_1
      integer, intent(out) :: index_2
 
      if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
        index_1 = 1
        index_2 = 1
      else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
        index_1 = 2
        index_2 = 1
      else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
        index_1 = 3
        index_2 = 1
      else if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
        index_1 = 1
        index_2 = 2
      else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
        index_1 = 2
        index_2 = 2
      else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
        index_1 = 3
        index_2 = 2
      else if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
        index_1 = 1
        index_2 = 3
      else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
        index_1 = 2
        index_2 = 3
      else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
        index_1 = 3
        index_2 = 3
      else
        index_1 = 0
        index_2 = 0
      end if
 
    end subroutine get_surface_indices
 
  end subroutine surface_grid_points_mapping
 
end module maxwell_boundary_op_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: