/*\file InternalFront.c
 *\brief: build pressureload
 */

#include "./InternalFront.h"

void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){

	bool*   elementonwater=NULL;
	int*    elements=NULL;
	int*    connectivity=NULL;
	int*    elementconnectivity=NULL;
	int*    front=NULL;
	double* front2=NULL;
	bool    found;
	int     numberofelements,numberofsegments;
	int     N,M;
	int     i,j,ii,jj,id;

	/*Boot module: */
	MODULEBOOT();

	/*checks on arguments on the matlab side: */
	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InternalFrontUsage);

	/*Fetch required fields*/
	FetchData(&numberofelements,mxGetAssignedField(MODEL,0,"numberofelements"));
	if(numberofelements<=0) _error_("No elements found in the model");
	FetchData(&elements,&M,&N,mxGetAssignedField(MODEL,0,"elements"));
	if(M!=numberofelements || N!=3) _error_("Field 'elements' should be of size [md.numberofelements 3]");
	FetchData(&elementonwater,&M,&N,mxGetAssignedField(MODEL,0,"elementonwater"));
	if(M!=numberofelements || N!=1) _error_("Field 'elementonwater' should be of size [md.numberofelements 1]");
	FetchData(&elementconnectivity,&M,&N,mxGetAssignedField(MODEL,0,"elementconnectivity"));
	if(M!=numberofelements || N!=3) _error_("Field 'elementconnectivity' should be of size [md.numberofelements 3]");

	/*Allocate and initialize all variables*/
	numberofsegments=0;
	front=xNew<int>(3*numberofelements*4);

	/*Loop over all elements on water*/
	for(i=0;i<numberofelements;i++){

		/*Skip if on water*/
		if(!elementonwater[i]) continue;

		/*Loop over all three adjacent elements*/
		for(j=0;j<3;j++){

			/*Skip if adjacent element does not exist or is on water*/
			id=elementconnectivity[i*3+j];
			if(id==0) continue;
			if(elementonwater[id-1]) continue;

			/*We have an ice front to add here, let's go!*/
			for(ii=0;ii<3;ii++){

				found=false;
				for(jj=0;jj<3;jj++){
					if(elements[(id-1)*3+ii]==elements[i*3+jj]){
						found=true;
						break;
					}
				}

				/*OK, we just found the node of id, which is not in i. We have the segment*/
				if(!found){
					front[numberofsegments*4+0]=elements[(id-1)*3+(ii+1)%3];
					front[numberofsegments*4+1]=elements[(id-1)*3+(ii+2)%3];
					front[numberofsegments*4+2]=id;
					front[numberofsegments*4+3]=IceEnum;
					numberofsegments+=1;
					break;
				}
			}

			/*In debugging mode, check that the segment has been found*/
			_assert_(!found);
		}
	}

	/*Now that we know how many segments there is we can allocate the final matrix*/
	if(numberofsegments){
		front2=xNew<double>(4*numberofsegments);
		for(i=0;i<4*numberofsegments;i++) front2[i]=(double)front[i];
	}
	xfree((void**)&front);

	/*write output datasets: */
	WriteData(FRONT,front2,numberofsegments,4);

	/*end module: */
	MODULEEND();
}

void InternalFrontUsage(void) {
	_pprintLine_("");
	_pprintLine_("   usage: icefront = " << __FUNCT__ << "(md);");
	_pprintLine_("");
}