Changeset 51

Timestamp:

04/24/09 15:36:04 (16 years ago)

Author:

Mathieu Morlighem

Message:

removed acceleratedtriaelem

Location:

issm/trunk/src/m/solutions/ice

Files:

• DeviatoricStressCompute.m (modified) (2 diffs)
• GetNumberOfActiveElements.m (deleted)
• ModelProcessorDiagnosticBaseVert.m (modified) (1 diff)
• ModelProcessorDiagnosticHoriz.m (modified) (3 diffs)
• ModelProcessorDiagnosticStokes.m (modified) (1 diff)
• ModelProcessorDiagnosticVert.m (modified) (1 diff)
• ModelProcessorMelting.m (modified) (1 diff)
• ModelProcessorSlopeCompute.m (modified) (1 diff)
• ModelProcessorThermal.m (modified) (1 diff)
• PressureElemCompute.m (modified) (2 diffs)
• StrainRateCompute.m (modified) (2 diffs)
• StressBedCompute.m (modified) (2 diffs)
• StressCompute.m (modified) (2 diffs)
• StressSurfaceCompute.m (modified) (2 diffs)
• SystemMatrices.m (modified) (3 diffs)
• ViscousHeatingCompute.m (modified) (1 diff)

issm/trunk/src/m/solutions/ice/DeviatoricStressCompute.m

@@ -14 +14 @@
 loads=m.loads;
 gridset=m.gridset;
-
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
+numberofelements=length(elements);
 
 if strcmpi(type,'2d')
         %initialize vectors
         deviatoricstress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[]);
-        deviatoricstress1=zeros((n2-n1)+1,3);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); 
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1);
+        deviatoricstress1=zeros(numberofelements,3);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); 
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1);
 
         %Go through all elements and call the deviatoricstress routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         deviatoricstressvector=DeviatoricStress(elements(n).element,grids,materials,inputs)';
@@ -58 +56 @@
         %initialize vectors
         deviatoricstress=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[]);
-        deviatoricstress1=zeros((n2-n1)+1,6);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
-        A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
+        deviatoricstress1=zeros(numberofelements,6);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
+        A3=zeros(numberofelements,1); Vx3=zeros(numberofelements,1); Vy3=zeros(numberofelements,1); Vz3=zeros(numberofelements,1);
 
         %Go through all elements and call the deviatoricstress routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         deviatoricstressvector=DeviatoricStress(elements(n).element,grids,materials,inputs)';

issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m

@@ -44 +44 @@
 elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos));
 elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos));
-elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
 elements(pos)=SetStructureField(elements(pos),'element','collapse',1);
 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);

issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m

@@ -52 +52 @@
                 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
                 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
-                elements(pos)=SetStructureField(elements(pos),'element','acceleration',md.acceleration);
                 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
 
@@ -87 +86 @@
                 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
                 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
-                elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
                 elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
                 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
@@ -125 +123 @@
         mygrids(md.elements(el3pos,:))=1;
         mygrids(md.elements(el6pos,:))=1;
-end
-
-%Deal with acceleration MacAyeal's element. 
-if md.acceleration & strcmpi(md.type,'2d'), 
-        %build accelerated element.
-        elements(end+1).element=acceleratedtriaelem;
-        elements(end).element.type='acceleratedtriaelem';
-        elements(end).element.x=md.x;
-        elements(end).element.y=md.y;
-        elements(end).element.index=md.elements;
-        elements(end).element.nods=md.numberofgrids;
-        elements(end).element.nel=md.numberofelements;
-        elements(end).element.thickness=md.thickness;
-        elements(end).element.surface=md.surface;
-        elements(end).element.bed=md.bed;
-        elements(end).element.thickness_el=(md.thickness(md.elements))*[1;1;1]/3;
-        elements(end).element.B_bar=(md.B(md.elements))*[1;1;1]/3;
-        elements(end).element.glen_coeff=md.n;
-        elements(end).element.friction_type=md.drag_type;
-        elements(end).element.drag=md.drag;
-        elements(end).element.p=md.p;
-        elements(end).element.q=md.q;
-        elements(end).element.meanvel=md.meanvel;
-        elements(end).element.epsvel=md.epsvel;
-        if isempty(md.segmentonneumann_diag),
-                elements(end).element.index_icefront=[];        
-        else
-                elements(end).element.index_icefront=md.segmentonneumann_diag(:,1:2);
-        end
-        elements(end).element.gridoniceshelf=md.gridoniceshelf;
-        elements(end).element.elementonicesheet=md.elementonicesheet;
-        elements(end).element.matid=1; %point to first matid element.
 end
 

issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticStokes.m

@@ -49 +49 @@
 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','meanvel',md.meanvel);
 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','epsvel',md.epsvel);
-elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','collapse',0);
 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','matid',stokesnewnumber);

issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m

@@ -46 +46 @@
 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
-elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
 elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);

issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m

@@ -44 +44 @@
 elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos));
 elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos));
-elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
 elements(pos)=SetStructureField(elements(pos),'element','collapse',1);
 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);

issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m

@@ -51 +51 @@
         elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos));
         elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos));
-        elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
         elements(pos)=SetStructureField(elements(pos),'element','collapse',1);
         elements(pos)=SetStructureField(elements(pos),'element','matid',pos);

issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m

@@ -46 +46 @@
 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
-elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);%acceleration 0 since no acceleration is posssible
 elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);

issm/trunk/src/m/solutions/ice/PressureElemCompute.m

@@ -18 +18 @@
 gridset=m.gridset;
 
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
-
 %initialization
 pressure=zeros((n2-n1)+1,1);
@@ -31 +28 @@
 
 %Go through elements and build pressure defined as P=-1/3*tr(stress)
-for n=n1:n2,
-        if ~isempty(elements(n).element),
-                stress_tensor=Stress(elements(n).element,grids,materials,inputs);       
-                trace=stress_tensor(1)+stress_tensor(2)+stress_tensor(3);
-                pressure(n)=-1/3*trace;
-        end
+for n=1:length(elements)
+        stress_tensor=Stress(elements(n).element,grids,materials,inputs);       
+        trace=stress_tensor(1)+stress_tensor(2)+stress_tensor(3);
+        pressure(n)=-1/3*trace;
 end

issm/trunk/src/m/solutions/ice/StrainRateCompute.m

@@ -17 +17 @@
 loads=m.loads;
 gridset=m.gridset;
-
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
+numberofelements=length(elements);
 
 if strcmpi(type,'2d')
         %initialize vectors
         strainrate=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[]);
-        strainrate1=zeros((n2-n1)+1,3);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1);
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1);
+        strainrate1=zeros(numberofelements,3);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1);
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1);
 
         %Go through all elements and call the strainrate routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         strainratevector=StrainRate(elements(n).element,grids,materials,inputs)';
@@ -63 +61 @@
         %initialize vectors
         strainrate=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[]);
-        strainrate1=zeros((n2-n1)+1,6);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
-        A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
+        strainrate1=zeros(numberofelements,6);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
+        A3=zeros(numberofelements,1); Vx3=zeros(numberofelements,1); Vy3=zeros(numberofelements,1); Vz3=zeros(numberofelements,1);
 
         %Go through all elements and call the strainrate routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         strainratevector=StrainRate(elements(n).element,grids,materials,inputs)';

issm/trunk/src/m/solutions/ice/StressBedCompute.m

@@ -17 +17 @@
 loads=m.loads;
 gridset=m.gridset;
-
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
+numberofelements=length(elements);
 
 %initialization
-stress_bed=zeros((n2-n1)+1,1);
+stress_bed=zeros(numberofelements,1);
 
 if strcmpi(type,'2d')
@@ -31 +29 @@
 %initialize vectors
 stress_bed=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'stress_n','stress_nn','normal_x',[],'normal_y',[],'normal_z',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[],'effectivevalue',[]);
-stress_bed1=zeros((n2-n1)+1,6);
-normal1=zeros((n2-n1)+1,3);
-stress_n1=zeros((n2-n1)+1,3);
-stress_nn1=zeros((n2-n1)+1,1);
-A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
-A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
+stress_bed1=zeros(numberofelements,6);
+normal1=zeros(numberofelements,3);
+stress_n1=zeros(numberofelements,3);
+stress_nn1=zeros(numberofelements,1);
+A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
+A3=zeros(numberofelements,1); Vx3=zeros(numberofelements,1); Vy3=zeros(numberofelements,1); Vz3=zeros(numberofelements,1);
 
 %Go through all elements and call the stress_bed routine, then compute eigen values and vector
-for n=n1:n2,
+for n=1:length(elements),
         if ~isempty(elements(n).element),
                 [stress_vector,normal]=StressBed(elements(n).element,grids,materials,inputs);

issm/trunk/src/m/solutions/ice/StressCompute.m

@@ -17 +17 @@
 loads=m.loads;
 gridset=m.gridset;
+numberofelements=length(elements);
 
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
 if strcmpi(type,'2d')
 
         %initialize vectors
         stress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[]);
