Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16789)
@@ -11,6 +11,6 @@
 /*{{{*/
 #include "../../datastructures/datastructures.h"
+#include "../../toolkits/toolkits.h"
 #include "../Update.h"
-
 class DataSet;
 class Parameters;
@@ -26,6 +26,6 @@
 template <class doublematrix> class Matrix;
 template <class doubletype> class Vector;
-
-#include "../../toolkits/toolkits.h"
+class ElementMatrix;
+class ElementVector;
 /*}}}*/
 
@@ -43,6 +43,8 @@
 		virtual void   AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
+		virtual ElementMatrix* CreateKMatrix(void)=0;
 		virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
 		virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
+		virtual ElementVector* CreatePVector(void)=0;
 		virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
 		virtual void   DeleteMaterials(void)=0;
@@ -110,4 +112,5 @@
 	   virtual Element*   SpawnBasalElement(void)=0;
 		virtual IssmDouble TMeltingPoint(IssmDouble pressure)=0;
+		virtual void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
 		virtual void   ResetCoordinateSystem()=0;
 		virtual int    VelocityInterpolation()=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16789)
@@ -2652,4 +2652,63 @@
 	/*clean-up*/
 	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Penta::ReduceMatrices) {{{*/
+void Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
+
+	int analysis_type;
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	if(pe){
+		if(analysis_type==StressbalanceAnalysisEnum){
+			if(this->element_type==MINIcondensedEnum){
+				int approximation;
+				inputs->GetInputValue(&approximation,ApproximationEnum);
+				if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+					//Do nothing, condensation already done in PVectorCoupling
+				}
+				else{
+					int indices[3]={18,19,20};
+					pe->StaticCondensation(Ke,3,&indices[0]);
+				}
+			}
+			else if(this->element_type==P1bubblecondensedEnum){
+				int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+				int offset = 0;
+				for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+				int* indices=xNew<int>(size);
+				for(int i=0;i<size;i++) indices[i] = offset+i;
+				pe->StaticCondensation(Ke,size,indices);
+				xDelete<int>(indices);
+				delete Ke;
+			}
+		}
+	}
+
+	if(Ke){
+		if(analysis_type==StressbalanceAnalysisEnum){
+			int approximation;
+			inputs->GetInputValue(&approximation,ApproximationEnum);
+			if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+				//Do nothing condensatino already done for Stokes part
+			}
+			else{
+				if(this->element_type==MINIcondensedEnum){
+					int indices[3]={18,19,20};
+					Ke->StaticCondensation(3,&indices[0]);
+				}
+				else if(this->element_type==P1bubblecondensedEnum){
+					int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+					int offset = 0;
+					for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+					int* indices=xNew<int>(size);
+					for(int i=0;i<size;i++) indices[i] = offset+i;
+					Ke->StaticCondensation(size,indices);
+					xDelete<int>(indices);
+				}
+			}
+		}
+	}
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16789)
@@ -115,4 +115,5 @@
 		void   ResultToVector(Vector<IssmDouble>* vector,int output_enum);
 		void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
+		void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
 		void   ResetCoordinateSystem(void);
 		void   SmbGradients();
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16789)
@@ -78,6 +78,8 @@
 		void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
 		void        CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){_error_("not implemented yet");};
+		ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");};
 		void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
 		void        CreatePVector(Vector<IssmDouble>* pf){_error_("not implemented yet");};
+		ElementVector* CreatePVector(void){_error_("not implemented yet");};
 		ElementVector* CreatePVectorL2Projection(void);
 		void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
@@ -147,4 +149,5 @@
 		void        ResultToVector(Vector<IssmDouble>* vector,int output_enum){_error_("not implemented");};
 		void        ResetCoordinateSystem(void){_error_("not implemented yet");};
+		void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
 		void	      SmbGradients(){_error_("not implemented yet");};
 		IssmDouble  TMeltingPoint(IssmDouble pressure){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16789)
