linear_medium.F90

Go to the documentation of this file.

!! Copyright (C) 2020 F. Bonafé, H. Appel, R. Jestädt
!!
!! 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 linear_medium_oct_m
  use algorithm_oct_m
  use calc_mode_par_oct_m
  use cgal_polyhedra_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use grid_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use iso_c_binding
  use io_oct_m
  use math_oct_m
  use linear_medium_to_em_field_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use output_low_oct_m
  use output_linear_medium_oct_m
  use parser_oct_m
  use propagator_exp_mid_2step_oct_m
  use propagator_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use regridding_oct_m
  use space_oct_m
  use system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                    &
    linear_medium_t,           &
    linear_medium_init,        &
    get_medium_box_points_map, &
    medium_box_init
 
  integer, parameter :: &
    MEDIUM_PARALLELEPIPED = 1,         &
    medium_box_file       = 2
 
  type, extends(system_t) :: linear_medium_t
    type(space_t)         :: space
    real(real64)          :: ep_factor
    real(real64)          :: mu_factor
    real(real64)          :: sigma_e_factor
    real(real64)          :: sigma_m_factor
    integer               :: edge_profile
    type(output_t)        :: outp
    type(grid_t)          :: gr
    type(multicomm_t)     :: mc
    type(single_medium_box_t) :: medium_box
 
  contains
    procedure :: init_interaction => linear_medium_init_interaction
    procedure :: init_interaction_as_partner => linear_medium_init_interaction_as_partner
    procedure :: initialize => linear_medium_initialize
    procedure :: do_algorithmic_operation => linear_medium_do_algorithmic_operation
    procedure :: is_tolerance_reached => linear_medium_is_tolerance_reached
    procedure :: copy_quantities_to_interaction => linear_medium_copy_quantities_to_interaction
    procedure :: restart_write_data => linear_medium_restart_write_data
    procedure :: restart_read_data => linear_medium_restart_read_data
    procedure :: update_kinetic_energy => linear_medium_update_kinetic_energy
    final :: linear_medium_finalize
  end type linear_medium_t
 
  interface linear_medium_t
    procedure linear_medium_constructor
  end interface linear_medium_t
 
