/*! \file  ContourToMeshx.c
 */

#include "./ContourToMeshx.h"

#undef __FUNCT__ 
#define __FUNCT__  "ContourToMeshx"

int ContourToMeshx( Vec* pin_nod,Vec* pin_elem, double* index, double* x, double* y,Contour** contours,int numcontours,char* interptype,int nel,int nods, int edgevalue) {

	int noerr=1;
	int i;
	int m,n;

	/*Contour:*/
	Contour* contouri=NULL;
	int      numgrids;
	double*  xc=NULL;
	double*  yc=NULL;
	double*  in_nod_serial;
	double   value;

	/*output: */
	Vec in_nod=NULL;
	Vec in_elem=NULL;

	in_nod=NewVec(nods);
	in_elem=NewVec(nel);

	/*Loop through all contours: */
	for (i=0;i<numcontours;i++){
		#ifdef _DEBUG_
			printf("Handling contour %i/%i\n",i,numcontours);
		#endif
		contouri=*(contours+i);
		numgrids=contouri->nods;
		xc=contouri->x;
		yc=contouri->y;
		IsInPoly(in_nod,xc,yc,numgrids,x,y,nods,edgevalue);
	}

	/*Get in_nod serialised for next operation: */
	VecToMPISerial(&in_nod_serial,in_nod);

	/*Take care of the case where an element interpolation has been requested: */
	if ((strcmp(interptype,"element")==0) || (strcmp(interptype,"element and node")==0)){
		for (n=0;n<nel;n++){
			if ( (in_nod_serial[ (int)*(index+3*n+0) -1] == 1) && (in_nod_serial[ (int)*(index+3*n+1) -1] == 1) && (in_nod_serial[ (int)*(index+3*n+2) -1] == 1) ){
				value=1; VecSetValues(in_elem,1,&n,&value,INSERT_VALUES);
			}
		}
	}

	/*Assemble vectors: */
	VecAssemblyBegin(in_nod);
	VecAssemblyEnd(in_nod);
	VecAssemblyBegin(in_elem);
	VecAssemblyEnd(in_elem);

	/*Assign output pointers: */
	*pin_nod=in_nod;
	*pin_elem=in_elem;

	/*Free ressources:*/
	xfree((void**)&in_nod_serial);

	return noerr;
}
