Index: /issm/trunk/src/c/modules/MaxAbsVxx/MaxAbsVxx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxAbsVxx/MaxAbsVxx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxAbsVxx/MaxAbsVxx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxAbsVxx( double* pmaxabsvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxAbsVyx/MaxAbsVyx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxAbsVyx/MaxAbsVyx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxAbsVyx/MaxAbsVyx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxAbsVyx( double* pmaxabsvy, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxAbsVzx/MaxAbsVzx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxAbsVzx/MaxAbsVzx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxAbsVzx/MaxAbsVzx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxAbsVzx( double* pmaxabsvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxVelx/MaxVelx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxVelx/MaxVelx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxVelx/MaxVelx.cpp	(revision 5414)
@@ -9,5 +9,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxVelx( double* pmaxvel, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxVxx/MaxVxx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxVxx/MaxVxx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxVxx/MaxVxx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxVxx( double* pmaxvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxVyx/MaxVyx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxVyx/MaxVyx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxVyx/MaxVyx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxVyx( double* pmaxvy, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MaxVzx/MaxVzx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MaxVzx/MaxVzx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MaxVzx/MaxVzx.cpp	(revision 5414)
@@ -9,5 +9,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MaxVzx( double* pmaxvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MinVelx/MinVelx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MinVelx/MinVelx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MinVelx/MinVelx.cpp	(revision 5414)
@@ -9,5 +9,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MinVelx( double* pminvel, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MinVxx/MinVxx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MinVxx/MinVxx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MinVxx/MinVxx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MinVxx( double* pminvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MinVyx/MinVyx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MinVyx/MinVyx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MinVyx/MinVyx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MinVyx( double* pminvy, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/MinVzx/MinVzx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MinVzx/MinVzx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/MinVzx/MinVzx.cpp	(revision 5414)
@@ -8,5 +8,4 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
-#include "../SurfaceAreax/SurfaceAreax.h"
 
 void MinVzx( double* pminvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){
Index: /issm/trunk/src/c/modules/SurfaceAreax/SurfaceAreax.cpp
===================================================================
--- /issm/trunk/src/c/modules/SurfaceAreax/SurfaceAreax.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/SurfaceAreax/SurfaceAreax.cpp	(revision 5414)
@@ -9,4 +9,5 @@
 #include "../../toolkits/toolkits.h"
 #include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h"
 
 void SurfaceAreax( double* pS, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
@@ -31,5 +32,8 @@
 	S=S_sum;
 
+	/*add surface area to element inputs:*/
+	InputUpdateFromConstantx(elements,nodes,vertices,loads,materials,parameters,S,SurfaceAreaEnum);
+
 	/*Assign output pointers: */
-	*pS=S;
+	if(pS) *pS=S;
 }
Index: /issm/trunk/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp
===================================================================
--- /issm/trunk/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp	(revision 5413)
+++ /issm/trunk/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp	(revision 5414)
@@ -14,5 +14,4 @@
 
 	/*Intermediary*/
-	double S;
 	int i;
 	Element* element=NULL;
@@ -22,12 +21,6 @@
 	double J_sum;
 
-	/*Compute surface area: */
-	SurfaceAreax(&S,elements,nodes,vertices, loads,materials,parameters);
-
-	/*add surface area to elements :*/
-	for(i=0;i<elements->Size();i++){
-		Element* element=(Element*)elements->GetObjectByOffset(i);
-		element->InputUpdateFromVector(&S,SurfaceAreaEnum,ConstantEnum);
-	}
+	/*Compute surface area and add to elements inputs */
+	SurfaceAreax(NULL,elements,nodes,vertices, loads,materials,parameters);
 
 	/*Compute Misfit: */
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5413)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5414)
@@ -6136,4 +6136,5 @@
 				name==AccumulationRateEnum ||
 				name==GeothermalFluxEnum ||
+				name==SurfaceAreaEnum||
 				name==PressureEnum ||
 				name==VxEnum ||
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5413)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5414)
@@ -7184,4 +7184,5 @@
 				name==MeltingRateEnum ||
 				name==AccumulationRateEnum ||
+				name==SurfaceAreaEnum||
 				name==ControlParameterEnum ||
 				name==VxEnum ||
Index: /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp	(revision 5413)
+++ /issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp	(revision 5414)
@@ -24,14 +24,18 @@
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
-	/*set analysis type to compute velocity: */
+	/*compute thickness */
 	_printf_("%s\n","      computing thickness");
 	femmodel->SetCurrentConfiguration(BalancedthicknessAnalysisEnum);
 	solver_linear(femmodel);
 
-	/*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
+	/*Call SurfaceAreax, because some it might be needed by PVector*/
+	SurfaceAreax(NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+	/*compute adjoint*/
 	_printf_("%s\n","      computing adjoint");
 	femmodel->SetCurrentConfiguration(BalancedthicknessAnalysisEnum,AdjointBalancedthicknessAnalysisEnum);
 	solver_adjoint_linear(femmodel);
 	
+	/*Save results*/
 	if(solution_type==AdjointSolutionEnum && !control_analysis){
 		if(verbose)_printf_("saving results:\n");
Index: /issm/trunk/src/c/solutions/adjointdiagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/adjointdiagnostic_core.cpp	(revision 5413)
+++ /issm/trunk/src/c/solutions/adjointdiagnostic_core.cpp	(revision 5414)
@@ -27,5 +27,5 @@
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
-	/*set analysis type to compute velocity: */
+	/*Compute velocities*/
 	_printf_("%s\n","      computing velocities");
 	if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
@@ -33,5 +33,8 @@
 	solver_diagnostic_nonlinear(femmodel,conserve_loads); 
 
-	/*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
+	/*Call SurfaceAreax, because some it might be needed by PVector*/
+	SurfaceAreax(NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+	/*Compute adjoint*/
 	_printf_("%s\n","      computing adjoint");
 	if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
@@ -39,4 +42,5 @@
 	solver_adjoint_linear(femmodel);
 	
+	/*Save results*/
 	if(solution_type==AdjointSolutionEnum && !control_analysis){
 		if(verbose)_printf_("saving results:\n");