contains
 
  ! ---------------------------------------------------------
  function linear_medium_constructor(namespace) result(sys)
    class(linear_medium_t), pointer    :: sys
    type(namespace_t),           intent(in) :: namespace
 
    push_sub(linear_medium_constructor)
 
    allocate(sys)
 
    call linear_medium_init(sys, namespace)
 
    pop_sub(linear_medium_constructor)
  end function linear_medium_constructor
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine linear_medium_init(this, namespace)
    class(linear_medium_t), target, intent(inout) :: this
    type(namespace_t),            intent(in)      :: namespace
 
    integer :: nlines, ncols,n_points, n_global_points
    integer(int64) :: index_range(4)
    integer, allocatable :: tmp(:)
    type(block_t) :: blk
 
    push_sub(linear_medium_init)
 
    call profiling_in('MEDIUM_BOX_INIT')
 
    this%namespace = namespace
    this%space = space_t(this%namespace)
    if (this%space%is_periodic()) then
      call messages_not_implemented('Linear medium for periodic systems', namespace=namespace)
    end if
    call grid_init_stage_1(this%gr, this%namespace, this%space)
    ! store the ranges for these two indices (serves as initial guess
    ! for parallelization strategy)
    index_range(1) = this%gr%np_global  ! Number of points in mesh
    index_range(2) = 1                      ! Number of states
    index_range(3) = 1                      ! Number of k-points
    index_range(4) = 100000                 ! Some large number
 
    ! create index and domain communicators
    call multicomm_init(this%mc, this%namespace, mpi_world, calc_mode_par, mpi_world%size, &
      index_range, (/ 5000, 1, 1, 1 /))
    call grid_init_stage_2(this%gr, this%namespace, this%space, this%mc)
 
    call medium_box_init(this%medium_box, this%namespace)
 
    !%Variable LinearMediumProperties
    !%Type block
    !%Section Maxwell::Medium
    !%Description
    !% Defines electromagnetic parameters for a linear medium box.
    !%
    !% Example:
    !%
    !% <tt>%LinearMediumProperties
    !% <br>&nbsp;&nbsp;   epsilon_factor | mu_factor | sigma_e | sigma_m
    !% <br>%</tt>
    !%
    !% Permittivity factor, permeability factor, electric conductivity and magnetic conductivity of the medium box.
    !%End
 
    if (parse_block(namespace, 'LinearMediumProperties', blk) == 0) then
 
      nlines = parse_block_n(blk)
      if (nlines /=  1) then
        call messages_input_error(namespace, 'LinearMediumProperties', 'should consist of one line', nlines)
      end if
      ncols = parse_block_cols(blk, 0)
      if (ncols /= 4) then
        call messages_input_error(namespace, 'LinearMediumProperties', 'should consist of four columns')
      end if
      call parse_block_float(blk, 0, 0, this%ep_factor)
      call parse_block_float(blk, 0, 1, this%mu_factor)
      call parse_block_float(blk, 0, 2, this%sigma_e_factor)
      call parse_block_float(blk, 0, 3, this%sigma_m_factor)
      write(message(1),'(a,es9.2)') 'Box epsilon factor: ', this%ep_factor
      write(message(2),'(a,es9.2)') 'Box mu factor:      ', this%mu_factor
      write(message(3),'(a,es9.2)') 'Box electric sigma: ', this%sigma_e_factor
      write(message(4),'(a,es9.2)') 'Box magnetic sigma: ', this%sigma_m_factor
      write(message(5),'(a)') ""
      call messages_info(5, namespace=namespace)
      call parse_block_end(blk)
 
      call messages_print_with_emphasis(namespace=namespace)
    else
      message(1) = 'You must specify the properties of your linear medium through the LinearMediumProperties block.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    !%Variable LinearMediumEdgeProfile
    !%Type integer
    !%Section Maxwell::Medium
    !%Description
    !% Defines the type of numerical approximation used for the derivatives at the edges of the medium box.
    !% Default is edged. When the box shape is read from file, only the edged profile is supported.
    !%
    !%Option edged 1
    !% Medium box edges are considered steep for derivatives.
    !%Option smooth 2
    !% Medium box edged and softened for derivatives.
    !%End
    call parse_variable(namespace, 'LinearMediumEdgeProfile', option__linearmediumedgeprofile__edged, this%edge_profile)
    if (this%edge_profile == option__linearmediumedgeprofile__edged) then
      write(message(1),'(a,a)')   'Box shape:          ', 'edged'
    else if (this%edge_profile == option__linearmediumedgeprofile__smooth) then
      write(message(1),'(a,a)')   'Box shape:          ', 'smooth'
    end if
    call messages_info(1, namespace=namespace)
 
    allocate(this%supported_interactions(0))
    this%supported_interactions_as_partner = [linear_medium_to_em_field]
    call this%quantities%add(quantity_t("permitivity", updated_on_demand = .false., always_available = .true.))
    call this%quantities%add(quantity_t("permeability", updated_on_demand = .false., always_available = .true.))
    call this%quantities%add(quantity_t("E conductivity", updated_on_demand = .false., always_available = .true.))
    call this%quantities%add(quantity_t("M conductivity", updated_on_demand = .false., always_available = .true.))
 
    call output_linear_medium_init(this%outp, this%namespace, this%space)
 
    call get_medium_box_points_map(this%medium_box, this%gr)
    call get_linear_medium_em_properties(this, this%medium_box, this%gr)
    call output_linear_medium(this%outp, this%namespace, this%space, this%gr, &
      this%outp%iter_dir, this%medium_box%points_map, this%medium_box%ep, &
      this%medium_box%mu, this%medium_box%c)
 
    if (this%medium_box%check_medium_points) then
      safe_allocate(tmp(1:this%gr%np))
      n_global_points = 0
      write(message(1),'(a, a, a)')   'Check of points inside surface of medium ', trim(this%medium_box%filename), ":"
      call messages_info(1, namespace=this%namespace)
      call get_points_map_from_file(this%medium_box%filename, this%gr, n_points, n_global_points, tmp, 0.99_real64)
      write(message(1),'(a, I8)')'Number of points inside medium (normal coordinates):', &
        this%medium_box%global_points_number
      write(message(2),'(a, I8)')'Number of points inside medium (rescaled coordinates):', n_global_points
      write(message(3), '(a)') ""
      call messages_info(3, namespace=this%namespace)
      safe_deallocate_a(tmp)
    end if
 
    call profiling_out('MEDIUM_BOX_INIT')
 
    pop_sub(linear_medium_init)
  end subroutine linear_medium_init
 
  ! ---------------------------------------------------------
  subroutine linear_medium_init_interaction(this, interaction)
    class(linear_medium_t), target, intent(inout) :: this
    class(interaction_t),           intent(inout) :: interaction
 
    push_sub(linear_medium_init_interaction)
 
    select type (interaction)
    class default
      message(1) = "Trying to initialize an unsupported interaction by a linear medium."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(linear_medium_init_interaction)
  end subroutine linear_medium_init_interaction
 
  ! ---------------------------------------------------------
  subroutine linear_medium_init_interaction_as_partner(partner, interaction)
    class(linear_medium_t),   intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    type(regridding_t), pointer :: regridding
    type(single_medium_box_t), pointer :: medium_box_grid
 
    push_sub(linear_medium_init_interaction_as_partner)
 
    select type (interaction)
    type is (linear_medium_to_em_field_t)
      regridding => regridding_t(interaction%system_gr, &
        partner%gr, partner%space, partner%namespace)
 
      ! create a medium box with the size of the partner grid and map all quantities to it
      medium_box_grid => partner%medium_box%to_grid(partner%gr)
 
      ! now do the grid transfer to the medium box stored in the interaction which has the size of the system grid
      call regridding%do_transfer(interaction%medium_box%aux_ep, medium_box_grid%aux_ep)
      call regridding%do_transfer(interaction%medium_box%aux_mu, medium_box_grid%aux_mu)
      call regridding%do_transfer(interaction%medium_box%ep, medium_box_grid%ep)
      call regridding%do_transfer(interaction%medium_box%mu, medium_box_grid%mu)
      call regridding%do_transfer(interaction%medium_box%c, medium_box_grid%c)
      call regridding%do_transfer(interaction%medium_box%sigma_e, medium_box_grid%sigma_e)
      call regridding%do_transfer(interaction%medium_box%sigma_m, medium_box_grid%sigma_m)
 
      safe_deallocate_p(regridding)
      safe_deallocate_p(medium_box_grid)
    class default
      message(1) = "Trying to initialize an unsupported interaction by a linear medium."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(linear_medium_init_interaction_as_partner)
  end subroutine linear_medium_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine linear_medium_initialize(this)
    class(linear_medium_t), intent(inout) :: this
 
    push_sub(linear_medium_initialize)
 
    pop_sub(linear_medium_initialize)
  end subroutine linear_medium_initialize
 
  ! ---------------------------------------------------------
  logical function linear_medium_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(linear_medium_t),         intent(inout) :: this
    class(algorithmic_operation_t), intent(in)    :: operation
    character(len=:), allocatable,  intent(out)   :: updated_quantities(:)
 
    push_sub(linear_medium_do_algorithmic_operation)
 
    ! Currently there are no linear medium specific algorithmic operations
    done = .false.
 
    pop_sub(linear_medium_do_algorithmic_operation)
  end function linear_medium_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  logical function linear_medium_is_tolerance_reached(this, tol) result(converged)
    class(linear_medium_t),   intent(in)    :: this
    real(real64),              intent(in)    :: tol
 
    push_sub(linear_medium_is_tolerance_reached)
 
    ! this routine is never called at present, no reason to be here
    assert(.false.)
    converged = .false.
 
    pop_sub(linear_medium_is_tolerance_reached)
  end function linear_medium_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine linear_medium_copy_quantities_to_interaction(partner, interaction)
    class(linear_medium_t),          intent(inout) :: partner
    class(interaction_surrogate_t),       intent(inout) :: interaction
 
    push_sub(linear_medium_copy_quantities_to_interaction)
 
    select type (interaction)
    type is (linear_medium_to_em_field_t)
      ! No need to copy anything
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(linear_medium_copy_quantities_to_interaction)
  end subroutine linear_medium_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine linear_medium_restart_write_data(this)
    class(linear_medium_t), intent(inout) :: this
 
    integer :: restart_file_unit
 
    push_sub(linear_medium_restart_write_data)
 
    message(1) = "Linear medium system "//trim(this%namespace%get())//" will only write a dummy restart file."
    call messages_warning(1, namespace=this%namespace)
 
    call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
    restart_file_unit = io_open('restart/'//td_dir// 'restart_linear_medium', this%namespace, action='write')
    write(restart_file_unit, *) 'Linear medium restart file'
    call io_close(restart_file_unit)
 
    message(1) = "Successfully wrote restart data for system "//trim(this%namespace%get())
    call messages_info(1, namespace=this%namespace)
 
    pop_sub(linear_medium_restart_write_data)
  end subroutine linear_medium_restart_write_data
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function linear_medium_restart_read_data(this)
    class(linear_medium_t), intent(inout) :: this
 
    integer :: restart_file_unit
 
    push_sub(linear_medium_restart_read_data)
 
    call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
    restart_file_unit = io_open('restart/'//td_dir// 'restart_linear_medium', this%namespace, action='read', die=.false.)
    if (restart_file_unit /= -1) then
      call io_close(restart_file_unit)
      linear_medium_restart_read_data = .true.
    else
      ! could not open file
      linear_medium_restart_read_data = .false.
    end if
 
    if (linear_medium_restart_read_data) then
      message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(linear_medium_restart_read_data)
  end function linear_medium_restart_read_data
 
  subroutine linear_medium_update_kinetic_energy(this)
    class(linear_medium_t), intent(inout) :: this
 
    push_sub(linear_medium_update_kinetic_energy)
 
    this%kinetic_energy = m_zero
 
    pop_sub(linear_medium_update_kinetic_energy)
 
  end subroutine linear_medium_update_kinetic_energy
 
  ! ---------------------------------------------------------
  subroutine linear_medium_finalize(this)
    type(linear_medium_t), intent(inout) :: this
 
    push_sub(linear_medium_finalize)
 
    call single_medium_box_end(this%medium_box)
    call system_end(this)
    call multicomm_end(this%mc)
    call grid_end(this%gr)
 
    pop_sub(linear_medium_finalize)
  end subroutine linear_medium_finalize
 
 
  ! Specific routines for this system:
 
  ! ---------------------------------------------------------
  subroutine get_medium_box_points_map(medium_box, gr)
    type(single_medium_box_t),   intent(inout) :: medium_box
    type(grid_t),                intent(in)    :: gr
 
    integer :: ip, ip_in, ip_bd, idim
    integer, allocatable :: tmp_points_map(:), tmp_bdry_map(:)
    real(real64)   :: bounds(2,gr%box%dim), xx(gr%box%dim)
    logical :: inside
 
    push_sub(get_medium_box_points_map)
 
    call profiling_in('GET_MEDIUM_BOX_POINTS_MAP')
 
    safe_allocate(tmp_points_map(1:gr%np))
    safe_allocate(tmp_bdry_map(1:gr%np))
    tmp_points_map = 0
    tmp_bdry_map = 0
 
    if (medium_box%box_shape == medium_box_file) then
 
      call get_points_map_from_file(medium_box%filename, gr, medium_box%points_number,&
        medium_box%global_points_number, tmp_points_map)
 
    else
 
      do idim = 1, 3
        bounds(1,idim) = medium_box%center(idim) - medium_box%lsize(idim)/m_two
        bounds(2,idim) = medium_box%center(idim) + medium_box%lsize(idim)/m_two
      end do
      ip_in = 0
      ip_bd = 0
      do ip = 1, gr%np
        xx(1:3) = gr%x(ip, 1:3)
        inside = check_point_in_bounds(xx, bounds(:,:))
        if (check_point_in_bounds(xx, bounds(:,:))) then
          ip_in = ip_in + 1
          tmp_points_map(ip_in) = ip
        end if
        if (check_point_on_bounds(xx, bounds(:,:))) then
          ip_bd = ip_bd + 1
          tmp_bdry_map(ip_bd) = ip
        end if
      end do
 
      medium_box%points_number = ip_in
      medium_box%bdry_number = ip_bd
 
      safe_allocate(medium_box%bdry_map(1:medium_box%bdry_number))
      medium_box%bdry_map = 0
      medium_box%bdry_map = tmp_bdry_map(1:medium_box%bdry_number)
 
    end if
    call single_medium_box_allocate(medium_box, medium_box%points_number)
    medium_box%points_map(:) = tmp_points_map(1:medium_box%points_number)
 
    ! TODO: add some check that medium boxes do not overlap (now they are different systems)
    safe_deallocate_a(tmp_points_map)
    safe_deallocate_a(tmp_bdry_map)
 
    call profiling_out('GET_MEDIUM_BOX_POINTS_MAP')
 
    pop_sub(get_medium_box_points_map)
  contains
 
    logical pure function check_point_in_bounds(xx, bounds) result (check)
      real(real64), intent(in) :: xx(:)
      real(real64), intent(in) :: bounds(:,:)
 
      check = .false.
      if ((xx(1) >= bounds(1,1)) .and. (xx(1) <= bounds(2,1)) .and. &
        (xx(2) >= bounds(1,2)) .and. (xx(2) <= bounds(2,2)) .and. &
        (xx(3) >= bounds(1,3)) .and. (xx(3) <= bounds(2,3)) ) then
        check = .true.
      end if
 
    end function check_point_in_bounds
 
    logical pure function check_point_on_bounds(xx, bounds) result (check)
      real(real64), intent(in) :: xx(:)
      real(real64), intent(in) :: bounds(:,:)
 
      check = .false.
      if (is_close(xx(1), bounds(1,1)) .and. (xx(2) >= bounds(1,2) .and. xx(3) >= bounds(1,3)) &
        .and. (xx(2) <= bounds(2,2) .and. xx(3) <= bounds(2,3)) .or. &
        is_close(xx(2), bounds(1,2)) .and. (xx(1) >= bounds(1,1) .and. xx(3) >= bounds(1,3)) &
        .and. (xx(1) <= bounds(2,1) .and. xx(3) <= bounds(2,3)) .or. &
        is_close(xx(3), bounds(1,3)) .and. (xx(1) >= bounds(1,1) .and. xx(2) >= bounds(1,2)) &
        .and. (xx(1) <= bounds(2,1) .and. xx(2) <= bounds(2,2)) .or. &
        is_close(xx(1), bounds(2,1)) .and. (xx(2) >= bounds(1,2) .and. xx(3) >= bounds(1,3)) &
        .and. (xx(2) <= bounds(2,2) .and. xx(3) <= bounds(2,3)) .or. &
        is_close(xx(2), bounds(2,2)) .and. (xx(1) >= bounds(1,1) .and. xx(3) >= bounds(1,3)) &
        .and. (xx(1) <= bounds(2,1) .and. xx(3) <= bounds(2,3)) .or. &
        is_close(xx(3), bounds(2,3)) .and. (xx(1) >= bounds(1,1) .and. xx(2) >= bounds(1,2)) &
        .and. (xx(1) <= bounds(2,1) .and. xx(2) <= bounds(2,2)) ) then
        check = .true.
      end if
 
    end function check_point_on_bounds
 
  end subroutine get_medium_box_points_map
 
  ! ----------------------------------------------------------
  subroutine get_linear_medium_em_properties(this, medium_box, gr)
    type(linear_medium_t), intent(in)          :: this
    type(single_medium_box_t),   intent(inout) :: medium_box
    type(grid_t),                intent(in)    :: gr
 
    integer :: ip, ip_in, ip_bd, ipp
    real(real64)   :: xx(gr%box%dim), xxp(gr%box%dim), dd, dd_max, dd_min
    real(real64), allocatable  :: tmp(:), tmp_grad(:,:)
 
    push_sub(get_linear_medium_em_properties)
 
    call profiling_in('GET_LIN_MEDIUM_EM_PROP')
 
    safe_allocate(tmp(1:gr%np_part))
    safe_allocate(tmp_grad(1:gr%np_part, 1:gr%box%dim))
    dd_max = max(2*gr%spacing(1), 2*gr%spacing(2), 2*gr%spacing(3))
 
    do ip_in = 1, medium_box%points_number
      ip = medium_box%points_map(ip_in)
      if (this%edge_profile == option__linearmediumedgeprofile__smooth) then
        assert(allocated(medium_box%bdry_map))
 
        xx(1:3) = gr%x(ip,1:3)
        dd_min = m_huge
 
        do ip_bd = 1, medium_box%bdry_number
          ipp = medium_box%bdry_map(ip_bd)
          xxp(1:3) = gr%x(ipp,1:3)
          dd = norm2(xx(1:3) - xxp(1:3))
          if (dd < dd_min) dd_min = dd
        end do
 
        medium_box%ep(ip_in) = p_ep + ((p_ep * this%ep_factor - p_ep)  &
          * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
        medium_box%mu(ip_in) = p_mu + ((p_mu * this%mu_factor - p_mu) &
          * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
        medium_box%c(ip_in) = m_one/sqrt(medium_box%ep(ip_in)*medium_box%mu(ip_in))
        medium_box%sigma_e(ip_in) = this%sigma_e_factor &
          * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max)))
        medium_box%sigma_m(ip_in) = this%sigma_m_factor &
          * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max)))
 
      else if (this%edge_profile == option__linearmediumedgeprofile__edged) then
 
        medium_box%ep(ip_in) = p_ep * this%ep_factor
        medium_box%mu(ip_in) = p_mu * this%mu_factor
        medium_box%c(ip_in) = m_one/sqrt(medium_box%ep(ip_in)*medium_box%mu(ip_in))
        medium_box%sigma_e(ip_in) = this%sigma_e_factor
        medium_box%sigma_m(ip_in) = this%sigma_m_factor
 
      end if
    end do
 
    tmp(:) = p_ep
    do ip_in = 1, medium_box%points_number
      ip = medium_box%points_map(ip_in)
      tmp(ip)= medium_box%ep(ip_in)
    end do
    call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
    do ip_in = 1, medium_box%points_number
      ip = medium_box%points_map(ip_in)
      medium_box%aux_ep(ip_in, :) = tmp_grad(ip, :)/(m_four * medium_box%ep(ip_in))
    end do
 
    tmp(:) = p_mu
    do ip_in = 1, medium_box%points_number
      ip = medium_box%points_map(ip_in)
      tmp(ip) = medium_box%mu(ip_in)
    end do
    call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
    do ip_in = 1, medium_box%points_number
      ip = medium_box%points_map(ip_in)
      medium_box%aux_mu(ip_in, :) = tmp_grad(ip, :)/(m_four * medium_box%mu(ip_in))
    end do
 
    !TODO: add print information about the medium box
 
    safe_deallocate_a(tmp)
    safe_deallocate_a(tmp_grad)
 
    call profiling_out('GET_LIN_MEDIUM_EM_PROP')
 
    pop_sub(get_linear_medium_em_properties)
  end subroutine get_linear_medium_em_properties
 
  ! ----------------------------------------------------------
  subroutine get_points_map_from_file(filename, mesh, n_points, global_points_number, tmp_map, scale_factor)
    character(len=256),       intent(in)    :: filename
    class(mesh_t),            intent(in)    :: mesh
    integer,                  intent(out)   :: n_points
    integer,                  intent(out)   :: global_points_number
    integer,                  intent(inout) :: tmp_map(:)
    real(real64), optional,          intent(in)    :: scale_factor
 
    integer :: ip_in, ip
    real(real64)   :: xx(3)
    type(cgal_polyhedra_t) :: cgal_poly
 
    push_sub(get_points_map_from_file)
 
    call profiling_in('GET_POINTS_MAP_FROM_FILE')
 
    call cgal_polyhedron_init(cgal_poly, trim(filename), verbose = .false.)
 
    ip_in = 0
    do ip = 1, mesh%np
      if (present(scale_factor)) then
        xx(1:3) = scale_factor * mesh%x(ip, 1:3)
      else
        xx(1:3) = mesh%x(ip, 1:3)
      end if
      if (cgal_polyhedron_point_inside(cgal_poly, xx(1), xx(2), xx(3))) then
        ip_in = ip_in + 1
        tmp_map(ip_in) = ip
      end if
    end do
    n_points = ip_in
    call cgal_polyhedron_end(cgal_poly)
 
    call mpi_world%allreduce(ip_in, global_points_number, 1, mpi_integer, mpi_sum)
 
    call profiling_out('GET_POINTS_MAP_FROM_FILE')
 
    pop_sub(get_points_map_from_file)
 
  end subroutine get_points_map_from_file
 
  ! ----------------------------------------------------------
  subroutine medium_box_init(medium_box, namespace)
    type(single_medium_box_t), intent(inout) :: medium_box
    type(namespace_t),      intent(in)       :: namespace
 
    integer :: nlines, ncols, idim
    type(block_t) :: blk
 
    push_sub(medium_box_init)
 
    !%Variable LinearMediumBoxShape
    !%Type integer
    !%Section Maxwell::Medium
    !%Description
    !% This variable defines the shape of the linear medium box.
    !% The default is <tt>medium_parallelepiped</tt>.
    !%Option medium_parallelepiped 1
    !% The medium box will be a parallelepiped whose center and dimensions are taken from
    !% the variable <tt>LinearMediumBoxSize</tt>.
    !%Option medium_box_file 2
    !% The simulation box will be read from an external file in OFF format, defined by the variable <tt>LinearMediumBoxFile</tt>.
    !%End
    call parse_variable(namespace, 'LinearMediumBoxShape', medium_parallelepiped, medium_box%box_shape)
 
    if (.not. varinfo_valid_option('LinearMediumBoxShape', medium_box%box_shape)) then
      call messages_input_error(namespace, 'LinearMediumBoxShape')
    end if
 
    select case (medium_box%box_shape)
    case (medium_parallelepiped)
 
      !%Variable LinearMediumBoxSize
      !%Type block
      !%Section Maxwell::Medium
      !%Description
      !% Defines center and size of a parallelepiped linear medium box.
      !%
      !% Example:
      !%
      !% <tt>%LinearMediumBoxSize
      !% <br>&nbsp;&nbsp;   center_x | center_y | center_z | x_length | y_length | z_length
      !% <br>%</tt>
      !%End
 
      if (parse_block(namespace, 'LinearMediumBoxSize', blk) == 0) then
        call messages_print_with_emphasis(msg=trim('Linear Medium box center and size:'), namespace=namespace)
 
        nlines = parse_block_n(blk)
        if (nlines /=  1) then
          call messages_input_error(namespace, 'LinearMediumBoxSize', 'should consist of one line')
        end if
        ncols = parse_block_cols(blk, 0)
        if (ncols /= 6) then
          call messages_input_error(namespace, 'LinearMediumBoxSize', 'should consist of six columns')
        end if
        do idim = 1, 3
          call parse_block_float(blk, 0, idim-1, medium_box%center(idim))
          call parse_block_float(blk, 0, idim+2, medium_box%lsize(idim))
        end do
        write(message(1),'(a,es9.2,a,es9.2,a,es9.2)') 'Box center:         ', medium_box%center(1), ' | ',&
          medium_box%center(2), ' | ', medium_box%center(3)
        write(message(2),'(a,es9.2,a,es9.2,a,es9.2)') 'Box size  :         ', medium_box%lsize(1), ' | ', &
          medium_box%lsize(2), ' | ', medium_box%lsize(3)
        write(message(3),'(a)') ""
        call messages_info(3)
        call parse_block_end(blk)
 
        call messages_print_with_emphasis(namespace=namespace)
      else
        message(1) = "For parallelepiped box shapes, you must provide a LinearMediumBoxSize block."
        call messages_fatal(1, namespace=namespace)
      end if
 
    case (medium_box_file)
      !%Variable LinearMediumBoxFile
      !%Type string
      !%Section Maxwell::Medium
      !%Description
      !% File in OFF format with the shape of the linear medium.
      !%End
      if (parse_is_defined(namespace, 'LinearMediumBoxFile')) then
        call parse_variable(namespace, 'LinearMediumBoxFile', 'mediumboxfile', medium_box%filename)
      else
        message(1) = "When using box_file as the box shape, you must provide a filename through the LinearMediumBoxFile variable."
        call messages_fatal(1, namespace=namespace)
      end if
 
      !%Variable CheckPointsMediumFromFile
      !%Type logical
      !%Default no
      !%Section Maxwell::Medium
      !%Description
      !% Whether to re-calculate the points map by artificially shrinking the coordinate system by a factor of
      !% 0.99 to check if the points inside the medium surface are properly detected. This works for only one
      !% medium surface which is centered in the origin of the coordinate system.
      !%End
      call parse_variable(namespace, 'CheckPointsMediumFromFile', .false., medium_box%check_medium_points)
 
    end select
 
    pop_sub(medium_box_init)
 
  end subroutine medium_box_init
 
end module linear_medium_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: