modelextract

PURPOSE ^

MODELEXTRACT - extract a model according to an Argus contour or flag list

SYNOPSIS ^

function md2=modelextract(md1,area)

DESCRIPTION ^

MODELEXTRACT - extract a model according to an Argus contour or flag list

   This routine extracts a submodel from a bigger model with respect to a given contour
   md must be followed by the corresponding exp file or flags list
   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
   If user wants every element outside the domain to be 
   modelextractd, add '~' to the name of the domain file (ex: '~Pattyn.exp');
   an empty string '' will be considered as an empty domain
   a string 'all' will be considered as the entire domain

   Usage:
      md2=modelextract(md1,area);

   Examples:
      md2=modelextract(md,'Domain.exp');
      md2=modelextract(md,md.elementoniceshelf);

   See also: EXTRUDE, COLLAPSE

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function md2=modelextract(md1,area)
0002 %MODELEXTRACT - extract a model according to an Argus contour or flag list
0003 %
0004 %   This routine extracts a submodel from a bigger model with respect to a given contour
0005 %   md must be followed by the corresponding exp file or flags list
0006 %   It can either be a domain file (argus type, .exp extension), or an array of element flags.
0007 %   If user wants every element outside the domain to be
0008 %   modelextractd, add '~' to the name of the domain file (ex: '~Pattyn.exp');
0009 %   an empty string '' will be considered as an empty domain
0010 %   a string 'all' will be considered as the entire domain
0011 %
0012 %   Usage:
0013 %      md2=modelextract(md1,area);
0014 %
0015 %   Examples:
0016 %      md2=modelextract(md,'Domain.exp');
0017 %      md2=modelextract(md,md.elementoniceshelf);
0018 %
0019 %   See also: EXTRUDE, COLLAPSE
0020 
0021 %some checks on list of arguments
0022 if ((nargin~=2) | (nargout~=1)),
0023     error('modelextract error message');
0024 end
0025 package='Ice';
0026 
0027 %first
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         %make sure that we actually don't want the elements outside the domain outline!
0037         if strcmpi(area(1),'~'),
0038             area=area(2:length(area));
0039             invert=1;
0040         else
0041             invert=0;
0042         end
0043         %ok, flag_elem elements
0044         A=expread(area,1);
0045         flag_elem=ArgusContourToMesh(md1.elements(:,1:3),md1.x,md1.y,A,'element',1);
0046         %check that the elements do not cross the contour -> flag as 0
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 %get the flags of grids
0072 flag_grid=ismember([1:md1.numberofgrids]',sort(unique(md1.elements(pos_elem,:)))');
0073 pos_grid=find(flag_grid);
0074 
0075 %keep track of some fields
0076 numberofgrids1=md1.numberofgrids;
0077 numberofelements1=md1.numberofelements;
0078 numberofgrids2=length(pos_grid);
0079 numberofelements2=length(pos_elem);
0080 
0081 %Create Pelem and Pgrid (transform old grids in new grids and same thing for the elements)
0082 Pelem=zeros(numberofelements1,1);
0083 Pelem(pos_elem)=[1:numberofelements2]';
0084 Pgrid=zeros(numberofgrids1,1);
0085 Pgrid(pos_grid)=[1:numberofgrids2]';
0086 
0087 %renumber the elements (some grid won't exist anymore)
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 %OK, now create the new model !
0101 
0102     %take every fields from model
0103     md2=md1;
0104 
0105     %Mesh
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     %Initial 2d mesh
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     %Elements
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     %Nodes
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     %Penalties
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     %Projections
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     %Extrusion
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     %Physical parameters
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     %Geometrical parameters
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     %Boundary conditions
0190     md2.gridonboundary=md1.gridonboundary(pos_grid);
0191 
0192     %Diagnostic
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     %Thermal
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     %Transient
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     %Observations
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     %Output parameters
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     %Results fields
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 %Now take care of the boudaries (Dirichlet).
0261 
0262 %Catch the elements that have not been extracted
0263 orphans_elem=find(~flag_elem);
0264 orphans_grid=unique(md1.elements(orphans_elem,:))';
0265 
0266 %Figure out which grid are on the boundary between md2 and md1
0267 gridstoflag1=intersect(orphans_grid,pos_grid);
0268 gridstoflag2=Pgrid(gridstoflag1);
0269 
0270 %Same thing for the elements
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 %These grids are "on boundary" and hence "on dirichlet"
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 %Keep track of pos_grid and pos_elem
0306 md2.extractedgrids=pos_grid;
0307 md2.extractedelements=pos_elem;
0308 
0309 %get segments
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         %ok, [grid1,grid2] are a segment only if they belong to one element
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 %function
0334 
0335 function bool=cross_extract(x1,y1,x2,y2);
0336 %return 1 if the polygon defined by (x1,y1) crosses the polygon defined by (x2,y2)
0337     x1 = x1(:); y1 = y1(:);
0338     x2 = x2(:); y2 = y2(:);
0339     l1 = length(x1);
0340     l2 = length(x2);
0341 
0342     % Indexing
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     % Calculate matrix of intersections
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

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003