@@ -2349,4 +2349,43 @@
 	/*clean-up*/
 	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::ReduceMatrices {{{*/
+void Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
+
+	/*Static condensation if requested*/
+	if(pe){
+		if(this->element_type==MINIcondensedEnum){
+			int indices[2]={6,7};
+			pe->StaticCondensation(Ke,2,&indices[0]);
+		}
+		else if(this->element_type==P1bubblecondensedEnum){
+			int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+			int offset = 0;
+			for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+			int* indices=xNew<int>(size);
+			for(int i=0;i<size;i++) indices[i] = offset+i;
+			pe->StaticCondensation(Ke,size,indices);
+			xDelete<int>(indices);
+		}
+	}
+
+	if(Ke){
+		if(this->element_type==MINIcondensedEnum){
+			int indices[2]={6,7};
+			Ke->StaticCondensation(2,&indices[0]);
+		}
+		else if(this->element_type==P1bubblecondensedEnum){
+			int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+			int offset = 0;
+			for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+			int* indices=xNew<int>(size);
+			for(int i=0;i<size;i++) indices[i] = offset+i;
+			Ke->StaticCondensation(size,indices);
+			xDelete<int>(indices);
+		}
+	}
+
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16788)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16789)
@@ -121,4 +121,5 @@
 		void        ResultInterpolation(int* pinterpolation,int output_enum);
 		void        ResultToVector(Vector<IssmDouble>* vector,int output_enum);
+		void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
 		void        ResetCoordinateSystem(void);
 		void	      SmbGradients();
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 16788)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 16789)
@@ -12,5 +12,5 @@
 	/*intermediary: */
 	int      i,M,N;
-	int      configuration_type;
+	int      configuration_type,analysisenum;
 	Element *element = NULL;
 	Load    *load    = NULL;
@@ -26,6 +26,8 @@
 	if(VerboseModule()) _printf0_("   Generating matrices");
 
-	/*retrive parameters: */
+	/*retrieve parameters: */
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&analysisenum,AnalysisTypeEnum);
+	Analysis* analysis = EnumToAnalysis(analysisenum);
 
 	/*First, we might need to do a dry run to get kmax if penalties are employed*/
@@ -62,23 +64,26 @@
 	}
 
-	/*Fill stiffness matrix from elements and loads */
+	/*Fill stiffness matrix and load vector from elements*/
 	for (i=0;i<femmodel->elements->Size();i++){
 		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->CreateKMatrix(Kff,Kfs);
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		//ElementVector* pe = analysis->CreatePVector(element);
+		//ElementMatrix* Ke = analysis->CreateKMatrix(element);
+		ElementVector* pe = element->CreatePVector();
+		ElementMatrix* Ke = element->CreateKMatrix();
+		element->ReduceMatrices(Ke,pe);
+		if(Ke) Ke->AddToGlobal(Kff,Kfs);
+		if(pe) pe->AddToGlobal(pf);
+		delete Ke;
+		delete pe;
 	}
 
-	for (i=0;i<femmodel->loads->Size();i++){
+	/*Fill stiffness matrix and load vector from loads*/
+	for(i=0;i<femmodel->loads->Size();i++){
 		load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
-		if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
-	}
-
-	/*Fill right hand side vector, from elements and loads */
-	for (i=0;i<femmodel->elements->Size();i++){
-		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->CreatePVector(pf);
-	}
-	for (i=0;i<femmodel->loads->Size();i++){
-		load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
-		if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
+		if(load->InAnalysis(configuration_type)){
+			load->CreateKMatrix(Kff,Kfs);
+			load->CreatePVector(pf);
+		}
 	}
 
@@ -87,9 +92,8 @@
 		for (i=0;i<femmodel->loads->Size();i++){
 			load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
-			if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
-		}
-		for (i=0;i<femmodel->loads->Size();i++){
-			load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
-			if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
+			if(load->InAnalysis(configuration_type)){
+				load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
+				load->PenaltyCreatePVector(pf,kmax);
+			}
 		}
 	}
@@ -109,5 +113,6 @@
 	//Kfs->AllocationInfo();
 
-	/*Assign output pointers: */
+	/*cleanu up and assign output pointers: */
+	delete analysis;
 	if(pKff) *pKff=Kff;
 	else      delete Kff;
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 16788)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 16789)
@@ -6,4 +6,5 @@
 
 #include "../../classes/classes.h"
+#include "../../analyses/analyses.h"
 
 /* local prototypes: */