Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16935)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16936)
@@ -828,4 +828,6 @@
 		case SSAHOApproximationEnum: 
 			return CreateKMatrixSSAHO(element);
+		case HOFSApproximationEnum: 
+			return CreateKMatrixHOFS(element);
 		case NoneApproximationEnum:
 			return NULL;
@@ -3065,4 +3067,23 @@
 
 /*Coupling (Tiling)*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixFS(element);
+	int indices[3]={18,19,20};
+	Ke1->StaticCondensation(3,&indices[0]);
+	int init = element->FiniteElement();
+	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+	ElementMatrix* Ke2=CreateKMatrixHO(element);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	delete Ke3;
+	return Ke;
+}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
 
@@ -3077,4 +3098,70 @@
 	delete Ke2;
 	delete Ke3;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/
+
+	/*Constants*/
+	int numnodes       = 3*6+1;
+	int numdofp        = 2*6;
+	int numdofs        = 4*6 + 3;
+	int numdoftotal    = (2+4)*6+ 3;
+
+	/*Intermediaries*/
+	int   i,j,init;
+	int*   cs_list   = xNew<int>(6*3+1);
+	int*   cs_list2  = xNew<int>(6*2+1);
+	Node  **node_list = xNew<Node*>(6*3+1);
+
+	/*Some parameters needed*/
+	init = element->FiniteElement();
+
+	/*prepare node list*/
+	for(i=0;i<6+1;i++){
+		node_list[i+6] = element->GetNode(i);
+		cs_list[i+6]   = XYZEnum;
+		cs_list2[i]    = XYZEnum;
+	}
+	for(i=0;i<6;i++){
+		node_list[i]       = element->GetNode(i);
+		node_list[i+2*6+1] = element->GetNode(i+6*1);
+		cs_list[i]         = XYEnum;
+		cs_list[i+2*6+1]   = PressureEnum;
+		cs_list2[i+6+1]    = PressureEnum;
+	}
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=element->NewElementMatrix(FSvelocityEnum);
+	ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/*Compute HO Matrix with P1 element type\n");*/
+	Ke1=CreateKMatrixFS(element); element->TransformInvStiffnessMatrixCoord(Ke1,node_list,2*6+1,cs_list2);
+	int indices[3]={18,19,20};
+	Ke1->StaticCondensation(3,&indices[0]);
+	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+	Ke2=CreateKMatrixHO(element); element->TransformInvStiffnessMatrixCoord(Ke2,XYEnum);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	/*Compute FS Matrix and condense it \n");*/
+
+	for(i=0;i<numdofs;i++) for(j=0;j<6;j++){
+		Ke->values[i*numdoftotal+numdofs+2*j+0]+=Ke1->values[i*numdofs+3*j+0];
+		Ke->values[i*numdoftotal+numdofs+2*j+1]+=Ke1->values[i*numdofs+3*j+1];
+	}
+	for(i=0;i<numdofp;i++) for(j=0;j<6;j++){
+		Ke->values[(i+numdofs)*numdoftotal+3*j+0]+=Ke2->values[i*numdofp+2*j+0];
+		Ke->values[(i+numdofs)*numdoftotal+3*j+1]+=Ke2->values[i*numdofp+2*j+1];
+	}
+
+	/*Transform Coordinate System*/ //Do not transform, already done in the matrices
+	element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
+	/*clean-up and return*/
+	xDelete<int>(cs_list);
+	xDelete<int>(cs_list2);
+	xDelete<Node*>(node_list);
+	delete Ke1;
+	delete Ke2;
 	return Ke;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16935)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16936)
@@ -71,5 +71,7 @@
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
 		/*Coupling*/
+		ElementMatrix* CreateKMatrixHOFS(Element* element);
 		ElementMatrix* CreateKMatrixSSAHO(Element* element);
+		ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
 		ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
 		ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 16935)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 16936)
@@ -296,4 +296,10 @@
 	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
 }/*}}}*/
+void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
+	::TransformInvStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum);
+}/*}}}*/
+void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
+	::TransformInvStiffnessMatrixCoord(Ke,nodes_list,numnodes,transformenum_list);
+}/*}}}*/
 void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/
 	::TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),transformenum);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16935)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16936)
@@ -73,4 +73,6 @@
 		void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
 		void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
+		void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list);
 		void       TransformLoadVectorCoord(ElementVector* pe,int transformenum);
 		void       TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list);
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16935)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16936)
@@ -7054,8 +7054,8 @@
 	/*Compute HO Matrix with P1 element type\n");*/
 	this->element_type=P1Enum;
-	Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
+	Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,XYEnum);
 	this->element_type=init;
 	/*Compute FS Matrix and condense it \n");*/
-	Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,2*NUMVERTICES+1,cs_list2);
+	Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,node_list,2*NUMVERTICES+1,cs_list2);
 	int indices[3]={18,19,20};
 	Ke2->StaticCondensation(3,&indices[0]);