Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17490)
@@ -738,4 +738,27 @@
 			tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3];
 			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			numnodes         = 5;
+			tetra_node_ids   = xNew<int>(numnodes);
+			tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0];
+			tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1];
+			tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2];
+			tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3];
+			tetra_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+index+1;
+			break;
+		case P2Enum:
+			numnodes        = 10;
+			tetra_node_ids   = xNew<int>(numnodes);
+			tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0];
+			tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1];
+			tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2];
+			tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3];
+			tetra_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+0]+1;
+			tetra_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+1]+1;
+			tetra_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+2]+1;
+			tetra_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+3]+1;
+			tetra_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+4]+1;
+			tetra_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+5]+1;
+			break;
 		default:
 			_error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.cpp	(revision 17490)
@@ -18,4 +18,5 @@
 #define NUMNODESP0  1
 #define NUMNODESP1  4
+#define NUMNODESP1b 5
 #define NUMNODESP2  10
 
@@ -72,4 +73,13 @@
 			basis[2]=gauss->coord3;
 			basis[3]=gauss->coord4;
+			return;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			/*Corner nodes*/
+			basis[0]=gauss->coord1;
+			basis[1]=gauss->coord2;
+			basis[2]=gauss->coord3;
+			basis[3]=gauss->coord4;
+			/*bubble*/
+			basis[4]=256.*gauss->coord1*gauss->coord2*gauss->coord3*gauss->coord4;
 			return;
 		case P2Enum:
@@ -171,4 +181,25 @@
 			dbasis[NUMNODESP1*2+3] = 1.;
 			return;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			dbasis[NUMNODESP1b*0+0] = -1.;
+			dbasis[NUMNODESP1b*1+0] = -1.;
+			dbasis[NUMNODESP1b*2+0] = -1.;
+
+			dbasis[NUMNODESP1b*0+1] = 1.;
+			dbasis[NUMNODESP1b*1+1] = 0.;
+			dbasis[NUMNODESP1b*2+1] = 0.;
+
+			dbasis[NUMNODESP1b*0+2] = 0.;
+			dbasis[NUMNODESP1b*1+2] = 1.;
+			dbasis[NUMNODESP1b*2+2] = 0.;
+
+			dbasis[NUMNODESP1b*0+3] = 0.;
+			dbasis[NUMNODESP1b*1+3] = 0.;
+			dbasis[NUMNODESP1b*2+3] = 1.;
+
+			dbasis[NUMNODESP1b*0+4] = 256.*(-gauss->coord2*gauss->coord3*gauss->coord4+gauss->coord1*gauss->coord3*gauss->coord4);
+			dbasis[NUMNODESP1b*1+4] = 256.*(-gauss->coord2*gauss->coord3*gauss->coord4+gauss->coord1*gauss->coord2*gauss->coord4);
+			dbasis[NUMNODESP1b*2+4] = 256.*(-gauss->coord2*gauss->coord3*gauss->coord4+gauss->coord1*gauss->coord2*gauss->coord3);
+			return;
 		case P2Enum:
 			dbasis[NUMNODESP2*0+0] = -4.*gauss->coord1+1.;
@@ -372,8 +403,10 @@
 
 	switch(finiteelement){
-		case P0Enum:   return NUMNODESP0;
-		case P1Enum:   return NUMNODESP1;
-		case P1DGEnum: return NUMNODESP1;
-		case P2Enum:   return NUMNODESP2;
+		case P0Enum:                return NUMNODESP0;
+		case P1Enum:                return NUMNODESP1;
+		case P1DGEnum:              return NUMNODESP1;
+		case P1bubbleEnum:          return NUMNODESP1b;
+		case P1bubblecondensedEnum: return NUMNODESP1b;
+		case P2Enum:                return NUMNODESP2;
 		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
 	}
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp	(revision 17490)
@@ -184,4 +184,30 @@
 			}
 			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			switch(iv){
+				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+				case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
+				case 4: coord1=1./4.; coord2=1./4.; coord3=1./4.; coord4=1./4.; break;
+				default: _error_("node index should be in [0 4]");
+			}
+			break;
+		case P2Enum:
+			switch(iv){
+				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+				case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
+
+				case 4: coord1=0.; coord2=.5; coord3=.5; coord4=0.; break;
+				case 5: coord1=.5; coord2=0.; coord3=.5; coord4=0.; break;
+				case 6: coord1=.5; coord2=.5; coord3=0.; coord4=0.; break;
+				case 7: coord1=.5; coord2=0.; coord3=0.; coord4=.5; break;
+				case 8: coord1=0.; coord2=.5; coord3=0.; coord4=.5; break;
+				case 9: coord1=0.; coord2=0.; coord3=.5; coord4=.5; break;
+				default: _error_("node index should be in [0 9]");
+			}
+			break;
 		default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
 	}
Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp	(revision 17490)
@@ -54,4 +54,5 @@
 				case Mesh2DhorizontalEnum: elementnbv = 3; break;
 				case Mesh2DverticalEnum:   elementnbv = 3; break;
+				case Mesh3DtetrasEnum:     elementnbv = 4; break;
 				case Mesh3DEnum:           elementnbv = 6; break;
 				default: _error_("mesh type not supported yet");
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp	(revision 17490)
@@ -40,4 +40,16 @@
 		elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
 		elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
+	}
+	else if(iomodel->meshtype==Mesh3DtetrasEnum){
+		elementnbv = 4;
+		elementnbe = 6;
+		elementedges         = xNew<int>(elementnbe*2);
+		elementedges_markers = xNew<int>(elementnbe);
+		elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
+		elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
+		elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
+		elementedges[2*3+0] = 0; elementedges[2*3+1] = 3; elementedges_markers[3] = 1;
+		elementedges[2*4+0] = 1; elementedges[2*4+1] = 3; elementedges_markers[4] = 1;
+		elementedges[2*5+0] = 2; elementedges[2*5+1] = 3; elementedges_markers[5] = 1;
 	}
 	else if(iomodel->meshtype==Mesh3DEnum){
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp	(revision 17490)
@@ -24,4 +24,7 @@
 		elementnbe = 3;
 	}
+	else if(iomodel->meshtype==Mesh3DtetrasEnum){
+		elementnbe = 6;
+	}
 	else if(iomodel->meshtype==Mesh3DEnum){
 		elementnbe = 9;
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17489)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17490)
@@ -64,6 +64,6 @@
 
 	/*Display size*/
+	Kff->GetSize(&M,&N);
 	if(VerboseModule()){
-		Kff->GetSize(&M,&N);
 		_printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
 	}
@@ -120,4 +120,6 @@
 	//Kff->AllocationInfo();
 	//Kfs->AllocationInfo();
+	//IssmDouble* k = Kff->ToSerial();
+	//printarray(k,N,N);
 
 	/*cleanu up and assign output pointers: */