Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18473)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18474)
@@ -69,6 +69,7 @@
 		ElementMatrix* CreateJacobianMatrixFS(Element* element);
 		ElementMatrix* CreateKMatrixFS(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousLATH(Element* element);
+		ElementMatrix* CreateKMatrixFSViscousLACR(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
 		ElementMatrix* CreateKMatrixFSViscous(Element* element);
 		ElementMatrix* CreateKMatrixFSFriction(Element* element);
@@ -76,6 +77,7 @@
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);
 		ElementVector* CreatePVectorFSViscousLATH(Element* element);
+		ElementVector* CreatePVectorFSViscousLACR(Element* element);
 		ElementVector* CreatePVectorFSViscousXTH(Element* element);
 		ElementVector* CreatePVectorFSShelf(Element* element);
 		ElementVector* CreatePVectorFSFront(Element* element);
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18474)
@@ -111,7 +111,7 @@
 
 	/*XTH LATH parameters*/
 	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
-	if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum){
+	if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
 		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
 		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
@@ -246,7 +246,7 @@
 	}
 	/*LATH parameters*/
 	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
-	if(fe_FS==LATaylorHoodEnum){
+	if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
 		iomodel->FetchDataToInput(elements,PressureEnum,0.);
 		InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
 	}
@@ -427,15 +427,16 @@
 		else if(isFS){  iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
 			/*Deduce velocity interpolation from finite element*/
 			switch(finiteelement){
-				case P1P1Enum            : finiteelement = P1Enum;       break;
-				case P1P1GLSEnum         : finiteelement = P1Enum;       break;
-				case MINIcondensedEnum   : finiteelement = P1bubbleEnum; break;
-				case MINIEnum            : finiteelement = P1bubbleEnum; break;
-				case TaylorHoodEnum      : finiteelement = P2Enum;       break;
-				case XTaylorHoodEnum     : finiteelement = P2Enum;       break;
-				case LATaylorHoodEnum    : finiteelement = P2Enum;       break;
-				case OneLayerP4zEnum     : finiteelement = P2xP4Enum;    break;
-				case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
+				case P1P1Enum              : finiteelement = P1Enum;       break;
+				case P1P1GLSEnum           : finiteelement = P1Enum;       break;
+				case MINIcondensedEnum     : finiteelement = P1bubbleEnum; break;
+				case MINIEnum              : finiteelement = P1bubbleEnum; break;
+				case TaylorHoodEnum        : finiteelement = P2Enum;       break;
+				case XTaylorHoodEnum       : finiteelement = P2Enum;       break;
+				case LATaylorHoodEnum      : finiteelement = P2Enum;       break;
+				case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
+				case OneLayerP4zEnum       : finiteelement = P2xP4Enum;    break;
+				case CrouzeixRaviartEnum   : finiteelement = P2bubbleEnum; break;
 				default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
 			}
 		}
@@ -881,7 +882,7 @@
 		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
 		if (fe_FS==XTaylorHoodEnum)
 		 solutionsequence_la_theta(femmodel);
-		else if (fe_FS==LATaylorHoodEnum)
+		else if (fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
 		 solutionsequence_la(femmodel);
 		else if(newton>0)
 		 solutionsequence_newton(femmodel);
@@ -2942,6 +2943,8 @@
 	 Ke1=CreateKMatrixFSViscousXTH(element);
 	else if(fe_FS==LATaylorHoodEnum)
 	 Ke1=CreateKMatrixFSViscousLATH(element);
+	else if(fe_FS==LACrouzeixRaviartEnum)
+	 Ke1=CreateKMatrixFSViscousLACR(element);
 	else
 	 Ke1=CreateKMatrixFSViscous(element);
 
@@ -3085,6 +3088,136 @@
 	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;
+	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() && 0){ 
+		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::CreateKMatrixFSViscousXTH(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -3450,6 +3583,13 @@
 		delete pe2;
 		return pe3;
 	}
+	else if(fe_FS==LACrouzeixRaviartEnum){
+		ElementVector* pe2=CreatePVectorFSViscousLACR(element);
+		ElementVector* pe3 = new ElementVector(pe,pe2);
+		delete pe;
+		delete pe2;
+		return pe3;
+	}
 	return pe;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
@@ -3680,6 +3820,19 @@
 		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);
