/*!\file MassFluxx
 * \brief: compute mass flux along a profile.
 */

#include "./MassFluxx.h"

#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool process_units){

	int i,j;
	extern int num_procs;
	extern int my_rank;
	
	Element* element=NULL;
	int element_id;
	
	/*output: */
	double mass_flux=0;
	double all_mass_flux=0;

	int  counter;

	/*all segments: */
	double** array=NULL;
	int      M;
	int*     mdims_array=NULL;
	int*     ndims_array=NULL;

	/*our segments of interest: */
	double*  segments=NULL;
	int      num_segments;

	/*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */
	parameters->FindParam(&array,&M,&mdims_array,&ndims_array,QmuMassFluxSegmentsEnum);

	/*Retrieve index of segments being used for MassFlux computation: */
	parameters->FindParam(&counter,IndexEnum);

	/*retrieve segments from array: */
	segments=array[counter-1]; //matlab to "C" indexing
	num_segments=mdims_array[counter-1];

	/*Go through segments, and then elements, and figure out which elements belong to a segment. 
	 * When we find one, use the element to compute the mass flux on the segment: */
	for(i=0;i<num_segments;i++){
		element_id=(int)*(segments+5*i+4);
		for(j=0;j<elements->Size();j++){
			element=(Element*)elements->GetObjectByOffset(j);
			if (element->Id()==element_id){
				/*We found the element which owns this segment, use it to compute the mass flux: */
				mass_flux+=element->MassFlux(segments+5*i+0,process_units);
				break;
			}
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
	mass_flux=all_mass_flux;
	#endif

	/*Free ressources:*/
	for(j=0;j<M;j++){
		double* matrix=array[j];
		xDelete<double>(matrix);
	}
	xDelete<int>(mdims_array);
	xDelete<int>(ndims_array);
	xDelete<double*>(array);
	
	/*Assign output pointers: */
	*pmass_flux=mass_flux;
}
