sll_m_particle_mesh_coupling_spline_cl_3d_feec.F90

Go to the documentation of this file.

 
 
module sll_m_particle_mesh_coupling_spline_cl_3d_feec
 
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_errors.h"
#include "sll_working_precision.h"
 
  use sll_m_gauss_legendre_integration, only : &
       sll_f_gauss_legendre_points_and_weights
 
  use sll_m_low_level_bsplines, only: &
       sll_s_uniform_bsplines_eval_basis, &
       sll_s_uniform_bsplines_eval_deriv
 
  use sll_m_particle_mesh_coupling_base_3d, only: &
       sll_c_particle_mesh_coupling_3d
 
  use sll_m_splines_pp, only: &
       sll_t_spline_pp_1d, &
       sll_t_spline_pp_3d, &
       sll_s_spline_pp_init_3d, &
       sll_s_spline_pp_free_3d, &
       sll_f_spline_pp_horner_1d, &
       sll_s_spline_pp_horner_m_1d, &
       sll_f_spline_pp_horner_3d
 
  implicit none
 
  public :: &
       sll_t_particle_mesh_coupling_spline_cl_3d_feec, &
       sll_s_uniform_bsplines_eval_basis_clamped
 
  private
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  type, extends(sll_c_particle_mesh_coupling_3d) :: sll_t_particle_mesh_coupling_spline_cl_3d_feec
 
     sll_real64 :: rdelta_x(3)  
     ! Information about the particles
     sll_int32  :: no_particles 
     sll_int32  :: n_quad_points 
 
     sll_real64, allocatable :: spline_val(:,:) 
     sll_real64, allocatable :: spline_0(:,:) 
     sll_real64, allocatable :: spline1_0(:,:) 
     sll_real64, allocatable :: spline_1(:,:) 
     sll_real64, allocatable :: spline_2d(:,:) 
     sll_real64, allocatable :: spline1_2d(:,:) 
     sll_real64, allocatable :: spline1_3d(:,:,:) 
     sll_real64, allocatable :: spline2_3d(:,:,:) 
     sll_real64, allocatable :: spline_val_more(:,:) 
     sll_real64, allocatable :: j1d(:) 
 
     sll_real64, allocatable :: quad_xw(:,:) 
 
     sll_int32, allocatable :: index1d(:,:) 
 
     sll_int32              :: n_total0 
     sll_int32              :: n_total1 
 
     ! For 1d quadrature in 3d space
     sll_int32 :: n_quad_points_line 
     sll_real64, allocatable :: quad_xw_line(:,:) 
 
   contains
 
     procedure :: add_charge => add_charge_single_spline_cl_3d_feec
     procedure :: add_particle_mass => add_particle_mass_spline_cl_3d_feec
     procedure :: add_particle_mass_od => add_particle_mass_od_spline_cl_3d_feec
     procedure :: evaluate => evaluate_field_single_spline_cl_3d_feec
     procedure :: evaluate_multiple => evaluate_multiple_spline_cl_3d_feec
     procedure :: add_current => add_current_cl_3d
     procedure :: add_current_evaluate
     procedure :: add_current_update_v_component1 => add_current_update_v_component1_spline_3d_feec
     procedure :: add_current_update_v_component2 => add_current_update_v_component2_spline_3d_feec
     procedure :: add_current_update_v_component3 => add_current_update_v_component3_spline_3d_feec
 
     procedure :: init => init_spline_cl_3d_feec
     procedure :: free => free_spline_cl_3d_feec
  end type sll_t_particle_mesh_coupling_spline_cl_3d_feec
 
  type :: vector
     sll_real64, allocatable :: vals(:)
  end type vector
 
