Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4989)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4990)
@@ -703,5 +703,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		J=tria->CostFunction(process_units);
-		delete tria;
+		delete tria->matice; delete tria;
 
 		/*delete B average*/
@@ -717,5 +717,5 @@
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 		J=tria->CostFunction(process_units);
-		delete tria;
+		delete tria->matice; delete tria;
 
 		/*delete B average*/
@@ -1577,5 +1577,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		J=tria->Misfit(process_units);
-		delete tria;
+		delete tria->matice; delete tria;
 		return J;
 	}
@@ -1584,5 +1584,5 @@
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 		J=tria->Misfit(process_units);
-		delete tria;
+		delete tria->matice; delete tria;
 		return J;
 	}
@@ -1712,5 +1712,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		S=tria->SurfaceArea();
-		delete tria;
+		delete tria->matice; delete tria;
 		return S;
 	}
@@ -1719,5 +1719,5 @@
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 		S=tria->SurfaceArea();
-		delete tria;
+		delete tria->matice; delete tria;
 		return S;
 	}
@@ -1982,5 +1982,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreateKMatrix(Kgg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -2020,5 +2020,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreateKMatrix(Kgg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -2151,5 +2151,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreateKMatrix(Kgg);
-		delete tria;
+		delete tria->matice; delete tria;
 
 		/*Delete B averaged*/
@@ -2251,5 +2251,5 @@
 			tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 			tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
-			delete tria;
+			delete tria->matice; delete tria;
 		}
 
@@ -2689,5 +2689,5 @@
 		tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
 		tria->CreateKMatrixDiagnosticSurfaceVert(Kgg);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 
@@ -2780,5 +2780,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreateKMatrixMelting(Kgg);
-		delete tria;
+		delete tria->matice; delete tria;
 		return;
 	}
@@ -2812,5 +2812,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreateKMatrix(Kgg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -2846,5 +2846,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreateKMatrix(Kgg);
-	delete tria;
+	delete tria->matice; delete tria;
 	return;
 
@@ -3092,5 +3092,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreateKMatrixThermal(Kgg);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 }
@@ -3123,5 +3123,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreatePVector(pg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -3159,5 +3159,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreatePVector(pg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -3245,5 +3245,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreatePVector(pg);
-		delete tria;
+		delete tria->matice; delete tria;
 		return;
 	}
@@ -3355,5 +3355,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		tria->CreatePVectorAdjointHoriz(p_g);
-		delete tria;
+		delete tria->matice; delete tria;
 		return;
 	}
@@ -3362,5 +3362,5 @@
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 		tria->CreatePVectorAdjointHoriz(p_g);
-		delete tria;
+		delete tria->matice; delete tria;
 		return;
 	}
@@ -3793,5 +3793,5 @@
 	tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 	tria->CreatePVectorAdjointStokes(p_g);
-	delete tria;
+	delete tria->matice; delete tria;
 	return;
 }
@@ -3862,5 +3862,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
 		tria->CreatePVectorDiagnosticBaseVert(pg);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 
@@ -3961,5 +3961,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreatePVector(pg);
-	delete tria;
+	delete tria->matice; delete tria;
 
 	/*Delete Vx and Vy averaged*/
@@ -3993,5 +3993,5 @@
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 	tria->CreatePVector(pg);
-	delete tria;
+	delete tria->matice; delete tria;
 	return;
 }
@@ -4166,5 +4166,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreatePVectorThermalShelf(pg);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 
@@ -4174,5 +4174,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreatePVectorThermalSheet(pg);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 
@@ -4702,5 +4702,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		tria->GradjB(gradient);
-		delete tria;
+		delete tria->matice; delete tria;
 
 		/*delete B average*/
@@ -4714,5 +4714,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
 		tria->GradjB(gradient);
-		delete tria;
+		delete tria->matice; delete tria;
 
 		/*delete B average*/
@@ -4758,5 +4758,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->GradjDrag(gradient);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 	else if (analysis_type==AdjointStokesAnalysisEnum){
@@ -4765,5 +4765,5 @@
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->GradjDragStokes(gradient);
-		delete tria;
+		delete tria->matice; delete tria;
 	}
 	else ISSMERROR("analysis %s not supported yet",EnumAsString(analysis_type));
@@ -5595,10 +5595,15 @@
 	tria->results=tria_results;
 	tria->parameters=tria_parameters;
-	tria->element_type=P1Enum; //Only P1 CG for now
+	tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
 	this->SpawnTriaHook(dynamic_cast<TriaHook*>(tria),&indices[0]);
+
+	/*Spawn matice*/
+	tria->matice=NULL;
+	tria->matice=(Matice*)this->matice->copy();
+	delete tria->matice->inputs;
+	tria->matice->inputs=(Inputs*)this->matice->inputs->SpawnTriaInputs(indices);
 
 	/*recover nodes, matice and matpar: */
 	tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
-	tria->matice=(Matice*)tria->hmatice->delivers();
 	tria->matpar=(Matpar*)tria->hmatpar->delivers();
 	
Index: /issm/trunk/src/c/objects/Elements/PentaHook.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaHook.cpp	(revision 4989)
+++ /issm/trunk/src/c/objects/Elements/PentaHook.cpp	(revision 4990)
@@ -92,5 +92,6 @@
 		}
 	}
-	triahook->hmatice=(Hook*)this->hmatice->copy();
+	// do not spawn hmatice. matice will be taken care of by Penta
+	triahook->hmatice=NULL;
 	triahook->hmatpar=(Hook*)this->hmatpar->copy();
 }
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 4989)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 4990)
@@ -871,9 +871,8 @@
 	double A1,A2,A3,z;
 
-	A1=gauss[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
-	A2=gauss[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
-	A3=gauss[2]; //third area coordinate value  In term of xi and eta: A3=y/SQRT3;
-	z=gauss[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
-
+	A1=gauss[0]; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
+	A2=gauss[1]; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
+	A3=gauss[2]; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value  In term of xi and eta: A3=y/SQRT3;
+	z=gauss[3];  ISSMASSERT(z>=-1 &&  z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
 
 	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
Index: /issm/trunk/src/c/objects/Materials/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 4989)
+++ /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 4990)
@@ -220,4 +220,7 @@
 	Matice* matice=NULL;
 
+	/*Initialize output*/
+	matice=new Matice();
+
 	/*copy fields: */
 	matice->mid=this->mid;