Changeset 13692
- Timestamp:
- 10/16/12 09:56:05 (12 years ago)
- Location:
- issm/trunk-jpl/src/m/classes/model
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/model/model.m
r13646 r13692 5 5 6 6 classdef model 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 7 properties (SetAccess=public) %Model fields 8 % {{{ 9 %Careful here: no other class should be used as default value this is a bug of matlab 10 mesh = 0; 11 mask = 0; 12 13 geometry = 0; 14 constants = 0; 15 surfaceforcings = 0; 16 basalforcings = 0; 17 materials = 0; 18 friction = 0; 19 flowequation = 0; 20 timestepping = 0; 21 initialization = 0; 22 rifts = 0; 23 24 debug = 0; 25 verbose = 0; 26 settings = 0; 27 solver = 0; 28 cluster = 0; 29 30 balancethickness = 0; 31 diagnostic = 0; 32 groundingline = 0; 33 hydrology = 0; 34 prognostic = 0; 35 thermal = 0; 36 steadystate = 0; 37 transient = 0; 38 39 autodiff = 0; 40 flaim = 0; 41 inversion = 0; 42 qmu = 0; 43 44 results = 0; 45 radaroverlay = 0; 46 miscellaneous = 0; 47 private = 0; 48 49 %}}} 50 end 51 methods (Static) 52 function md = loadobj(md) % {{{ 53 % This function is directly called by matlab when a model object is 54 % loaded. If the input is a struct it is an old version of model and 55 % old fields must be recovered (make sure they are in the deprecated 56 % model properties) 57 58 if verLessThan('matlab','7.9'), 59 disp('Warning: your matlab version is old and there is a risk that load does not work correctly'); 60 disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it'); 61 62 % This is a Matlab bug: all the fields of md have their default value 63 % Example of error message: 64 % Warning: Error loading an object of class 'model': 65 % Undefined function or method 'exist' for input arguments of type 'cell' 66 % 67 % This has been fixed in MATLAB 7.9 (R2009b) and later versions 68 end 69 70 if isstruct(md) 71 disp('Recovering model object from a previous version'); 72 md = structtomodel(model,md); 73 end 74 75 %2012 August 4th 76 if isa(md.materials,'materials'), 77 disp('Recovering old materials'); 78 if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z), 79 md.materials=matice(md.materials); 80 else 81 md.materials=matdamageice(md.materials); 82 end 83 end 84 85 end% }}} 86 end 87 methods 88 function md = model(varargin) % {{{ 89 90 switch nargin 91 case 0 92 md=setdefaultparameters(md); 93 otherwise 94 error('model constructor error message: 0 of 1 argument only in input.'); 95 end 96 end 97 %}}} 98 function md = checkmessage(md,string) % {{{ 99 if(nargout~=1) error('wrong usage, model must be an output'); end 100 disp(['model not consistent: ' string]); 101 md.private.isconsistent=false; 102 end 103 %}}} 104 function md = collapse(md)% {{{ 105 %COLLAPSE - collapses a 3d mesh into a 2d mesh 106 % 107 % This routine collapses a 3d model into a 2d model 108 % and collapses all the fileds of the 3d model by 109 % taking their depth-averaged values 110 % 111 % Usage: 112 % md=collapse(md) 113 % 114 % See also: EXTRUDE, MODELEXTRACT 115 116 %Check that the model is really a 3d model 117 if ~md.mesh.dimension==3, 118 error('collapse error message: only 3d mesh can be collapsed') 119 end 120 121 %Start with changing alle the fields from the 3d mesh 122 123 %drag is limited to nodes that are on the bedrock. 124 md.friction.coefficient=project2d(md,md.friction.coefficient,1); 125 126 %p and q (same deal, except for element that are on the bedrock: ) 127 md.friction.p=project2d(md,md.friction.p,1); 128 md.friction.q=project2d(md,md.friction.q,1); 129 130 %observations 131 if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; 132 if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; 133 if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; 134 if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; 135 if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; 136 if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; 137 if ~isnan(md.surfaceforcings.mass_balance), 138 md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); 139 end; 140 if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; 141 142 %results 143 if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; 144 if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; 145 if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; 146 if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; 147 if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; 148 149 %bedinfo and surface info 150 md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); 151 md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); 152 md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); 153 md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); 154 155 %elementstype 156 if ~isnan(md.flowequation.element_equation) 157 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); 158 md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); 159 md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); 160 md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); 161 md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); 162 end 163 164 %boundary conditions 165 md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); 166 md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); 167 md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); 168 md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); 169 md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); 170 md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); 171 172 %Extrusion of Neumann BC 173 if ~isnan(md.diagnostic.icefront), 174 numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); 175 md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer 176 end 177 178 %materials 179 md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); 180 md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); 181 if isa(md.materials,'matdamageice') 182 md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); 183 end 184 185 %special for thermal modeling: 186 md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); 187 md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux 188 189 %update of connectivity matrix 190 md.mesh.average_vertex_connectivity=25; 191 192 %Collapse the mesh 193 nodes2d=md.mesh.numberofvertices2d; 194 elements2d=md.mesh.numberofelements2d; 195 196 %parameters 197 md.geometry.surface=project2d(md,md.geometry.surface,1); 198 md.geometry.thickness=project2d(md,md.geometry.thickness,1); 199 md.geometry.bed=project2d(md,md.geometry.bed,1); 200 md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); 201 md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); 202 md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); 203 md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); 204 md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); 205 md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); 206 md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); 207 md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); 208 md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); 209 210 %lat long 211 if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end 212 if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end 213 214 %Initialize with the 2d mesh 215 md.mesh.x=md.mesh.x2d; 216 md.mesh.y=md.mesh.y2d; 217 md.mesh.z=zeros(size(md.mesh.x2d)); 218 md.mesh.numberofvertices=md.mesh.numberofvertices2d; 219 md.mesh.numberofelements=md.mesh.numberofelements2d; 220 md.mesh.elements=md.mesh.elements2d; 221 222 %Keep a trace of lower and upper nodes 223 md.mesh.lowervertex=NaN; 224 md.mesh.uppervertex=NaN; 225 md.mesh.lowerelements=NaN; 226 md.mesh.upperelements=NaN; 227 228 %Remove old mesh 229 md.mesh.x2d=NaN; 230 md.mesh.y2d=NaN; 231 md.mesh.elements2d=NaN; 232 md.mesh.numberofelements2d=md.mesh.numberofelements; 233 md.mesh.numberofvertices2d=md.mesh.numberofvertices; 234 md.mesh.numberoflayers=0; 235 236 %Update mesh type 237 md.mesh.dimension=2; 238 end % }}} 239 function md2 = extract(md,area) % {{{ 240 %extract - extract a model according to an Argus contour or flag list 241 % 242 % This routine extracts a submodel from a bigger model with respect to a given contour 243 % md must be followed by the corresponding exp file or flags list 244 % It can either be a domain file (argus type, .exp extension), or an array of element flags. 245 % If user wants every element outside the domain to be 246 % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); 247 % an empty string '' will be considered as an empty domain 248 % a string 'all' will be considered as the entire domain 249 % add an argument 0 if you do not want the elements to be checked (faster) 250 % 251 % Usage: 252 % md2=extract(md,area); 253 % 254 % Examples: 255 % md2=extract(md,'Domain.exp'); 256 % md2=extract(md,md.mask.elementonfloatingice); 257 % 258 % See also: EXTRUDE, COLLAPSE 259 260 %copy model 261 md1=md; 262 263 %some checks 264 if ((nargin~=2) | (nargout~=1)), 265 help extract 266 error('extract error message: bad usage'); 267 end 268 269 %get check option 270 if (nargin==3 & varargin{1}==0), 271 checkoutline=0; 272 else 273 checkoutline=1; 274 end 275 276 %get elements that are inside area 277 flag_elem=FlagElements(md1,area); 278 if ~any(flag_elem), 279 error('extracted model is empty'); 280 end 281 282 %kick out all elements with 3 dirichlets 283 spc_elem=find(~flag_elem); 284 spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); 285 flag=ones(md1.mesh.numberofvertices,1); 286 flag(spc_node)=0; 287 pos=find(sum(flag(md1.mesh.elements),2)==0); 288 flag_elem(pos)=0; 289 290 %extracted elements and nodes lists 291 pos_elem=find(flag_elem); 292 pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); 293 294 %keep track of some fields 295 numberofvertices1=md1.mesh.numberofvertices; 296 numberofelements1=md1.mesh.numberofelements; 297 numberofvertices2=length(pos_node); 298 numberofelements2=length(pos_elem); 299 flag_node=zeros(numberofvertices1,1); 300 flag_node(pos_node)=1; 301 302 %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) 303 Pelem=zeros(numberofelements1,1); 304 Pelem(pos_elem)=[1:numberofelements2]'; 305 Pnode=zeros(numberofvertices1,1); 306 Pnode(pos_node)=[1:numberofvertices2]'; 307 308 %renumber the elements (some node won't exist anymore) 309 elements_1=md1.mesh.elements; 310 elements_2=elements_1(pos_elem,:); 311 elements_2(:,1)=Pnode(elements_2(:,1)); 312 elements_2(:,2)=Pnode(elements_2(:,2)); 313 elements_2(:,3)=Pnode(elements_2(:,3)); 314 if md1.mesh.dimension==3, 315 elements_2(:,4)=Pnode(elements_2(:,4)); 316 elements_2(:,5)=Pnode(elements_2(:,5)); 317 elements_2(:,6)=Pnode(elements_2(:,6)); 318 end 319 320 %OK, now create the new model ! 321 322 %take every fields from model 323 md2=md1; 324 325 %automatically modify fields 326 327 %loop over model fields 328 model_fields=fields(md1); 329 for i=1:length(model_fields), 330 %get field 331 field=md1.(model_fields{i}); 332 fieldsize=size(field); 333 if isobject(field), %recursive call 334 object_fields=fields(md1.(model_fields{i})); 335 for j=1:length(object_fields), 336 %get field 337 field=md1.(model_fields{i}).(object_fields{j}); 338 fieldsize=size(field); 339 %size = number of nodes * n 340 if fieldsize(1)==numberofvertices1 341 md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); 342 elseif (fieldsize(1)==numberofvertices1+1) 343 md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; 344 %size = number of elements * n 345 elseif fieldsize(1)==numberofelements1 346 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); 347 end 348 end 349 else 350 %size = number of nodes * n 351 if fieldsize(1)==numberofvertices1 352 md2.(model_fields{i})=field(pos_node,:); 353 elseif (fieldsize(1)==numberofvertices1+1) 354 md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; 355 %size = number of elements * n 356 elseif fieldsize(1)==numberofelements1 357 md2.(model_fields{i})=field(pos_elem,:); 358 end 359 end 360 end 361 362 %modify some specific fields 363 364 %Mesh 365 md2.mesh.numberofelements=numberofelements2; 366 md2.mesh.numberofvertices=numberofvertices2; 367 md2.mesh.elements=elements_2; 368 369 %mesh.uppervertex mesh.lowervertex 370 if md1.mesh.dimension==3 371 md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); 372 pos=find(~isnan(md2.mesh.uppervertex)); 373 md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); 374 375 md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); 376 pos=find(~isnan(md2.mesh.lowervertex)); 377 md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); 378 379 md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); 380 pos=find(~isnan(md2.mesh.upperelements)); 381 md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); 382 383 md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); 384 pos=find(~isnan(md2.mesh.lowerelements)); 385 md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); 386 end 387 388 %Initial 2d mesh 389 if md1.mesh.dimension==3 390 flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); 391 pos_elem_2d=find(flag_elem_2d); 392 flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); 393 pos_node_2d=find(flag_node_2d); 394 395 md2.mesh.numberofelements2d=length(pos_elem_2d); 396 md2.mesh.numberofvertices2d=length(pos_node_2d); 397 md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); 398 md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); 399 md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); 400 md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); 401 402 md2.mesh.x2d=md1.mesh.x(pos_node_2d); 403 md2.mesh.y2d=md1.mesh.y(pos_node_2d); 404 end 405 406 %Edges 407 if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... 408 %renumber first two columns 409 pos=find(md2.mesh.edges(:,4)~=-1); 410 md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); 411 md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); 412 md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); 413 md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); 414 %remove edges when the 2 vertices are not in the domain. 415 md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); 416 %Replace all zeros by -1 in the last two columns; 417 pos=find(md2.mesh.edges(:,3)==0); 418 md2.mesh.edges(pos,3)=-1; 419 pos=find(md2.mesh.edges(:,4)==0); 420 md2.mesh.edges(pos,4)=-1; 421 %Invert -1 on the third column with last column (Also invert first two columns!!) 422 pos=find(md2.mesh.edges(:,3)==-1); 423 md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); 424 md2.mesh.edges(pos,4)=-1; 425 values=md2.mesh.edges(pos,2); 426 md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); 427 md2.mesh.edges(pos,1)=values; 428 %Finally remove edges that do not belong to any element 429 pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); 430 md2.mesh.edges(pos,:)=[]; 431 end 432 433 %Penalties 434 if ~isnan(md2.diagnostic.vertex_pairing), 435 for i=1:size(md1.diagnostic.vertex_pairing,1); 436 md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); 437 end 438 md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); 439 end 440 if ~isnan(md2.prognostic.vertex_pairing), 441 for i=1:size(md1.prognostic.vertex_pairing,1); 442 md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); 443 end 444 md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); 445 end 446 447 %recreate segments 448 if md1.mesh.dimension==2 449 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 450 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 451 md2.mesh.segments=contourenvelope(md2); 452 md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 453 else 454 %First do the connectivity for the contourenvelope in 2d 455 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); 456 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); 457 md2.mesh.segments=contourenvelope(md2); 458 md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 459 md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); 460 %Then do it for 3d as usual 461 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 462 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 463 end 464 465 %Boundary conditions: Dirichlets on new boundary 466 %Catch the elements that have not been extracted 467 orphans_elem=find(~flag_elem); 468 orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; 469 %Figure out which node are on the boundary between md2 and md1 470 nodestoflag1=intersect(orphans_node,pos_node); 471 nodestoflag2=Pnode(nodestoflag1); 472 if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, 473 if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 474 md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 475 md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); 476 else 477 md2.diagnostic.spcvx(nodestoflag2)=NaN; 478 md2.diagnostic.spcvy(nodestoflag2)=NaN; 479 disp(' ') 480 disp('!! extract warning: spc values should be checked !!') 481 disp(' ') 482 end 483 %put 0 for vz 484 md2.diagnostic.spcvz(nodestoflag2)=0; 485 end 486 if ~isnan(md1.thermal.spctemperature), 487 md2.thermal.spctemperature(nodestoflag2,1)=1; 488 end 489 490 %Diagnostic 491 if ~isnan(md2.diagnostic.icefront) 492 md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); 493 md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); 494 md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); 495 if md1.mesh.dimension==3 496 md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); 497 md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); 498 end 499 md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); 500 end 501 502 %Results fields 503 if isstruct(md1.results), 504 md2.results=struct(); 505 solutionfields=fields(md1.results); 506 for i=1:length(solutionfields), 507 %get subfields 508 solutionsubfields=fields(md1.results.(solutionfields{i})); 509 for j=1:length(solutionsubfields), 510 field=md1.results.(solutionfields{i}).(solutionsubfields{j}); 511 if length(field)==numberofvertices1, 512 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); 513 elseif length(field)==numberofelements1, 514 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); 515 else 516 md2.results.(solutionfields{i}).(solutionsubfields{j})=field; 517 end 518 end 519 end 520 end 521 522 %Keep track of pos_node and pos_elem 523 md2.mesh.extractedvertices=pos_node; 524 md2.mesh.extractedelements=pos_elem; 525 end % }}} 526 function md = extrude(md,varargin) % {{{ 527 %EXTRUDE - vertically extrude a 2d mesh 528 % 529 % vertically extrude a 2d mesh and create corresponding 3d mesh. 530 % The vertical distribution can: 531 % - follow a polynomial law 532 % - follow two polynomial laws, one for the lower part and one for the upper part of the mesh 533 % - be discribed by a list of coefficients (between 0 and 1) 534 % 535 % 536 % Usage: 537 % md=extrude(md,numlayers,extrusionexponent); 538 % md=extrude(md,numlayers,lowerexponent,upperexponent); 539 % md=extrude(md,listofcoefficients); 540 % 541 % Example: 542 % md=extrude(md,8,3); 543 % md=extrude(md,8,3,2); 544 % md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); 545 % 546 % See also: MODELEXTRACT, COLLAPSE 547 548 %some checks on list of arguments 549 if ((nargin>4) | (nargin<2) | (nargout~=1)), 550 help extrude; 551 error('extrude error message'); 552 end 553 554 %Extrude the mesh 555 if nargin==2, %list of coefficients 556 clist=varargin{1}; 557 if any(clist<0) | any(clist>1), 558 error('extrusioncoefficients must be between 0 and 1'); 559 end 560 extrusionlist=sort(unique([clist(:);0;1])); 561 numlayers=length(extrusionlist); 562 elseif nargin==3, %one polynomial law 563 if varargin{2}<=0, 564 help extrude; 565 error('extrusionexponent must be >=0'); 566 end 567 numlayers=varargin{1}; 568 extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; 569 elseif nargin==4, %two polynomial laws 570 numlayers=varargin{1}; 571 lowerexp=varargin{2}; 572 upperexp=varargin{3}; 573 574 if varargin{2}<=0 | varargin{3}<=0, 575 help extrude; 576 error('lower and upper extrusionexponents must be >=0'); 577 end 578 579 lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; 580 upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; 581 extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); 582 583 end 584 585 if numlayers<2, 586 error('number of layers should be at least 2'); 587 end 588 if md.mesh.dimension==3, 589 error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); 590 end 591 592 %Initialize with the 2d mesh 593 x3d=[]; 594 y3d=[]; 595 z3d=[]; %the lower node is on the bed 596 thickness3d=md.geometry.thickness; %thickness and bed for these nodes 597 bed3d=md.geometry.bed; 598 599 %Create the new layers 600 for i=1:numlayers, 601 x3d=[x3d; md.mesh.x]; 602 y3d=[y3d; md.mesh.y]; 603 %nodes are distributed between bed and surface accordingly to the given exponent 604 z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; 605 end 606 number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh 607 608 %Extrude elements 609 elements3d=[]; 610 for i=1:numlayers-1, 611 elements3d=[elements3d;[md.mesh.elements+(i-1)*md.mesh.numberofvertices md.mesh.elements+i*md.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part 612 end 613 number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh 614 615 %Keep a trace of lower and upper nodes 616 mesh.lowervertex=NaN*ones(number_nodes3d,1); 617 mesh.uppervertex=NaN*ones(number_nodes3d,1); 618 mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; 619 mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; 620 md.mesh.lowervertex=mesh.lowervertex; 621 md.mesh.uppervertex=mesh.uppervertex; 622 623 %same for lower and upper elements 624 mesh.lowerelements=NaN*ones(number_el3d,1); 625 mesh.upperelements=NaN*ones(number_el3d,1); 626 mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; 627 mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; 628 md.mesh.lowerelements=mesh.lowerelements; 629 md.mesh.upperelements=mesh.upperelements; 630 631 %Save old mesh 632 md.mesh.x2d=md.mesh.x; 633 md.mesh.y2d=md.mesh.y; 634 md.mesh.elements2d=md.mesh.elements; 635 md.mesh.numberofelements2d=md.mesh.numberofelements; 636 md.mesh.numberofvertices2d=md.mesh.numberofvertices; 637 638 %Update mesh type 639 md.mesh.dimension=3; 640 641 %Build global 3d mesh 642 md.mesh.elements=elements3d; 643 md.mesh.x=x3d; 644 md.mesh.y=y3d; 645 md.mesh.z=z3d; 646 md.mesh.numberofelements=number_el3d; 647 md.mesh.numberofvertices=number_nodes3d; 648 md.mesh.numberoflayers=numlayers; 649 650 %Ok, now deal with the other fields from the 2d mesh: 651 652 %lat long 653 md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); 654 md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); 655 656 %drag coefficient is limited to nodes that are on the bedrock. 657 md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1); 658 659 %p and q (same deal, except for element that are on the bedrock: ) 660 md.friction.p=project3d(md,'vector',md.friction.p,'type','element'); 661 md.friction.q=project3d(md,'vector',md.friction.q,'type','element'); 662 663 %observations 664 md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node'); 665 md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node'); 666 md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node'); 667 md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node'); 668 md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node'); 669 md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node'); 670 md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node'); 671 672 %results 673 if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; 674 if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; 675 if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; 676 if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; 677 if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; 678 if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; 679 680 %bedinfo and surface info 681 md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); 682 md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); 683 md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); 684 md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); 685 686 %elementstype 687 if ~isnan(md.flowequation.element_equation) 688 oldelements_type=md.flowequation.element_equation; 689 md.flowequation.element_equation=zeros(number_el3d,1); 690 md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); 691 end 692 693 %verticestype 694 if ~isnan(md.flowequation.vertex_equation) 695 oldvertices_type=md.flowequation.vertex_equation; 696 md.flowequation.vertex_equation=zeros(number_nodes3d,1); 697 md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); 698 end 699 md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node'); 700 md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node'); 701 md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node'); 702 703 %boundary conditions 704 md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node'); 705 md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node'); 706 md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node'); 707 md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); 708 md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node'); 709 md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node'); 710 md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node'); 711 712 %in 3d, pressureload: [node1 node2 node3 node4 element] 713 pressureload_layer1=[md.diagnostic.icefront(:,1:2) md.diagnostic.icefront(:,2)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,1)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,3:4)]; %Add two columns on the first layer 714 pressureload=[]; 715 for i=1:numlayers-1, 716 pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; 717 end 718 md.diagnostic.icefront=pressureload; 719 720 %connectivity 721 md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); 722 md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; 723 for i=2:numlayers-1, 724 md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... 725 =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; 726 end 727 md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0; 728 729 %materials 730 md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); 731 md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); 732 if isa(md.materials,'matdamageice') 733 md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); 734 end 735 736 %parameters 737 md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node'); 738 md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node'); 739 md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node'); 740 md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node'); 741 md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node'); 742 md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); 743 md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element'); 744 md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node'); 745 md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element'); 746 md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node'); 747 md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element'); 748 md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'type','node'); 749 if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; 750 if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; 751 if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; 752 if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end 753 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end 754 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end 755 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end 756 757 %Put lithostatic pressure if there is an existing pressure 758 if ~isnan(md.initialization.pressure), 759 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); 760 end 761 762 %special for thermal modeling: 763 md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); 764 if ~isnan(md.basalforcings.geothermalflux) 765 md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux 766 end 767 768 %increase connectivity if less than 25: 769 if md.mesh.average_vertex_connectivity<=25, 770 md.mesh.average_vertex_connectivity=100; 771 end 772 772 end % }}} 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 773 function md = structtomodel(md,structmd) % {{{ 774 775 if ~isstruct(structmd) error('input model is not a structure'); end 776 777 %loaded model is a struct, initialize output and recover all fields 778 md = structtoobj(model,structmd); 779 780 %Old field now classes 781 if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end 782 if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end 783 784 %Field name change 785 if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end 786 if isfield(structmd,'p'), md.friction.p=structmd.p; end 787 if isfield(structmd,'q'), md.friction.q=structmd.p; end 788 if isfield(structmd,'melting'), md.basalforcings.melting_rate=structmd.melting; end 789 if isfield(structmd,'melting_rate'), md.basalforcings.melting_rate=structmd.melting_rate; end 790 if isfield(structmd,'accumulation'), md.surfaceforcings.mass_balance=structmd.accumulation; end 791 if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end 792 if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end 793 if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end 794 if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end 795 if isfield(structmd,'gridonbed'), md.mesh.vertexonbed=structmd.gridonbed; end 796 if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end 797 if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end 798 if isfield(structmd,'gridoniceshelf'), md.mask.vertexonfloatingice=structmd.gridoniceshelf; end 799 if isfield(structmd,'gridonicesheet'), md.mask.vertexongroundedice=structmd.gridonicesheet; end 800 if isfield(structmd,'gridonwater'), md.mask.vertexonwater=structmd.gridonwater; end 801 if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end 802 if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.solver=structmd.petscoptions; end 803 if isfield(structmd,'g'), md.constants.g=structmd.g; end 804 if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end 805 if isfield(structmd,'surface_mass_balance'), md.surfaceforcings.mass_balance=structmd.surface_mass_balance; end 806 if isfield(structmd,'basal_melting_rate'), md.basalforcings.melting_rate=structmd.basal_melting_rate; end 807 if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end 808 if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end 809 if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end 810 if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end 811 if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end 812 if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end 813 if isfield(structmd,'riftproperties'), %old implementation 814 md.rifts=rifts(); 815 md.rifts.riftproperties=structmd.riftproperties; 816 md.rifts.riftstruct=structmd.rifts; 817 md.rifts.riftproperties=structmd.riftinfo; 818 end 819 if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end 820 if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end 821 if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end 822 if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end 823 if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end 824 if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end 825 if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end 826 if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end 827 if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end 828 if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end 829 if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end 830 if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end 831 if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end 832 if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end 833 if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end 834 if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end 835 if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end 836 if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end 837 if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end 838 if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end 839 if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end 840 if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end 841 if isfield(structmd,'spcthickness'), md.prognostic.spcthickness=structmd.spcthickness; end 842 if isfield(structmd,'artificial_diffusivity'), md.prognostic.stabilization=structmd.artificial_diffusivity; end 843 if isfield(structmd,'hydrostatic_adjustment'), md.prognostic.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end 844 if isfield(structmd,'penalties'), md.prognostic.vertex_pairing=structmd.penalties; end 845 if isfield(structmd,'penalty_offset'), md.prognostic.penalty_factor=structmd.penalty_offset; end 846 if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end 847 if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end 848 if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end 849 if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end 850 if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end 851 if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end 852 if isfield(structmd,'elementonwater'), md.mask.elementonwater=structmd.elementonwater; end 853 if isfield(structmd,'nodeoniceshelf'), md.mask.vertexonfloatingice=structmd.nodeoniceshelf; end 854 if isfield(structmd,'nodeonicesheet'), md.mask.vertexongroundedice=structmd.nodeonicesheet; end 855 if isfield(structmd,'nodeonwater'), md.mask.vertexonwater=structmd.nodeonwater; end 856 if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end 857 if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end 858 if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end 859 if isfield(structmd,'ismacayealpattyn'), md.flowequation.ismacayealpattyn=structmd.ismacayealpattyn; end 860 if isfield(structmd,'ishutter'), md.flowequation.ishutter=structmd.ishutter; end 861 if isfield(structmd,'isstokes'), md.flowequation.isstokes=structmd.isstokes; end 862 if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end 863 if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end 864 if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end 865 if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end 866 if isfield(structmd,'isdiagnostic'), md.transient.isdiagnostic=structmd.isdiagnostic; end 867 if isfield(structmd,'isprognostic'), md.transient.isprognostic=structmd.isprognostic; end 868 if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end 869 if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end 870 if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end 871 if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end 872 if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end 873 if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end 874 if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end 875 if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end 876 if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end 877 if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end 878 if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end 879 if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end 880 if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end 881 if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end 882 if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end 883 if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end 884 if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end 885 if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end 886 if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end 887 if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end 888 if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end 889 if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end 890 if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end 891 if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end 892 if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end 893 if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end 894 if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end 895 if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end 896 if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end 897 if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end 898 if isfield(structmd,'long'), md.mesh.long=structmd.long; end 899 if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end 900 if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end 901 if isfield(structmd,'dim'), md.mesh.dimension=structmd.dim; end 902 if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end 903 if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end 904 if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end 905 if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end 906 if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end 907 if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end 908 if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end 909 if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end 910 if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end 911 if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end 912 if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end 913 if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end 914 if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end 915 if isfield(structmd,'elementonbed'), md.mesh.elementonbed=structmd.elementonbed; end 916 if isfield(structmd,'elementonsurface'), md.mesh.elementonsurface=structmd.elementonsurface; end 917 if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end 918 if isfield(structmd,'nodeonbed'), md.mesh.vertexonbed=structmd.nodeonbed; end 919 if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end 920 if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end 921 if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end 922 if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end 923 if isfield(structmd,'edges'), md.mesh.edges=structmd.edges; end 924 if isfield(structmd,'y'), md.mesh.y=structmd.y; end 925 if isfield(structmd,'x'), md.mesh.x=structmd.x; end 926 if isfield(structmd,'z'), md.mesh.z=structmd.z; end 927 if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end 928 if isfield(structmd,'pressureload'), md.diagnostic.icefront=structmd.pressureload; end 929 if isfield(structmd,'diagnostic_ref'), md.diagnostic.referential=structmd.diagnostic_ref; end 930 if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end 931 if isfield(structmd,'part'); md.qmu.partition=structmd.part; end 932 933 %Field changes 934 if (isfield(structmd,'type') & ischar(structmd.type)), 935 if strcmpi(structmd.type,'2d'), md.mesh.dimension=2; end 936 if strcmpi(structmd.type,'3d'), md.mesh.dimension=3; end 937 end 938 if isnumeric(md.verbose), 939 md.verbose=verbose; 940 end 941 if size(md.diagnostic.icefront,2)==3 || size(md.diagnostic.icefront,2)==5, 942 front=md.diagnostic.icefront; 943 md.diagnostic.icefront=[front 1*md.mask.elementonfloatingice(front(:,end))]; 944 end 945 if isfield(structmd,'spcvelocity'), 946 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); 947 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); 948 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); 949 pos=find(structmd.spcvelocity(:,1)); md.diagnostic.spcvx(pos)=structmd.spcvelocity(pos,4); 950 pos=find(structmd.spcvelocity(:,2)); md.diagnostic.spcvy(pos)=structmd.spcvelocity(pos,5); 951 pos=find(structmd.spcvelocity(:,3)); md.diagnostic.spcvz(pos)=structmd.spcvelocity(pos,6); 952 end 953 if isfield(structmd,'spcvx'), 954 md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); 955 pos=find(~isnan(structmd.spcvx)); md.diagnostic.spcvx(pos)=structmd.spcvx(pos); 956 end 957 if isfield(structmd,'spcvy'), 958 md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); 959 pos=find(~isnan(structmd.spcvy)); md.diagnostic.spcvy(pos)=structmd.spcvy(pos); 960 end 961 if isfield(structmd,'spcvz'), 962 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); 963 pos=find(~isnan(structmd.spcvz)); md.diagnostic.spcvz(pos)=structmd.spcvz(pos); 964 end 965 if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]), 966 pos=find(structmd.pressureload(:,end)==120); md.diagnostic.icefront(pos,end)=0; 967 pos=find(structmd.pressureload(:,end)==118); md.diagnostic.icefront(pos,end)=1; 968 pos=find(structmd.pressureload(:,end)==119); md.diagnostic.icefront(pos,end)=2; 969 end 970 if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50, 971 pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0; 972 pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1; 973 pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2; 974 pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3; 975 pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4; 976 pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5; 977 pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6; 978 pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7; 979 end 980 if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50, 981 pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0; 982 pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1; 983 pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2; 984 pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3; 985 pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4; 986 pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5; 987 pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6; 988 pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7; 989 end 990 if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law), 991 if (structmd.rheology_law==272), md.materials.rheology_law='None'; end 992 if (structmd.rheology_law==368), md.materials.rheology_law='Paterson'; end 993 if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end 994 end 995 if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration), 996 if (structmd.groundingline_migration==272), md.groundingline.migration='None'; end 997 if (structmd.groundingline_migration==273), md.groundingline.migration='AgressiveMigration'; end 998 if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end 999 end 1000 if isfield(structmd,'control_type') & isnumeric(structmd.control_type), 1001 if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end 1002 if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end 1003 if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end 1004 end 1005 if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]), 1006 pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101; 1007 pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102; 1008 pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103; 1009 pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104; 1010 pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105; 1011 pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201; 1012 pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501; 1013 pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502; 1014 pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503; 1015 end 1016 1017 if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2, 1018 md.thermal.stabilization=2; 1019 md.prognostic.stabilization=1; 1020 md.balancethickness.stabilization=1; 1021 end 1022 if isnumeric(md.prognostic.hydrostatic_adjustment) 1023 if md.prognostic.hydrostatic_adjustment==269, 1024 md.prognostic.hydrostatic_adjustment='Incremental'; 1025 else 1026 md.prognostic.hydrostatic_adjustment='Absolute'; 1027 end 1028 end 1029 1030 %New fields 1031 if ~isfield(structmd,'upperelements'); 1032 md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d; 1033 md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN; 1034 end 1035 if ~isfield(structmd,'lowerelements'); 1036 md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d; 1037 md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN; 1038 end 1039 1040 if ~isfield(structmd,'diagnostic_ref'); 1041 md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); 1042 end 1043 1044 end% }}} 1045 function md = setdefaultparameters(md) % {{{ 1046 1047 %initialize subclasses 1048 md.mesh = mesh(); 1049 md.mask = mask(); 1050 md.constants = constants(); 1051 md.geometry = geometry(); 1052 md.initialization = initialization(); 1053 md.surfaceforcings = surfaceforcings(); 1054 md.basalforcings = basalforcings(); 1055 md.friction = friction(); 1056 md.rifts = rifts(); 1057 md.timestepping = timestepping(); 1058 md.groundingline = groundingline(); 1059 md.materials = matice(); 1060 md.flowequation = flowequation(); 1061 md.debug = debug(); 1062 md.verbose = verbose('solution',true,'qmu',true,'control',true); 1063 md.settings = settings(); 1064 md.solver = solver(); 1065 if ismumps(), 1066 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),mumpsoptions()); 1067 else 1068 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),iluasmoptions()); 1069 end 1070 md.cluster = generic(); 1071 md.balancethickness = balancethickness(); 1072 md.diagnostic = diagnostic(); 1073 md.hydrology = hydrology(); 1074 md.prognostic = prognostic(); 1075 md.thermal = thermal(); 1076 md.steadystate = steadystate(); 1077 md.transient = transient(); 1078 md.autodiff = autodiff(); 1079 md.flaim = flaim(); 1080 md.inversion = inversion(); 1081 md.qmu = qmu(); 1082 md.radaroverlay = radaroverlay(); 1083 md.results = struct(); 1084 md.miscellaneous = miscellaneous(); 1085 md.private = private(); 1086 end 1087 %}}} 1088 function disp(obj) % {{{ 1089 disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(obj.mesh) ']'],'mesh properties')); 1090 disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(obj.mask) ']'],'defines grounded and floating elements')); 1091 disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(obj.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); 1092 disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(obj.constants) ']'],'physical constants')); 1093 disp(sprintf('%19s: %-22s -- %s','surfaceforcings' ,['[1x1 ' class(obj.surfaceforcings) ']'],'surface forcings')); 1094 disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings')); 1095 disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(obj.materials) ']'],'material properties')); 1096 disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties')); 1097 disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(obj.flowequation) ']'],'flow equations')); 1098 disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(obj.timestepping) ']'],'time stepping for transient models')); 1099 disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(obj.initialization) ']'],'initial guess/state')); 1100 disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(obj.rifts) ']'],'rifts properties')); 1101 disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(obj.debug) ']'],'debugging tools (valgrind, gprof)')); 1102 disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve')); 1103 disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(obj.settings) ']'],'settings properties')); 1104 disp(sprintf('%19s: %-22s -- %s','solver' ,['[1x1 ' class(obj.solver) ']'],'PETSc options for each solution')); 1105 disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)')); 1106 disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution')); 1107 disp(sprintf('%19s: %-22s -- %s','diagnostic' ,['[1x1 ' class(obj.diagnostic) ']'],'parameters for diagnostic solution')); 1108 disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution')); 1109 disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution')); 1110 disp(sprintf('%19s: %-22s -- %s','prognostic' ,['[1x1 ' class(obj.prognostic) ']'],'parameters for prognostic solution')); 1111 disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution')); 1112 disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution')); 1113 disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution')); 1114 disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters')); 1115 disp(sprintf('%19s: %-22s -- %s','flaim' ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters')); 1116 disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(obj.inversion) ']'],'parameters for inverse methods')); 1117 disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(obj.qmu) ']'],'dakota properties')); 1118 disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(obj.results) ']'],'model results')); 1119 disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay')); 1120 disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields')); 1121 end % }}} 1122 end 1123 1123 end -
issm/trunk-jpl/src/m/classes/model/model.py
r13470 r13692 170 170 # }}} 171 171 172 """ 173 function md2 = extract(md,area) % {{{ 174 %extract - extract a model according to an Argus contour or flag list 175 % 176 % This routine extracts a submodel from a bigger model with respect to a given contour 177 % md must be followed by the corresponding exp file or flags list 178 % It can either be a domain file (argus type, .exp extension), or an array of element flags. 179 % If user wants every element outside the domain to be 180 % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); 181 % an empty string '' will be considered as an empty domain 182 % a string 'all' will be considered as the entire domain 183 % add an argument 0 if you do not want the elements to be checked (faster) 184 % 185 % Usage: 186 % md2=extract(md,area); 187 % 188 % Examples: 189 % md2=extract(md,'Domain.exp'); 190 % md2=extract(md,md.mask.elementonfloatingice); 191 % 192 % See also: EXTRUDE, COLLAPSE 193 194 %copy model 195 md1=md; 196 197 %some checks 198 if ((nargin~=2) | (nargout~=1)), 199 help extract 200 error('extract error message: bad usage'); 201 end 202 203 %get check option 204 if (nargin==3 & varargin{1}==0), 205 checkoutline=0; 206 else 207 checkoutline=1; 208 end 209 210 %get elements that are inside area 211 flag_elem=FlagElements(md1,area); 212 if ~any(flag_elem), 213 error('extracted model is empty'); 214 end 215 216 %kick out all elements with 3 dirichlets 217 spc_elem=find(~flag_elem); 218 spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); 219 flag=ones(md1.mesh.numberofvertices,1); 220 flag(spc_node)=0; 221 pos=find(sum(flag(md1.mesh.elements),2)==0); 222 flag_elem(pos)=0; 223 224 %extracted elements and nodes lists 225 pos_elem=find(flag_elem); 226 pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); 227 228 %keep track of some fields 229 numberofvertices1=md1.mesh.numberofvertices; 230 numberofelements1=md1.mesh.numberofelements; 231 numberofvertices2=length(pos_node); 232 numberofelements2=length(pos_elem); 233 flag_node=zeros(numberofvertices1,1); 234 flag_node(pos_node)=1; 235 236 %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) 237 Pelem=zeros(numberofelements1,1); 238 Pelem(pos_elem)=[1:numberofelements2]'; 239 Pnode=zeros(numberofvertices1,1); 240 Pnode(pos_node)=[1:numberofvertices2]'; 241 242 %renumber the elements (some node won't exist anymore) 243 elements_1=md1.mesh.elements; 244 elements_2=elements_1(pos_elem,:); 245 elements_2(:,1)=Pnode(elements_2(:,1)); 246 elements_2(:,2)=Pnode(elements_2(:,2)); 247 elements_2(:,3)=Pnode(elements_2(:,3)); 248 if md1.mesh.dimension==3, 249 elements_2(:,4)=Pnode(elements_2(:,4)); 250 elements_2(:,5)=Pnode(elements_2(:,5)); 251 elements_2(:,6)=Pnode(elements_2(:,6)); 252 end 253 254 %OK, now create the new model ! 255 256 %take every fields from model 257 md2=md1; 258 259 %automatically modify fields 260 261 %loop over model fields 262 model_fields=fields(md1); 263 for i=1:length(model_fields), 264 %get field 265 field=md1.(model_fields{i}); 266 fieldsize=size(field); 267 if isobject(field), %recursive call 268 object_fields=fields(md1.(model_fields{i})); 269 for j=1:length(object_fields), 270 %get field 271 field=md1.(model_fields{i}).(object_fields{j}); 272 fieldsize=size(field); 273 %size = number of nodes * n 274 if fieldsize(1)==numberofvertices1 275 md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); 276 elseif (fieldsize(1)==numberofvertices1+1) 277 md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; 278 %size = number of elements * n 279 elseif fieldsize(1)==numberofelements1 280 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); 281 end 282 end 283 else 284 %size = number of nodes * n 285 if fieldsize(1)==numberofvertices1 286 md2.(model_fields{i})=field(pos_node,:); 287 elseif (fieldsize(1)==numberofvertices1+1) 288 md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; 289 %size = number of elements * n 290 elseif fieldsize(1)==numberofelements1 291 md2.(model_fields{i})=field(pos_elem,:); 292 end 293 end 294 end 295 296 %modify some specific fields 297 298 %Mesh 299 md2.mesh.numberofelements=numberofelements2; 300 md2.mesh.numberofvertices=numberofvertices2; 301 md2.mesh.elements=elements_2; 302 303 %mesh.uppervertex mesh.lowervertex 304 if md1.mesh.dimension==3 305 md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); 306 pos=find(~isnan(md2.mesh.uppervertex)); 307 md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); 308 309 md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); 310 pos=find(~isnan(md2.mesh.lowervertex)); 311 md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); 312 313 md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); 314 pos=find(~isnan(md2.mesh.upperelements)); 315 md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); 316 317 md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); 318 pos=find(~isnan(md2.mesh.lowerelements)); 319 md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); 320 end 321 322 %Initial 2d mesh 323 if md1.mesh.dimension==3 324 flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); 325 pos_elem_2d=find(flag_elem_2d); 326 flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); 327 pos_node_2d=find(flag_node_2d); 328 329 md2.mesh.numberofelements2d=length(pos_elem_2d); 330 md2.mesh.numberofvertices2d=length(pos_node_2d); 331 md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); 332 md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); 333 md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); 334 md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); 335 336 md2.mesh.x2d=md1.mesh.x(pos_node_2d); 337 md2.mesh.y2d=md1.mesh.y(pos_node_2d); 338 end 339 340 %Edges 341 if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... 342 %renumber first two columns 343 pos=find(md2.mesh.edges(:,4)~=-1); 344 md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); 345 md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); 346 md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); 347 md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); 348 %remove edges when the 2 vertices are not in the domain. 349 md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); 350 %Replace all zeros by -1 in the last two columns; 351 pos=find(md2.mesh.edges(:,3)==0); 352 md2.mesh.edges(pos,3)=-1; 353 pos=find(md2.mesh.edges(:,4)==0); 354 md2.mesh.edges(pos,4)=-1; 355 %Invert -1 on the third column with last column (Also invert first two columns!!) 356 pos=find(md2.mesh.edges(:,3)==-1); 357 md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); 358 md2.mesh.edges(pos,4)=-1; 359 values=md2.mesh.edges(pos,2); 360 md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); 361 md2.mesh.edges(pos,1)=values; 362 %Finally remove edges that do not belong to any element 363 pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); 364 md2.mesh.edges(pos,:)=[]; 365 end 366 367 %Penalties 368 if ~isnan(md2.diagnostic.vertex_pairing), 369 for i=1:size(md1.diagnostic.vertex_pairing,1); 370 md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); 371 end 372 md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); 373 end 374 if ~isnan(md2.prognostic.vertex_pairing), 375 for i=1:size(md1.prognostic.vertex_pairing,1); 376 md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); 377 end 378 md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); 379 end 380 381 %recreate segments 382 if md1.mesh.dimension==2 383 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 384 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 385 md2.mesh.segments=contourenvelope(md2); 386 md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 387 else 388 %First do the connectivity for the contourenvelope in 2d 389 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); 390 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); 391 md2.mesh.segments=contourenvelope(md2); 392 md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 393 md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); 394 %Then do it for 3d as usual 395 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 396 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 397 end 398 399 %Boundary conditions: Dirichlets on new boundary 400 %Catch the elements that have not been extracted 401 orphans_elem=find(~flag_elem); 402 orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; 403 %Figure out which node are on the boundary between md2 and md1 404 nodestoflag1=intersect(orphans_node,pos_node); 405 nodestoflag2=Pnode(nodestoflag1); 406 if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, 407 if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 408 md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 409 md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); 410 else 411 md2.diagnostic.spcvx(nodestoflag2)=NaN; 412 md2.diagnostic.spcvy(nodestoflag2)=NaN; 413 disp(' ') 414 disp('!! extract warning: spc values should be checked !!') 415 disp(' ') 416 end 417 %put 0 for vz 418 md2.diagnostic.spcvz(nodestoflag2)=0; 419 end 420 if ~isnan(md1.thermal.spctemperature), 421 md2.thermal.spctemperature(nodestoflag2,1)=1; 422 end 423 424 %Diagnostic 425 if ~isnan(md2.diagnostic.icefront) 426 md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); 427 md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); 428 md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); 429 if md1.mesh.dimension==3 430 md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); 431 md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); 432 end 433 md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); 434 end 435 436 %Results fields 437 if isstruct(md1.results), 438 md2.results=struct(); 439 solutionfields=fields(md1.results); 440 for i=1:length(solutionfields), 441 %get subfields 442 solutionsubfields=fields(md1.results.(solutionfields{i})); 443 for j=1:length(solutionsubfields), 444 field=md1.results.(solutionfields{i}).(solutionsubfields{j}); 445 if length(field)==numberofvertices1, 446 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); 447 elseif length(field)==numberofelements1, 448 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); 449 else 450 md2.results.(solutionfields{i}).(solutionsubfields{j})=field; 451 end 452 end 453 end 454 end 455 456 %Keep track of pos_node and pos_elem 457 md2.mesh.extractedvertices=pos_node; 458 md2.mesh.extractedelements=pos_elem; 459 end % }}} 460 """ 461 172 462 def extrude(md,*args): # {{{ 173 463 """
Note:
See TracChangeset
for help on using the changeset viewer.