Index: /issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp	(revision 27243)
+++ /issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp	(revision 27244)
@@ -4,4 +4,6 @@
 #include "../classes/classes.h"
 #include "../classes/Inputs/TransientInput.h"
+#include "../classes/Inputs/TriaInput.h"
+#include "../classes/gauss/Gauss.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
@@ -100,39 +102,6 @@
 void           MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
-	/*retrieve thickness, ice and ocean levelsets from MmemasstransportAnalysis inputs in our element:*/
-
-	IssmDouble height,ice,ocean;
-	int       *doflist = NULL;
-
-	/*Fetch number of nodes and initialize values*/
-	int         numnodes = element->GetNumberOfNodes();
-	int         numdof   = numnodes*3;
-	IssmDouble* values   = xNew<IssmDouble>(numdof);
-
-	/*Get dof list and inputs */
-	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	Input* h_input=element->GetInput(MmemasstransportThicknessEnum); _assert_(h_input);
-	Input* i_input=element->GetInput(MmemasstransportMaskIceLevelsetEnum); _assert_(i_input);
-	Input* o_input=element->GetInput(MmemasstransportMaskOceanLevelsetEnum); _assert_(o_input);
-
-	/*Ok, we have the thickness in inputs, fill in solution */
-	Gauss* gauss=element->NewGauss();
-	for(int i=0;i<numnodes;i++){
-		gauss->GaussVertex(i);
-		h_input->GetInputValue(&height,gauss);
-		i_input->GetInputValue(&ice,gauss);
-		o_input->GetInputValue(&ocean,gauss);
-		values[3*i+0]=height;
-		values[3*i+1]=ice;
-		values[3*i+2]=ocean;
-	}
-
-	/*Add value to global vector*/
-	solution->SetValues(numdof,doflist,values,INS_VAL);
-
-	/*Free ressources:*/
-	delete gauss;
-	xDelete<int>(doflist);
-	xDelete<IssmDouble>(values);
+	/*do nothing:*/
+	return;
 }/*}}}*/
 void           MmemasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
@@ -141,48 +110,29 @@
 void           MmemasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
-	int         i,domaintype;
-	int*        doflist=NULL;
+	/*Fetch number of nodes and dof for this finite element*/
+	IssmDouble time;
+	IssmDouble ice[3];
+	IssmDouble ocean[3];
+	IssmDouble height;
+	int numnodes = element->GetNumberOfNodes();
+	
+	element->parameters->FindParam(&time,TimeEnum);
 
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-	int numdof   = numnodes*3;
+	TriaInput* h_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportThicknessEnum,time)); _assert_(h_input);
+	TriaInput* i_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskIceLevelsetEnum,time)); _assert_(i_input);
+	TriaInput* o_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskOceanLevelsetEnum,time)); _assert_(o_input);
 
-	/*Fetch dof list and allocate solution vectors*/
-	element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numdof);
-	IssmDouble* height        = xNew<IssmDouble>(numnodes);
-	IssmDouble* ice        = xNew<IssmDouble>(numnodes);
-	IssmDouble* ocean        = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Retrieve h:*/
-	for(i=0;i<numnodes;i++){
-		height[i]=values[3*i+0];
-		ice[i]=values[3*i+1];
-		ocean[i]=values[3*i+2];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(height[i])) _error_("NaN found in ice thickness solution vector");
-		if(xIsInf<IssmDouble>(height[i])) _error_("Inf found in ice thickness solution vector");
-		if(xIsNan<IssmDouble>(ice[i])) _error_("NaN found in ice mask solution vector");
-		if(xIsInf<IssmDouble>(ice[i])) _error_("Inf found in ice mask solution vector");
-		if(xIsNan<IssmDouble>(ocean[i])) _error_("NaN found in ocean mask solution vector");
-		if(xIsInf<IssmDouble>(ocean[i])) _error_("Inf found in ocean mask solution vector");
-
+	Gauss* gauss=element->NewGauss();
+	for(int iv=0;iv<3;iv++){
+		gauss->GaussVertex(iv);
+		i_input->GetInputValue(&ice[iv],gauss);
+		o_input->GetInputValue(&ocean[iv],gauss);
 	}
+	h_input->GetInputAverage(&height);
 
 	/*Add thickness ice and ocean levelsets as inputs to the tria element: */
-	element->AddInput(ThicknessEnum,height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0! 
+	element->AddInput(ThicknessEnum,&height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0! 
 	element->AddInput(MaskIceLevelsetEnum,ice,P1Enum);
 	element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum);
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(height);
-	xDelete<IssmDouble>(ice);
-	xDelete<IssmDouble>(ocean);
-	xDelete<IssmDouble>(values);
-	xDelete<int>(doflist);
 
 }/*}}}*/
