#ifndef GALERKIN_INTEGRALS_HPP
#define GALERKIN_INTEGRALS_HPP

#include "green.hpp"
#include "element.hpp"
#include "quadrature.hpp"

#include "vector.h" /* dot, cross */

#include <cmath>   /* sqrt */
#include <cstdlib> /* malloc */

/* ------------------------------------------------------------------------ */
/* Compute regular double L and integral over two linear TRIA elements */
template <typename kType>
void dint_tritri_lin(const gauss_t *g3,
                     const double *nodesx,
                     const accelerator_t *acceleratorx,
                     const double *nodesy,
                     const accelerator_t *acceleratory,
                     kType const& k,
                     complex_scalar *l,
                     complex_scalar *m
                    )
{
    int ix, iy, j, s, qx, qy;
    double rx[3], jacx, jacy;
    double r[3], nx[3], ny[3];
	complex_scalar g, dg, ltemp, mtemp;

    for (s = 0; s < 3*3; s++)
        l[s] = m[s] = complex_scalar(0.0, 0.0);

    jacx = sqrt(dot(acceleratorx->n0, acceleratorx->n0));
    jacy = sqrt(dot(acceleratory->n0, acceleratory->n0));
	
	for (j = 0; j < 3; ++j)
	{
		nx[j] = acceleratorx->n0[j] / jacx;
		ny[j] = acceleratory->n0[j] / jacy;
	}

    /* for each x gaussian integration point */
    for (ix = 0; ix < g3->num; ix++)
    {
        /* computing integration x location */
        for (j = 0; j < 3; j++)
        {
            rx[j] = 0.0;
            for (s = 0; s < 3; s++)
                rx[j] += g3->N[ix+s*g3->num]*nodesx[s+3*j];
        }
        /* for each y gaussian integration point */
        for (iy = 0; iy < g3->num; iy++)
        {
            /* computing integration y location */
            for (j = 0; j < 3; j++)
            {
                r[j] = -rx[j];
                for (s = 0; s < 3; s++)
                    r[j] += g3->N[iy+s*g3->num]*nodesy[s+3*j];
            }
            
			green(r, k, g, ny, dg);

            ltemp = g3->w[ix] * g3->w[iy] * dg;
            mtemp = g3->w[ix] * g3->w[iy] * g;

            for (qx = 0; qx < 3; qx++)
            {
                for (qy = 0; qy < 3; qy++)
                {
                    l[qx+3*qy] += g3->N[ix+qx*g3->num] * g3->N[iy+qy*g3->num] * ltemp;
                    m[qx+3*qy] += g3->N[ix+qx*g3->num] * g3->N[iy+qy*g3->num] * mtemp;
                }
            }
        }
    }

    /* Perform multiplication */
    for (s = 0; s < 3*3; s++)
    {
        l[s] *= jacx*jacy;
        m[s] *= jacx*jacy;
    }
}

/* ------------------------------------------------------------------------ */
/* Compute singular double B integral over a linear TRIA element            */
template <typename kType>
void dint_tri_lin_sing(const gauss_t *g3,
                       const double *nodes,
                       const accelerator_t *accelerator,
                       const gauss_t *g4,
                       kType const& k,
                       complex_scalar *l,
                       complex_scalar *m
                      )
{
    int ix, iy, j, s, qx, qy;
    double rx[3], ry[3], jac;
    double N[3], nx[3];
    complex_scalar g, dg;
    double *wprime, *xiprime, *etaprime;
    int gnum;

	wprime = new double [3*g4->num];
    xiprime = new double [3*g4->num];
    etaprime = new double [3*g4->num];
	
    for (s = 0; s < 3*3; s++)
        l[s] = m[s] = complex_scalar(0.0, 0.0);

    /* Jacobian over the triangle */
    jac = sqrt(dot(accelerator->n0, accelerator->n0));
	
	/* Compute the normal vector */
	for (j = 0; j < 3; ++j)
		nx[j] = accelerator->n0[j] / jac;

    /* for each x gaussian integration point */
    for (ix = 0; ix < g3->num; ix++)
    {
        /* computing integration x location */
        for (j = 0; j < 3; j++)
        {
            rx[j] = 0.0;
            for (s = 0; s < 3; s++)
                rx[j] += g3->N[ix+s*g3->num]*nodes[s+3*j];
        }
        
		/* get singular quadrature points */
        sing_quadr_face<TriaElem>(*g4, g3->xi[ix], g3->xi[ix+g3->num], &gnum, xiprime, etaprime, wprime);
		
        /* for each gaussian integration point */
        for (iy = 0; iy < gnum; iy++)
        {
			/* Evaluate shape function values */
            shapefun<TriaElem>(xiprime[iy], etaprime[iy], N);

            /* for each coordinate direction */
            for (j = 0; j < 3; j++)
            {
                ry[j] = -rx[j];
                /* computing integration location and its derivatives */
                for (s = 0; s < 3; s++)
                    ry[j] += N[s]*nodes[s+3*j];
            }
            
            green(ry, k, g, nx, dg);

            for (qx = 0; qx < 3; qx++)
            {
                for (qy = 0; qy < 3; qy++)
                {
                    l[qx+3*qy] += g3->N[ix+qx*g3->num] * g3->w[ix] * N[qy] * wprime[iy] * dg;
                    m[qx+3*qy] += g3->N[ix+qx*g3->num] * g3->w[ix] * N[qy] * wprime[iy] * g;
                }
            }
        }
    }

    for (s = 0; s < 3*3; s++)
    {
        l[s] *= jac*jac;
        m[s] *= jac*jac;
    }

    delete [] wprime;
    delete [] xiprime;
    delete [] etaprime;
}

#endif // GALERKIN_INTEGRALS_HPP