Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26500)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26501)
@@ -788,8 +788,13 @@
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.);
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.);
+	/*MLHO*/
 	if(isMLHO){
-		/*itapopo FIXME applying the same initialization for shear vx and shear vy for now*/
 		iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.);
 		iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.);
+		/*3D MLHO also need to have VxBase and VyBase for reconstruting Vx and Vy*/
+		if (iomodel->domaintype==Domain3DEnum) {
+			iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.);
+			iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.);
+		}
 	}
    if(iomodel->domaintype==Domain2DhorizontalEnum){
@@ -2754,14 +2759,31 @@
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
 
-	//_error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line.");
-
 	/* Check if ice in element */
 	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	int      domaintype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Domain3DEnum: 
+			if(!element->IsOnBase()) return NULL;
+			basalelement = element->SpawnBasalElement(true);
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixMLHOViscous(element);
-	ElementMatrix* Ke2=CreateKMatrixMLHOFriction(element);
+	ElementMatrix* Ke1=CreateKMatrixMLHOViscous(basalelement);
+	ElementMatrix* Ke2=CreateKMatrixMLHOFriction(basalelement);
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
 	/*clean-up and return*/
+	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
 	delete Ke1;
 	delete Ke2;
@@ -2849,50 +2871,9 @@
 	xDelete<IssmDouble>(basis);
 	return Ke;
-
-	//itapopo OLD - testing above
-	/*Intermediaries*/
-	//int      domaintype;
-   Element* basalelement;
-
-   /*Get basal element*/
-   element->FindParam(&domaintype,DomainTypeEnum);
-   switch(domaintype){
-      case Domain2DhorizontalEnum:
-         basalelement = element;
-         break;
-      case Domain3DEnum: case Domain2DverticalEnum:
-          _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-			//basalelement = element->SpawnBasalElement();
-         break;
-      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-   }
-
-	//Element* basalelement = element->SpawnBasalElement();
-	//ElementMatrix* Ke    = basalelement->NewElementMatrix(MLHOApproximationEnum);
-	ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33
-
-	/*Fetch number of nodes and dof for this finite element*/
-	//int numnodes = basalelement->GetNumberOfNodes();
-
-	for(int i=0;i<numnodes;i++){
-      for(int j=0;j<numnodes;j++){
-         Ke->values[(4*i+0)*2*2*numnodes+4*j+0] = KeSSA->values[2*i*2*numnodes+2*j]; //K11
-         Ke->values[(4*i+2)*2*2*numnodes+4*j+2] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K33
-      }  
-   }
-
-	KeSSA->Echo();
-	Ke->Echo();
-	_error_("mesh ");
-
-	/*Transform Coordinate System*/
-	//basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum);
-
-	/*clean-up and return*/
-	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
-	delete KeSSA;
-	return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
 
 	/*Intermediaries*/
@@ -2900,27 +2881,11 @@
 	IssmDouble *xyz_list = NULL;
 	int      domaintype;
-	Element* basalelement;
-
-	/*Get basal element*/
-	element->FindParam(&domaintype,DomainTypeEnum);
-	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			basalelement = element;
-			break;
-		case Domain3DEnum: case Domain2DverticalEnum:
-			_error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-			//basalelement = element->GetBasalElement()->SpawnBasalElement(true);
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
-
-	/*Get element on base*/
-	//Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true);
+	int dim=2;
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
+	int numnodes = element->GetNumberOfNodes();
 
 	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
+	ElementMatrix* Ke     = element->NewElementMatrix(MLHOApproximationEnum);
 	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes); // like SSA
@@ -2930,6 +2895,4 @@
 	Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input* surface_input  = element->GetInput(SurfaceEnum);   _assert_(surface_input);
-	//Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input); //vertically integrated vx
-	//Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input); //vertically integrated vy
 	Input* vxbase_input   = element->GetInput(VxBaseEnum);    _assert_(vxbase_input); 
 	Input* vybase_input   = element->GetInput(VyBaseEnum);    _assert_(vybase_input);
