![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
- Timestamp:
- 08/20/12 17:39:30 (13 years ago)
- Location:
- issm/branches/trunk-jpl-damage
- Files:
-
- 2 edited
-
. (modified) (2 props)
-
src/m/classes/model/model.m (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-jpl-damage
- Property svn:ignore
-
old new 1 projects 1 2 autom4te.cache 2 3 aclocal.m4
-
- Property svn:mergeinfo changed
/issm/trunk-jpl merged: 12948-13099
- Property svn:ignore
-
issm/branches/trunk-jpl-damage/src/m/classes/model/model.m
r12878 r13101 91 91 end 92 92 %}}} 93 function md = collapse(md)% {{{ 94 %COLLAPSE - collapses a 3d mesh into a 2d mesh 95 % 96 % This routine collapses a 3d model into a 2d model 97 % and collapses all the fileds of the 3d model by 98 % taking their depth-averaged values 99 % 100 % Usage: 101 % md=collapse(md) 102 % 103 % See also: EXTRUDE, MODELEXTRACT 104 105 %Check that the model is really a 3d model 106 if ~md.mesh.dimension==3, 107 error('collapse error message: only 3d mesh can be collapsed') 108 end 109 110 %Start with changing alle the fields from the 3d mesh 111 112 %drag is limited to nodes that are on the bedrock. 113 md.friction.coefficient=project2d(md,md.friction.coefficient,1); 114 115 %p and q (same deal, except for element that are on the bedrock: ) 116 md.friction.p=project2d(md,md.friction.p,1); 117 md.friction.q=project2d(md,md.friction.q,1); 118 119 %observations 120 if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; 121 if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; 122 if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; 123 if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; 124 if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; 125 if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; 126 if ~isnan(md.surfaceforcings.mass_balance), 127 md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); 128 end; 129 if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; 130 131 %results 132 if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; 133 if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; 134 if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; 135 if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; 136 if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; 137 138 %bedinfo and surface info 139 md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); 140 md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); 141 md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); 142 md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); 143 144 %elementstype 145 if ~isnan(md.flowequation.element_equation) 146 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); 147 md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); 148 md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); 149 md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); 150 md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); 151 end 152 153 %boundary conditions 154 md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); 155 md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); 156 md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); 157 md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); 158 md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); 159 md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); 160 161 %Extrusion of Neumann BC 162 if ~isnan(md.diagnostic.icefront), 163 numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); 164 md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer 165 end 166 167 %materials 168 md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); 169 md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); 170 171 %special for thermal modeling: 172 md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); 173 md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux 174 175 %update of connectivity matrix 176 md.mesh.average_vertex_connectivity=25; 177 178 %Collapse the mesh 179 nodes2d=md.mesh.numberofvertices2d; 180 elements2d=md.mesh.numberofelements2d; 181 182 %parameters 183 md.geometry.surface=project2d(md,md.geometry.surface,1); 184 md.geometry.thickness=project2d(md,md.geometry.thickness,1); 185 md.geometry.bed=project2d(md,md.geometry.bed,1); 186 md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); 187 md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); 188 md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); 189 md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); 190 md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); 191 md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); 192 md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); 193 md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); 194 md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); 195 196 %lat long 197 if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end 198 if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end 199 200 %Initialize with the 2d mesh 201 md.mesh.x=md.mesh.x2d; 202 md.mesh.y=md.mesh.y2d; 203 md.mesh.z=zeros(size(md.mesh.x2d)); 204 md.mesh.numberofvertices=md.mesh.numberofvertices2d; 205 md.mesh.numberofelements=md.mesh.numberofelements2d; 206 md.mesh.elements=md.mesh.elements2d; 207 208 %Keep a trace of lower and upper nodes 209 md.mesh.lowervertex=NaN; 210 md.mesh.uppervertex=NaN; 211 md.mesh.lowerelements=NaN; 212 md.mesh.upperelements=NaN; 213 214 %Remove old mesh 215 md.mesh.x2d=NaN; 216 md.mesh.y2d=NaN; 217 md.mesh.elements2d=NaN; 218 md.mesh.numberofelements2d=md.mesh.numberofelements; 219 md.mesh.numberofvertices2d=md.mesh.numberofvertices; 220 md.mesh.numberoflayers=0; 221 222 %Update mesh type 223 md.mesh.dimension=2; 224 end % }}} 225 function md2 = extract(md,area) % {{{ 226 %extract - extract a model according to an Argus contour or flag list 227 % 228 % This routine extracts a submodel from a bigger model with respect to a given contour 229 % md must be followed by the corresponding exp file or flags list 230 % It can either be a domain file (argus type, .exp extension), or an array of element flags. 231 % If user wants every element outside the domain to be 232 % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); 233 % an empty string '' will be considered as an empty domain 234 % a string 'all' will be considered as the entire domain 235 % add an argument 0 if you do not want the elements to be checked (faster) 236 % 237 % Usage: 238 % md2=extract(md,area); 239 % 240 % Examples: 241 % md2=extract(md,'Domain.exp'); 242 % md2=extract(md,md.mask.elementonfloatingice); 243 % 244 % See also: EXTRUDE, COLLAPSE 245 246 %copy model 247 md1=md; 248 249 %some checks 250 if ((nargin~=2) | (nargout~=1)), 251 help extract 252 error('extract error message: bad usage'); 253 end 254 255 %get check option 256 if (nargin==3 & varargin{1}==0), 257 checkoutline=0; 258 else 259 checkoutline=1; 260 end 261 262 %get elements that are inside area 263 flag_elem=FlagElements(md1,area); 264 if ~any(flag_elem), 265 error('extracted model is empty'); 266 end 267 268 %kick out all elements with 3 dirichlets 269 spc_elem=find(~flag_elem); 270 spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); 271 flag=ones(md1.mesh.numberofvertices,1); 272 flag(spc_node)=0; 273 pos=find(sum(flag(md1.mesh.elements),2)==0); 274 flag_elem(pos)=0; 275 276 %extracted elements and nodes lists 277 pos_elem=find(flag_elem); 278 pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); 279 280 %keep track of some fields 281 numberofvertices1=md1.mesh.numberofvertices; 282 numberofelements1=md1.mesh.numberofelements; 283 numberofvertices2=length(pos_node); 284 numberofelements2=length(pos_elem); 285 flag_node=zeros(numberofvertices1,1); 286 flag_node(pos_node)=1; 287 288 %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) 289 Pelem=zeros(numberofelements1,1); 290 Pelem(pos_elem)=[1:numberofelements2]'; 291 Pnode=zeros(numberofvertices1,1); 292 Pnode(pos_node)=[1:numberofvertices2]'; 293 294 %renumber the elements (some node won't exist anymore) 295 elements_1=md1.mesh.elements; 296 elements_2=elements_1(pos_elem,:); 297 elements_2(:,1)=Pnode(elements_2(:,1)); 298 elements_2(:,2)=Pnode(elements_2(:,2)); 299 elements_2(:,3)=Pnode(elements_2(:,3)); 300 if md1.mesh.dimension==3, 301 elements_2(:,4)=Pnode(elements_2(:,4)); 302 elements_2(:,5)=Pnode(elements_2(:,5)); 303 elements_2(:,6)=Pnode(elements_2(:,6)); 304 end 305 306 %OK, now create the new model ! 307 308 %take every fields from model 309 md2=md1; 310 311 %automatically modify fields 312 313 %loop over model fields 314 model_fields=fields(md1); 315 for i=1:length(model_fields), 316 %get field 317 field=md1.(model_fields{i}); 318 fieldsize=size(field); 319 if isobject(field), %recursive call 320 object_fields=fields(md1.(model_fields{i})); 321 for j=1:length(object_fields), 322 %get field 323 field=md1.(model_fields{i}).(object_fields{j}); 324 fieldsize=size(field); 325 %size = number of nodes * n 326 if fieldsize(1)==numberofvertices1 327 md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); 328 elseif (fieldsize(1)==numberofvertices1+1) 329 md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; 330 %size = number of elements * n 331 elseif fieldsize(1)==numberofelements1 332 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); 333 end 334 end 335 else 336 %size = number of nodes * n 337 if fieldsize(1)==numberofvertices1 338 md2.(model_fields{i})=field(pos_node,:); 339 elseif (fieldsize(1)==numberofvertices1+1) 340 md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; 341 %size = number of elements * n 342 elseif fieldsize(1)==numberofelements1 343 md2.(model_fields{i})=field(pos_elem,:); 344 end 345 end 346 end 347 348 %modify some specific fields 349 350 %Mesh 351 md2.mesh.numberofelements=numberofelements2; 352 md2.mesh.numberofvertices=numberofvertices2; 353 md2.mesh.elements=elements_2; 354 355 %mesh.uppervertex mesh.lowervertex 356 if md1.mesh.dimension==3 357 md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); 358 pos=find(~isnan(md2.mesh.uppervertex)); 359 md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); 360 361 md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); 362 pos=find(~isnan(md2.mesh.lowervertex)); 363 md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); 364 365 md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); 366 pos=find(~isnan(md2.mesh.upperelements)); 367 md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); 368 369 md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); 370 pos=find(~isnan(md2.mesh.lowerelements)); 371 md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); 372 end 373 374 %Initial 2d mesh 375 if md1.mesh.dimension==3 376 flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); 377 pos_elem_2d=find(flag_elem_2d); 378 flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); 379 pos_node_2d=find(flag_node_2d); 380 381 md2.mesh.numberofelements2d=length(pos_elem_2d); 382 md2.mesh.numberofvertices2d=length(pos_node_2d); 383 md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); 384 md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); 385 md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); 386 md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); 387 388 md2.mesh.x2d=md1.mesh.x(pos_node_2d); 389 md2.mesh.y2d=md1.mesh.y(pos_node_2d); 390 end 391 392 %Edges 393 if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... 394 %renumber first two columns 395 pos=find(md2.mesh.edges(:,4)~=-1); 396 md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); 397 md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); 398 md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); 399 md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); 400 %remove edges when the 2 vertices are not in the domain. 401 md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); 402 %Replace all zeros by -1 in the last two columns; 403 pos=find(md2.mesh.edges(:,3)==0); 404 md2.mesh.edges(pos,3)=-1; 405 pos=find(md2.mesh.edges(:,4)==0); 406 md2.mesh.edges(pos,4)=-1; 407 %Invert -1 on the third column with last column (Also invert first two columns!!) 408 pos=find(md2.mesh.edges(:,3)==-1); 409 md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); 410 md2.mesh.edges(pos,4)=-1; 411 values=md2.mesh.edges(pos,2); 412 md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); 413 md2.mesh.edges(pos,1)=values; 414 %Finally remove edges that do not belong to any element 415 pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); 416 md2.mesh.edges(pos,:)=[]; 417 end 418 419 %Penalties 420 if ~isnan(md2.diagnostic.vertex_pairing), 421 for i=1:size(md1.diagnostic.vertex_pairing,1); 422 md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); 423 end 424 md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); 425 end 426 if ~isnan(md2.prognostic.vertex_pairing), 427 for i=1:size(md1.prognostic.vertex_pairing,1); 428 md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); 429 end 430 md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); 431 end 432 433 %recreate segments 434 if md1.mesh.dimension==2 435 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 436 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 437 md2.mesh.segments=contourenvelope(md2); 438 md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 439 else 440 %First do the connectivity for the contourenvelope in 2d 441 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); 442 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); 443 md2.mesh.segments=contourenvelope(md2); 444 md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; 445 md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); 446 %Then do it for 3d as usual 447 md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); 448 md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); 449 end 450 451 %Boundary conditions: Dirichlets on new boundary 452 %Catch the elements that have not been extracted 453 orphans_elem=find(~flag_elem); 454 orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; 455 %Figure out which node are on the boundary between md2 and md1 456 nodestoflag1=intersect(orphans_node,pos_node); 457 nodestoflag2=Pnode(nodestoflag1); 458 if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, 459 if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 460 md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 461 md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); 462 else 463 md2.diagnostic.spcvx(nodestoflag2)=NaN; 464 md2.diagnostic.spcvy(nodestoflag2)=NaN; 465 disp(' ') 466 disp('!! extract warning: spc values should be checked !!') 467 disp(' ') 468 end 469 %put 0 for vz 470 md2.diagnostic.spcvz(nodestoflag2)=0; 471 end 472 if ~isnan(md1.thermal.spctemperature), 473 md2.thermal.spctemperature(nodestoflag2,1)=1; 474 end 475 476 %Diagnostic 477 if ~isnan(md2.diagnostic.icefront) 478 md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); 479 md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); 480 md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); 481 if md1.mesh.dimension==3 482 md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); 483 md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); 484 end 485 md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); 486 end 487 488 %Results fields 489 if isstruct(md1.results), 490 md2.results=struct(); 491 solutionfields=fields(md1.results); 492 for i=1:length(solutionfields), 493 %get subfields 494 solutionsubfields=fields(md1.results.(solutionfields{i})); 495 for j=1:length(solutionsubfields), 496 field=md1.results.(solutionfields{i}).(solutionsubfields{j}); 497 if length(field)==numberofvertices1, 498 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); 499 elseif length(field)==numberofelements1, 500 md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); 501 else 502 md2.results.(solutionfields{i}).(solutionsubfields{j})=field; 503 end 504 end 505 end 506 end 507 508 %Keep track of pos_node and pos_elem 509 md2.mesh.extractedvertices=pos_node; 510 md2.mesh.extractedelements=pos_elem; 511 end % }}} 512 function md = extrude(md,varargin) % {{{ 513 %EXTRUDE - vertically extrude a 2d mesh 514 % 515 % vertically extrude a 2d mesh and create corresponding 3d mesh. 516 % The vertical distribution can: 517 % - follow a polynomial law 518 % - follow two polynomial laws, one for the lower part and one for the upper part of the mesh 519 % - be discribed by a list of coefficients (between 0 and 1) 520 % 521 % 522 % Usage: 523 % md=extrude(md,numlayers,extrusionexponent); 524 % md=extrude(md,numlayers,lowerexponent,upperexponent); 525 % md=extrude(md,listofcoefficients); 526 % 527 % Example: 528 % md=extrude(md,8,3); 529 % md=extrude(md,8,3,2); 530 % md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); 531 % 532 % See also: MODELEXTRACT, COLLAPSE 533 534 %some checks on list of arguments 535 if ((nargin>4) | (nargin<2) | (nargout~=1)), 536 help extrude; 537 error('extrude error message'); 538 end 539 540 %Extrude the mesh 541 if nargin==2, %list of coefficients 542 list=varargin{1}; 543 if any(list<0) | any(list>1), 544 error('extrusioncoefficients must be between 0 and 1'); 545 end 546 extrusionlist=sort(unique([list(:);0;1])); 547 numlayers=length(extrusionlist); 548 elseif nargin==3, %one polynomial law 549 if varargin{2}<=0, 550 help extrude; 551 error('extrusionexponent must be >=0'); 552 end 553 numlayers=varargin{1}; 554 extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; 555 elseif nargin==4, %two polynomial laws 556 numlayers=varargin{1}; 557 lowerexp=varargin{2}; 558 upperexp=varargin{3}; 559 560 if varargin{2}<=0 | varargin{3}<=0, 561 help extrude; 562 error('lower and upper extrusionexponents must be >=0'); 563 end 564 565 lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; 566 upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; 567 extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); 568 569 end 570 571 if numlayers<2, 572 error('number of layers should be at least 2'); 573 end 574 if md.mesh.dimension==3, 575 error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); 576 end 577 578 %Initialize with the 2d mesh 579 x3d=[]; 580 y3d=[]; 581 z3d=[]; %the lower node is on the bed 582 thickness3d=md.geometry.thickness; %thickness and bed for these nodes 583 bed3d=md.geometry.bed; 584 585 %Create the new layers 586 for i=1:numlayers, 587 x3d=[x3d; md.mesh.x]; 588 y3d=[y3d; md.mesh.y]; 589 %nodes are distributed between bed and surface accordingly to the given exponent 590 z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; 591 end 592 number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh 593 594 %Extrude elements 595 elements3d=[]; 596 for i=1:numlayers-1, 597 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 598 end 599 number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh 600 601 %Keep a trace of lower and upper nodes 602 mesh.lowervertex=NaN*ones(number_nodes3d,1); 603 mesh.uppervertex=NaN*ones(number_nodes3d,1); 604 mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; 605 mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; 606 md.mesh.lowervertex=mesh.lowervertex; 607 md.mesh.uppervertex=mesh.uppervertex; 608 609 %same for lower and upper elements 610 mesh.lowerelements=NaN*ones(number_el3d,1); 611 mesh.upperelements=NaN*ones(number_el3d,1); 612 mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; 613 mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; 614 md.mesh.lowerelements=mesh.lowerelements; 615 md.mesh.upperelements=mesh.upperelements; 616 617 %Save old mesh 618 md.mesh.x2d=md.mesh.x; 619 md.mesh.y2d=md.mesh.y; 620 md.mesh.elements2d=md.mesh.elements; 621 md.mesh.numberofelements2d=md.mesh.numberofelements; 622 md.mesh.numberofvertices2d=md.mesh.numberofvertices; 623 624 %Update mesh type 625 md.mesh.dimension=3; 626 627 %Build global 3d mesh 628 md.mesh.elements=elements3d; 629 md.mesh.x=x3d; 630 md.mesh.y=y3d; 631 md.mesh.z=z3d; 632 md.mesh.numberofelements=number_el3d; 633 md.mesh.numberofvertices=number_nodes3d; 634 md.mesh.numberoflayers=numlayers; 635 636 %Ok, now deal with the other fields from the 2d mesh: 637 638 %lat long 639 md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); 640 md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); 641 642 %drag coefficient is limited to nodes that are on the bedrock. 643 md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1); 644 645 %p and q (same deal, except for element that are on the bedrock: ) 646 md.friction.p=project3d(md,'vector',md.friction.p,'type','element'); 647 md.friction.q=project3d(md,'vector',md.friction.q,'type','element'); 648 649 %observations 650 md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node'); 651 md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node'); 652 md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node'); 653 md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node'); 654 md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node'); 655 md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node'); 656 md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node'); 657 658 %results 659 if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; 660 if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; 661 if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; 662 if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; 663 if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; 664 if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; 665 666 %bedinfo and surface info 667 md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); 668 md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); 669 md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); 670 md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); 671 672 %elementstype 673 if ~isnan(md.flowequation.element_equation) 674 oldelements_type=md.flowequation.element_equation; 675 md.flowequation.element_equation=zeros(number_el3d,1); 676 md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); 677 end 678 679 %verticestype 680 if ~isnan(md.flowequation.vertex_equation) 681 oldvertices_type=md.flowequation.vertex_equation; 682 md.flowequation.vertex_equation=zeros(number_nodes3d,1); 683 md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); 684 end 685 md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node'); 686 md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node'); 687 md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node'); 688 689 %boundary conditions 690 md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node'); 691 md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node'); 692 md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node'); 693 md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); 694 md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node'); 695 md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node'); 696 md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node'); 697 698 %in 3d, pressureload: [node1 node2 node3 node4 element] 699 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 700 pressureload=[]; 701 for i=1:numlayers-1, 702 pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; 703 end 704 md.diagnostic.icefront=pressureload; 705 706 %connectivity 707 md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); 708 md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; 709 for i=2:numlayers-1, 710 md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... 711 =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; 712 end 713 md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0; 714 715 %materials 716 md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); 717 md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); 718 719 %parameters 720 md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node'); 721 md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node'); 722 md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node'); 723 md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node'); 724 md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node'); 725 md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); 726 md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element'); 727 md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node'); 728 md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element'); 729 md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node'); 730 md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element'); 731 md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'type','node'); 732 if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; 733 if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; 734 if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; 735 if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end 736 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end 737 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end 738 if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end 739 740 %Put lithostatic pressure if there is an existing pressure 741 if ~isnan(md.initialization.pressure), 742 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); 743 end 744 745 %special for thermal modeling: 746 md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); 747 if ~isnan(md.basalforcings.geothermalflux) 748 md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux 749 end 750 751 %increase connectivity if less than 25: 752 if md.mesh.average_vertex_connectivity<=25, 753 md.mesh.average_vertex_connectivity=100; 754 end 755 end % }}} 93 756 function md = structtomodel(md,structmd) % {{{ 94 757 … … 385 1048 md.settings = settings(); 386 1049 md.solver = solver(); 387 if ismumps ,388 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum ,mumpsoptions);1050 if ismumps(), 1051 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),mumpsoptions()); 389 1052 else 390 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum ,iluasmoptions);1053 md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum(),iluasmoptions()); 391 1054 end 392 1055 md.cluster = generic();
Note:
See TracChangeset
for help on using the changeset viewer.
