Index: sm/trunk/src/m/classes/public/area.m
===================================================================
--- /issm/trunk/src/m/classes/public/area.m	(revision 1235)
+++ 	(revision )
@@ -1,22 +1,0 @@
-function areas=area(md)
-%AREA - compute areas for every element
-%
-%   Usage:
-%      areas=area(md)
-
-areas=zeros(md.numberofelements,1);
-
-%load some variables (it is much faster if the variab;es are loaded from md once for all)
-index=md.elements;
-x=md.x; y=md.y; z=md.z;
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
-
-%compute the volume of each element
-if strcmpi(md.type,'3d')
-	%V=area(triangle)*1/3(z1+z2+z3)
-	thickness=mean(z(index(:,4:6)))-mean(z(index(:,1:3)));
-	areas=abs(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
-else
-	%compute the surface of the triangle
-	areas=abs(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
-end
Index: /issm/trunk/src/m/classes/public/averaging.m
===================================================================
--- /issm/trunk/src/m/classes/public/averaging.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/averaging.m	(revision 1236)
@@ -22,5 +22,4 @@
 %initialization
 weights=zeros(md.numberofgrids,1);
-areas=zeros(md.numberofelements,1);
 data=data(:);
 
@@ -34,12 +33,11 @@
 if strcmpi(md.type,'3d')
 	rep=6;
+	areas=Getareas(index,md.x,md.y,md.z);
 else
 	rep=3;
+	areas=Getareas(index,md.x,md.y);
 end
 summation=1/rep*ones(rep,1);
 linesize=rep*numberofelements;
-
-%compute the volume of each element
-areas=area(md);
 
 %update weights that holds the volume of all the element holding the grid i
Index: /issm/trunk/src/m/classes/public/mechanicalproperties.m
===================================================================
--- /issm/trunk/src/m/classes/public/mechanicalproperties.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/mechanicalproperties.m	(revision 1236)
@@ -28,7 +28,4 @@
 index=md.elements;
 summation=[1;1;1];
-alpha=zeros(numberofelements,3);
-beta=zeros(numberofelements,3);
-%gamma=zeros(numberofelements,3);
 directionsstress=zeros(numberofelements,4);
 directionsstrain=zeros(numberofelements,4);
@@ -37,11 +34,5 @@
 
 %compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
-x1=md.x(index(:,1)); x2=md.x(index(:,2)); x3=md.x(index(:,3));
-y1=md.y(index(:,1)); y2=md.y(index(:,2)); y3=md.y(index(:,3));
-invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
-alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
-beta =[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
-%gamma=[invdet.*(x2.*y3-x3.*y2) invdet.*(y1.*x3-y3.*x1) invdet.*(x1.*y2-x2.*y1)];
-clear invdet x1 x2 x3 y1 y2 y3
+[alpha beta]=GetNodalFunctionsCoeff(index,md.x,md.y);
 
 %compute shear
Index: /issm/trunk/src/m/classes/public/mesh/meshadaptation.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshadaptation.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/mesh/meshadaptation.m	(revision 1236)
@@ -22,16 +22,11 @@
 disp(sprintf('      metric computation') )
 
-%load some variables (it is much faster if the variab;es are loaded from md once for all) 
+%initialization
+index=md.elements;
+numberofgrids=md.numberofgrids;
 numberofelements=md.numberofelements;
-numberofgrids=md.numberofgrids;
-index=md.elements;
-x=md.x; y=md.y; z=md.z;
-
-%initialization
-alpha=zeros(md.numberofelements,3);
-beta=zeros(md.numberofelements,3);
-gradx=zeros(md.numberofgrids,1);
-grady=zeros(md.numberofgrids,1);
-metric=zeros(md.numberofelements,1);
+gradx=zeros(numberofgrids,1);
+grady=zeros(numberofgrids,1);
+metric=zeros(numberofelements,1);
 
 %build some usefull variables
@@ -40,17 +35,12 @@
 summation=1/3*ones(3,1);
 linesize=3*numberofelements;
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
 
-%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
-alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
-beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
+%get areas and  nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
+[alpha beta]=GetNodalFunctionsCoeff(index,md.x,md.y);
+areas=Getareas(index,md.x,md.y);
 
 %Compute gradient for each element
 grad_elx=sum(field(index).*alpha,2); 
 grad_ely=sum(field(index).*beta,2);
-
-%compute the volume of each element
-areas=area(md);
 
 %update weights that holds the volume of all the element holding the grid i
Index: /issm/trunk/src/m/classes/public/mesh/meshexprefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshexprefine.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/mesh/meshexprefine.m	(revision 1236)
@@ -38,7 +38,7 @@
 end
 
-%Read domainname file into a matlab array (x,y):
+%Read domainame file into a matlab array (x,y):
 refinearea=ContourToMesh(md.elements,md.x,md.y,expread(domainname,1),'element',1);
-aires=area(md);
+aires=Getareas(md.elements,md.x,md.y);
 
 %flags areas within the domain
Index: /issm/trunk/src/m/classes/public/removeholes.m
===================================================================
--- /issm/trunk/src/m/classes/public/removeholes.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/removeholes.m	(revision 1236)
@@ -24,5 +24,4 @@
 end
 
-
 %Ok, retrieve and write domain outline without holes, to disk.
 domainoutline_string=md.domainoutline;
@@ -33,5 +32,5 @@
 %Now create new model with mesh based on DomainOutlineTemp: 
 %get average resolution
-resolution=mean(sqrt(2*area(md)));
+resolution=mean(sqrt(2*area(md.elements,md.x,md.y)));
 md2=model;
 md2=mesh(md2,'DomainOutlineTemp.exp',resolution);
Index: /issm/trunk/src/m/classes/public/shear2d.m
===================================================================
--- /issm/trunk/src/m/classes/public/shear2d.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/shear2d.m	(revision 1236)
@@ -7,14 +7,5 @@
 %      [sx,sy,sxy,s]=shear2d(md);
 
-alpha=zeros(md.numberofelements,3); beta=zeros(md.numberofelements,3); 
-gamma=zeros(md.numberofelements,3); area=zeros(md.numberofelements,1);
-
-for n=1:md.numberofelements
-	X=inv([md.x(md.elements(n,:)) md.y(md.elements(n,:)) ones(3,1)]);
-	alpha(n,:)=X(1,:);
-	beta(n,:)=X(2,:);
-end
-clear X;
-
+[alpha beta]=GetNodalFunctionsCoeff(md.elements,md.x,md.y); 
 
 summation=[1;1;1];
Index: /issm/trunk/src/m/classes/public/slope.m
===================================================================
--- /issm/trunk/src/m/classes/public/slope.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/slope.m	(revision 1236)
@@ -18,16 +18,6 @@
 end
 
-%initialization
-alpha=zeros(numberofelements,3);
-beta=zeros(numberofelements,3);
-
-%build some usefull variables
-summation=1/3*ones(3,1);
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
-
 %compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
-alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
-beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
+[alpha beta]=GetNodalFunctionsCoeff(index,x,y);
 
 summation=[1;1;1];
Index: /issm/trunk/src/m/classes/public/thicknessevolution.m
===================================================================
--- /issm/trunk/src/m/classes/public/thicknessevolution.m	(revision 1235)
+++ /issm/trunk/src/m/classes/public/thicknessevolution.m	(revision 1236)
@@ -13,30 +13,14 @@
 end
 
-%load some variables (it is much faster if the variab;es are loaded from md once for all) 
-numberofelements=md.numberofelements;
+%load some variables 
 H=md.thickness;
 vx=md.vx;
 vy=md.vy;
-index=md.elements;
-x=md.x; y=md.y;
-
-%initialization
-alpha=zeros(md.numberofelements,3);
-beta=zeros(md.numberofelements,3);
-gradx=zeros(md.numberofgrids,1);
-grady=zeros(md.numberofgrids,1);
-
-%build some usefull variables
-line=index(:);
-summation=1/3*ones(3,1);
-linesize=3*numberofelements;
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
 
 %compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
-alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
-beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
+[alpha beta]=GetNodalFunctionsCoeff(md.elements,md.x,md.y); 
 
 %compute dhdt=div(Hu)
+summation=1/3*ones(3,1);
 dhdt=(vx(index)*summation).*sum( H(index).*alpha,2) + (vy(index)*summation).*sum(H(index).*beta,2) ...
 	+ ( H(index)*summation).*sum(vx(index).*alpha,2) + ( H(index)*summation).*sum(vy(index).*beta,2);