Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18559)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18560)
@@ -2942,8 +2942,6 @@
 	if(fe_FS==XTaylorHoodEnum)
 	 Ke1=CreateKMatrixFSViscousXTH(element);
-	else if(fe_FS==LATaylorHoodEnum)
-	 Ke1=CreateKMatrixFSViscousLATH(element);
-	else if(fe_FS==LACrouzeixRaviartEnum)
-	 Ke1=CreateKMatrixFSViscousLACR(element);
+	else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
+	 Ke1=CreateKMatrixFSViscousLA(element);
 	else
 	 Ke1=CreateKMatrixFSViscous(element);
@@ -2959,5 +2957,5 @@
 	return Ke;
 }/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLATH(Element* element){/*{{{*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLA(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -2977,134 +2975,4 @@
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->GetNumberOfNodes(P1Enum);
-	int lnumnodes = element->GetNumberOfNodes(P2Enum);
-	int numdof    = vnumnodes*dim;
-	int pnumdof   = pnumnodes;
-	int lnumdof   = lnumnodes;
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(vnumnodes);
-	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
-	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
-
-	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke       = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B        = xNew<IssmDouble>(epssize*numdof);
-	IssmDouble*    Bprime   = xNew<IssmDouble>(epssize*numdof);
-	IssmDouble*    BtBUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
-	IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
-	IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
-	IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input* vx_input = element->GetInput(VxEnum);     _assert_(vx_input);
-	Input* vy_input = element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input = NULL;
-	if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss = element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBFSvel(B,element,dim,xyz_list,gauss);
-		this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss);
-
-		element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-		for(i=0;i<epssize;i++)   D[i*epssize+i] = 2*viscosity*gauss->weight*Jdet;
-
-		TripleMultiply(B,epssize,numdof,1,
-					D,epssize,epssize,0,
-					Bprime,epssize,numdof,0,
-					&Ke->values[0],1);
-
-		this->GetBFSUzawa(BU,element,dim,xyz_list,gauss);
-		this->GetBFSprimeUzawa(BprimeU,element,dim,xyz_list,gauss);
-
-		DU = gauss->weight*Jdet*sqrt(r);
-
-		TripleMultiply(BU,1,pnumdof,1,
-					&DU,1,1,0,
-					BprimeU,1,numdof,0,
-					BtBUzawa,1);
-	}
-
-	if(element->IsOnBase()){ 
-		element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
-		element->GetVerticesCoordinatesBase(&xyz_list_base);
-		element->NormalBase(&normal[0],xyz_list_base);
-
-		IssmDouble* Dlambda  = xNewZeroInit<IssmDouble>(dim*dim);
-		IssmDouble* C        = xNewZeroInit<IssmDouble>(dim*lnumdof);
-		IssmDouble* Cprime   = xNewZeroInit<IssmDouble>(dim*numdof);
-		IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
-
-		delete gauss;
-		gauss = element->NewGaussBase(5);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
-
-			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-			this->GetCFS(C,element,dim,xyz_list,gauss);
-			this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
-			for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
-			TripleMultiply(C,dim,lnumdof,1,
-						Dlambda,dim,dim,0,
-						Cprime,dim,numdof,0,
-						CtCUzawa,1);
-		}
-
-		/*The sigma naugmentation should not be transformed*/
-		MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
-					CtCUzawa,lnumdof,numdof,0,
-					&Ke->values[0],1);
-
-		/*Delete base part*/
-		xDelete<IssmDouble>(Dlambda);
-		xDelete<IssmDouble>(C);
-		xDelete<IssmDouble>(Cprime);
-		xDelete<IssmDouble>(CtCUzawa);
-		xDelete<IssmDouble>(xyz_list_base);
-	}
-
-	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,cs_list);
-
-	/*The pressure augmentation should not be transformed*/
-	MatrixMultiply(BtBUzawa,pnumdof,numdof,1,
-				BtBUzawa,pnumdof,numdof,0,
-				&Ke->values[0],1);
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(BprimeU);
-	xDelete<IssmDouble>(BU);
-	xDelete<IssmDouble>(BtBUzawa);
-	xDelete<int>(cs_list);
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLACR(Element* element){/*{{{*/
-
-	/*Intermediaries*/
-	int         i,dim,epssize;
-	IssmDouble  r,rl,Jdet,viscosity,DU,DUl;
-	IssmDouble	normal[3];
-	IssmDouble *xyz_list = NULL;
-	IssmDouble *xyz_list_base = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-	element->FindParam(&r,AugmentedLagrangianREnum);
-	if(dim==2) epssize = 3;
-	else       epssize = 6;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->GetNumberOfNodes(P1DGEnum);
 	int lnumnodes = element->GetNumberOfNodes(P2Enum);
 	int numdof    = vnumnodes*dim;
@@ -3808,23 +3676,10 @@
 		delete pe4;
 	}
