0001 function md2=modelextract(md1,area)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 if ((nargin~=2) | (nargout~=1)),
0023 error('modelextract error message');
0024 end
0025 package='Ice';
0026
0027
0028 if ischar(area),
0029 if isempty(area),
0030 flag_elem=zeros(md1.numberofelements,1);
0031 invert=0;
0032 elseif strcmpi(area,'all')
0033 flag_elem=ones(md1.numberofelements,1);
0034 invert=0;
0035 else
0036
0037 if strcmpi(area(1),'~'),
0038 area=area(2:length(area));
0039 invert=1;
0040 else
0041 invert=0;
0042 end
0043
0044 A=expread(area,1);
0045 flag_elem=ArgusContourToMesh(md1.elements(:,1:3),md1.x,md1.y,A,'element',1);
0046
0047 pos=find(flag_elem);
0048 for contour=1:length(A)
0049 for i=1: length(pos)
0050 test=cross_extract(md1.x(md1.elements(pos(i),:)),md1.y(md1.elements(pos(i),:)),A(contour).x,A(contour).y);
0051 if test
0052 flag_elem(pos(i))=0;
0053 end
0054 end
0055 end
0056 end
0057 if invert,
0058 flag_elem=~flag_elem;
0059 end
0060 elseif isfloat(area),
0061 if size(area,1)~=md1.numberofelements,
0062 setelementstypeusage();
0063 error('Flags must be of same size as number of elements in model');
0064 end
0065 flag_elem=area;
0066 else
0067 error('Invalide option');
0068 end
0069 pos_elem=find(flag_elem);
0070
0071
0072 flag_grid=ismember([1:md1.numberofgrids]',sort(unique(md1.elements(pos_elem,:)))');
0073 pos_grid=find(flag_grid);
0074
0075
0076 numberofgrids1=md1.numberofgrids;
0077 numberofelements1=md1.numberofelements;
0078 numberofgrids2=length(pos_grid);
0079 numberofelements2=length(pos_elem);
0080
0081
0082 Pelem=zeros(numberofelements1,1);
0083 Pelem(pos_elem)=[1:numberofelements2]';
0084 Pgrid=zeros(numberofgrids1,1);
0085 Pgrid(pos_grid)=[1:numberofgrids2]';
0086
0087
0088 elements_1=md1.elements;
0089 elements_2=elements_1(pos_elem,:);
0090
0091 elements_2(:,1)=Pgrid(elements_2(:,1));
0092 elements_2(:,2)=Pgrid(elements_2(:,2));
0093 elements_2(:,3)=Pgrid(elements_2(:,3));
0094 if strcmpi(md1.type,'3d'),
0095 elements_2(:,4)=Pgrid(elements_2(:,4));
0096 elements_2(:,5)=Pgrid(elements_2(:,5));
0097 elements_2(:,6)=Pgrid(elements_2(:,6));
0098 end
0099
0100
0101
0102
0103 md2=md1;
0104
0105
0106 md2.numberofelements=numberofelements2;
0107 md2.numberofgrids=numberofgrids2;
0108 md2.elements=elements_2;
0109 if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
0110 md2.x=md1.x(pos_grid);
0111 md2.y=md1.y(pos_grid);
0112 md2.z=md1.z(pos_grid);
0113 if ~isnan(md2.bed_slopex), md2.bed_slopex=md1.bed_slopex(pos_grid); end;
0114 if ~isnan(md2.bed_slopey), md2.bed_slopex=md1.bed_slopey(pos_grid); end;
0115 if ~isnan(md2.surface_slopex), md2.surface_slopex=md1.surface_slopex(pos_grid); end;
0116 if ~isnan(md2.surface_slopey), md2.surface_slopex=md1.surface_slopey(pos_grid); end;
0117
0118
0119 if strcmpi(md1.type,'3d')
0120 flag_elem_2d=flag_elem(1:md1.numberofelements2d);
0121 pos_elem_2d=find(flag_elem_2d);
0122 flag_grid_2d=flag_grid(1:md1.numberofgrids2d);
0123 pos_grid_2d=find(flag_grid_2d);
0124
0125 md2.numberofelements2d=length(pos_elem_2d);
0126 md2.numberofgrids2d=length(pos_grid_2d);
0127 md2.elements2d=md1.elements2d(pos_elem_2d,:);
0128 md2.elements2d(:,1)=Pgrid(md2.elements2d(:,1));
0129 md2.elements2d(:,2)=Pgrid(md2.elements2d(:,2));
0130 md2.elements2d(:,3)=Pgrid(md2.elements2d(:,3));
0131
0132 if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d,:); end;
0133 md2.x2d=md1.x(pos_grid_2d);
0134 md2.y2d=md1.y(pos_grid_2d);
0135 md2.z2d=md1.z(pos_grid_2d);
0136 end
0137
0138
0139 if ~isnan(md2.elementonhutter), md2.elementonhutter=md1.elementonhutter(pos_elem); end;
0140 if ~isnan(md2.elementonmacayeal), md2.elementonmacayeal=md1.elementonmacayeal(pos_elem); end;
0141 if ~isnan(md2.elementonpattyn), md2.elementonpattyn=md1.elementonpattyn(pos_elem); end;
0142 if ~isnan(md2.elementonstokes), md2.elementonstokes=md1.elementonstokes(pos_elem); end;
0143
0144
0145 if ~isnan(md2.gridonhutter), md2.gridonhutter=md1.gridonhutter(pos_grid); end;
0146 if ~isnan(md2.gridonmacayeal), md2.gridonmacayeal=md1.gridonmacayeal(pos_grid); end;
0147 if ~isnan(md2.gridonpattyn), md2.gridonpattyn=md1.gridonpattyn(pos_grid); end;
0148 if ~isnan(md2.gridonstokes), md2.gridonstokes=md1.gridonstokes(pos_grid); end;
0149
0150
0151 if ~isnan(md2.penalties),
0152 for i=1:size(md1.penalties,1);
0153 md2.penalties(i,:)=Pgrid(md1.penalties(i,:));
0154 end
0155 md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
0156 end
0157 md2.rifts=NaN;
0158 md2.numrifts=0;
0159
0160
0161 if ~isnan(md2.uppergrids), md2.uppergrids=Pgrid(md1.uppergrids(pos_grid)); end;
0162 if ~isnan(md2.lowergrids), md2.lowergrids=Pgrid(md1.lowergrids(pos_grid)); end;
0163 if ~isnan(md2.deadgrids), md2.deadgrids=md1.deadgrids(pos_grid);end;
0164
0165
0166 md2.elementonbed=md1.elementonbed(pos_elem);
0167 md2.elementonsurface=md1.elementonsurface(pos_elem);
0168 md2.gridonbed=md1.gridonbed(pos_grid);
0169 md2.gridonsurface=md1.gridonsurface(pos_grid);
0170 if ~isnan(md2.firn_layer), md2.firn_layer=md1.firn_layer(pos_grid); end;
0171
0172
0173 if ~isnan(md2.drag), md2.drag=md1.drag(pos_grid);end;
0174 if ~isnan(md2.p), md2.p=md1.p(pos_elem);end;
0175 if ~isnan(md2.q), md2.q=md1.q(pos_elem);end;
0176 if ~isnan(md2.n), md2.n=md1.n(pos_elem);end;
0177 if ~isnan(md2.B), md2.B=md1.B(pos_grid);end;
0178
0179
0180 if ~isnan(md2.elementoniceshelf), md2.elementoniceshelf=md1.elementoniceshelf(pos_elem);end;
0181 if ~isnan(md2.elementonicesheet), md2.elementonicesheet=md1.elementonicesheet(pos_elem);end;
0182 if ~isnan(md2.gridoniceshelf), md2.gridoniceshelf=md1.gridoniceshelf(pos_grid);end;
0183 if ~isnan(md2.gridonicesheet), md2.gridonicesheet=md1.gridonicesheet(pos_grid);end;
0184 if ~isnan(md2.surface), md2.surface=md1.surface(pos_grid); end;
0185 if ~isnan(md2.thickness), md2.thickness=md1.thickness(pos_grid); end;
0186 if ~isnan(md2.new_thickness), md2.new_thickness=md1.new_thickness(pos_grid); end;
0187 if ~isnan(md2.bed), md2.bed=md1.bed(pos_grid); end;
0188
0189
0190 md2.gridonboundary=md1.gridonboundary(pos_grid);
0191
0192
0193 if ~isnan(md2.segmentonneumann_diag)
0194 md2.segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1));
0195 md2.segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2));
0196 md2.segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end));
0197 if strcmpi(md1.type,'3d')
0198 md2.segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3));
0199 md2.segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4));
0200 end
0201 md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);
0202 end
0203 md2.neumannvalues_diag=NaN;
0204 if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag=md1.gridondirichlet_diag(pos_grid); end;
0205 if ~isnan(md2.dirichletvalues_diag),
0206 md2.dirichletvalues_diag=[];
0207 for i=1:size(md1.dirichletvalues_diag,2)
0208 md2.dirichletvalues_diag(:,i)=md1.dirichletvalues_diag(pos_grid,i);
0209 end
0210 end
0211
0212
0213 if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal=md1.gridondirichlet_thermal(pos_grid); end;
0214 if ~isnan(md2.dirichletvalues_thermal), md2.dirichletvalues_thermal=md1.dirichletvalues_thermal(pos_grid); end;
0215
0216
0217 if ~isnan(md2.segmentonneumann_prog)
0218 md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1));
0219 md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2));
0220 md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end));
0221 md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:);
0222 end
0223 md2.neumannvalues_prog=NaN;
0224 if ~isnan(md2.segmentonneumann_prog2)
0225 md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1));
0226 md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2));
0227 md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end));
0228 md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:);
0229 end
0230 md2.neumannvalues_prog=NaN;
0231 if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog=md1.gridondirichlet_prog(pos_grid); end;
0232 if ~isnan(md2.dirichletvalues_prog), md2.dirichletvalues_prog=md1.dirichletvalues_prog(pos_grid); end;
0233
0234
0235 if ~isnan(md2.vx_obs), md2.vx_obs=md1.vx_obs(pos_grid); end;
0236 if ~isnan(md2.vy_obs), md2.vy_obs=md1.vy_obs(pos_grid); end;
0237 if ~isnan(md2.vel_obs), md2.vel_obs=md1.vel_obs(pos_grid); end;
0238 if ~isnan(md2.accumulation), md2.accumulation=md1.accumulation(pos_grid); end;
0239 if ~isnan(md2.geothermalflux), md2.geothermalflux=md1.geothermalflux(pos_grid); end;
0240 if ~isnan(md2.observed_temperature), md2.observed_temperature=md1.observed_temperature(pos_grid); end;
0241
0242
0243 md2.viscousheating=NaN;
0244 md2.pressure_elem=NaN;
0245 md2.stress=NaN;
0246 md2.stress_surface=NaN;
0247 md2.stress_bed=NaN;
0248 md2.deviatoricstress=NaN;
0249 md2.strainrate=NaN;
0250
0251
0252 if ~isnan(md2.vx), md2.vx=md1.vx(pos_grid); end;
0253 if ~isnan(md2.vy), md2.vy=md1.vy(pos_grid); end;
0254 if ~isnan(md2.vz), md2.vz=md1.vz(pos_grid); end;
0255 if ~isnan(md2.vel), md2.vel=md1.vel(pos_grid); end;
0256 if ~isnan(md2.temperature), md2.temperature=md1.temperature(pos_grid); end;
0257 if ~isnan(md2.melting), md2.melting=md1.melting(pos_grid); end;
0258 if ~isnan(md2.pressure), md2.pressure=md1.pressure(pos_grid); end;
0259
0260
0261
0262
0263 orphans_elem=find(~flag_elem);
0264 orphans_grid=unique(md1.elements(orphans_elem,:))';
0265
0266
0267 gridstoflag1=intersect(orphans_grid,pos_grid);
0268 gridstoflag2=Pgrid(gridstoflag1);
0269
0270
0271 flaglist=ismember(md1.elements,gridstoflag1);
0272 elementstoflag1=find((flaglist(:,1) | flaglist(:,2) | flaglist(:,3)) & ~flag_elem);
0273 flaglist=ismember(md2.elements,gridstoflag2);
0274 elementstoflag2=find(flaglist(:,1) | flaglist(:,2) | flaglist(:,3));
0275
0276
0277 test=0;
0278 if ~isnan(md2.gridonboundary), md2.gridonboundary(gridstoflag2)=1; end;
0279 if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag(gridstoflag2)=1; end;
0280 if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog(gridstoflag2)=1; end;
0281 if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal(gridstoflag2)=1; end;
0282 if ~isnan(md2.dirichletvalues_diag),
0283 test=1;
0284 if strcmpi(md1.type,'3d')
0285 if ~isnan(md2.vx) & ~isnan(md2.vy)
0286 md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx(gridstoflag2);
0287 md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy(gridstoflag2);
0288 end
0289 else
0290 if ~isnan(md2.vx_obs) & ~isnan(md2.vy_obs)
0291 md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2);
0292 md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
0293 end
0294 end
0295 end
0296 if ~isnan(md2.dirichletvalues_prog), test=1; end;
0297 if ~isnan(md2.dirichletvalues_thermal), test=1; end;
0298
0299 if test
0300 disp(' ')
0301 disp('!! modelextract warning: dirichlet values should be checked !!')
0302 disp(' ')
0303 end
0304
0305
0306 md2.extractedgrids=pos_grid;
0307 md2.extractedelements=pos_elem;
0308
0309
0310 elementsonboundary=find(sum(md2.gridonboundary(md2.elements),2)>=2);
0311 md2.segments=[];
0312 for i=1:length(elementsonboundary),
0313 el=elementsonboundary(i);
0314 for j=1:3,
0315 grid1=md2.elements(el,j);
0316 if (j==3),
0317 grid2=md2.elements(el,1);
0318 else
0319 grid2=md2.elements(el,j+1);
0320 end
0321
0322 pos=find( ((md2.elements(:,1)==grid1) | (md2.elements(:,2)==grid1) | (md2.elements(:,3)==grid1)) & ...
0323 ((md2.elements(:,1)==grid2) | (md2.elements(:,2)==grid2) | (md2.elements(:,3)==grid2))...
0324 );
0325 if length(pos)==1,
0326 md2.segments=[md2.segments; [grid1 grid2 pos]];
0327 end
0328 end
0329 end
0330 md2.segmentmarkers=ones(size(md2.segments,1));
0331
0332
0333 end
0334
0335 function bool=cross_extract(x1,y1,x2,y2);
0336
0337 x1 = x1(:); y1 = y1(:);
0338 x2 = x2(:); y2 = y2(:);
0339 l1 = length(x1);
0340 l2 = length(x2);
0341
0342
0343 [i11,i12] = meshgrid(1:l1,1:l2);
0344 [i21,i22] = meshgrid([2:l1 1],[2:l2 1]);
0345 i11 = i11(:); i12 = i12(:);
0346 i21 = i21(:); i22 = i22(:);
0347
0348
0349 S = iscross([x1(i11) x1(i21)]',[y1(i11) y1(i21)]',...
0350 [x2(i12) x2(i22)]',[y2(i12) y2(i22)]')';
0351 S(find(S==0.5))=0;
0352
0353 bool=any(any(S));
0354 end