/*!\file RiftConstraintsState.cpp
 * \brief: manage penalties for rifts 
 */

#include "./ConstraintsStateLocal.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../include/include.h"
#include "../../io/io.h"
#include "../../shared/shared.h"

#define _ZIGZAGCOUNTER_

/*current module: */
/*RiftIsPresent(Loads* loads,int configuration_type){{{*/
int RiftIsPresent(Loads* loads,int configuration_type){


	int i;
	
	int found=0;
	int mpi_found=0;

	/*go though loads, and figure out if one of the loads is a Riftfront: */
	for (i=0;i<loads->Size();i++){
		Load* load=(Load*)loads->GetObjectByOffset(i);
		if(load->InAnalysis(configuration_type)){
			if(RiftfrontEnum==loads->GetEnum(i)){
				found=1;
				break;
			}
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
	found=mpi_found;
	#endif

	return found;
}
/*}}}*/
/*RiftConstraintsState(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int configuration_type){{{*/
void RiftConstraintsState(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int configuration_type){

	int num_unstable_constraints=0;
	int converged=0;
	int potential;

	RiftConstrain(&num_unstable_constraints,loads,configuration_type);
	if(num_unstable_constraints==0)converged=1;
	
	if(RiftIsFrozen(loads,configuration_type)){
		converged=1;
		num_unstable_constraints=0;
	}
	else if(num_unstable_constraints<=min_mechanical_constraints){
		if(VerboseModule()) _pprintLine_("   freezing constraints");
		RiftFreezeConstraints(loads,configuration_type);
	}

	/*Assign output pointers: */
	*pconverged=converged;
	*pnum_unstable_constraints=num_unstable_constraints;
}
/*}}}*/
/*RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int configuration_type){{{*/
void RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int configuration_type){

	int			i;
	
	/* generic object pointer: */
	Riftfront* riftfront=NULL;
	Load*      load=NULL;

	int unstable;
	int sum_num_unstable_constraints;
	int num_unstable_constraints=0;	
		
	/*Enforce constraints: */
	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){

			load=(Load*)loads->GetObjectByOffset(i);
			if(load->InAnalysis(configuration_type)){

				riftfront=(Riftfront*)load;

				riftfront->Constrain(&unstable);

				num_unstable_constraints+=unstable;
			}
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
	num_unstable_constraints=sum_num_unstable_constraints;
	#endif
	
	/*Assign output pointers: */
	*pnum_unstable_constraints=num_unstable_constraints;

}
/*}}}*/
/*RiftIsFrozen(Loads* loads,int configuration_type){{{*/
int RiftIsFrozen(Loads* loads,int configuration_type){

	int			i;
	
	/* generic object pointer: */
	Load*      load=NULL;
	Riftfront* riftfront=NULL;
	int found=0;
	int mpi_found=0;

	/*Enforce constraints: */
	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){
			
			load=(Load*)loads->GetObjectByOffset(i);
			if(load->InAnalysis(configuration_type)){

				riftfront=(Riftfront*)load;
				if (riftfront->IsFrozen()){
					found=1;
					break;
				}
			}
		}
	}
	
	/*Is there just one found? that would mean we have frozen! : */
	#ifdef _HAVE_MPI_
	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
	found=mpi_found;
	#endif

	return found;
}
/*}}}*/
/*RiftFreezeConstraints(Loads* loads,int configuration_type){{{*/
void RiftFreezeConstraints(Loads* loads,int configuration_type){

	int			i;
	
	/* generic object pointer: */
	Load*      load=NULL;
	Riftfront* riftfront=NULL;

	/*Enforce constraints: */
	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){

			load=(Load*)loads->GetObjectByOffset(i);
			if(load->InAnalysis(configuration_type)){
				
				riftfront=(Riftfront*)load;
				riftfront->FreezeConstraints();
			}

		}
	}

}
/*}}}*/

