Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25944)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25945)
@@ -5048,5 +5048,4 @@
 	/*coefficient of Nitsche's method*/
 	IssmDouble gamma, hx, hy, hz;
-	IssmDouble pnu, qnv;
 
 	/*Get problem dimension*/
@@ -5065,7 +5064,4 @@
 	IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
 	IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
-	IssmDouble* vdbasisn = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* nsigma = xNew<IssmDouble>(numdof);
-	IssmDouble  bed_normal[3], bed_tangent[6];
 
 	/*Prepare coordinate system list*/
@@ -5100,6 +5096,4 @@
 	}
 
-	element->NormalBase(&bed_normal[0],xyz_list_base);
-	element->TangentBase(&bed_tangent[0], &bed_normal[0]);
 	element->ElementSizes(&hx,&hy,&hz);
 	element->FindParam(&gamma,FeFSNitscheGammaEnum);
@@ -5126,25 +5120,13 @@
 		element->NodalFunctionsPressure(pbasis,gauss);
 		element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-
-		/*Frictions*/
+		
 		for(int i=0;i<vnumnodes;i++){
 			for(int j=0;j<vnumnodes;j++){
-				for (int d=0; d<2; d++) {
-					for (int p=0; p<dim; p++) {
-						for (int q=0; q<dim; q++) {
-							Ke->values[(dim*i+p)*numdof+dim*j+q] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j]*bed_tangent[d*3+p]*bed_tangent[d*3+q];
-						}
-					}
+				for (int d=0; d<dim; d++) {
+						Ke->values[(dim*i+d)*numdof+dim*j+d] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
 				}
 			}
 		}
 		/* -------- Nitsche terms -------- */
-		/* n\cdot\nabla\phi*/
-		for(int i=0;i<vnumnodes;i++){
-			vdbasisn[i] = 0.0;
-			for (int d=0; d<dim; d++) {
-				vdbasisn[i] += bed_normal[d]*vdbasis[d*vnumnodes+i];
-			}
-		}
 		/* Boundary terms for integration by parts.
 			The coefficient matrix of n*sigma*n-gamma*n*u is stored in the following order:
@@ -5153,50 +5135,26 @@
 		*/
 		for(int i=0;i<vnumnodes;i++){
-			for(int d=0;d<dim;d++){
-				nsigma[d*vnumnodes+i] = 1.0 * viscosity * bed_normal[d]*vdbasisn[i];
+			for(int j=0;j<vnumnodes;j++){
+				/* gamma*(n*u)*(n*v) */
+				Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += gauss->weight*Jdet * (gamma*viscosity) * vbasis[j] * vbasis[i];
+				/* -sigma(v)*(n*u) */
+				Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += (- gauss->weight*Jdet) * 2.0 * viscosity * (-vbasis[j]) * vdbasis[(dim-1)*vnumnodes+i];
+				/* -sigma(u)*(n*v) */
+				Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += (- gauss->weight*Jdet) * 2.0 * viscosity * vdbasis[(dim-1)*vnumnodes+j] * (-vbasis[i]);
 			}
 		}
-		for(int i=0;i<pnumnodes;i++){
-			nsigma[dim*vnumnodes+i] = -FSreconditioning*pbasis[i];
-		}
-
-		/*velocity components A11 */
-		for(int i=0;i<vnumnodes;i++){
-			for (int p=0; p<dim; p++) {
-				for(int j=0;j<vnumnodes;j++){
-					for (int q=0; q<dim; q++) {
-						/* gamma*(n*u)*(n*v) */
-						Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
-						/* sigma(u)*(n*v) */
-						Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j];
-						/* sigma(v)*(n*u) */
-						Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* bed_normal[p] * vbasis[i] * nsigma[q*vnumnodes+j];
-					}
-				}
+		/* pressure x velocity  component A12, +p*(n*v) */
+		for(int k=0;k<pnumnodes;k++){
+			for(int i=0;i<vnumnodes;i++){
+				Ke->values[(dim*i+dim-1)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * FSreconditioning*pbasis[k] * (-vbasis[i]);
 			}
 		}
-
-		/* pressure x velocity  component A12 */
+		/* velocity x pressure component A21, +(n*u)*q */
 		for(int k=0;k<pnumnodes;k++){
-			pnu = nsigma[dim*vnumnodes+k];
-			for(int i=0;i<vnumnodes;i++){
-				for (int p=0;p<dim;p++) {
-					Ke->values[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * pnu * (- bed_normal[p] * vbasis[i]);
-				}
+			for(int j=0;j<vnumnodes;j++){
+				Ke->values[(dim*vnumnodes+k)*numdof+dim*j+dim-1] += gauss->weight*Jdet * (-vbasis[j]) * FSreconditioning*pbasis[k];
 			}
 		}
-		/* velocity x pressure component A21 */
-		for(int k=0;k<pnumnodes;k++){
-			qnv = nsigma[dim*vnumnodes+k];
-			for(int j=0;j<vnumnodes;j++){
-				for (int q=0;q<dim;q++) {
-					Ke->values[(dim*vnumnodes+k)*numdof+dim*j+q] += gauss->weight*Jdet * ( - bed_normal[q] * vbasis[j]) * qnv;
-				}
-			}
-		}
-	}
-
-	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,cs_list);
+	}
 
 	/*Clean up and return*/
@@ -5208,7 +5166,4 @@
 	xDelete<IssmDouble>(vdbasis);
 	xDelete<IssmDouble>(pbasis);
-	xDelete<IssmDouble>(vdbasisn);
-	xDelete<IssmDouble>(nsigma);
-	xDelete<int>(cs_list);
 	return Ke;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25944)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25945)
@@ -3327,4 +3327,5 @@
 	int          numindices;
 	int         *indices = NULL;
+	bool		isNitsche;
 	IssmDouble   slopex,slopey,groundedice;
 	IssmDouble   xz_plane[6];
@@ -3333,4 +3334,6 @@
 	/*For FS only: we want the CS to be tangential to the bedrock*/
 	this->Element::GetInputValue(&approximation,ApproximationEnum);
+	this->FindParam(&isNitsche,FlowequationIsNitscheEnum);
+
 	if(!IsOnBase() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
 
@@ -3357,23 +3360,24 @@
 		xz_plane[2]=slopex;   xz_plane[5]=1.;
 
-		if(groundedice>=0){
-			if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
-				this->nodes[indices[i]]->DofInSSet(2); //vz
+		if (!isNitsche) {
+			if(groundedice>=0){
+				if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
+					this->nodes[indices[i]]->DofInSSet(2); //vz
+				}
+				else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
+					this->nodes[indices[i]]->DofInSSet(4); //vz
+				}
+				else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
 			}
-			else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
-				this->nodes[indices[i]]->DofInSSet(4); //vz
+			else{
+				if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
+					this->nodes[indices[i]]->DofInFSet(2); //vz
+				}
+				else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
+					this->nodes[indices[i]]->DofInFSet(4); //vz
+				}
+				else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
 			}
-			else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
-		}
-		else{
-			if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
-				this->nodes[indices[i]]->DofInFSet(2); //vz
-			}
-			else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
-				this->nodes[indices[i]]->DofInFSet(4); //vz
-			}
-			else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
-		}
-
+		}
 		XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
 	}
Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 25944)
+++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 25945)
@@ -47,7 +47,5 @@
 		bedslope_core(femmodel);
 		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
-		if (!isNitsche) {
-			ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-		}
+		ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 
 		/*We need basal melt rates for the shelf dampening*/