@@ -2940,15 +2903,12 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss      = element->NewGauss(2);
-	Gauss* gauss_base = basalelement->NewGauss();
 	while(gauss->next()){
-		gauss->SynchronizeGaussBase(gauss_base);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss_base);
-		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
-		basalelement->NodalFunctions(basis, gauss_base);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
 
 		thickness_input->GetInputValue(&thickness, gauss);
 		n_input->GetInputValue(&n,gauss);
-		int dim=2;//itapopo
 		element->material->ViscosityMLHO(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input);
 
@@ -3021,6 +2981,4 @@
 	/*Clean up and return*/
 	delete gauss;
-	delete gauss_base;
-	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(dbasis);
@@ -3196,16 +3154,20 @@
 	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
 
-	element->FindParam(&dim,DomainDimensionEnum);
 	element->FindParam(&domaintype,DomainTypeEnum);
 	rho_ice =element->FindParam(MaterialsRhoIceEnum);
 	g       =element->FindParam(ConstantsGEnum);
-	if(dim==2){
-		element->GetInputListOnVertices(thickness,ThicknessEnum);
-		for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
-	}
-	else{
-		element->GetVerticesCoordinates(&xyz_list);
-		element->GetInputListOnVertices(surface,SurfaceEnum);
-		for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			element->GetInputListOnVertices(thickness,ThicknessEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+			dim=2;
+			break;
+		case Domain3DEnum:
+			element->GetVerticesCoordinates(&xyz_list);
+			element->GetInputListOnVertices(surface,SurfaceEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+			dim=2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
 	element->AddInput(PressureEnum,pressure,P1Enum);
@@ -3217,12 +3179,10 @@
 	switch(domaintype){
 		case Domain2DhorizontalEnum:
-			// itapopo FIXME use basalelement or element only here
 			basalelement = element;
 			break;
 		case Domain3DEnum:
-			_error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-			//if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
-			//basalelement=element->SpawnBasalElement();
-			//break;
+			if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
+			basalelement=element->SpawnBasalElement();
+			break;
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
@@ -3230,9 +3190,8 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
-	int numdof   = numnodes*4; //4 DOFs per node
+	int numdof   = numnodes*dim*2; //2xdim DOFs per node
 
 	/*Fetch dof list and allocate solution vectors*/
 	basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); 
-	//basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numdof);
 	IssmDouble* vbx       = xNew<IssmDouble>(numnodes);
@@ -3247,16 +3206,12 @@
 	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
 	IssmDouble* n			 = xNew<IssmDouble>(numnodes);
+	IssmDouble* H			 = xNew<IssmDouble>(numnodes);
+	IssmDouble* s			 = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
-	//std::cout<<"MLHO  ********  Element ID="<<element->Id()<<"\n";
-	//for(i=0;i<numdof;i++){
-	//	std::cout<<values[i]*31536000<<"\n";
-	//}
-
 	/*Transform solution in Cartesian Space*/
-	basalelement->TransformSolutionCoord(&values[0],XYEnum);
-	//basalelement->FindParam(&domaintype,DomainTypeEnum);
+	if(dim==2) basalelement->TransformSolutionCoord(&values[0],XYEnum);
 
    /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -3292,17 +3247,27 @@
 	/*Compute the vertically averaged velocities on each node*/
 	basalelement->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); 
-   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
-		vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2);
-		vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
-	}
-
-	/*Get Vz and compute vel (vertically averaged vel)*/
-	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); 
-	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 
-
-	/*Add vx and vy as inputs to the tria element (vertically averaged velocities): */
-	element->AddBasalInput(VxEnum,vx,element->GetElementType());
-	element->AddBasalInput(VyEnum,vy,element->GetElementType());
-	element->AddBasalInput(VelEnum,vel,element->GetElementType()); 
+
+	/* Reconstruct vx, vy and vz solutions for 3D problem
+	 * Add vx and vy as inputs to the tria element (vertically averaged velocities): */
+
+   switch(domaintype){
+      case Domain2DhorizontalEnum:
+			for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
+				vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2);
+				vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
+				vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]); 
+			}
+			element->AddBasalInput(VxEnum,vx,element->GetElementType());
+			element->AddBasalInput(VyEnum,vy,element->GetElementType());
+			element->AddBasalInput(VelEnum,vel,element->GetElementType()); 
+         break;
+      case Domain3DEnum:
+		   basalelement->GetInputListOnNodes(&H[0],ThicknessEnum,0.);
+		   basalelement->GetInputListOnNodes(&s[0],SurfaceEnum,0.);
+			element->Recover3DMLHOInput(VxEnum, numnodes, vbx, vshx, n, H, s);
+			element->Recover3DMLHOInput(VyEnum, numnodes, vby, vshy, n, H, s);
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
 
 	/*Free ressources:*/
