#include "mex.h"
#include <math.h>

/* The computational routine */
void assembLoop(double* bDet, double *tri, double *f, double *b, mwSize n)
{
    double S4[3][3] = { {2, 1, 1},
                        {1, 2, 1},
                        {1, 1, 2}};

    mwSize i,j,m;
    double temp, myDet; 
                        
    for (i = 0; i < n; i++) {
      myDet = fabs(bDet[i]);
      for (m = 0; m < 3; m++) {
        temp = 0;
        for (j = 0; j < 3; j++) {
          temp += f[(int)tri[i + n*j]-1] * S4[m][j] / 24;
        }
        b[(int)tri[i + n*m]-1] += myDet * temp;
      }
    }
}

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    double* bDet;
    double* tri;
    double* f;
    size_t columnsF,rowsF,columnsT,rowsT;
    double* b;

    /* Inputs parsen und analysisieren */
    bDet = mxGetPr(prhs[0]);
    tri  =  mxGetPr(prhs[1]);
    f    = mxGetPr(prhs[2]);

    rowsF    = (int) mxGetM(prhs[2]);
    columnsF = (int) mxGetN(prhs[2]);

    rowsT    = (int) mxGetM(prhs[1]);
    columnsT = (int) mxGetN(prhs[1]);

    /* Output Erzeugen */
    plhs[0] = mxCreateDoubleMatrix((mwSize)rowsF, (mwSize)columnsF, mxREAL);
    b       = mxGetPr(plhs[0]);

    /* call the computational routine */
    assembLoop(bDet,tri,f,b,(mwSize)rowsT);
}