+		pe = new ElementVector(petemp,pe4);
+		delete pe1;
+		delete pe2;
+		delete pe3;
+		delete petemp;
+		delete pe4;
+	}
 	else{
 		ElementVector* pe1=CreatePVectorFSViscous(element);
 		ElementVector* pe2=CreatePVectorFSShelf(element);
@@ -3986,6 +4139,89 @@
 		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);
 
Index: ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 18474)
@@ -16,8 +16,15 @@
 void UzawaPressureAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	/*Update elements: */
-	int finiteelement = P1Enum;
+	int finiteelement;
 	int counter=0;
+	int fe_FS;
+
+	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS=LATaylorHoodEnum) finiteelement = P1Enum;
+	else if(fe_FS=LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
+	else _error_("solution not supported yet");
+
 	for(int i=0;i<iomodel->numberofelements;i++){
 		if(iomodel->my_elements[i]){
 			Element* element=(Element*)elements->GetObjectByOffset(counter);
@@ -34,7 +41,14 @@
 }/*}}}*/
 void UzawaPressureAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
-	int finiteelement = P1Enum;
+	int finiteelement;
+	int fe_FS;
+
+	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS=LATaylorHoodEnum) finiteelement = P1Enum;
+	else if(fe_FS=LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
+	else _error_("solution not supported yet");
+
 	::CreateNodes(nodes,iomodel,UzawaPressureAnalysisEnum,finiteelement);
 }/*}}}*/
 void UzawaPressureAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -175,7 +189,10 @@
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes      = element->GetNumberOfNodes();
-	int numnodessigma = element->GetNumberOfNodes(P2Enum);
+	int numnodessigma;
+	if(element->element_type==P1Enum) numnodessigma=element->GetNumberOfNodes(P2Enum);
+	else if(element->element_type==P1DGEnum) numnodessigma=element->GetNumberOfNodes(P2Enum);
+	else _error_("finite element not supported yet");
 
 	element->FindParam(&dim,DomainDimensionEnum);
 
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18473)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18474)
@@ -608,6 +608,7 @@
 	XTaylorHoodEnum,
 	OneLayerP4zEnum,
 	CrouzeixRaviartEnum,
+	LACrouzeixRaviartEnum,
 	/*}}}*/
 	/*Results{{{*/
 	SaveResultsEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18474)
@@ -597,6 +597,7 @@
 		case XTaylorHoodEnum : return "XTaylorHood";
 		case OneLayerP4zEnum : return "OneLayerP4z";
 		case CrouzeixRaviartEnum : return "CrouzeixRaviart";
+		case LACrouzeixRaviartEnum : return "LACrouzeixRaviart";
 		case SaveResultsEnum : return "SaveResults";
 		case BoolExternalResultEnum : return "BoolExternalResult";
 		case DoubleExternalResultEnum : return "DoubleExternalResult";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18474)
@@ -609,6 +609,7 @@
 	      else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
 	      else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
 	      else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
+	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
 	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
 	      else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum;
 	      else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
 	      else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
-	      else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+	      if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
+	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
 	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp	(revision 18474)
@@ -532,7 +532,48 @@
 				}
 			}
 			break;
+		case LACrouzeixRaviartEnum:
+			_assert_(approximation==FSApproximationEnum);
+			/*P2b velocity*/
+			EdgesPartitioning(&my_edges,iomodel);
+			for(i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			for(i=0;i<iomodel->numberofedges;i++){
+				if(my_edges[i]){
+					nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
+				}
+			}
+			id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
+			if(iomodel->meshelementtype==PentaEnum){
+				FacesPartitioning(&my_faces,iomodel);
+				for(i=0;i<iomodel->numberoffaces;i++){
+					if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
+						if(my_faces[i]){
+							node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
+							nodes->AddObject(node);
+						}
+					}
+					else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
+						/*Nothing*/
+					}
+					else{
+						_error_("not supported");
+					}
+				}
+				id0 = id0+iomodel->numberoffaces;
+			}
+			for(i=0;i<iomodel->numberofelements;i++){
+				if(iomodel->my_elements[i]){
+					nodes->AddObject(new Node(id0+i+1,id0-iomodel->nodecounter+i,lid++,0,iomodel,analysis,FSvelocityEnum));
+				}
+			}
 