contains
 
 
  !---------------------------------------------------------------------------!
  subroutine add_charge_single_spline_cl_3d_feec(self, position, marker_charge, degree, rho_dofs)
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(3) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_int32,                                intent( in )    :: degree(3) 
    sll_real64,                               intent( inout ) :: rho_dofs(:) 
    !local variables
    sll_int32 :: box(3)
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k
 
 
    call convert_x_to_xbox( self, position, xi, box )
 
    self%spline_val = 0.0_f64
    if(degree(1)==self%spline_degree(1)) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    else if(degree(1)==self%spline_degree(1)-1) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    end if
 
    do j=2,3
       call sll_s_uniform_bsplines_eval_basis( degree(j), xi(j), self%spline_val(1:degree(j)+1,j) )
    end do
 
    box(2:3) = box(2:3)-degree(2:3)
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, degree(3)+1       
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+degree(1))*self%n_cells(2)
       do j = 1, degree(2)+1
          index3d(2) =  index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+degree(1))
          do i = 1, degree(1)+1
             index1d = index3d(2) + box(1)+i-1
             rho_dofs(index1d) = rho_dofs(index1d) + &
                  marker_charge * self%spline_val(i,1) *&
                  self%spline_val(j,2) * self%spline_val(k,3)
          end do
       end do
    end do
 
 
  end subroutine add_charge_single_spline_cl_3d_feec
 
 
  !---------------------------------------------------------------------------!
  subroutine add_particle_mass_spline_cl_3d_feec(self, position, marker_charge, degree, particle_mass )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(3) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_int32,                                intent( in )    :: degree(3) 
    sll_real64,                               intent( inout ) :: particle_mass(:,:) 
    !local variables
    sll_int32 :: box(3)
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k, col1, col2, col3, ind
 
 
    call convert_x_to_xbox( self, position, xi, box )
 
    self%spline_val = 0.0_f64
    if(degree(1)==self%spline_degree(1)) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    else if(degree(1)==self%spline_degree(1)-1) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    end if
 
    do j = 2, 3
       call sll_s_uniform_bsplines_eval_basis( degree(j), xi(j), self%spline_val(1:degree(j)+1,j) )
    end do
 
    !3d array
    do k = 1, degree(3)+1
       do j = 1, degree(2)+1
          do i = 1, degree(1)+1
             self%spline1_3d(i,j,k) = self%spline_val(i,1)*self%spline_val(j,2)*self%spline_val( k, 3)
          end do
       end do
    end do
 
 
    box(2:3) = box(2:3)-degree(2:3)
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, degree(3)+1       
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+degree(1))*self%n_cells(2)
       do j = 1, degree(2)+1
          index3d(2) =  index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+degree(1))
          do i = 1, degree(1)+1
             index1d = index3d(2) + box(1)+i-1
             ind = 1+(degree(1)+1-i)+(degree(2)+1-j)*(2*degree(1)+1)+ &
                  (degree(3)+1-k)*(2*degree(1)+1)*(2*degree(2)+1)
             do col3 = 1, degree(3)+1
                do col2 = 1, degree(2)+1
                   do col1 = 1, degree(1)+1                     
                      particle_mass(ind, index1d) = particle_mass( ind, index1d) +&
                           self%spline1_3d(i, j, k) * self%spline1_3d(col1, col2, col3)*&
                           marker_charge
                      ind = ind+1
                   end do
                   ind = ind+degree(1)
                end do
                ind = ind+degree(2) * (2*degree(1)+1)
             end do
          end do
       end do
    end do
 
 
  end subroutine add_particle_mass_spline_cl_3d_feec
 
 
  subroutine add_particle_mass_od_spline_cl_3d_feec(self, position, marker_charge, degree1, degree2, particle_mass )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(3) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_int32,                                intent( in )    :: degree1(3), degree2(3) 
    sll_real64,                               intent( inout ) :: particle_mass(:,:) 
    !local variables
    sll_int32 :: box(3)
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k, col1, col2, col3, ind
 
    call convert_x_to_xbox( self, position, xi, box )
 
    self%spline_0 = 0.0_f64
    self%spline1_0 = 0.0_f64
 
    if(degree1(1)==self%spline_degree(1)) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), degree1(1), xi(1), box(1), self%spline_0(1:degree1(1)+1,1) )
    else if(degree1(1)==self%spline_degree(1)-1) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), degree1(1), xi(1), box(1), self%spline_0(1:degree1(1)+1,1) )
    end if
 
    if(degree2(1)==self%spline_degree(1)) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), degree2(1), xi(1), box(1), self%spline1_0(1:degree2(1)+1,1) )
    else if(degree2(1)==self%spline_degree(1)-1) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), degree2(1), xi(1), box(1), self%spline1_0(1:degree2(1)+1,1) )
    end if
 
    do j = 2, 3
       call sll_s_uniform_bsplines_eval_basis( degree1(j), xi(j), self%spline_0(1:degree1(j)+1,j) )
       call sll_s_uniform_bsplines_eval_basis( degree2(j), xi(j), self%spline1_0(1:degree2(j)+1,j) )
    end do
 
    self%spline1_3d=0._f64
    !3d array
    do k =1, degree1(3)+1
       do j = 1, degree1(2)+1
          do i = 1, degree1(1)+1
             self%spline1_3d(i,j,k) = self%spline_0(i,1)*self%spline_0(j,2)*self%spline_0(k,3)
          end do
       end do
    end do
    self%spline2_3d=0._f64
    do k = 1, degree2(3)+1
       do j = 1, degree2(2)+1
          do i = 1, degree2(1)+1
             self%spline2_3d(i,j,k) = self%spline1_0(i,1)*self%spline1_0(j,2)*self%spline1_0(k,3)
          end do
       end do
    end do
 
    box(2:3) = box(2:3)-degree1(2:3)
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, degree1(3)+1
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+degree1(1))*self%n_cells(2)
       do j = 1, degree1(2)+1
          index3d(2) =  index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+degree1(1))
          do i = 1, degree1(1)+1
             index1d = index3d(2) + box(1)+i-1
             ind = 1+(degree1(1)+1-i)+(degree1(2)+1-j)*(degree1(1)+degree2(1)+1)+ &
                  (degree1(3)+1-k)*(degree1(1)+degree2(1)+1)*(degree1(2)+degree2(2)+1)
             do col3 = 1, degree2(3)+1
                do col2 = 1, degree2(2)+1
                   do col1 = 1, degree2(1)+1
                      particle_mass( ind, index1d) = particle_mass( ind, index1d) +&
                           self%spline1_3d(i, j, k) * self%spline2_3d(col1, col2, col3)*&
                           marker_charge
                      ind = ind+1
                   end do
                   ind = ind+degree1(1)
                end do
                ind = ind+degree1(2) * (degree1(1)+degree2(1)+1)
             end do
          end do
       end do
    end do
 
  end subroutine add_particle_mass_od_spline_cl_3d_feec
 
 
  !---------------------------------------------------------------------------!
  subroutine evaluate_field_single_spline_cl_3d_feec(self, position, degree, field_dofs, field_value)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent( inout ) :: self 
    sll_real64,                              intent( in )    :: position(3) 
    sll_int32 ,                              intent( in )    :: degree(3) 
    sll_real64,                              intent( in )    :: field_dofs(:) 
    sll_real64,                              intent( out )   :: field_value 
    !local variables
    sll_int32 :: i,j,k
    sll_int32 :: box(3)
    sll_int32 :: index3d(3)
    sll_int32 :: index1d
    sll_real64 :: xi(3)
 
    call convert_x_to_xbox( self, position, xi, box )
 
    self%spline_val = 0.0_f64
    if(degree(1)==self%spline_degree(1)) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    else if(degree(1)==self%spline_degree(1)-1) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), degree(1), xi(1), box(1), self%spline_val(1:degree(1)+1,1) )
    end if
 
    do j=2,3
       call sll_s_uniform_bsplines_eval_basis( degree(j), xi(j), self%spline_val(1:degree(j)+1,j) )
    end do
 
    field_value = 0.0_f64
    box(2:3) = box(2:3)-degree(2:3)
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, degree(3)+1       
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+degree(1))*self%n_cells(2)
       do j = 1, degree(2)+1
          index3d(2) =  index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+degree(1))
          do i = 1, degree(1)+1
             index1d = index3d(2) + box(1)+i-1
             field_value = field_value + field_dofs(index1d) * &
                  self%spline_val(i,1)*self%spline_val(j,2)*self%spline_val(k,3)
          end do
       end do
    end do
 
  end subroutine evaluate_field_single_spline_cl_3d_feec
 
  subroutine evaluate_multiple_spline_cl_3d_feec(self, position, components, field_dofs, field_value)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent( inout ) :: self 
    sll_real64,                              intent( in )    :: position(3) 
    sll_int32,                               intent(in)      :: components(:) 
    sll_real64,                              intent( in )    :: field_dofs(:,:) 
    sll_real64,                              intent(out)     :: field_value(:) 
 
 
    print*, 'evaluate_multiple_spline_cl_3d is not yet implemented'
 
  end subroutine evaluate_multiple_spline_cl_3d_feec
 
 
  !---------------------------------------------------------------------------!
  subroutine add_current_evaluate ( self, position_old, position_new, xdot, efield_dofs, j_dofs, efield_val)
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position_old(3) 
    sll_real64,                               intent( in )    :: position_new(3) 
    sll_real64,                               intent( in )    :: xdot(3) 
    sll_real64,                               intent( in )    :: efield_dofs(:) 
    sll_real64,                               intent( inout ) :: j_dofs(:) 
    sll_real64,                               intent( out   ) :: efield_val(3) 
    !local variables
    sll_real64 :: xold(3), xnew(3), xnewtilde(3), xbox(3)
    sll_int32 :: boxold(3), boxnew(3), sigma_counter(3), boxdiff(3), increment(3), box(3)
    sll_real64 :: sigma_r, sigma_l, sigma, sigma_next(3), field_value(3), weight
    sll_int32 :: j, q, bl, maxind, index_is
    type(vector) :: sigmas(3)
 
    call convert_x_to_xbox( self, position_old, xold, boxold )
    call convert_x_to_xbox( self, position_new, xnew, boxnew )
 
    ! Now we need to compute the normalized 1d line along which to integrate:
    boxdiff = boxnew-boxold
    xnewtilde = xnew + real(boxdiff,f64)
 
    sigma_l = 0.0_f64
    box = boxold ! We start in the old box
 
    sigma_counter = 0
 
    efield_val = 0.0_f64
 
    do j = 1, 3
       if (boxdiff(j) > 0 ) then
          allocate ( sigmas(j)%vals(boxdiff(j)+1) )
          do bl = 1, boxdiff(j)
             sigmas(j)%vals(bl) = (real(bl,f64)  - xold(j))/(xnewtilde(j)-xold(j))
          end do
          sigmas(j)%vals(boxdiff(j)+1) = 1.0_f64
          sigma_next(j) = sigmas(j)%vals(1)
          increment(j) = 1
       elseif (boxdiff(j) < 0 ) then
          allocate ( sigmas(j)%vals(-boxdiff(j)+1) )
          do bl = boxdiff(j)+1, 0
             sigmas(j)%vals(-bl+1) = (real(bl,f64)  - xold(j))/(xnewtilde(j)-xold(j))
          end do
          sigmas(j)%vals(-boxdiff(j)+1) = 1.0_f64
          sigma_next(j) = sigmas(j)%vals(1)
          increment(j) = -1
       else
          sigma_next(j) = 1.0_f64
          increment(j) = 0
       end if
    end do
 
    sigma_r = 0.0_f64
    do while ( sigma_r < 1.0_f64 )
       ! Identify index of next intersection
       index_is = minloc(sigma_next, dim=1)
       sigma_r = sigma_next(index_is)
 
       do q = 1, self%n_quad_points_line
          sigma = sigma_l + (sigma_r-sigma_l) * self%quad_xw_line(1,q)
          xbox = xold* (1.0_f64 - sigma) + xnewtilde * sigma  - real(sigma_counter*increment, f64)
          if (maxval(xbox) > 1.0_f64 ) then
             print*, 'box error 1', xbox, sigma, sigma_counter, increment
             xbox = min(xbox, 1._f64)
          elseif (minval(xbox) < 0.0_f64 ) then
             print*, 'box error 2', xbox, sigma, sigma_counter, increment
             xbox = max(xbox, 0._f64)
          end if
 
          weight = self%quad_xw_line(2,q)*(sigma_r-sigma_l)
 
          call point_add_eval (self, box, xbox, efield_dofs, xdot*weight, j_dofs, field_value )
          efield_val = efield_val + field_value * weight
       end do
       if (sigma_r < 1.0_f64 ) then
          ! Update the index
          sigma_counter(index_is) = sigma_counter(index_is)+1
          sigma_next(index_is) = sigmas(index_is)%vals(sigma_counter(index_is)+1)
          box(index_is) = box(index_is) + increment(index_is)
          sigma_l = sigma_r
       end if
    end do
 
  end subroutine add_current_evaluate
 
 
  subroutine point_add_eval ( self, box_in, xbox, field_dofs, weight, j_dofs, field_value )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_int32, intent(in) :: box_in(3)
    sll_real64, intent(in) :: xbox(3)
    sll_real64, intent(in) :: field_dofs(:)
    sll_real64, intent(in) :: weight(3)
    sll_real64, intent(inout) :: j_dofs(:)
    sll_real64, intent(out) :: field_value(3)
    !local variables
    sll_int32 :: i,j,k, index1d
    sll_int32 :: box(3), index3d(3), deg
    sll_real64 :: spval
 
    field_value = 0.0_f64
 
    deg=self%spline_degree(1)
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), deg, xbox(1), box_in(1), self%spline_0(1:deg+1,1) )
 
    deg=self%spline_degree(1)-1
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), deg, xbox(1), box_in(1), self%spline_1(1:deg+1,1) )
 
    do j = 2, 3
       call sll_s_uniform_bsplines_eval_basis( self%spline_degree(j), xbox(j), &
            self%spline_0(1:self%spline_degree(j)+1,j) )
       call sll_s_uniform_bsplines_eval_basis( self%spline_degree(j)-1, xbox(j), &
            self%spline_1(1:self%spline_degree(j),j) )
    end do
 
    box(1) = box_in(1)
    box(2:3) = box_in(2:3)-self%spline_degree(2:3)
 
    deg = self%spline_degree(1)-1
 
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 )      
    do k = 1, self%spline_degree(3)+1 
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+deg)*self%n_cells(2)
       do j = 1, self%spline_degree(2)+1
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+deg)
          do i = 1, self%spline_degree(1)
             index1d = index3d(2) + box(1)+i-1
 
             spval =  self%spline_1(i,1) * &
                  self%spline_0(j,2) * self%spline_0(k,3)
             field_value(1) = field_value(1) + field_dofs(index1d) * spval
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(1) * spval
          end do
       end do
    end do
 
    box(2) = box(2)+1
    deg = self%spline_degree(1)
    call box_index( self, box(2), 2 )   
    do k = 1, self%spline_degree(3)+1       
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+deg)*self%n_cells(2)
       do j = 1, self%spline_degree(2)
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+deg)
          do i = 1, self%spline_degree(1)+1
             index1d = index3d(2) + box(1)+i-1 + self%n_total1 
 
             spval = self%spline_0(i,1) * &
                  self%spline_1(j,2) * self%spline_0(k,3)
             field_value(2) = field_value(2) + field_dofs(index1d ) * spval
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(2) * spval
          end do
       end do
    end do
 
    box(2) = box(2)-1
    box(3) = box(3)+1
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, self%spline_degree(3)    
       index3d(3) = (self%index1d(k,3)-1)*(self%n_cells(1)+deg)*self%n_cells(2)
       do j = 1, self%spline_degree(2)+1
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*(self%n_cells(1)+deg)
          do i = 1, self%spline_degree(1)+1
             index1d = index3d(2) + box(1)+i-1 + self%n_total1+self%n_total0
 
             spval = self%spline_0(i,1) * &
                  self%spline_0(j,2) * self%spline_1(k,3)
             field_value(3) = field_value(3) + field_dofs(index1d ) * spval
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(3) * spval
          end do
       end do
    end do
 
 
  end subroutine point_add_eval
 
 
  subroutine add_current_cl_3d  ( self, position_old, position_new, xdot, j_dofs) 
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in    ) :: position_old(3) 
    sll_real64,                               intent( in    ) :: position_new(3) 
    sll_real64,                               intent( in    ) :: xdot(3) 
    sll_real64,                               intent( inout ) :: j_dofs(:) 
    !local variables
    sll_real64 :: xold(3), xnew(3), xnewtilde(3), xbox(3)
    sll_int32 :: boxold(3), boxnew(3), sigma_counter(3), boxdiff(3), increment(3), box(3)
    sll_real64 :: sigma_r, sigma_l, sigma, sigma_next(3), field_value(3), weight
    sll_int32 :: j, q, bl, maxind, index_is
    type(vector) :: sigmas(3)
 
    call convert_x_to_xbox( self, position_old, xold, boxold )
    call convert_x_to_xbox( self, position_new, xnew, boxnew )
 
    ! Now we need to compute the normalized 1d line along which to integrate:
    boxdiff = boxnew-boxold
    xnewtilde = xnew + real(boxdiff,f64)
    sigma_l = 0.0_f64
    box = boxold ! We start in the old box
    sigma_counter = 0
 
    do j = 1, 3
       if (boxdiff(j) > 0 ) then
          allocate ( sigmas(j)%vals(boxdiff(j)+1) )
          do bl = 1, boxdiff(j)
             sigmas(j)%vals(bl) = (real(bl,f64)  - xold(j))/(xnewtilde(j)-xold(j))
          end do
          sigmas(j)%vals(boxdiff(j)+1) = 1.0_f64
          sigma_next(j) = sigmas(j)%vals(1)
          increment(j) = 1
       elseif (boxdiff(j) < 0 ) then
          allocate ( sigmas(j)%vals(-boxdiff(j)+1) )
          do bl = boxdiff(j)+1, 0
             sigmas(j)%vals(-bl+1) = (real(bl,f64)  - xold(j))/(xnewtilde(j)-xold(j))
          end do
          sigmas(j)%vals(-boxdiff(j)+1) = 1.0_f64
          sigma_next(j) = sigmas(j)%vals(1)
          increment(j) = -1
       else
          sigma_next(j) = 1.0_f64
          increment(j) = 0
       end if
    end do
 
    sigma_r = 0.0_f64
    do while ( sigma_r < 1.0_f64 )
       ! Identify index of next intersection
       index_is = minloc(sigma_next, dim=1)
       sigma_r = sigma_next(index_is)
 
       do q = 1, self%n_quad_points_line
          sigma = sigma_l + (sigma_r-sigma_l) * self%quad_xw_line(1,q)
          xbox = xold* (1.0_f64 - sigma) + xnewtilde * sigma  - real(sigma_counter*increment, f64)
          if (maxval(xbox) > 1.0_f64 ) then
             print*, 'box error1', xbox, sigma, sigma_counter, increment
             xbox = min(xbox, 1._f64)
          elseif (minval(xbox) < 0.0_f64 ) then
             print*, 'box error2', xbox, sigma, sigma_counter, increment
             xbox = max(xbox, 0._f64)
          end if
 
          weight = self%quad_xw_line(2,q)*(sigma_r-sigma_l)
          call integrate_spline_cl_3d(self, box, xbox, xdot*weight, j_dofs )
       end do
       if (sigma_r < 1.0_f64 ) then
          ! Update the index
          sigma_counter(index_is) = sigma_counter(index_is)+1
          sigma_next(index_is) = sigmas(index_is)%vals(sigma_counter(index_is)+1)
          box(index_is) = box(index_is) + increment(index_is)
          sigma_l = sigma_r
       end if
    end do
  end subroutine add_current_cl_3d
 
 
  subroutine integrate_spline_cl_3d( self, box_in, xbox, weight, j_dofs )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_int32,  intent( in    ) :: box_in(3)
    sll_real64, intent( in    ) :: xbox(3)
    sll_real64, intent( in    ) :: weight(3)
    sll_real64, intent( inout ) :: j_dofs(:)
    !local variables
    sll_int32 :: i,j,k, index1d
    sll_int32 :: box(3), index3d(3), deg
 
    deg=self%spline_degree(1)
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), deg, xbox(1), box_in(1), self%spline_0(1:deg+1,1)  )
 
    deg=self%spline_degree(1)-1
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), deg, xbox(1), box_in(1), self%spline_1(1:deg+1,1) )
 
 
    do j = 2, 3
       call sll_s_uniform_bsplines_eval_basis( self%spline_degree(j), xbox(j), &
            self%spline_0(1:self%spline_degree(j)+1,j) )
       call sll_s_uniform_bsplines_eval_basis( self%spline_degree(j)-1, xbox(j), &
            self%spline_1(1:self%spline_degree(j),j) )
    end do
 
    box(1) = box_in(1)
    box(2:3) = box_in(2:3)-self%spline_degree(2:3)
 
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, self%spline_degree(3)+1 
       index3d(3) = (self%index1d(k,3)-1)*(self%n_dofs(1)-1)*self%n_cells(2)
       do j = 1, self%spline_degree(2)+1
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*(self%n_dofs(1)-1)
          do i = 1, self%spline_degree(1)
             index1d = index3d(2) + box(1)+i-1
 
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(1) * self%spline_1(i,1) * &
                  self%spline_0(j,2) * self%spline_0(k,3)
          end do
       end do
    end do
 
    box(2) = box(2)+1
    call box_index( self, box(2), 2 )   
    do k = 1, self%spline_degree(3)+1       
       index3d(3) = (self%index1d(k,3)-1)*self%n_dofs(1)*self%n_cells(2)
       do j = 1, self%spline_degree(2)
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*self%n_dofs(1)
          do i = 1, self%spline_degree(1)+1
             index1d = index3d(2) + box(1)+i-1 + self%n_total1
 
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(2) * self%spline_0(i,1) * &
                  self%spline_1(j,2) * self%spline_0(k,3)
          end do
       end do
    end do
 
    box(2) = box(2)-1
    box(3) = box(3)+1
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 )
    do k = 1, self%spline_degree(3)    
       index3d(3) = (self%index1d(k,3)-1)*self%n_dofs(1)*self%n_cells(2)
       do j = 1, self%spline_degree(2)+1
          index3d(2) = index3d(3) + (self%index1d(j,2)-1)*self%n_dofs(1)
          do i = 1, self%spline_degree(1)+1
             index1d = index3d(2) + box(1)+i-1 + self%n_total1+self%n_total0
             j_dofs(index1d) = j_dofs(index1d) + &
                  weight(3) * self%spline_0(i,1) * &
                  self%spline_0(j,2) * self%spline_1(k,3)
          end do
       end do
    end do
 
 
  end subroutine integrate_spline_cl_3d
 
 
  subroutine add_current_update_v_component1_spline_3d_feec (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_real64, intent(in)    :: position_old(3) 
    sll_real64, intent(in)    :: position_new 
    sll_real64, intent(in)    :: marker_charge 
    sll_real64, intent(in)    :: qoverm 
    sll_real64, intent(in)    :: bfield_dofs(:) 
    sll_real64, intent(inout) :: vi(3) 
    sll_real64, intent(inout) :: j_dofs(:) 
    !local variables
    sll_int32 :: box(3), boxnew, boxold
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k
    sll_real64 :: xnew
    sll_real64 :: vt(2), vtt2, vtt3
    sll_int32  :: start1, start2
    sll_int32  :: local_size
    sll_int32  :: startjk, stride(2)
    sll_real64 :: splinejk
 
    start1 = self%n_total0+self%n_total1
    start2 = self%n_total0
    stride = 1
 
    call convert_x_to_xbox( self, position_old, xi, box )
    call convert_x_to_xbox_1d( self, 1, position_new, xnew, boxnew )
    boxold = box(1)
 
    ! Achtung wir brauchen nun splines von beidem Grad
    self%spline_0 = 0.0_f64
    self%spline_1 = 0.0_f64
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_0%spline2, self%spline_0(1:self%spline_degree(2)+1,2), self%spline_degree(2), xi(2))
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_1%spline2, self%spline_1(1:self%spline_degree(2),2), self%spline_degree(2)-1, xi(2))
 
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_0%spline3, self%spline_0(1:self%spline_degree(3)+1,3), self%spline_degree(3), xi(3))
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_1%spline3, self%spline_1(1:self%spline_degree(3),3), self%spline_degree(3)-1, xi(3))
 
    box(2:3) = box(2:3) - self%spline_degree(2:3)
    local_size = abs(boxnew-boxold)+self%spline_degree(1)
 
    call box_index( self, box(2), 2 )   
    call box_index( self, box(3), 3 ) 
    do k = 1, self%spline_degree(3)+1      
       vtt2 = 0.0_f64
       vtt3 = 0.0_f64
       index3d(3) = self%index1d(k,3)
       do j = 1, self%spline_degree(2)+1
          index3d(2) = self%index1d(j,2)
 
          startjk  = 1 + (index3d(2)-1)*(self%n_dofs(1)-1) + &
               (index3d(3)-1)*(self%n_dofs(1)-1)*self%n_cells(2)
          splinejk = self%spline_0(j, 2) * self%spline_0(k,3)
 
          vt = 0.0_f64
          self%j1d = 0.0_f64
          call  add_current_1d(self, 1, xi(1), boxold, xnew, boxnew, marker_charge, bfield_dofs, start1+startjk, start2+startjk, stride, vt, self%j1d)
          if (j>1) then
             vtt2 = vtt2 + vt(1)*self%spline_1(j-1, 2)
          end if
          vtt3 = vtt3 + vt(2)*self%spline_0(j, 2)
 
          do i = 1, self%spline_degree(1) !local_size?
             index3d(1) = box(1)+i-1
             index1d = startjk+index3d(1)-1
             j_dofs(index1d) = j_dofs(index1d) + self%j1d(index3d(1)) * &
                  splinejk 
          end do
       end do
       vi(2) = vi(2) - qoverm*vtt2*self%spline_0(k, 3)
       if ( k> 1) then
          vi(3) = vi(3) + qoverm*vtt3*self%spline_1(k-1, 3)
       end if
    end do
 
  end subroutine add_current_update_v_component1_spline_3d_feec
 
 
  subroutine add_current_update_v_component2_spline_3d_feec (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_real64, intent(in)    :: position_old(3) 
    sll_real64, intent(in)    :: position_new 
    sll_real64, intent(in)    :: marker_charge 
    sll_real64, intent(in)    :: qoverm 
    sll_real64, intent(in)    :: bfield_dofs(:) 
    sll_real64, intent(inout) :: vi(3) 
    sll_real64, intent(inout) :: j_dofs(:) 
    !local variables
    sll_int32 :: box(3), boxnew, boxold
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k
    sll_real64 :: xnew
    sll_real64 :: vt(2), vtt2, vtt3
    sll_int32  :: start1, start2
    sll_int32  :: local_size
    sll_int32 :: stride(2), startjk0, startjk1
    sll_real64 :: splineik
 
    start1 = self%n_total0+self%n_total1
    start2 = 0
    stride(1) = self%n_dofs(1)-1
    stride(2) = self%n_dofs(1)
 
    call convert_x_to_xbox( self, position_old, xi, box )
    call convert_x_to_xbox_1d( self, 2, position_new, xnew, boxnew )
    boxold = box(2)
 
 
    ! Achtung wir brauchen nun splines von beidem Grad
    self%spline_0 = 0.0_f64
    self%spline_1 = 0.0_f64
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), &
         self%spline_degree(1), xi(1), box(1), &
         self%spline_0(1:self%spline_degree(1)+1,1) )
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), &
         self%spline_degree(1)-1, xi(1), box(1), &
         self%spline_1(1:self%spline_degree(1),1) )
 
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_0%spline3, self%spline_0(1:self%spline_degree(3)+1,3), self%spline_degree(3), xi(3))
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_1%spline3, self%spline_1(1:self%spline_degree(3),3), self%spline_degree(3)-1, xi(3))
 
    box(3) = box(3) - self%spline_degree(3)
    local_size = abs(boxnew-boxold)+self%spline_degree(2)
    ! Define the range of the first component
    if (boxold<boxnew) then
       box(2) = boxold-self%spline_degree(2)+1
    else
       box(2) = boxnew-self%spline_degree(2)+1
    end if
 
    call box_index( self, box(3), 3 ) 
    do k = 1, self%spline_degree(3)+1      
       vtt2 = 0.0_f64
       vtt3 = 0.0_f64
       index3d(3) = self%index1d(k,3)
       do j = 1, self%spline_degree(1)+1
          index3d(1) = box(1)+j-1
 
          startjk0  = index3d(1) + &
               (index3d(3)-1)*self%n_dofs(1)*self%n_cells(2)
          startjk1  = index3d(1)-1 + &
               (index3d(3)-1)*(self%n_dofs(1)-1)*self%n_cells(2)
          splineik = self%spline_0(j, 1) * self%spline_0(k,3) 
 
          vt = 0.0_f64
          self%j1d = 0.0_f64
          call  add_current_1d(self, 2, xi(2), boxold-1, xnew, boxnew-1, marker_charge, bfield_dofs,  start1+startjk1, start2+startjk0, stride, vt, self%j1d)
 
          if (j>1) then
             vtt2 = vtt2 + vt(1)*self%spline_1(j-1, 1)
          end if
          vtt3 = vtt3 + vt(2)*self%spline_0(j, 1)
 
          do i=1,local_size
             index3d(2) = modulo(box(2)+i-2,self%n_cells(2))+1
             index1d = startjk0+(index3d(2)-1)*stride(2)
             j_dofs(index1d) = j_dofs(index1d) + self%j1d(index3d(2)) * &
                  splineik
 
          end do
       end do
       vi(1) = vi(1) + qoverm*vtt2*self%spline_0(k, 3)
       if ( k> 1) then
          vi(3) = vi(3) - qoverm*vtt3*self%spline_1(k-1, 3)
       end if
    end do
 
 
  end subroutine add_current_update_v_component2_spline_3d_feec
 
 
  subroutine add_current_update_v_component3_spline_3d_feec (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_real64, intent(in)    :: position_old(3) 
    sll_real64, intent(in)    :: position_new 
    sll_real64, intent(in)    :: marker_charge 
    sll_real64, intent(in)    :: qoverm 
    sll_real64, intent(in)    :: bfield_dofs(:) 
    sll_real64, intent(inout) :: vi(3) 
    sll_real64, intent(inout) :: j_dofs(:) 
    !local variables
    sll_int32 :: box(3), boxnew, boxold
    sll_real64 :: xi(3)
    sll_int32  :: index3d(3)
    sll_int32  :: index1d
    sll_int32  :: i,j,k
    sll_real64 :: xnew
    sll_real64 :: vt(2), vtt2, vtt3
    sll_int32  :: start1, start2
    sll_int32  :: local_size
    sll_int32  :: stride(2), startjk0, startjk1
    sll_real64 :: splineij
 
    start1 = self%n_total0
    start2 = 0
    stride(1) = (self%n_dofs(1)-1)*self%n_cells(2)
    stride(2) = self%n_dofs(1)*self%n_cells(2)
 
    call convert_x_to_xbox( self, position_old, xi, box )
    call convert_x_to_xbox_1d( self, 3, position_new, xnew, boxnew )
    boxold = box(3)
 
    ! Achtung wir brauchen nun splines von beidem Grad
    self%spline_0 = 0.0_f64
    self%spline_1 = 0.0_f64
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_0%spline1, self%n_cells(1), &
         self%spline_degree(1), xi(1), box(1), &
         self%spline_0(1:self%spline_degree(1)+1,1) )
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, self%n_cells(1), &
         self%spline_degree(1)-1, xi(1), box(1), &
         self%spline_1(1:self%spline_degree(1),1) )
 
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_0%spline2, self%spline_0(1:self%spline_degree(2)+1,2), self%spline_degree(2), xi(2))
    call sll_s_spline_pp_horner_m_1d(self%spline_pp_1%spline2, self%spline_1(1:self%spline_degree(2),2), self%spline_degree(2)-1, xi(2))
 
    box(2) = box(2) -  self%spline_degree(2)
    local_size = abs(boxnew-boxold)+self%spline_degree(3)
    ! Define the range of the third component
    if (boxold<boxnew) then
       box(3) = boxold-self%spline_degree(3)+1
    else
       box(3) = boxnew-self%spline_degree(3)+1
    end if
 
    call box_index( self, box(2), 2 )
    do k = 1, self%spline_degree(2)+1      
       vtt2 = 0.0_f64
       vtt3 = 0.0_f64
       index3d(2) = self%index1d(k,2)
       do j = 1, self%spline_degree(1)+1
          index3d(1) = box(1)+j-1
 
          startjk0  = index3d(1) + (index3d(2)-1)*self%n_dofs(1)
          startjk1  = index3d(1)-1 + (index3d(2)-1)*(self%n_dofs(1)-1)
          splineij = self%spline_0(j, 1) * self%spline_0(k,2)
 
          vt = 0.0_f64
          self%j1d = 0.0_f64
          call  add_current_1d(self, 3, xi(3), boxold-1, xnew, boxnew-1, marker_charge, bfield_dofs, start1+startjk1, start2+startjk0, stride, vt, self%j1d)
          if (j>1) then
             vtt2 = vtt2 + vt(1)*self%spline_1(j-1, 1)
          end if
          vtt3 = vtt3 + vt(2)*self%spline_0(j, 1)
 
          do i = 1, local_size
             index3d(3) = modulo(box(3)+i-2,self%n_cells(3))+1
             index1d = startjk0 + (index3d(3)-1)*stride(2)
 
             j_dofs(index1d) = j_dofs(index1d) + self%j1d(index3d(3)) * &
                  splineij  
          end do
 
       end do
       vi(1) = vi(1) - qoverm*vtt2*self%spline_0(k, 2)
       if ( k> 1) then
          vi(2) = vi(2) + qoverm*vtt3*self%spline_1(k-1, 2)
       end if
    end do
 
  end subroutine add_current_update_v_component3_spline_3d_feec
 
 
  subroutine add_current_1d (self, component, r_old, index_old,  r_new, index_new, marker_charge, bfield_dofs, start1, start2, stride, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_int32,  intent(in)    :: component
    sll_real64, intent(in)    :: r_old
    sll_int32,  intent(in)    :: index_old
    sll_real64, intent(in)    :: r_new
    sll_int32,  intent(in)    :: index_new
    sll_real64, intent(in)    :: marker_charge 
    sll_real64, intent(in)    :: bfield_dofs(self%n_total0+self%n_total1*2) 
    sll_int32,                              intent(in)    :: start1 
    sll_int32,                              intent(in)    :: start2 
    sll_int32,                              intent(in)    :: stride(2) 
    sll_real64, intent(inout) :: vi(2) 
    sll_real64, intent(inout) :: j_dofs(:) 
    ! local variables
    sll_int32  ::  ind
 
    if (index_old == index_new) then
       if (r_old < r_new) then
          call update_jv(self, component, r_old, r_new, index_old, marker_charge, &
               1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       else
          call update_jv(self, component, r_new, r_old, index_old, marker_charge,  &
               -1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       end if
    elseif (index_old < index_new) then
       call update_jv (self, component, r_old, 1.0_f64, index_old, marker_charge, &
            1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       call update_jv (self, component, 0.0_f64, r_new, index_new, marker_charge, &
            1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       do ind = index_old+1, index_new-1
          call update_jv (self, component, 0.0_f64, 1.0_f64, ind, marker_charge, &
               1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       end do
    else
       call update_jv (self, component, r_new, 1.0_f64, index_new, marker_charge, &
            -1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       call update_jv (self, component, 0.0_f64, r_old, index_old, marker_charge, &
            -1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       do ind = index_new+1, index_old-1
          call update_jv (self, component, 0.0_f64, 1.0_f64, ind, marker_charge, &
               -1.0_f64, bfield_dofs, start1, start2, stride, vi, j_dofs)
       end do
    end if
 
 
  end subroutine add_current_1d
 
 
  subroutine update_jv(self, component, lower, upper, index, marker_charge, sign, bfield_dofs, start1, start2, stride, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_int32,                              intent(in)    :: component
    sll_real64,                             intent(in)    :: lower
    sll_real64,                             intent(in)    :: upper
    sll_int32,                              intent(in)    :: index
    sll_real64,                             intent(in)    :: marker_charge
    sll_real64,                             intent(in)    :: sign
    sll_real64,                             intent(inout) :: vi(2)
    sll_real64,                             intent(in)    :: bfield_dofs(self%n_total0+self%n_total1*2)
    sll_int32,                              intent(in)    :: start1
    sll_int32,                              intent(in)    :: start2
    sll_int32,                              intent(in)    :: stride(2)
    sll_real64,                             intent(inout) :: j_dofs(:)
    !Local variables
    sll_int32  :: ind, i_grid, i_mod, j
    sll_real64 :: c1, c2
    sll_int32 :: spline_degree
    sll_int32 :: n_cells
 
    self%spline_val(:,1) = 0._f64
    self%spline_val_more(:,1) = 0._f64
 
    n_cells = self%n_cells(component)
    spline_degree = self%spline_degree(component)-1
 
    c1 =  0.5_f64*(upper-lower)
    c2 =  0.5_f64*(upper+lower)
 
    if( component == 1 ) then
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, n_cells, &
            spline_degree, c1*self%quad_xw(1,1)+c2, index, &
            self%spline_val(1:spline_degree+1,1) )
       self%spline_val(:,1) = self%spline_val(:,1) * (self%quad_xw(2,1)*c1)
       do j=2,self%n_quad_points
          call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp_1%spline1, n_cells, &
               spline_degree, c1*self%quad_xw(1,j)+c2, index, &
               self%spline_val_more(1:spline_degree+1,1))
          self%spline_val(:,1) = self%spline_val(:,1) + self%spline_val_more(:,1) * (self%quad_xw(2,j)*c1)
       end do
       self%spline_val(:,1) = self%spline_val(:,1) * (sign*self%delta_x(component))
 
       ind = 1
       do i_grid = index-1, index-1+spline_degree
          j_dofs(i_grid+1) = j_dofs(i_grid+1) + &
               (marker_charge*self%spline_val(ind,1))
          vi(1) = vi(1) + self%spline_val(ind,1)*bfield_dofs(i_grid*stride(1)+start1)
          vi(2) = vi(2) + self%spline_val(ind,1)*bfield_dofs(i_grid*stride(2)+start2)
          ind = ind + 1
       end do
    else
 
       call sll_s_uniform_bsplines_eval_basis(spline_degree, c1*self%quad_xw(1,1)+c2, &
            self%spline_val(1:spline_degree,1))
       self%spline_val(:,1) = self%spline_val(:,1) * (self%quad_xw(2,1)*c1)
       do j=2,self%n_quad_points
          call sll_s_uniform_bsplines_eval_basis(spline_degree, c1*self%quad_xw(1,j)+c2, &
               self%spline_val_more(1:spline_degree,1))
          self%spline_val(:,1) = self%spline_val(:,1) + self%spline_val_more(:,1) * (self%quad_xw(2,j)*c1)
       end do
       self%spline_val(:,1) = self%spline_val(:,1) * (sign*self%delta_x(component))
 
       ind = 1
       do i_grid = index - spline_degree, index
          i_mod = modulo(i_grid, n_cells )
          j_dofs(i_mod+1) = j_dofs(i_mod+1) + &
               (marker_charge*self%spline_val(ind,1))
          vi(1) = vi(1) + self%spline_val(ind,1)*bfield_dofs(i_mod*stride(1)+start1)
          vi(2) = vi(2) + self%spline_val(ind,1)*bfield_dofs(i_mod*stride(2)+start2)
          ind = ind + 1
       end do
    end if
 
  end subroutine update_jv
 
  
  function convert_index_3d_to_1d( index3d, n_cells, degree ) result( index1d )
    sll_int32, intent( in ) :: index3d(3)
    sll_int32, intent( in ) :: n_cells(3)
    sll_int32 :: degree
    sll_int32 :: index1d
 
    index1d = index3d(1) + (index3d(2)-1)*(n_cells(1)+degree) + (index3d(3)-1)*(n_cells(1)+degree)*n_cells(2)
 
  end function convert_index_3d_to_1d
 
  
  subroutine convert_x_to_xbox( self, position, xi, box )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(3) 
    sll_real64,                               intent( out )    :: xi(3) 
    sll_int32,                                intent( out )    :: box(3) 
 
    xi = (position - self%domain(:,1)) * self%rdelta_x
    box = floor( xi ) + 1
    xi = xi - real(box-1, f64)
 
    if( box(1) == self%n_cells(1) + 1 ) then
       if( xi(1) == 0._f64)then 
          xi(1) = 1._f64
          box(1) = self%n_cells(1)
       else
          print*, 'box,x, xbox', box, position(1), xi(1)
          sll_error('convert_x_to_xbox', 'box1 value to high' ) 
       end if
    else if( box(1) > self%n_cells(1) + 1 ) then
       print*, 'box,x, xbox', box, position(1), xi(1)
       sll_error('convert_x_to_xbox', 'box1 value to high' ) 
    else if( box(1) <= 0 ) then
       print*, 'box,x, xbox', box, position(1), xi(1)
       sll_error('convert_x_to_xbox', 'box1 value to low' ) 
    end if
 
  end subroutine convert_x_to_xbox
 
  
  subroutine convert_x_to_xbox_1d( self, component, position, xi, box )
    class( sll_t_particle_mesh_coupling_spline_cl_3d_feec ), intent(inout)   :: self
    sll_int32,                                intent( in )    :: component
    sll_real64,                               intent( in )    :: position 
    sll_real64,                               intent( out )    :: xi 
    sll_int32,                                intent( out )    :: box 
 
    xi = (position - self%domain(component,1))/self%delta_x(component)!* self%rdelta_x ! destroys perfect gauss conservation for symplectic splitting with fft solver
    box = floor( xi ) + 1
    xi = xi - real(box-1, f64)
 
    if( component ==1 .and. box == self%n_cells(1) + 1 ) then
       if( xi == 0._f64)then 
          xi = 1._f64
          box = self%n_cells(1)
       else
          print*, 'box,x, xbox', box, position, xi
          sll_error('convert_x_to_xbox', 'box1 value to high' ) 
       end if
    end if
 
  end subroutine convert_x_to_xbox_1d
 
  
  subroutine box_index( self, box, comp )
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent(inout) :: self
    sll_int32, intent(in) :: box
    sll_int32, intent(in) :: comp
    !local variables
    sll_int32 :: j
 
    do j = 1, self%spline_degree(comp)+1
       self%index1d(j,comp) = modulo(box+j-2,self%n_cells(comp))+1
    end do
 
  end subroutine box_index
 
  subroutine sll_s_uniform_bsplines_eval_basis_clamped( spline, n_cells, degree, xi, box, spline_val )
    type(sll_t_spline_pp_1d), intent( in    ) :: spline
    sll_int32,                intent( in    ) :: n_cells
    sll_int32,                intent( in    ) :: degree
    sll_real64,               intent( in    ) :: xi
    sll_int32,                intent( in    ) :: box
    sll_real64,               intent(   out ) :: spline_val(:)
    !local variables
    sll_int32 :: i
 
    if(box >= 1 .and. box <= degree-1)then
       do i=1, degree+1
          spline_val(i) = sll_f_spline_pp_horner_1d( degree, spline%poly_coeffs_boundary_left(:,:,box), xi, i)
       end do
    else if (box >= degree .and. box <= n_cells-degree+1)then
       call sll_s_uniform_bsplines_eval_basis( degree, xi, spline_val )
    else if(box >= n_cells-degree+2 .and. box <= n_cells)then
       do i=1, degree+1
          spline_val(i) = sll_f_spline_pp_horner_1d( degree, spline%poly_coeffs_boundary_right(:,:,box-1-n_cells+degree), xi, i)
       end do
    else
       print*, 'box, xbox', box, xi
       sll_error( "uniform_bsplines_eval_basis_clamped", "box out of range" )
    end if
 
  end subroutine sll_s_uniform_bsplines_eval_basis_clamped
 
 
  subroutine init_spline_cl_3d_feec ( self, n_cells, domain, spline_degree, boundary, no_particles )
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent( out ) :: self 
    sll_int32,                               intent(in) :: n_cells(3) 
    sll_real64,                              intent(in) :: domain(3,2) 
    sll_int32,                               intent(in) :: spline_degree(3) 
    sll_int32,                               intent(in) :: boundary(3) 
    sll_int32, optional,                     intent(in) :: no_particles 
    !local variables
    sll_int32 :: maxspan, j, maxind
 
 
    ! Store cells information
    self%domain = domain
    self%n_cells = n_cells
    self%n_dofs = n_cells + spline_degree
    self%n_total = product(n_cells)
    self%n_total0 = (n_cells(1)+spline_degree(1))*n_cells(2)*n_cells(3)
    self%n_total1 = (n_cells(1)+spline_degree(1)-1)*n_cells(2)*n_cells(3)
    self%delta_x = (self%domain(:,2)-self%domain(:,1))/real(n_cells, f64)
    self%rdelta_x = 1.0_f64/self%delta_x
 
    ! Store basis function information
    if( present(no_particles) )then
       self%no_particles = no_particles
    else
       self%no_particles = 1
    end if
 
    ! Initialize information on the spline
    self%spline_degree = spline_degree
    maxspan = maxval(spline_degree + 1)
 
    self%n_quad_points = (maxval(self%spline_degree)+2)/2
    allocate( self%quad_xw(2,self%n_quad_points) )
    ! normalized Gauss Legendre points and weights
    self%quad_xw = sll_f_gauss_legendre_points_and_weights(self%n_quad_points)
 
    ! For the line integral
    self%n_quad_points_line = ceiling(real(sum(self%spline_degree), f64)*0.5_f64)
    allocate( self%quad_xw_line(2,self%n_quad_points_line) )
    ! normalized Gauss Legendre points and weights
    self%quad_xw_line = sll_f_gauss_legendre_points_and_weights(self%n_quad_points_line)
 
    self%quad_xw_line(1,:) = 0.5_f64*(self%quad_xw_line(1,:)+1.0_f64)
    self%quad_xw_line(2,:) = 0.5_f64*self%quad_xw_line(2,:)
 
    allocate( self%spline_val(maxspan,3) )
    allocate( self%spline_0(maxspan,3) )
    allocate( self%spline1_0(maxspan,3) )
    allocate( self%spline_1(maxspan,3) )
    allocate( self%spline_2d(maxspan, maxspan) )
    allocate( self%spline1_2d(1:maxspan, 1:maxspan) )
    allocate( self%spline1_3d(1:maxspan, 1:maxspan, 1:maxspan) )
    allocate( self%spline2_3d(1:maxspan, 1:maxspan, 1:maxspan) )
    allocate( self%index1d(maxspan, 3) )
 
    allocate( self%spline_val_more(maxspan,3) )
    allocate( self%j1d( maxval(self%n_dofs)+maxval(self%spline_degree) ))
 
    call sll_s_spline_pp_init_3d(self%spline_pp_0, spline_degree, n_cells, boundary)
    call sll_s_spline_pp_init_3d(self%spline_pp_1, spline_degree-1, n_cells, boundary)
 
  end subroutine init_spline_cl_3d_feec
 
 
  subroutine free_spline_cl_3d_feec(self)
    class(sll_t_particle_mesh_coupling_spline_cl_3d_feec), intent( inout ) :: self 
 
    deallocate( self%spline_val)
    deallocate( self%spline_0)
    deallocate( self%spline1_0 )
    deallocate( self%spline_1 )
    deallocate( self%spline_2d )
    deallocate( self%spline1_2d )
    deallocate( self%spline1_3d )
    deallocate( self%spline2_3d )
 
    call sll_s_spline_pp_free_3d( self%spline_pp_0 )
    call sll_s_spline_pp_free_3d( self%spline_pp_1 )
 
  end subroutine free_spline_cl_3d_feec
 
 
end module sll_m_particle_mesh_coupling_spline_cl_3d_feec