Changeset 1236

Timestamp:

07/06/09 11:12:27 (15 years ago)

Author:

Mathieu Morlighem

Message:

use GetNodalFunctionsCoeff and GetArea now

Location:

issm/trunk/src/m/classes/public

Files:

• area.m (deleted)
• averaging.m (modified) (2 diffs)
• mechanicalproperties.m (modified) (2 diffs)
• mesh/meshadaptation.m (modified) (2 diffs)
• mesh/meshexprefine.m (modified) (1 diff)
• removeholes.m (modified) (2 diffs)
• shear2d.m (modified) (1 diff)
• slope.m (modified) (1 diff)
• thicknessevolution.m (modified) (1 diff)

issm/trunk/src/m/classes/public/averaging.m

@@ -22 +22 @@
 %initialization
 weights=zeros(md.numberofgrids,1);
-areas=zeros(md.numberofelements,1);
 data=data(:);
 
@@ -34 +33 @@
 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

issm/trunk/src/m/classes/public/mechanicalproperties.m

@@ -28 +28 @@
 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 +34 @@
 
 %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

issm/trunk/src/m/classes/public/mesh/meshadaptation.m

@@ -22 +22 @@
 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 +35 @@
 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

issm/trunk/src/m/classes/public/mesh/meshexprefine.m

@@ -38 +38 @@
 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

issm/trunk/src/m/classes/public/removeholes.m

@@ -24 +24 @@
 end
 
-
 %Ok, retrieve and write domain outline without holes, to disk.
 domainoutline_string=md.domainoutline;
@@ -33 +32 @@
 %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);

issm/trunk/src/m/classes/public/shear2d.m

@@ -7 +7 @@
 %      [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];

issm/trunk/src/m/classes/public/slope.m

@@ -18 +18 @@
 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];

issm/trunk/src/m/classes/public/thicknessevolution.m

@@ -13 +13 @@
 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);