#ifndef BEMHG_GALERKIN_HPP
#define BEMHG_GALERKIN_HPP

#include "elem_traits.hpp"

#include "galerkin_integrals.hpp"
#include "vector.h"
#include "mesh.h"

/* ------------------------------------------------------------------------ */
/* Compute surface system matrices of a bem model */
template <typename kType>
        void matrix_surf_galerkin(
		int nnodes,
        double const *nodes,
        int nelements,
        double const *elements,
        gauss_t const *g3,
        gauss_t const *g4,
        double const *dist,
        kType const &k,
        double *Ar,
        double *Ai,
        double *Br,
        double *Bi)
{
    enum {MAXNVERT = 4, NDIM = 3};
    
	/* Intermediate storage */
	int elemx[MAXNVERT];					/* node indices of x element */
    double nodx[MAXNVERT*NDIM];				/* nodal coordinates of x element */
	complex_scalar a[MAXNVERT*MAXNVERT];
	complex_scalar b[MAXNVERT*MAXNVERT];
	
    /* Compute element centres */
    accelerator_t *accelerators = new accelerator_t[nelements];
    init_accelerators(nnodes, nodes, nelements, elements, accelerators);
    
    /* Clear output matrices */
    for (int n = 0; n < nnodes; ++n)
        for (int j = 0; j < nnodes; ++j)
            Ar[n+nnodes*j] = Ai[n+nnodes*j] = Br[n+nnodes*j] = Bi[n+nnodes*j] = 0.0;
    
	/* Integration for each element (x) as reference */
    for (int ex = 0; ex < nelements; ex++)
    {
        int nvertx = (int)elements[ex]; 		/* number of vertices of x elem */
        
      
        
        /* Collect element vertex nodes and coordinates */
        for (int s = 0; s < nvertx; s++)
        {
            elemx[s] = (int)elements[ex+(s+1)*nelements];
            for (int j = 0; j < 3; j++)
                nodx[s+nvertx*j] = nodes[elemx[s]+j*nnodes];
		}
        
        /* singular integration */
        dint_tri_lin_sing<kType>(&g3[0], nodx, &accelerators[ex], &g4[0], k, a, b);
        
        /* Fill in the element matrix into the system matrix */
        for (int ix = 0; ix < nvertx; ix++)
        {
            for (int iy = 0; iy < nvertx; iy++)
            {
				Ar[elemx[ix]+nnodes*elemx[iy]] += a[ix+nvertx*iy].real();
				Ai[elemx[ix]+nnodes*elemx[iy]] += a[ix+nvertx*iy].imag();
				Br[elemx[ix]+nnodes*elemx[iy]] += b[ix+nvertx*iy].real();
				Bi[elemx[ix]+nnodes*elemx[iy]] += b[ix+nvertx*iy].imag();
            }
        }

        /* Integration for the half of elements (y) */
        for (int ey = ex+1; ey < nelements; ey++)
        {
			int nverty = (int)elements[ey]; /* number of vertices of y elem */
			int elemy[MAXNVERT];					/* node indices of x element */
			double nody[MAXNVERT*NDIM];				/* nodal coordinates of x element */
			/* Collect element vertex nodes and coordinates */
			for (int s = 0; s < nverty; s++)
			{
				elemy[s] = (int)elements[ey+(s+1)*nelements];
				for (int j = 0; j < 3; j++)
					nody[s+nverty*j] = nodes[elemy[s]+j*nnodes];
			}
            int gs = gauss_division(accelerators[ex].center, accelerators[ey].center, dist);
           
			/* regular integration */
			dint_tritri_lin<kType>(&g3[gs], nodx, &accelerators[ex], nody, &accelerators[ey], k, a, b);
			
			/* Fill in the element matrix */
            for (int ix = 0; ix < nvertx; ix++)
            {
                for (int iy = 0; iy < nverty; iy++)
                {
					Ar[elemx[ix]+nnodes*elemy[iy]] += a[ix+nvertx*iy].real();
					Ar[elemy[iy]+nnodes*elemx[ix]] += a[ix+nvertx*iy].real();
					Ai[elemx[ix]+nnodes*elemy[iy]] += a[ix+nvertx*iy].imag();
					Ai[elemy[iy]+nnodes*elemx[ix]] += a[ix+nvertx*iy].imag();
					Br[elemx[ix]+nnodes*elemy[iy]] += b[ix+nvertx*iy].real();
					Br[elemy[iy]+nnodes*elemx[ix]] += b[ix+nvertx*iy].real();
					Bi[elemx[ix]+nnodes*elemy[iy]] += b[ix+nvertx*iy].imag();
					Bi[elemy[iy]+nnodes*elemx[ix]] += b[ix+nvertx*iy].imag();
                }
            }
        }
    }
    delete [] accelerators;
}

#endif // BEMHG_GALERKIN_HPP