Changeset 1755
- Timestamp:
- 08/18/09 15:16:26 (15 years ago)
- Location:
- issm/trunk/src/m/classes
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/classes/@model/model.m
r1666 r1755 133 133 %Boundary conditions 134 134 md.gridonboundary=NaN; 135 %Diagnostic 136 md.segmentonneumann_diag=NaN; 137 md.neumannvalues_diag=NaN; 138 md.gridondirichlet_diag=NaN; 139 md.dirichletvalues_diag=NaN; 140 141 %Diagnostic stokes 142 md.segmentonneumann_diag_stokes=NaN; 143 144 %Thermal 135 md.pressureload=NaN; 136 md.pressureload_stokes=NaN; 137 md.spcvelocity=NaN; 138 md.spctemperature=NaN; 139 md.spcthickness=NaN; 145 140 md.min_thermal_constraints=0; 146 md.gridondirichlet_thermal=NaN;147 md.dirichletvalues_thermal=NaN;148 149 %Transient150 md.segmentonneumann_prog=NaN;151 md.neumannvalues_prog=NaN;152 md.segmentonneumann_prog2=NaN;153 md.neumannvalues_prog2=NaN;154 md.gridondirichlet_prog=NaN;155 md.dirichletvalues_prog=NaN;156 141 157 142 %Observations -
issm/trunk/src/m/classes/public/BasinConstrain.m
r1338 r1755 48 48 49 49 %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. 50 md. gridondirichlet_diag(gridnotondomain)=1;51 md. dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);52 md. dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);50 md.spcvelocity(gridnotondomain,1:2)=1; 51 md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); 52 md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); 53 53 md.elementonwater(elementnotondomain)=1; 54 54 … … 57 57 numpos=unique(md.elements(pos,:)); 58 58 grids=setdiff(1:1:md.numberofgrids,numpos); 59 md. gridondirichlet_diag(grids)=1;60 md. dirichletvalues_diag(grids,1)=md.vx_obs(grids);61 md. dirichletvalues_diag(grids,2)=md.vy_obs(grids);59 md.spcvelocity(grids,1:2)=1; 60 md.spcvelocity(grids,4)=md.vx_obs(grids); 61 md.spcvelocity(grids,5)=md.vy_obs(grids); 62 62 63 63 %make sure icefronts that are completely spc'd are taken out: 64 free_segments=find(sum(md. gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);65 md. segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);64 free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); 65 md.pressureload=md.pressureload(free_segments,:); -
issm/trunk/src/m/classes/public/BasinConstrain2.m
r1317 r1755 48 48 49 49 %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. 50 md. gridondirichlet_diag(gridnotondomain)=1;51 md. dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);52 md. dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);50 md.spcvelocity(gridnotondomain,1:2)=1; 51 md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); 52 md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); 53 53 md.elementonwater(elementnotondomain)=1; 54 54 … … 57 57 numpos=unique(md.elements(pos,:)); 58 58 grids=setdiff(1:1:md.numberofgrids,numpos); 59 md. gridondirichlet_diag(grids)=1;60 md. dirichletvalues_diag(grids,1)=md.vx_obs(grids);61 md. dirichletvalues_diag(grids,2)=md.vy_obs(grids);59 md.spcvelocity(grids,1:2)=1; 60 md.spcvelocity(grids,4)=md.vx_obs(grids); 61 md.spcvelocity(grids,5)=md.vy_obs(grids); 62 62 63 63 64 64 %make sure icefronts that are completely spc'd are taken out: 65 free_segments=find(sum(md. gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);66 md. segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);65 free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); 66 md.pressureload=md.pressureload(free_segments,:); -
issm/trunk/src/m/classes/public/BasinConstrainShelf.m
r1317 r1755 48 48 49 49 %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. 50 md. gridondirichlet_diag(gridnotondomain)=1;51 md. dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);52 md. dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);50 md.spcvelocity(gridnotondomain,1:2)=1; 51 md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); 52 md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); 53 53 md.elementonwater(elementnotondomain)=1; 54 54 … … 57 57 numpos=unique(md.elements(pos,:)); 58 58 grids=setdiff(1:1:md.numberofgrids,numpos); 59 md. gridondirichlet_diag(grids)=1;60 md. dirichletvalues_diag(grids,1)=md.vx_obs(grids);61 md. dirichletvalues_diag(grids,2)=md.vy_obs(grids);59 md.spcvelocity(grids,1:2)=1; 60 md.spcvelocity(grids,4)=md.vx_obs(grids); 61 md.spcvelocity(grids,5)=md.vy_obs(grids); 62 62 63 63 64 64 %make sure any grid with NaN velocity is spc'd: 65 65 pos=find(isnan(md.vel_obs_raw)); 66 md. gridondirichlet_diag(pos)=1;66 md.spcvelocity(pos,1:2)=1; 67 67 %we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand. 68 md. dirichletvalues_diag(pos,1)=md.vx_obs(pos);69 md. dirichletvalues_diag(pos,2)=md.vy_obs(pos);68 md.spcvelocity(pos,4)=md.vx_obs(pos); 69 md.spcvelocity(pos,5)=md.vy_obs(pos); 70 70 71 71 %iceshelves: any grid on icesheet is spc'd 72 72 pos=find(md.gridonicesheet); 73 md. gridondirichlet_diag(pos)=1;74 md. dirichletvalues_diag(pos,1)=md.vx_obs(pos);75 md. dirichletvalues_diag(pos,2)=md.vy_obs(pos);73 md.spcvelocity(pos,1:2)=1; 74 md.spcvelocity(pos,4)=md.vx_obs(pos); 75 md.spcvelocity(pos,5)=md.vy_obs(pos); 76 76 77 77 %make sure icefronts that are completely spc'd are taken out: 78 free_segments=find(sum(md. gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);79 md. segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);78 free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); 79 md.pressureload=md.pressureload(free_segments,:); -
issm/trunk/src/m/classes/public/collapse.m
r1 r1755 61 61 62 62 %boundary conditions 63 md.gridondirichlet_diag=project2d(md,md.gridondirichlet_diag,md.numlayers); 64 dirichletvalues_diag2d(:,1)=project2d(md,md.dirichletvalues_diag(:,1),md.numlayers); 65 dirichletvalues_diag2d(:,2)=project2d(md,md.dirichletvalues_diag(:,2),md.numlayers); 66 md.dirichletvalues_diag=dirichletvalues_diag2d; 63 md.spcvelocity=project2d(md,md.spcvelocity,md.numlayers); 64 md.spcthickness=project2d(md,md.spcthickness,md.numlayers); 65 md.spctemperature=project2d(md,md.spctemperature,md.numlayers); 67 66 68 67 %Extrusion of Neumann BC 69 68 %in 2d, segmentonnumann is: [grid1 grid2 element] 70 numberofneumann2d=size(md.segmentonneumann_diag,1)/md.numlayers; 71 md.segmentonneumann_diag=[md.segmentonneumann_diag(1:numberofneumann2d,1:2) md.segmentonneumann_diag(1:numberofneumann2d,5)]; %Add two columns on the first layer 72 73 %Prognostic 74 md.gridondirichlet_prog=project2d(md,md.gridondirichlet_prog,md.numlayers); 75 md.dirichletvalues_prog=project2d(md,md.dirichletvalues_prog,md.numlayers); 76 %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,5))]; 77 %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,5))]; 69 numberofneumann2d=size(md.pressureload,1)/md.numlayers; 70 md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5)]; %Add two columns on the first layer 78 71 79 72 %materials … … 85 78 md.observed_temperature=DepthAverage(md,md.observed_temperature); 86 79 md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux 87 md.gridondirichlet_thermal=project2d(md,md.gridondirichlet_thermal,md.numlayers); %surface temperature88 md.dirichletvalues_thermal=project2d(md,md.dirichletvalues_thermal,md.numlayers); %surface temperature89 90 %NaN the values that are not on an spc'd temperature.91 pos=find(~md.gridondirichlet_thermal);92 md.dirichletvalues_thermal(pos)=NaN;93 80 94 81 %update of connectivity matrix -
issm/trunk/src/m/classes/public/display/displaybc.m
r1277 r1755 19 19 20 20 disp(sprintf('\n diagnostic:')); 21 fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list'); 22 fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]'); 23 fielddisplay(md,'segmentonneumann_diag','segments on ice front list'); 24 fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]'); 21 fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)'); 22 fielddisplay(md,'pressureload','segments on ice front list'); 25 23 26 24 disp(sprintf('\n prognostic:')); 27 fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list'); 28 fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]'); 29 fielddisplay(md,'segmentonneumann_prog','segments on ice front list'); 25 fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)'); 30 26 fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]'); 31 27 32 28 disp(sprintf('\n thermal:')); 33 fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list'); 34 fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]'); 29 fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)'); 35 30 fielddisplay(md,'melting','melting rate [m/a]'); -
issm/trunk/src/m/classes/public/display/displaydiagnostic.m
r1666 r1755 18 18 19 19 disp(sprintf('\n boundary conditions:')); 20 fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list'); 21 fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]'); 22 fielddisplay(md,'segmentonneumann_diag','segments on ice front list'); 23 fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]'); 20 fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)'); 21 fielddisplay(md,'pressureload','segments on ice front list'); 24 22 25 23 disp(sprintf('\n %s','Penalties:')); -
issm/trunk/src/m/classes/public/display/displayprognostic.m
r1277 r1755 17 17 18 18 disp(sprintf('\n boundary conditions:')); 19 fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list'); 20 fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]'); 19 fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)'); 21 20 fielddisplay(md,'segmentonneumann_prog','segments on ice front list'); 22 21 fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]'); -
issm/trunk/src/m/classes/public/display/displaythermal.m
r1277 r1755 12 12 13 13 disp(sprintf('\n boundary conditions:')); 14 fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list'); 15 fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]'); 14 fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)'); 16 15 fielddisplay(md,'melting','melting rate [m/a]'); 17 16 -
issm/trunk/src/m/classes/public/extrude.m
r1651 r1755 145 145 146 146 %boundary conditions 147 md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node'); 148 md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node'); 147 md.spcvelocity=project3d(md,md.spcvelovity,'node'); 148 md.spctemperature=project3d(md,md.spcvelovity,'node'); 149 md.spcthickness=project3d(md,md.spcvelovity,'node'); 149 150 150 151 %Extrusion of Neumann BC 151 152 %in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element] 152 old segmentonneumann_diag=md.segmentonneumann_diag;153 segmentonneumann_diag_layer1=[oldsegmentonneumann_diag(:,1:2) oldsegmentonneumann_diag(:,2)+md.numberofgrids2d oldsegmentonneumann_diag(:,1)+md.numberofgrids2d oldsegmentonneumann_diag(:,3)]; %Add two columns on the first layer154 segmentonneumann_diag=[];153 oldpressureload=md.pressureload; 154 pressureload_layer1=[oldpressureload(:,1:2) oldpressureload(:,2)+md.numberofgrids2d oldpressureload(:,1)+md.numberofgrids2d oldpressureload(:,3)]; %Add two columns on the first layer 155 pressureload=[]; 155 156 for i=1:numlayers-1, 156 segmentonneumann_diag=[segmentonneumann_diag ;segmentonneumann_diag_layer1(:,1:4)+(i-1)*md.numberofgrids2d segmentonneumann_diag_layer1(:,5)+(i-1)*md.numberofelements2d ];157 pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.numberofgrids2d pressureload_layer1(:,5)+(i-1)*md.numberofelements2d ]; 157 158 end 158 159 159 160 %plug into md 160 md.segmentonneumann_diag=segmentonneumann_diag; 161 md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node'); 162 md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node'); 163 %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))]; 164 %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))]; 161 md.pressureload=pressureload; 165 162 166 163 %materials … … 184 181 md.observed_temperature=project3d(md,md.observed_temperature,'node'); 185 182 md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux 186 md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers); %surface temperature187 md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers); %surface temperature188 189 %NaN the values that are not on an spc'd temperature.190 pos=find(~md.gridondirichlet_thermal);191 md.dirichletvalues_thermal(pos)=NaN;192 183 193 184 %increase connectivity if less than 25: 194 185 if md.connectivity<25, 195 md.connectivity= 48;186 md.connectivity=70; 196 187 end 197 188 -
issm/trunk/src/m/classes/public/ismodelselfconsistent.m
r1737 r1755 66 66 67 67 %NO NAN 68 fields={'numberofelements','numberofgrids','x','y','z',' segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...68 fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',... 69 69 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',... 70 70 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; … … 79 79 80 80 %FIELDS >= 0 81 fields={'numberofelements','numberofgrids','elements',' segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...81 fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',... 82 82 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',... 83 83 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; … … 96 96 97 97 %FIELDS ~=0 98 fields={'numberofelements','numberofgrids','elements',' segmentonneumann_diag','drag_type',...98 fields={'numberofelements','numberofgrids','elements','drag_type',... 99 99 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',... 100 100 'sparsity','deltaH','DeltaH','timeacc','timedec'}; … … 118 118 119 119 %SIZE NUMBEROFGRIDS 120 fields={'x','y','z','B','drag',' gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};120 fields={'x','y','z','B','drag','scpvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'}; 121 121 for i=1:length(fields), 122 122 if length(md.(fields{i}))~=md.numberofgrids, … … 194 194 195 195 %SINGULAR 196 if ~any( md.gridondirichlet_diag),196 if ~any(sum(md.spcvelocity(:,1:3),2)==3), 197 197 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) 198 bool=0;return;199 end200 201 %DIRICHLET VALUES202 if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),203 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);204 198 bool=0;return; 205 199 end … … 208 202 if any(md.thickness<=0), 209 203 pos=find(md.thickness<=0); 210 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),204 if any(find(md.spcthickness(pos,1)==0)), 211 205 disp(['model ' md.name ' has some grids with 0 thickness']); 212 206 bool=0; return; … … 217 211 %PROGNOSTIC 218 212 if strcmp(md.analysis_type,'prognostic'), 219 220 %DIRICHLET VALUES221 if isempty(md.gridondirichlet_prog) | isempty(md.dirichletvalues_prog),222 disp(['model ' md.name ' is not well posed. Missing dirichlet values for prognostic run']);223 bool=0;return;224 end225 213 226 214 %VELOCITIES … … 339 327 340 328 %SINGULAR 341 if ~any( md.gridondirichlet_diag),329 if ~any(sum(md.spcvelocity(:,1:3),2)==3), 342 330 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) 343 bool=0;return;344 end345 346 %DIRICHLET VALUES347 if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),348 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);349 331 bool=0;return; 350 332 end … … 353 335 if any(md.thickness<=0), 354 336 pos=find(md.thickness<=0); 355 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),337 if any(find(md.spcthickness(pos,1)==0)), 356 338 disp(['model ' md.name ' has some grids with 0 thickness']); 357 339 bool=0; return; -
issm/trunk/src/m/classes/public/marshall.m
r1666 r1755 77 77 WriteData(fid,md.gridoniceshelf,'Mat','gridoniceshelf'); 78 78 79 WriteData(fid,md.segmentonneumann_diag,'Mat','segmentonneumann_diag'); 80 WriteData(fid,md.segmentonneumann_diag_stokes,'Mat','segmentonneumann_diag_stokes'); 81 WriteData(fid,md.neumannvalues_diag,'Mat','neumannvalues_diag'); 82 WriteData(fid,md.gridondirichlet_diag,'Mat','gridondirichlet_diag'); 83 WriteData(fid,md.dirichletvalues_diag,'Mat','dirichletvalues_diag'); 79 WriteData(fid,md.spcvelocity,'Mat','spcvelocity'); 80 WriteData(fid,md.spctemperature,'Mat','spctemperature'); 81 WriteData(fid,md.spcthickness,'Mat','spcthickness'); 84 82 85 WriteData(fid,md.segmentonneumann_prog,'Mat','segmentonneumann_prog'); 86 WriteData(fid,md.neumannvalues_prog,'Mat','neumannvalues_prog'); 87 WriteData(fid,md.gridondirichlet_prog,'Mat','gridondirichlet_prog'); 88 WriteData(fid,md.dirichletvalues_prog,'Mat','dirichletvalues_prog'); 89 WriteData(fid,md.segmentonneumann_prog2,'Mat','segmentonneumann_prog2'); 90 WriteData(fid,md.neumannvalues_prog2,'Mat','neumannvalues_prog2'); 83 WriteData(fid,md.pressureload,'Mat','pressureload'); 84 WriteData(fid,md.pressureload_stokes,'Mat','pressureload_stokes'); 91 85 92 WriteData(fid,md.gridondirichlet_thermal,'Mat','gridondirichlet_thermal');93 WriteData(fid,md.dirichletvalues_thermal,'Mat','dirichletvalues_thermal');94 86 WriteData(fid,md.geothermalflux,'Mat','geothermalflux'); 95 87 WriteData(fid,md.melting,'Mat','melting'); -
issm/trunk/src/m/classes/public/modelextract.m
r1730 r1755 188 188 gridstoflag1=intersect(orphans_grid,pos_grid); 189 189 gridstoflag2=Pgrid(gridstoflag1); 190 if ~isnan(md1. gridondirichlet_diag),191 md2. gridondirichlet_diag(gridstoflag2)=1;190 if ~isnan(md1.spcvelocity), 191 md2.spcvelocity(gridstoflag2,1:2)=1; 192 192 if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs) 193 md2. dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2);194 md2. dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);193 md2.spcvelocity(gridstoflag2,4)=md2.vx_obs(gridstoflag2); 194 md2.spcvelocity(gridstoflag2,5)=md2.vy_obs(gridstoflag2); 195 195 else 196 md2. dirichletvalues_diag=zeros(numberofgrids2,2);197 disp(' ') 198 disp('!! modelextract warning: dirichletvalues should be checked !!')199 disp(' ') 200 end 201 end 202 if ~isnan(md1. gridondirichlet_thermal),203 md2. gridondirichlet_thermal(gridstoflag2)=1;196 md2.spcvelocity(gridstoflag2,4:5)=zeros(length(gridstoflag2)),2); 197 disp(' ') 198 disp('!! modelextract warning: spc values should be checked !!') 199 disp(' ') 200 end 201 end 202 if ~isnan(md1.spctemperature), 203 md2.spctemperature(gridstoflag2,1)=1; 204 204 if ~isnan(md1.observed_temperature) 205 md2. dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2);205 md2.spctemperature(gridstoflag2,2)=md2.observed_temperature(gridstoflag2); 206 206 else 207 md2. dirichletvalues_thermal=zeros(numberofgrids2,2);208 disp(' ') 209 disp('!! modelextract warning: dirichletvalues should be checked !!')207 md2.spctemperature(gridstoflag2,2)=zeros(length(gridstoflag2),2); 208 disp(' ') 209 disp('!! modelextract warning: spc values should be checked !!') 210 210 disp(' ') 211 211 end … … 214 214 %Recreate border stokes 215 215 if md1.isstokes, 216 pos=find( md2.gridondirichlet_diag); %find all the grids on the boundary of the domain without icefront216 pos=find(sum(md2.spcvelocity(:,1:3),2)); %find all the grids on the boundary of the domain without icefront 217 217 md2.gridonstokes(pos)=0; %we have to constrain all the boundary of the domain without icefront 218 218 stokes_elements=find(md2.elements_type(:,2)==StokesFormulationEnum()); %find the elements on the stokes domain … … 223 223 224 224 %Diagnostic 225 if ~isnan(md2. segmentonneumann_diag)226 md2. segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1));227 md2. segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2));228 md2. segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end));225 if ~isnan(md2.pressureload) 226 md2.pressureload(:,1)=Pgrid(md1.pressureload(:,1)); 227 md2.pressureload(:,2)=Pgrid(md1.pressureload(:,2)); 228 md2.pressureload(:,end)=Pelem(md1.pressureload(:,end)); 229 229 if strcmpi(md1.type,'3d') 230 md2. segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3));231 md2. segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4));232 end 233 md2. segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);234 end 235 if ~isnan(md2. segmentonneumann_diag_stokes)236 md2. segmentonneumann_diag_stokes(:,1)=Pgrid(md1.segmentonneumann_diag_stokes(:,1));237 md2. segmentonneumann_diag_stokes(:,2)=Pgrid(md1.segmentonneumann_diag_stokes(:,2));238 md2. segmentonneumann_diag_stokes(:,end)=Pelem(md1.segmentonneumann_diag_stokes(:,end));230 md2.pressureload(:,3)=Pgrid(md1.pressureload(:,3)); 231 md2.pressureload(:,4)=Pgrid(md1.pressureload(:,4)); 232 end 233 md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2)),:); 234 end 235 if ~isnan(md2.pressureload_stokes) 236 md2.pressureload_stokes(:,1)=Pgrid(md1.pressureload_stokes(:,1)); 237 md2.pressureload_stokes(:,2)=Pgrid(md1.pressureload_stokes(:,2)); 238 md2.pressureload_stokes(:,end)=Pelem(md1.pressureload_stokes(:,end)); 239 239 if strcmpi(md1.type,'3d') 240 md2.segmentonneumann_diag_stokes(:,3)=Pgrid(md1.segmentonneumann_diag_stokes(:,3)); 241 md2.segmentonneumann_diag_stokes(:,4)=Pgrid(md1.segmentonneumann_diag_stokes(:,4)); 242 end 243 md2.segmentonneumann_diag_stokes=md2.segmentonneumann_diag_stokes(find(md2.segmentonneumann_diag_stokes(:,1) & md2.segmentonneumann_diag_stokes(:,2)),:); 244 end 245 md2.neumannvalues_diag=NaN; 246 247 %Transient 248 if ~isnan(md2.segmentonneumann_prog) 249 md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1)); 250 md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2)); 251 md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end)); 252 md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:); 253 end 254 md2.neumannvalues_prog=NaN; 255 if ~isnan(md2.segmentonneumann_prog2) 256 md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1)); 257 md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2)); 258 md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end)); 259 md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:); 240 md2.pressureload_stokes(:,3)=Pgrid(md1.pressureload_stokes(:,3)); 241 md2.pressureload_stokes(:,4)=Pgrid(md1.pressureload_stokes(:,4)); 242 end 243 md2.pressureload_stokes=md2.pressureload_stokes(find(md2.pressureload_stokes(:,1) & md2.pressureload_stokes(:,2)),:); 260 244 end 261 245 -
issm/trunk/src/m/classes/public/setelementstype.m
r1682 r1755 61 61 gridonstokes=ones(md.numberofgrids,1); 62 62 gridonstokes(md.elements(nonstokesflag,:))=0; %non stokes grids 63 pos=find( md.gridondirichlet_diag);%find all the grids on the boundary of the domain without icefront63 pos=find(sum(md.spcvelocity(:,1:3),2)==3); %find all the grids on the boundary of the domain without icefront 64 64 gridonstokes(pos)=0; %we have to constrain all the boundary of the domain without icefront 65 65
Note:
See TracChangeset
for help on using the changeset viewer.