/*
 * \file GrdLoads.cpp
 * \brief: Implementation of GrdLoads class
 */

/*Headers: {{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./GrdLoads.h"
#include "./SealevelGeometry.h"
using namespace std;
/*}}}*/

/*Object constructors and destructor*/
GrdLoads::GrdLoads(int nel,SealevelGeometry* slgeom){ /*{{{*/

	vloads=new Vector<IssmDouble>(nel);
	for (int i=0;i<SLGEOM_NUMLOADS;i++) vsubloads[i]=new Vector<IssmDouble>(slgeom->nbar[i]);

	vsealevelloads=new Vector<IssmDouble>(nel);
	vsealevelloads->Set(0);vsealevelloads->Assemble();

	vsubsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);

}; /*}}}*/
GrdLoads::~GrdLoads(){ /*{{{*/

	delete vloads;
	xDelete<IssmDouble>(loads);
	for(int i=0;i<SLGEOM_NUMLOADS;i++){
		delete vsubloads[i];
		xDelete<IssmDouble>(subloads[i]);
	}
	delete vsealevelloads;
	xDelete<IssmDouble>(sealevelloads);
	delete vsubsealevelloads;
	xDelete<IssmDouble>(subsealevelloads);
}; /*}}}*/

void GrdLoads::BroadcastLoads(void){ /*{{{*/

	/*Initialize barycentre vectors, now that we know their size: */
	vloads->Assemble();
	for (int i=0;i<SLGEOM_NUMLOADS;i++){
		vsubloads[i]->Assemble();
	}

	loads=vloads->ToMPISerial();
	for (int i=0;i<SLGEOM_NUMLOADS;i++){
		subloads[i]=vsubloads[i]->ToMPISerial();
	}

} /*}}}*/
void GrdLoads::AssembleSealevelLoads(void){ /*{{{*/

	vsealevelloads->Assemble();
	vsubsealevelloads->Assemble();

} /*}}}*/
void GrdLoads::BroadcastSealevelLoads(void){ /*{{{*/

	sealevelloads=vsealevelloads->ToMPISerial();
	subsealevelloads=vsubsealevelloads->ToMPISerial();

} /*}}}*/
void GrdLoads::SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom){

	IssmDouble area,re, S;
	int ylmindex, intj;
	IssmDouble deg2coeff_local[5];

	femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);

	for (int c=0;c<5;c++) deg2coeff_local[c]=0;

	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		//first, loads on centroids
		S=0;
		S+=loads[element->Sid()];

		if(sealevelloads) S+=sealevelloads[element->Sid()];
		if(S!=0){
			element->Element::GetInputValue(&area,AreaEnum);

			for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients: 2,0; 2,1cos; 2,1sin; 2,2cos; 2,2sin
				ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
				deg2coeff_local[c] += S*area/re/re*slgeom->Ylm[ylmindex];
			}
		}
		//add loads on subelement barycenters
		for (int i=0;i<SLGEOM_NUMLOADS;i++){
			if (slgeom->issubelement[i][element->lid]){
				intj=slgeom->subelementmapping[i][element->lid];
				S=0;
				S+=subloads[i][intj];
				if(i==SLGEOM_OCEAN && sealevelloads) S+=subsealevelloads[intj];
				if(S!=0){
					area=slgeom->area_subel[i][intj];
					for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients
						ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
						deg2coeff_local[c] += S*area/re/re*slgeom->Ylm_subel[i][ylmindex];
					}
				}
			}
		}
	}

	//sum each degree 2 coefficient across all cpus to get global total
	ISSM_MPI_Reduce (&deg2coeff_local[0],&deg2coeff[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&deg2coeff[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	ISSM_MPI_Reduce (&deg2coeff_local[1],&deg2coeff[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&deg2coeff[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
	ISSM_MPI_Reduce (&deg2coeff_local[2],&deg2coeff[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
	ISSM_MPI_Bcast(&deg2coeff[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

}