-	else if(fe_FS==LATaylorHoodEnum){
+	else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
 		ElementVector* pe1=CreatePVectorFSViscous(element);
 		ElementVector* pe2=CreatePVectorFSShelf(element);
 		ElementVector* pe3=CreatePVectorFSFront(element);
 		ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
-		ElementVector* pe4=CreatePVectorFSViscousLATH(element);
-		pe = new ElementVector(petemp,pe4);
-		delete pe1;
-		delete pe2;
-		delete pe3;
-		delete petemp;
-		delete pe4;
-	}
-	else if(fe_FS==LACrouzeixRaviartEnum){
-		ElementVector* pe1=CreatePVectorFSViscous(element);
-		ElementVector* pe2=CreatePVectorFSShelf(element);
-		ElementVector* pe3=CreatePVectorFSFront(element);
-		ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
-		ElementVector* pe4=CreatePVectorFSViscousLACR(element);
+		ElementVector* pe4=CreatePVectorFSViscousLA(element);
 		pe = new ElementVector(petemp,pe4);
 		delete pe1;
@@ -3917,4 +3772,87 @@
 }/*}}}*/
 #endif
+ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLA(Element* element){/*{{{*/
+
+	int         i,dim;
+	IssmDouble  Jdet,r,pressure;
+	IssmDouble  bed_normal[3];
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(numnodes);
+	if(dim==2) for(i=0;i<numnodes;i++) cs_list[i] = XYEnum;
+	else       for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
+
+	/*Initialize vectors*/
+	ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    dbasis  = xNew<IssmDouble>(3*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->FindParam(&r,AugmentedLagrangianREnum);
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/*Get pressure and sigmann*/
+	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+	Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
+
+	gauss=element->NewGauss(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		
+		pressure_input->GetInputValue(&pressure, gauss);
+		element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
+
+		for(i=0;i<numnodes;i++){
+			pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
+			pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
+			if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
+		}
+	}
+
+	if(element->IsOnBase()){ 
+		IssmDouble   sigmann;
+		IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
+
+		element->GetVerticesCoordinatesBase(&xyz_list_base);
+		element->NormalBase(&bed_normal[0],xyz_list_base);
+
+		delete gauss;
+		gauss=element->NewGaussBase(5);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+			gauss->GaussPoint(ig);
+
+			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+			element->NodalFunctionsVelocity(vbasis,gauss);
+			sigmann_input->GetInputValue(&sigmann, gauss);
+
+			for(i=0;i<numnodes;i++){
+				pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
+				pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
+				if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
+			}
+		}
+		xDelete<IssmDouble>(xyz_list_base);
+		xDelete<IssmDouble>(vbasis);
+	}
+
+	/*Transform coordinate system*/
+	//element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(dbasis);
+	return pe;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
 
@@ -4087,170 +4025,4 @@
 	xDelete<IssmDouble>(vdbasis);
 	xDelete<IssmDouble>(tbasis);
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLATH(Element* element){/*{{{*/
-
-	int         i,dim;
-	IssmDouble  Jdet,r,pressure;
-	IssmDouble  bed_normal[3];
-	IssmDouble *xyz_list      = NULL;
-	IssmDouble *xyz_list_base = NULL;
-	Gauss*      gauss         = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(numnodes);
-	if(dim==2) for(i=0;i<numnodes;i++) cs_list[i] = XYEnum;
-	else       for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
-
-	/*Initialize vectors*/
-	ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    dbasis  = xNew<IssmDouble>(3*numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->FindParam(&r,AugmentedLagrangianREnum);
-	element->GetVerticesCoordinates(&xyz_list);
-
-	/*Get pressure and sigmann*/
-	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
-
-	gauss=element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		
-		pressure_input->GetInputValue(&pressure, gauss);
-		element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
-
-		for(i=0;i<numnodes;i++){
-			pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
-			pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
-			if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
-		}
-	}
-
-	if(element->IsOnBase()){ 
-		IssmDouble   sigmann;
-		IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
-
-		element->GetVerticesCoordinatesBase(&xyz_list_base);
-		element->NormalBase(&bed_normal[0],xyz_list_base);
-
-		delete gauss;
-		gauss=element->NewGaussBase(5);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
-
-			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-			element->NodalFunctionsVelocity(vbasis,gauss);
-			sigmann_input->GetInputValue(&sigmann, gauss);
-
-			for(i=0;i<numnodes;i++){
-				pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
-				pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
-				if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
-			}
-		}
-		xDelete<IssmDouble>(xyz_list_base);
-		xDelete<IssmDouble>(vbasis);
-	}
-
-	/*Transform coordinate system*/
-	//element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<int>(cs_list);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(dbasis);
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLACR(Element* element){/*{{{*/
-
-	int         i,dim;
-	IssmDouble  Jdet,r,pressure;
-	IssmDouble  bed_normal[3];
-	IssmDouble *xyz_list      = NULL;
-	IssmDouble *xyz_list_base = NULL;
-	Gauss*      gauss         = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Prepare coordinate system list*/
-	int* cs_list = xNew<int>(numnodes);
-	if(dim==2) for(i=0;i<numnodes;i++) cs_list[i] = XYEnum;
-	else       for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
-
-	/*Initialize vectors*/
-	ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    dbasis  = xNew<IssmDouble>(3*numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->FindParam(&r,AugmentedLagrangianREnum);
-	element->GetVerticesCoordinates(&xyz_list);
-
-	/*Get pressure and sigmann*/
-	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
-
-	gauss=element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		
-		pressure_input->GetInputValue(&pressure, gauss);
-		element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
-
-		for(i=0;i<numnodes;i++){
-			pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
-			pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
-			if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
-		}
-	}
-
-	if(element->IsOnBase()){ 
-		IssmDouble   sigmann;
-		IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
-
-		element->GetVerticesCoordinatesBase(&xyz_list_base);
-		element->NormalBase(&bed_normal[0],xyz_list_base);
-
-		delete gauss;
-		gauss=element->NewGaussBase(5);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
-
-			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-			element->NodalFunctionsVelocity(vbasis,gauss);
-			sigmann_input->GetInputValue(&sigmann, gauss);
-
-			for(i=0;i<numnodes;i++){
-				pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
-				pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
-				if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
-			}
-		}
-		xDelete<IssmDouble>(xyz_list_base);
-		xDelete<IssmDouble>(vbasis);
-	}
-
-	/*Transform coordinate system*/
-	//element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<int>(cs_list);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(dbasis);
 	return pe;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18559)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18560)
@@ -69,6 +69,5 @@
 		ElementMatrix* CreateJacobianMatrixFS(Element* element);
 		ElementMatrix* CreateKMatrixFS(Element* element);
-		ElementMatrix* CreateKMatrixFSViscousLATH(Element* element);
-		ElementMatrix* CreateKMatrixFSViscousLACR(Element* element);
+		ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
 		ElementMatrix* CreateKMatrixFSViscous(Element* element);
@@ -77,6 +76,5 @@
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);
-		ElementVector* CreatePVectorFSViscousLATH(Element* element);
-		ElementVector* CreatePVectorFSViscousLACR(Element* element);
+		ElementVector* CreatePVectorFSViscousLA(Element* element);
 		ElementVector* CreatePVectorFSViscousXTH(Element* element);
 		ElementVector* CreatePVectorFSShelf(Element* element);