-        stress1=zeros((n2-n1)+1,3);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
+        stress1=zeros(numberofelements,3);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
 
         %Go through all elements and call the stress routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         stressvector=Stress(elements(n).element,grids,materials,inputs)';
@@ -64 +63 @@
         %initialize vectors
         stress=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[]);
-        stress1=zeros((n2-n1)+1,6);
-        A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-        A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
-        A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
+        stress1=zeros(numberofelements,6);
+        A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+        A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
+        A3=zeros(numberofelements,1); Vx3=zeros(numberofelements,1); Vy3=zeros(numberofelements,1); Vz3=zeros(numberofelements,1);
 
         %Go through all elements and call the stress routine, then compute eigen values and vector
-        for n=n1:n2,
+        for n=1:length(elements),
                 if ~isempty(elements(n).element),
                         stressvector=Stress(elements(n).element,grids,materials,inputs)';

issm/trunk/src/m/solutions/ice/StressSurfaceCompute.m

@@ -17 +17 @@
 loads=m.loads;
 gridset=m.gridset;
-
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
+numberofelements=length(elements);
 
 %initialization
-stress_surface=zeros((n2-n1)+1,1);
+stress_surface=zeros(numberofelements,1);
 
 if strcmpi(type,'2d')
@@ -31 +29 @@
 %initialize vectors
 stress_surface=struct('xx',[],'yy',[],'zz',[],'xy',[],'xz',[],'yz',[],'stress_n','stress_nn','normal_x',[],'normal_y',[],'normal_z',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'principalvalue3',[],'principalaxis3',[],'effectivevalue',[]);
-stress_surface1=zeros((n2-n1)+1,6);
-normal1=zeros((n2-n1)+1,3);
-stress_n1=zeros((n2-n1)+1,3);
-stress_nn1=zeros((n2-n1)+1,1);
-A1=zeros((n2-n1)+1,1); Vx1=zeros((n2-n1)+1,1); Vy1=zeros((n2-n1)+1,1); Vz1=zeros((n2-n1)+1,1);
-A2=zeros((n2-n1)+1,1); Vx2=zeros((n2-n1)+1,1); Vy2=zeros((n2-n1)+1,1); Vz2=zeros((n2-n1)+1,1);
-A3=zeros((n2-n1)+1,1); Vx3=zeros((n2-n1)+1,1); Vy3=zeros((n2-n1)+1,1); Vz3=zeros((n2-n1)+1,1);
+stress_surface1=zeros(numberofelements,6);
+normal1=zeros(numberofelements,3);
+stress_n1=zeros(numberofelements,3);
+stress_nn1=zeros(numberofelements,1);
+A1=zeros(numberofelements,1); Vx1=zeros(numberofelements,1); Vy1=zeros(numberofelements,1); Vz1=zeros(numberofelements,1);
+A2=zeros(numberofelements,1); Vx2=zeros(numberofelements,1); Vy2=zeros(numberofelements,1); Vz2=zeros(numberofelements,1);
+A3=zeros(numberofelements,1); Vx3=zeros(numberofelements,1); Vy3=zeros(numberofelements,1); Vz3=zeros(numberofelements,1);
 
 %Go through all elements and call the stress_surface routine, then compute eigen values and vector
-for n=n1:n2,
+for n=1:length(elements),
         if ~isempty(elements(n).element),
                 [stress_vector,normal]=StressSurface(elements(n).element,grids,materials,inputs);

issm/trunk/src/m/solutions/ice/SystemMatrices.m

@@ -13 +13 @@
 pg={};
 
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
-
 if kflag, 
         
@@ -22 +19 @@
 
         %Go through elements and build stiffness matrix 
-        for n=n1:n2,
+        for n=1:length(elements),
 
                 if ~isempty(elements(n).element),
@@ -42 +39 @@
 
         %Go through elements and build loads
-        for n=n1:n2,
+        for n=1:length(elements),
                 
                 if ~isempty(elements(n).element),

issm/trunk/src/m/solutions/ice/ViscousHeatingCompute.m

@@ -18 +18 @@
 gridset=m.gridset;
 
-%figure out active elements that will take part in the stiffness and load generation
-[n1,n2]=GetNumberOfActiveElements(elements);
-
 %initialize vector
-viscousheating=zeros((n2-n1)+1,1);
+viscousheating=zeros(numberofelements,1);
 
 %Go through all elements and call the vicous heating routine.
-for n=n1:n2,
+for n=1:length(elements),
                 
         if ~isempty(elements(n).element),