Changeset 13008
- Timestamp:
- 08/13/12 10:52:25 (13 years ago)
- Location:
- issm/trunk-jpl/src/m
- Files:
-
- 1 added
- 7 deleted
- 2 edited
- 15 moved
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/modellist.m
r12030 r13008 10 10 end 11 11 methods 12 function md_list=modelsextract(md,flags,minel,varargin) % {{{ 13 %modelsextract - extract several self contained models according to a list of element flags. 14 % 15 % The difference between this routine and the modelextract.m routine (without an 's') is that 16 % as many models are extracted as there are closed contours defined in area. 17 % This routine is needed for example when doing data assimilation of ice shelves in Antarctica. 18 % Many independent ice shelves are present, and we don't want data assimilation on one ice shelf 19 % to be hindered by another totally independent ice shelf. 20 % 21 % Usage: 22 % md_list=modelsextract(md,elementfalgs,minel); 23 % 24 % Examples: 25 % md_list=modelsextract(md,md.mask.elementonfloatingice,1000); 26 % 27 % See also: EXTRUDE, COLLAPSE, MODELEXTRACT 28 29 disp('selecting pools of elements'); 30 %go through flags and build as many independent element flags as there are groups of connected 1s 31 %in flags. 32 33 %2D or 3D? 34 if md.mesh.dimension==3, 35 numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out. 36 flags=project2d(md,flags,1); 37 else 38 numberofelements=md.mesh.numberofelements; 39 end 40 41 %recover extra arguments: 42 distance=0; 43 if nargin==4, 44 distance=varargin{1}; 45 end 46 47 flag_list=cell(0,1); 48 49 for i=1:size(flags,1), 50 51 if (flags(i)), 52 53 %ok, we are sure element i is part of a new pool. 54 pool=zeros(numberofelements,1); 55 pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags); 56 flag_list{end+1,1}=pool; 57 58 %speed up rest of computation by taking pool out of flags: 59 pos=find(pool);flags(pos)=0; 60 61 end 62 end 63 64 %go through flag_list and discard any pool of less than minel elements: 65 ex_pos=[]; 66 for i=1:length(flag_list), 67 if length(find(flag_list{i}))<minel, 68 ex_pos=[ex_pos; i]; 69 end 70 end 71 flag_list(ex_pos)=[]; 72 73 %now, if distance was specified, expand the flag_list by distance km: 74 if distance, 75 for i=1:length(flag_list), 76 flag_list{i}=PropagateFlagsUntilDistance(md,flag_list{i},distance); 77 end 78 end 79 80 %now, go use the pools of flags to extract models: 81 disp(['extracting ' num2str(size(flag_list,1)) ' models']); 82 models=cell(0,1); 83 84 for i=1:size(flag_list,1), 85 disp([' ' num2str(i) '/' num2str(size(flag_list,1))]); 86 if md.mesh.dimension==3, 87 flags2d=flag_list{i}; 88 realflags=project3d(md,flags2d,'element'); 89 else 90 realflags=flag_list{i}; 91 end 92 models{end+1,1}=modelextract(md,realflags); 93 end 94 95 %return model list 96 md_list=modellist(models); 97 98 end %end of this function }}} 99 function md_list=modelsextractfromdomains(md,directory) % {{{ 100 %modelsextractfromdomains- extract several self contained models according to a list of domains 101 % 102 % Usage: 103 % md_list=modelsextractfromdomains(md,'Basins/'); 104 % 105 % Examples: 106 % md_list=modelsextract(md,'Basins/'); 107 % 108 % See also: MODELSEXTRACTS, MODELEXTRACT 109 110 %go into directory and get list of files. 111 cd(directory); 112 basins=listfiles; 113 cd .. 114 115 models=cell(0,1); 116 for i=1:length(basins), 117 models{end+1,1}=modelextract(md,[directory '/' basins{i}]); 118 end 119 120 %return model list: 121 md_list=modellist(models); 122 123 end % }}} 12 124 function obj = modellist(varargin) % {{{ 13 125 -
issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m
r12996 r13008 7 7 % 8 8 9 10 11 9 new_flags=flags; 12 10 … … 24 22 sum_conn=sum(conn,2); 25 23 border_elements=flag_elements(find(sum_conn>=1)); 26 27 24 28 25 %average x and y over elements: … … 49 46 end 50 47 51 52 48 %check which of these new elements are more than distance away from the border elements 53 49 for i=1:length(new_elements), -
issm/trunk-jpl/src/m/model/mesh/bamg.m
r12365 r13008 339 339 error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)'); 340 340 end 341 end 342 343 function geom=processgeometry(geom,tol,outline); % {{{ 344 345 %Deal with edges 346 disp('Checking Edge crossing...'); 347 i=0; 348 while (i<size(geom.Edges,1)), 349 350 %edge counter 351 i=i+1; 352 353 %Get coordinates 354 x1=geom.Vertices(geom.Edges(i,1),1); 355 y1=geom.Vertices(geom.Edges(i,1),2); 356 x2=geom.Vertices(geom.Edges(i,2),1); 357 y2=geom.Vertices(geom.Edges(i,2),2); 358 color1=geom.Edges(i,3); 359 360 j=i; %test edges located AFTER i only 361 while (j<size(geom.Edges,1)), 362 363 %edge counter 364 j=j+1; 365 366 %Skip if the two edges already have a vertex in common 367 if any(ismember(geom.Edges(i,1:2),geom.Edges(j,1:2))), 368 continue 369 end 370 371 %Get coordinates 372 x3=geom.Vertices(geom.Edges(j,1),1); 373 y3=geom.Vertices(geom.Edges(j,1),2); 374 x4=geom.Vertices(geom.Edges(j,2),1); 375 y4=geom.Vertices(geom.Edges(j,2),2); 376 color2=geom.Edges(j,3); 377 378 %Check if the two edges are crossing one another 379 if SegIntersect([x1 y1; x2 y2],[x3 y3; x4 y4]), 380 381 %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) 382 x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 383 y=det([det([x1 y1; x2 y2]) y1-y2;det([x3 y3; x4 y4]) y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); 384 385 %Add vertex to the list of vertices 386 geom.Vertices(end+1,:)=[x y min(color1,color2)]; 387 id=size(geom.Vertices,1); 388 389 %Update edges i and j 390 edgei=geom.Edges(i,:); 391 edgej=geom.Edges(j,:); 392 geom.Edges(i,:) =[edgei(1) id edgei(3)]; 393 geom.Edges(end+1,:)=[id edgei(2) edgei(3)]; 394 geom.Edges(j,:) =[edgej(1) id edgej(3)]; 395 geom.Edges(end+1,:)=[id edgej(2) edgej(3)]; 396 397 %update current edge second tip 398 x2=x; y2=y; 399 end 400 end 401 402 end 403 404 %Check point outside 405 disp('Checking for points outside the domain...'); 406 i=0; 407 num=0; 408 while (i<size(geom.Vertices,1)), 409 410 %vertex counter 411 i=i+1; 412 413 %Get coordinates 414 x=geom.Vertices(i,1); 415 y=geom.Vertices(i,2); 416 color=geom.Vertices(i,3); 417 418 %Check that the point is inside the domain 419 if (color~=1 & ~ContourToNodes(x,y,outline(1),1)), 420 421 %Remove points from list of Vertices 422 num=num+1; 423 geom.Vertices(i,:)=[]; 424 425 %update edges 426 [posedges dummy]=find(geom.Edges==i); 427 geom.Edges(posedges,:)=[]; 428 posedges=find(geom.Edges>i); 429 geom.Edges(posedges)=geom.Edges(posedges)-1; 430 431 %update counter 432 i=i-1; 433 end 434 end 435 if num, 436 disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']); 437 end 438 439 %Check point spacing 440 if ~isnan(tol), 441 disp('Checking point spacing...'); 442 i=0; 443 while (i<size(geom.Vertices,1)), 444 445 %vertex counter 446 i=i+1; 447 448 %Get coordinates 449 x1=geom.Vertices(i,1); 450 y1=geom.Vertices(i,2); 451 452 j=i; %test edges located AFTER i only 453 while (j<size(geom.Vertices,1)), 454 455 %vertex counter 456 j=j+1; 457 458 %Get coordinates 459 x2=geom.Vertices(j,1); 460 y2=geom.Vertices(j,2); 461 462 %Check whether the two vertices are too close 463 if ((x2-x1)^2+(y2-y1)^2<tol^2) 464 465 %Remove points from list of Vertices 466 geom.Vertices(j,:)=[]; 467 468 %update edges 469 posedges=find(ismember(geom.Edges,j)); 470 geom.Edges(posedges)=i; 471 posedges=find(geom.Edges>j); 472 geom.Edges(posedges)=geom.Edges(posedges)-1; 473 474 %update counter 475 j=j-1; 476 477 end 478 end 479 end 480 end 481 %remove empty edges 482 geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[]; 483 end % }}}
Note:
See TracChangeset
for help on using the changeset viewer.