+			/*No pressure*/
+			break;
+
 		default:
 			_error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
 	}
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18474)
@@ -2050,6 +2050,17 @@
 			tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+2;
 			tria_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+3;
 			break;
+		case LACrouzeixRaviartEnum:
+			numnodes        = 7;
+			tria_node_ids   = xNew<int>(numnodes);
+			tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
+			tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
+			tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
+			tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
+			tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
+			tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
+			tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+index+1;
+			break;
 		default:
 			_error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
 	}
Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 18474)
@@ -965,6 +965,7 @@
 		case LATaylorHoodEnum:      return NUMNODESP2;
 		case OneLayerP4zEnum:       return NUMNODESP2xP4+NUMNODESP1;
 		case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
+		case LACrouzeixRaviartEnum: return NUMNODESP2b;
 		default:       _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -974,14 +975,15 @@
 int  PentaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
-		case P1P1Enum:           return P1Enum;
-		case P1P1GLSEnum:        return P1Enum;
-		case MINIcondensedEnum:  return P1bubbleEnum;
-		case MINIEnum:           return P1bubbleEnum;
-		case TaylorHoodEnum:     return P2Enum;
-		case LATaylorHoodEnum:   return P2Enum;
-		case OneLayerP4zEnum:    return P2xP4Enum;
-		case CrouzeixRaviartEnum:return P2bubbleEnum;
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1bubbleEnum;
+		case MINIEnum:              return P1bubbleEnum;
+		case TaylorHoodEnum:        return P2Enum;
+		case LATaylorHoodEnum:      return P2Enum;
+		case OneLayerP4zEnum:       return P2xP4Enum;
+		case CrouzeixRaviartEnum:   return P2bubbleEnum;
+		case LACrouzeixRaviartEnum: return P2bubbleEnum;
 		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -991,14 +993,15 @@
 int  PentaRef::PressureInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
-		case P1P1Enum:           return P1Enum;
-		case P1P1GLSEnum:        return P1Enum;
-		case MINIcondensedEnum:  return P1Enum;
-		case MINIEnum:           return P1Enum;
-		case TaylorHoodEnum:     return P1Enum;
-		case LATaylorHoodEnum:   return NoneEnum;
-		case OneLayerP4zEnum:    return P1Enum;
-		case CrouzeixRaviartEnum:return P1DGEnum;
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1Enum;
+		case MINIEnum:              return P1Enum;
+		case TaylorHoodEnum:        return P1Enum;
+		case LATaylorHoodEnum:      return NoneEnum;
+		case OneLayerP4zEnum:       return P1Enum;
+		case CrouzeixRaviartEnum:   return P1DGEnum;
+		case LACrouzeixRaviartEnum: return NoneEnum;
 		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18474)
@@ -2416,6 +2416,29 @@
 			penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+5;
 			penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+6;
 			break;
