Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4011)
+++ /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4012)
@@ -87,13 +87,9 @@
 
 	/*Compute slopes: */
-	slope_core(&surfaceslopex,&surfaceslopey,femmodel,SurfaceAnalysisEnum);
-	slope_core(&bedslopex,&bedslopey,femmodel,BedAnalysisEnum);
-		
-	/*Update: */
-	model->UpdateInputsFromVector(surfaceslopex,SurfaceSlopexEnum,VertexEnum);
-	model->UpdateInputsFromVector(surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
-	model->UpdateInputsFromVector(bedslopex,BedSlopexEnum,VertexEnum);
-	model->UpdateInputsFromVector(bedslopey,BedSlopeyEnum,VertexEnum);
-
+	slope_core(femmodel,SurfaceXAnalysisEnum);
+	slope_core(femmodel,SurfaceYAnalysisEnum);
+	slope_core(femmodel,BedSlopeXAnalysisEnum);
+	slope_core(femmodel,BedSlopeYAnalysisEnum);
+		
 	
 	if(ishutter){
Index: /issm/trunk/src/c/solutions/slope_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/slope_core.cpp	(revision 4011)
+++ /issm/trunk/src/c/solutions/slope_core.cpp	(revision 4012)
@@ -9,5 +9,5 @@
 #include "../modules/modules.h"
 
-void slope_core(Vec* pslopex,Vec* pslopey,FemModel* femmodel, int AnalysisEnum){
+void slope_core(FemModel* femmodel, int sub_analysis_type){
 
 	/*parameters: */
@@ -16,10 +16,7 @@
 	bool isstokes;
 	bool ishutter;
-	int xanalysis;
-	int yanalysis;
 
 	/*output: */
-	Vec slopex=NULL;
-	Vec slopey=NULL; 
+	Vec slope=NULL;
 
 	/*Recover some parameters: */
@@ -29,46 +26,33 @@
 	femmodel->parameters->FindParam(&ishutter,IsHutterEnum);
 
-	if(verbose)_printf_("%s\n","computing slope (x and y derivatives)...");
-
-	/*Specify type of computations: */
-	if(AnalysisEnum==SurfaceAnalysisEnum){
-		xanalysis=SurfaceXAnalysisEnum;
-		yanalysis=SurfaceYAnalysisEnum;
-	}
-	else if(AnalysisEnum==BedAnalysisEnum){
-		xanalysis=BedXAnalysisEnum;
-		yanalysis=BedYAnalysisEnum;
-	}
-	else ISSMERROR("%s%s%s"," analysis ",EnumAsString(AnalysisEnum)," not supported yet!");
-
+	if(verbose)_printf_("%s\n","computing slope...");
 
 	/*Early return possible? */
-	if(!ishutter && AnalysisEnum==SurfaceAnalysisEnum){
-		/*no need to compute Surface Slope except for Hutter: */
-		*pslopex=NULL;
-		*pslopey=NULL;
-		return;
-	}
-	if(!isstokes && AnalysisEnum==BedAnalysisEnum){
-		/*no need to compute Bed Slope except for full Stokes: */
-		*pslopex=NULL;
-		*pslopey=NULL;
-		return;
-	}
+	if(!ishutter && AnalysisEnum==SurfaceAnalysisEnum) return;
+	if(!isstokes && AnalysisEnum==BedAnalysisEnum) return;
+
+	/*Call on core computations: */
+	solver_linear(&slope,femmodel,SlopecomputeAnalysisEnum,sub_analysis_type);
 	
-	
-	/*Call on core computations: */
-	diagnostic_core_linear(&slopex,femmodel,SlopecomputeAnalysisEnum,xanalysis);
-	diagnostic_core_linear(&slopey,femmodel,SlopecomputeAnalysisEnum,yanalysis);
+	/*Plug solution into inputs: */
+	UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,slope,analysis_type, sub_analysis_type);
 
 	/*extrude if we are in 3D: */
 	if (dim==3){
-		if(verbose)_printf_("%s\n","extruding slopes in 3d...");
-		FieldExtrudex( slopex, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"slopex",0);
-		FieldExtrudex( slopey, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"slopey",0);
+		if(verbose)_printf_("%s\n","extruding slope in 3d...");
+		switch(sub_analysis_type){
+			case SurfaceXAnalysisEnum:
+				InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum,0);
+			case SurfaceYAnalysisEnum:
+				InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum,0);
+			case BedXAnalysisEnum:
+				InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum,0);
+			case BedYAnalysisEnum:
+				InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum,0);
+			default: ISSMERROR("%s%s%s"," sub_analysis_type ",EnumAsString(sub_analysis_type)," not supported yet!");
+		}
 	}
 
