//
// Created by jundafeng on 3/23/22.
//

#ifndef DGTD_KERNELS_CUH
  #define DGTD_KERNELS_CUH

  #include "Constants.h"
  #include "../dgtd-performance.hpp"
  
  typedef struct{
    // Plane Wave Excitation
    int ExcitationFlag; // Scattering Planewave Waveform: (0)Time Harmonic (1)Gauss (2)Neumann (3) Modulated Gauss
    fp_t_ts to; //Tdelay: Delay of excitation pulse (sec)
    fp_t_ts tau; //Tau: Parameter for pulse width (sec)
    fp_t_ts freq_m; //Freq: Central Modulation Frequency (MHz)
    fp_t_ts Emagnitude; //magnitude of the E field
    fp_t_ts theta; //angle of excitation
    fp_t_ts phi; //angle of excitation
    fp_t_ts Epol[3]; //Efield direction
    fp_t_ts ro[3]; //position where the excitation starts
    
    // Port Excitation
    int TimeDistributionFlag; // Port Waveform (0)TH (1)Gauss
    int PortFlag;
    fp_t_ts ro_port[3];
    fp_t_ts r1_port[3];
    fp_t_ts PortImpedance;
    fp_t_ts PortDirection[3];
    
  } ExcitationProp;

  // Ceiling funciton for X / Y.
  __host__ __device__ static inline int ceil_div(int x, int y){return (x - 1) / y + 1;}
  __host__ static inline size_t ceil_div(size_t x, size_t y){return (x - 1) / y + 1;}
  __host__ static inline size_t floor_div(size_t x, size_t y){return (x - 1) / y;}

  constexpr int BLOCK_SIZE_Matrix_EVAL = 128;

  #if defined(DGTD_USE_CUDA)
    #include <cuda.h>
    #include <cuda_runtime.h>
    #include <cooperative_groups.h>
    #include "cuda_utils.h"
    #include <cuComplex.h>
    
    //#include <cooperative_groups/memcpy_async.h>
    //#include <cuda/pipeline>
    namespace cg = cooperative_groups;
    #define WARP_SIZE 32
    __global__ void makeNeighField(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements);
    __global__ void checker(fp_t_ts* a, fp_t_ts* b, int length);
    __global__ void addCouplingResults(fp_t_ts* Field, fp_t_ts* auxiliar, int* neighs, int* neighsOffset, int elements);

    __global__ void makeNeighFieldRegular(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int polyOrder);

    __global__ void makeNeighFieldRegular_custom(int* MappingArray, fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int polyOrder, int numFaces);

    __global__ void addCouplingResultsRegular(fp_t_ts* Field, fp_t_ts* auxiliar, int elements);

    __global__ void addCouplingResultsRegular_custom(fp_t_ts* Field, fp_t_ts* auxiliar, int elements, int numFaces);

    // Testing
    __global__ void fill_ones_cols(fp_t_ts* __restrict__ M, int rows, int cols, int ld);
    

    __global__ void addExcitationH_PML( int* ExcitationFacesCnt, 
                                    int* ExcitationFacesOffset, 
                                    int* ExcitationFacesNum, 
                                    int excitationCNT,
                                    int excitationPMLCNT,
                                    int8_t* map,
                                    ExcitationProp exci_prop,
                                    int PolyFlag,
                                    fp_t_ts dt_muo, fp_t_ts t,
                                    fp_t_ts* nd_coords_tet, 
                                    fp_t_ts* nd_coords_face, 
                                    fp_t_ts* Z_face,
                                    fp_t_ts* InvMat,
                                    fp_t_ts* Field);

    __global__ void addExcitationE_PML( int* ExcitationFacesCnt,
                                    int* ExcitationFacesOffset,
                                    int* ExcitationFacesNum,
                                    int excitationCNT,
                                    int excitationPMLCNT,
                                    int8_t* map,
                                    ExcitationProp exci_prop,
                                    int PolyFlag,
                                    fp_t_ts dt_epsilon, fp_t_ts t,
                                    fp_t_ts* nd_coords_tet,
                                    fp_t_ts* nd_coords_face,
                                    fp_t_ts* Z_face,
                                    fp_t_ts* InvMat,
                                    fp_t_ts* Field);
                               
    __global__ void addExcitationH( int* ExcitationFacesCnt, 
                                      int* ExcitationFacesOffset, 
                                      int* ExcitationFacesNum, 
                                      int excitationCNT,
                                      int8_t* map,
                                      ExcitationProp exci_prop,
                                      int PolyFlag,
                                      fp_t_ts dt_muo, fp_t_ts t,
                                      fp_t_ts* nd_coords_tet, 
                                      fp_t_ts* nd_coords_face, 
                                      fp_t_ts* Z_face,
                                      fp_t_ts* InvMat,
                                      fp_t_ts* Field);
  
      __global__ void addExcitationE( int* ExcitationFacesCnt,
                                      int* ExcitationFacesOffset,
                                      int* ExcitationFacesNum,
                                      int excitationCNT,
                                      int8_t* map,
                                      ExcitationProp exci_prop,
                                      int PolyFlag,
                                      fp_t_ts dt_epsilon, fp_t_ts t,
                                      fp_t_ts* nd_coords_tet,
                                      fp_t_ts* nd_coords_face,
                                      fp_t_ts* Z_face,
                                      fp_t_ts* InvMat,
                                      fp_t_ts* Field);
                                    

    __global__ void CalculatePadeFreq(fp_t* a_k, fp_t* b_k, int M_global, int N_2, int* constantM, cuDoubleComplex* Hf);


    __global__ void fillOnes(fp_t_ts* arr, size_t N);

    #endif
#endif //DGTD_KERNELS_CU