/*diverse trials and errors: */
/*RiftIsMaterialStable(Loads* loads){{{*/
int RiftIsMaterialStable(Loads* loads){

	int i;
	
	Riftfront* riftfront=NULL;
	int found=0;
	int mpi_found=0;

	/*go though loads, and if non-linearity of the material has converged, let all penalties know: */
	for (i=0;i<loads->Size();i++){

		if(RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);

			if (riftfront->IsMaterialStable()){
				found=1;
				/*do not break! all penalties should get informed the non-linearity converged!*/
			}
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
	found=mpi_found;
	#endif

	return found;
}
/*}}}*/
/*RiftIsPreStable(Loads* loads){{{*/
int RiftIsPreStable(Loads* loads){


	int i;
	
	Riftfront* riftfront=NULL;
	int found=0;
	int mpi_found=0;

	/*go though loads, and figure out if one of the penpair loads is still not stable: */
	for (i=0;i<loads->Size();i++){

		if(RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);

			if (riftfront->PreStable()==0){
				found=1;
				break;
			}
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
	found=mpi_found;
	#endif

	if (found){
		/*We found an unstable constraint. : */
		return 0;
	}
	else{
		return 1;
	}
}
/*}}}*/
/*RiftSetPreStable(Loads* loads){{{*/
void RiftSetPreStable(Loads* loads){


	int i;
	
	Riftfront* riftfront=NULL;
	int found=0;
	int mpi_found=0;

	/*go though loads, and set loads to pre stable.:*/
	for (i=0;i<loads->Size();i++){

		if(RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
			riftfront->SetPreStable();
		}
	}
}
/*}}}*/
/*RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){{{*/
void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){

	int			i;
	
	/* generic object pointer: */
	Riftfront* riftfront=NULL;

	int unstable;
	int sum_num_unstable_constraints;
	int num_unstable_constraints=0;	
		
	/*Enforce constraints: */
	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);

			riftfront->PreConstrain(&unstable);

			num_unstable_constraints+=unstable;
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
	num_unstable_constraints=sum_num_unstable_constraints;
	#endif
	
	/*Assign output pointers: */
	*pnum_unstable_constraints=num_unstable_constraints;

}
/*}}}*/
/*RiftMaxPenetrationInInputs(Loads* loads){{{*/
void RiftMaxPenetrationInInputs(Loads* loads){

	int			i;
	
	/* generic object pointer: */
	Riftfront* riftfront=NULL;

	/*rift penetration: */
	IssmDouble max_penetration=0;
	IssmDouble mpi_max_penetration;
	IssmDouble penetration;

	/*Ok, we are going to find the node pairs which are not penetrating, even though they 
	 * are penalised. We will release only the one with has least <0 penetration. : */

	max_penetration=0;
	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);

			riftfront->MaxPenetration(&penetration);

			if (penetration>max_penetration)max_penetration=penetration;
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
	MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD);                
	max_penetration=mpi_max_penetration;
	#endif

	/*feed max_penetration to inputs: */
	for(i=0;i<loads->Size();i++){
		Load* load=(Load*)loads->GetObjectByOffset(i);
		load->InputUpdateFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum);
	}
}
/*}}}*/
/*RiftPotentialUnstableConstraints(Loads* loads){{{*/
int RiftPotentialUnstableConstraints(Loads* loads){

	int			i;
	
	/* generic object pointer: */
	Riftfront* riftfront=NULL;

	/*Ok, we are going to find the node pairs which are not penetrating, even though they 
	 * are penalised. We will release only the one with has least <0 penetration. : */
	int unstable=0;
	int sum_num_unstable_constraints=0;
	int num_unstable_constraints=0;

	for (i=0;i<loads->Size();i++){

		if (RiftfrontEnum==loads->GetEnum(i)){

			riftfront=(Riftfront*)loads->GetObjectByOffset(i);

			riftfront->PotentialUnstableConstraint(&unstable);

			num_unstable_constraints+=unstable;
		}
	}

	#ifdef _HAVE_MPI_
	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
	num_unstable_constraints=sum_num_unstable_constraints;
	#endif

	return num_unstable_constraints;
}
/*}}}*/