-	/*Assign output pointers:*/
-	*pslopex=slopex;
-	*pslopey=slopey;
+	/*Free ressources:*/
+	VecFree(&slope);
 }
Index: /issm/trunk/src/c/solutions/solutions.h
===================================================================
--- /issm/trunk/src/c/solutions/solutions.h	(revision 4011)
+++ /issm/trunk/src/c/solutions/solutions.h	(revision 4012)
@@ -32,8 +32,6 @@
 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,int analysis_type,int sub_analysis_type);
 void diagnostic_core_linear(Vec* ppg,FemModel* fem,int  analysis_type,int sub_analysis_type);
-void slope_core(Vec* pslopex,Vec* pslopey,FemModel* fem,int AnalysisEnum);
-
-
-
+void solver_linear(Vec* ppg,FemModel* fem,int  analysis_type,int sub_analysis_type);
+void slope_core(FemModel* fem,int sub_analysis_type);
 
 
Index: /issm/trunk/src/c/solutions/solver_linear.cpp
===================================================================
--- /issm/trunk/src/c/solutions/solver_linear.cpp	(revision 4012)
+++ /issm/trunk/src/c/solutions/solver_linear.cpp	(revision 4012)
@@ -0,0 +1,65 @@
+/*!\file: solver_linear.cpp
+ * \brief: numerical core of linear solutions
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../modules/modules.h"
+
+void solver_linear(Vec* pug,FemModel* fem,int analysis_type,int sub_analysis_type){
+
+	/*parameters:*/
+	int kflag,pflag;
+	int verbose=0;
+	char* solver_string=NULL;
+
+	/*output: */
+	Vec ug=NULL;
+	Vec uf=NULL; 
+	
+	/*intermediary: */
+	Mat Kgg=NULL;
+	Mat Kff=NULL;
+	Mat Kfs=NULL;
+	Vec pg=NULL;
+	Vec pf=NULL;
+
+	/*Recover parameters: */
+	kflag=1; pflag=1;
+	fem->parameters->FindParam(&verbose,VerboseEnum);
+	fem->parameters->FindParam(&solver_string,SolverStringEnum);
+
+	//*Generate system matrices
+	if (verbose) _printf_("   Generating matrices\n");
+	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+
+	if (verbose) _printf_("   Generating penalty matrices\n");
+	//*Generate penalty system matrices
+	PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+
+	/*!Reduce matrix from g to f size:*/
+	if (verbose) _printf_("   reducing matrix from g to f set\n");
+	Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg);
+
+	/*!Reduce load from g to f size: */
+	if (verbose) _printf_("   reducing load from g to f set\n");
+	Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);VecFree(&pg); MatFree(&Kfs);
+
+	/*Solve: */
+	if (verbose) _printf_("   solving\n");
+	Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf);
+
+	//Merge back to g set
+	if (verbose) _printf_("   merging solution from f to g set\n");
+	Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);VecFree(&uf);
+
+	//Update inputs using new solution:
+	UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type);
+
+	/*free ressources: */
+	xfree((void**)&solver_string);
+
+	/*Assign output pointers:*/
+	*pug=ug;
+}