@@ -3320,4 +3285,6 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(n);
+	xDelete<IssmDouble>(s);
+	xDelete<IssmDouble>(H);
 	xDelete<int>(doflist);
 	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
@@ -3333,5 +3300,6 @@
 	switch(domaintype){
 		case Domain2DhorizontalEnum: dim = 2; dofpernode = 4; break;
-		case Domain2DverticalEnum: case Domain3DEnum: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); break;
+		case Domain3DEnum: dim = 2; dofpernode = 4; break;
+		case Domain2DverticalEnum: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); break;
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
@@ -3361,10 +3329,10 @@
 		vxbase_input->GetInputValue(&vbx,gauss);	//base vx
 		vxshear_input->GetInputValue(&vshx,gauss);//shear vx
-		values[i*4+0]=vbx;  //base vx
-		values[i*4+1]=vshx; //shear vx
+		values[i*dofpernode+0]=vbx;  //base vx
+		values[i*dofpernode+1]=vshx; //shear vx
 		vybase_input->GetInputValue(&vby,gauss);	//base vy
 		vyshear_input->GetInputValue(&vshy,gauss);//shear vy
-		values[i*4+2]=vby; //base vy
-		values[i*4+3]=vshy;//shear vy  
+		values[i*dofpernode+2]=vby; //base vy
+		values[i*dofpernode+3]=vshy;//shear vy  
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26500)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26501)
@@ -343,4 +343,5 @@
 		virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
 		virtual int        PressureInterpolation()=0;
+      virtual void       Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){_error_("not implemented yet");};
 		virtual void       ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
 		virtual void       ResetFSBasalBoundaryCondition()=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26500)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26501)
@@ -179,4 +179,39 @@
 		else _error_("not implemented yet");
 	}
+
+}
+/*}}}*/
+void       Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/
+   /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H
+    * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs
+    * To project to penta the DOFs are doubled in size
+    *
+    */
+   _assert_(this->inputs);
+   if(!IsOnBase()) return;
+   else{
+      if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){
+         IssmDouble vel[NUMVERTICES];
+         IssmDouble* xyz_list;
+         Penta* penta = this;
+         _assert_(NUMVERTICES-2*numnodes==0);
+
+         for(;;){
+            penta->GetVerticesCoordinates(&xyz_list);
+            for(int i=0;i<numnodes;i++) {
+               vel[i] = vb[i] + vsh[i]*(1.0-pow((s[i]-xyz_list[i*3+2])/H[i], (n[i]+1)));
+               vel[i+NUMVERTICES2D] = vb[i] + vsh[i]*(1-pow((s[i]-xyz_list[(i+NUMVERTICES2D)*3+2])/H[i], (n[i]+1)));
+            }
+
+            /*Add to the bottom side of the element*/
+            penta->AddInput(targetVel_enum,&vel[0],P1Enum);
+            if (penta->IsOnSurface()) break;
+            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
+         }
+
+         xDelete<IssmDouble>(xyz_list);
+      }
+      else _error_("not implemented yet");
+   }
 
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26500)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26501)
@@ -168,4 +168,5 @@
 		void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		int            PressureInterpolation();
+      void           Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s);
 		void           ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
 		void           ResetFSBasalBoundaryCondition(void);