Index: /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp
===================================================================
--- /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp	(revision 4698)
+++ /issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp	(revision 4699)
@@ -1,1 +1,350 @@
-s:
+/*!\file RiftConstraints.cpp
+ * \brief: manage penalties for rifts 
+ */
+
+#include "./PenaltyConstraintsLocal.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../include/include.h"
+#include "../../shared/shared.h"
+
+#define _ZIGZAGCOUNTER_
+
+int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int analysis_type){
+
+	int num_unstable_constraints=0;
+	int converged=0;
+	int potential;
+	ISSMERROR(" check analysis_type across all routine!");
+
+	Constrain(&num_unstable_constraints,loads);
+	if(num_unstable_constraints==0)converged=1;
+	
+	
+	if(IsFrozen(loads)){
+		converged=1;
+		num_unstable_constraints=0;
+	}
+	else if(num_unstable_constraints<=min_mechanical_constraints){
+		_printf_("   freezing constraints\n");
+		FreezeConstraints(loads);
+	}
+
+	/*Assign output pointers: */
+	*pconverged=converged;
+	*pnum_unstable_constraints=num_unstable_constraints;
+}
+
+int IsMaterialStable(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 _PARALLEL_
+	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;
+}
+
+int RiftIsPresent(Loads* loads,int analysis_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(analysis_type)){
+			if(RiftfrontEnum==loads->GetEnum(i)){
+				found=1;
+				break;
+			}
+		}
+	}
+
+	#ifdef _PARALLEL_
+	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;
+}
+
+int IsPreStable(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 _PARALLEL_
+	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;
+	}
+}
+
+int SetPreStable(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();
+		}
+	}
+}
+
+int PreConstrain(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 _PARALLEL_
+	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;
+
+}
+
+int Constrain(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->Constrain(&unstable);
+
+			num_unstable_constraints+=unstable;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	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;
+
+}
+
+void FreezeConstraints(Loads* loads){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->FreezeConstraints();
+
+		}
+	}
+
+}
+
+int IsFrozen(Loads* loads){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+	int found=0;
+	int mpi_found=0;
+
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+			if (riftfront->IsFrozen()){
+				found=1;
+				break;
+			}
+		}
+	}
+	
+	/*Is there just one found? that would mean we have frozen! : */
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	return found;
+}
+
+int MaxPenetrationInInputs(Loads* loads){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*rift penetration: */
+	double max_penetration=0;
+	double mpi_max_penetration;
+	double penetration;
+
+	/*Ok, we are going to find the grid 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 _PARALLEL_
+	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);
+	}
+}
+
+int PotentialUnstableConstraints(Loads* loads){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*Ok, we are going to find the grid 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 _PARALLEL_
+	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;
+}