Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20017)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20018)
@@ -3648,5 +3648,8 @@
 				else{
 					for(int n=0;n<nl;n++){
-						Pn=legendre(Pn1,Pn2,reCast<int,IssmDouble>(cosalpha),n); Pn1=Pn2; Pn2=Pn;
+						if(n==0)Pn=1;
+						else if(n==1)Pn=cosalpha;
+						else Pn= ( (2*n-1)*cosalpha*Pn1 - (n-1)*Pn2 ) /n;
+						Pn2=Pn1; Pn1=Pn;
 						G_elastic += deltalove[n]*Pn;
 					}
@@ -3809,9 +3812,9 @@
 
 		/*initialize G_elastic:*/
-		G_elastic=xNew<IssmDouble>(gsize);
+		G_elastic=xNewZeroInit<IssmDouble>(gsize);
 	}
 	if(compute_G_rigid){
 		/*Initialize G_rigid and G_elastic:*/
-		G_rigid=xNew<IssmDouble>(gsize);
+		G_rigid=xNewZeroInit<IssmDouble>(gsize);
 	}
 
@@ -3846,5 +3849,8 @@
 				else{
 					for(int n=0;n<nl;n++){
-						Pn=legendre(Pn1,Pn2,reCast<int,IssmDouble>(cosalpha),n); Pn1=Pn2; Pn2=Pn;
+						if(n==0)Pn=1;
+						else if(n==1)Pn=cosalpha;
+						else Pn= ( (2*n-1)*cosalpha*Pn1 - (n-1)*Pn2 ) /n;
+						Pn2=Pn1; Pn1=Pn;
 						G_elastic[i] += deltalove[n]*Pn;
 					}