Index: /issm/trunk/src/m/solutions/ice/DeviatoricStressCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/DeviatoricStressCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/DeviatoricStressCompute.m (revision 51) @@ -14,17 +14,15 @@ 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,11 +56,11 @@ %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)'; Index: sm/trunk/src/m/solutions/ice/GetNumberOfActiveElements.m =================================================================== --- /issm/trunk/src/m/solutions/ice/GetNumberOfActiveElements.m (revision 50) +++ (revision ) @@ -1,22 +1,0 @@ -function [n1,n2]=GetNumberOfActiveElements(elements) -%GETNUMBEROFACTIVEELEMENTS - get the number of active elements -% -% Go through all elements, and if we find one with special characteristics, -% like the MacAyeal's accelerated element, then we need to updat nel. -% -% Usage: -% [n1,n2]=GetNumberOfActiveElements(elements) - -for n=length(elements):-1:1, %start from the end, where we usually put the special elements. - - if strcmp(elements(n).element.type,'acceleratedtriaelem'), - %we have an accelerated triaelem, stop the loop, and return n1 n2 - n1=n; - n2=n; - return; - end -end - -%If we are here, no special elements, return n1=1, n2=length(elements) -n1=1; -n2=length(elements); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticBaseVert.m (revision 51) @@ -44,5 +44,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m (revision 51) @@ -52,5 +52,4 @@ 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,5 +86,4 @@ 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,36 +123,4 @@ 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 Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticStokes.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticStokes.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticStokes.m (revision 51) @@ -49,5 +49,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticVert.m (revision 51) @@ -46,5 +46,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorMelting.m (revision 51) @@ -44,5 +44,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorSlopeCompute.m (revision 51) @@ -51,5 +51,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ModelProcessorThermal.m (revision 51) @@ -46,5 +46,4 @@ 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); Index: /issm/trunk/src/m/solutions/ice/PressureElemCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/PressureElemCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/PressureElemCompute.m (revision 51) @@ -18,7 +18,4 @@ 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,9 +28,7 @@ %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 Index: /issm/trunk/src/m/solutions/ice/StrainRateCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/StrainRateCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/StrainRateCompute.m (revision 51) @@ -17,17 +17,15 @@ 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,11 +61,11 @@ %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)'; Index: /issm/trunk/src/m/solutions/ice/StressBedCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/StressBedCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/StressBedCompute.m (revision 51) @@ -17,10 +17,8 @@ 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,14 +29,14 @@ %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); Index: /issm/trunk/src/m/solutions/ice/StressCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/StressCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/StressCompute.m (revision 51) @@ -17,17 +17,16 @@ 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,11 +63,11 @@ %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)'; Index: /issm/trunk/src/m/solutions/ice/StressSurfaceCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/StressSurfaceCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/StressSurfaceCompute.m (revision 51) @@ -17,10 +17,8 @@ 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,14 +29,14 @@ %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); Index: /issm/trunk/src/m/solutions/ice/SystemMatrices.m =================================================================== --- /issm/trunk/src/m/solutions/ice/SystemMatrices.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/SystemMatrices.m (revision 51) @@ -13,7 +13,4 @@ pg={}; -%figure out active elements that will take part in the stiffness and load generation -[n1,n2]=GetNumberOfActiveElements(elements); - if kflag, @@ -22,5 +19,5 @@ %Go through elements and build stiffness matrix - for n=n1:n2, + for n=1:length(elements), if ~isempty(elements(n).element), @@ -42,5 +39,5 @@ %Go through elements and build loads - for n=n1:n2, + for n=1:length(elements), if ~isempty(elements(n).element), Index: /issm/trunk/src/m/solutions/ice/ViscousHeatingCompute.m =================================================================== --- /issm/trunk/src/m/solutions/ice/ViscousHeatingCompute.m (revision 50) +++ /issm/trunk/src/m/solutions/ice/ViscousHeatingCompute.m (revision 51) @@ -18,12 +18,9 @@ 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),