Index: /issm/trunk-jpl/src/c/analyses/Analysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/Analysis.h	(revision 17008)
+++ /issm/trunk-jpl/src/c/analyses/Analysis.h	(revision 17009)
@@ -33,4 +33,5 @@
 
 		/*Finite element Analysis*/
+		virtual void           Core(FemModel* femmodel)=0;
 		virtual ElementVector* CreateDVector(Element* element)=0;
 		virtual ElementMatrix* CreateJacobianMatrix(Element* element)=0;
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17008)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17009)
@@ -4,4 +4,5 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
 
 /*Model processing*/
@@ -813,5 +814,40 @@
 /*Finite Element Analysis*/
 void           StressbalanceAnalysis::Core(FemModel* femmodel){/*{{{*/
-	_error_("not implemented");
+
+	/*Intermediaries*/
+	bool isSIA,isSSA,isL1L2,isHO,isFS;
+	bool conserve_loads = true;
+	int  meshtype,newton;
+
+	/* recover parameters:*/
+	femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum);
+	femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
+	femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum);
+	femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+
+	if((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
+		if(VerboseSolution()) _printf0_("   computing velocities\n");
+
+		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+		if(newton>0)
+		 solutionsequence_newton(femmodel);
+		else
+		 solutionsequence_nonlinear(femmodel,conserve_loads); 
+	}
+
+	if ((isSSA || isL1L2 || isHO) && isFS){
+		if(VerboseSolution()) _printf0_("   computing coupling between lower order models and FS\n");
+		solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
+	}
+
+	if (meshtype==Mesh3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
+		if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
+		femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
+		solutionsequence_linear(femmodel);
+	}
+
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreateDVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17008)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17009)
@@ -4,4 +4,5 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
 
 /*Model processing*/
@@ -136,5 +137,8 @@
 /*Finite Element Analysis*/
 void           StressbalanceSIAAnalysis::Core(FemModel* femmodel){/*{{{*/
-	_error_("not implemented");
+
+		if(VerboseSolution()) _printf0_("   computing SIA velocities\n");
+		femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
+		solutionsequence_linear(femmodel);
 }/*}}}*/
 ElementVector* StressbalanceSIAAnalysis::CreateDVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 17008)
+++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 17009)
@@ -6,4 +6,5 @@
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
+#include "../analyses/analyses.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
@@ -13,15 +14,12 @@
 
 	/*parameters: */
-	bool  dakota_analysis;
-	int   meshtype;
-	bool  isSIA,isSSA,isL1L2,isHO,isFS;
-	bool  conserve_loads    = true;
-	bool  save_results;
-	int   newton;
-	int   solution_type;
-	int   numoutputs        = 0;
-	char** requested_outputs = NULL;
-	int    i;
-
+	bool       dakota_analysis;
+	int        meshtype;
+	bool       isSIA,isSSA,isL1L2,isHO,isFS;
+	bool       save_results;
+	int        solution_type;
+	int        numoutputs        = 0;
+	char     **requested_outputs = NULL;
+	Analysis  *analysis          = NULL;
 
 	/* recover parameters:*/
@@ -32,5 +30,4 @@
 	femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
 	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
-	femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
 	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
@@ -47,5 +44,5 @@
 	}
 
-	/*Compute slopes: */
+	/*Compute slopes if necessary */
 	if(isSIA || (isFS && meshtype==Mesh2DverticalEnum)) surfaceslope_core(femmodel);
 	if(isFS){
@@ -55,33 +52,23 @@
 	}
 
+	/*Compute SIA velocities*/
 	if(isSIA){
-		if(VerboseSolution()) _printf0_("   computing SIA velocities\n");
 
 		/*Take the last velocity into account so that the velocity on the SSA domain is not zero*/
 		if(isSSA || isL1L2 || isHO ) ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum);
-		femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
-		solutionsequence_linear(femmodel);
+
+		analysis = new StressbalanceSIAAnalysis();
+		analysis->Core(femmodel);
+		delete analysis;
+
+		/*Reset velocities for other ice flow models*/
 		if(isSSA || isL1L2 || isHO) ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum);
 	}
 
-	if ((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
-		if(VerboseSolution()) _printf0_("   computing velocities\n");
-
-		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
-		if(newton>0)
-		 solutionsequence_newton(femmodel);
-		else
-		 solutionsequence_nonlinear(femmodel,conserve_loads); 
-	}
-
-	if ((isSSA || isL1L2 || isHO) && isFS){
-		if(VerboseSolution()) _printf0_("   computing coupling betweem lower order models and full-FS\n");
-		solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
-	}
-
-	if (meshtype==Mesh3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
-		if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
-		femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
-		solutionsequence_linear(femmodel);
+	/*Compute stressbalance for SSA L1L2 HO and FS*/
+	if(isSSA || isL1L2 || isHO || isFS){
+		analysis = new StressbalanceAnalysis();
+		analysis->Core(femmodel);
+		delete analysis;
 	}