+		case LACrouzeixRaviartEnum:
+			numnodes         = 19;
+			penta_node_ids   = xNew<int>(numnodes);
+			penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
+			penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
+			penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
+			penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
+			penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
+			penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
+			penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
+			penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
+			penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
+			penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
+			penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
+			penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
+			penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
+			penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
+			penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
+			penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
+			penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
+			penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
+			penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
+			break;
 		default:
 			_error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
 	}
Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18473)
+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 18474)
@@ -407,23 +407,24 @@
 int  TriaRef::NumberofNodes(int finiteelement){/*{{{*/
 
 	switch(finiteelement){
-		case NoneEnum:              return 0;
-		case P0Enum:                return NUMNODESP0;
-		case P1Enum:                return NUMNODESP1;
-		case P1DGEnum:              return NUMNODESP1;
-		case P1bubbleEnum:          return NUMNODESP1b;
-		case P1bubblecondensedEnum: return NUMNODESP1b;
-		case P2Enum:                return NUMNODESP2;
-		case P2bubbleEnum:          return NUMNODESP2b;
-		case P2bubblecondensedEnum: return NUMNODESP2b;
-		case P1P1Enum:              return NUMNODESP1*2;
-		case P1P1GLSEnum:           return NUMNODESP1*2;
-		case MINIcondensedEnum:     return NUMNODESP1b+NUMNODESP1;
-		case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
-		case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
-		case LATaylorHoodEnum:      return NUMNODESP2;
-		case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
-		case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
+		case NoneEnum:                return 0;
+		case P0Enum:                  return NUMNODESP0;
+		case P1Enum:                  return NUMNODESP1;
+		case P1DGEnum:                return NUMNODESP1;
+		case P1bubbleEnum:            return NUMNODESP1b;
+		case P1bubblecondensedEnum:   return NUMNODESP1b;
+		case P2Enum:                  return NUMNODESP2;
+		case P2bubbleEnum:            return NUMNODESP2b;
+		case P2bubblecondensedEnum:   return NUMNODESP2b;
+		case P1P1Enum:                return NUMNODESP1*2;
+		case P1P1GLSEnum:             return NUMNODESP1*2;
+		case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
+		case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
+		case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
+		case LATaylorHoodEnum:        return NUMNODESP2;
+		case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
+		case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
+		case LACrouzeixRaviartEnum:   return NUMNODESP2b;
 		default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
 	}
 
@@ -433,14 +434,15 @@
 int  TriaRef::VelocityInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
-		case P1P1Enum:           return P1Enum;
-		case P1P1GLSEnum:        return P1Enum;
-		case MINIcondensedEnum:  return P1bubbleEnum;
-		case MINIEnum:           return P1bubbleEnum;
-		case TaylorHoodEnum:     return P2Enum;
-		case LATaylorHoodEnum:   return P2Enum;
-		case XTaylorHoodEnum:    return P2Enum;
-		case CrouzeixRaviartEnum:return P2bubbleEnum;
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1bubbleEnum;
+		case MINIEnum:              return P1bubbleEnum;
+		case TaylorHoodEnum:        return P2Enum;
+		case LATaylorHoodEnum:      return P2Enum;
+		case XTaylorHoodEnum:       return P2Enum;
+		case CrouzeixRaviartEnum:   return P2bubbleEnum;
+		case LACrouzeixRaviartEnum: return P2bubbleEnum;
 		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
@@ -450,14 +452,15 @@
 int  TriaRef::PressureInterpolation(int fe_stokes){/*{{{*/
 
 	switch(fe_stokes){
-		case P1P1Enum:            return P1Enum;
-		case P1P1GLSEnum:         return P1Enum;
-		case MINIcondensedEnum:   return P1Enum;
-		case MINIEnum:            return P1Enum;
-		case TaylorHoodEnum:      return P1Enum;
-		case LATaylorHoodEnum:    return NoneEnum;
-		case XTaylorHoodEnum:     return P1Enum;
-		case CrouzeixRaviartEnum: return P1DGEnum;
+		case P1P1Enum:              return P1Enum;
+		case P1P1GLSEnum:           return P1Enum;
+		case MINIcondensedEnum:     return P1Enum;
+		case MINIEnum:              return P1Enum;
+		case TaylorHoodEnum:        return P1Enum;
+		case LATaylorHoodEnum:      return NoneEnum;
+		case XTaylorHoodEnum:       return P1Enum;
+		case CrouzeixRaviartEnum:   return P1DGEnum;
+		case LACrouzeixRaviartEnum: return NoneEnum;
 		default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
 	}
 
