Changeset 9714
- Timestamp:
- 09/09/11 07:28:59 (14 years ago)
- Location:
- issm/trunk
- Files:
-
- 60 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/classes/mesh.m
r9706 r9714 9 9 extractedvertices = modelfield('default',NaN,'marshall',false); 10 10 extractedelements = modelfield('default',NaN,'marshall',false); 11 vertexonboundary = modelfield('default',NaN,'marshall',false); 12 lat = modelfield('default',NaN,'marshall',false); 13 long = modelfield('default',NaN,'marshall',false); 14 hemisphere = modelfield('default',NaN,'marshall',false); 15 segments = modelfield('default',NaN,'marshall',false); 16 segmentmarkers = modelfield('default',NaN,'marshall',false); 11 17 end 12 18 methods -
issm/trunk/src/m/classes/model/model.m
r9706 r9714 43 43 %FIXME: all other fields should belong to other classes 44 44 45 %Mesh46 45 dim = modelfield('default',0,'marshall',true,'format','Integer'); 47 46 numberofelements = modelfield('default',0,'marshall',true,'format','Integer'); … … 55 54 edges = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); 56 55 57 %Initial 2d mesh58 56 numberofelements2d = modelfield('default',0,'marshall',true,'format','Integer'); 59 57 numberofnodes2d = modelfield('default',0,'marshall',true,'format','Integer'); … … 62 60 y2d = modelfield('default',NaN,'marshall',false); 63 61 64 %latlon of the coorinates65 lat = modelfield('default',NaN,'marshall',false);66 long = modelfield('default',NaN,'marshall',false);67 hemisphere = modelfield('default',NaN,'marshall',false);68 69 %Penalties70 segments = modelfield('default',NaN,'marshall',false);71 segmentmarkers = modelfield('default',NaN,'marshall',false);72 73 %Projections74 62 uppernodes = modelfield('default',NaN,'marshall',false); 75 63 upperelements = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',2); … … 77 65 lowernodes = modelfield('default',NaN,'marshall',false); 78 66 79 %Extrusion80 67 numlayers = modelfield('default',0,'marshall',true,'format','Integer'); 81 68 elementonbed = modelfield('default',NaN,'marshall',true,'format','BooleanMat','mattype',2); … … 83 70 nodeonbed = modelfield('default',NaN,'marshall',true,'format','BooleanMat','mattype',1); 84 71 nodeonsurface = modelfield('default',NaN,'marshall',true,'format','BooleanMat','mattype',1); 85 86 %Boundary conditions87 nodeonboundary = modelfield('default',NaN,'marshall',false);88 72 %}}} 89 73 end … … 234 218 if isfield(structmd,'gridonicesheet'), md.mask.vertexongroundedice=structmd.gridonicesheet; end 235 219 if isfield(structmd,'gridonwater'), md.mask.vertexonwater=structmd.gridonwater; end 236 if isfield(structmd,'gridonboundary'), md. nodeonboundary=structmd.gridonboundary; end220 if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end 237 221 if isfield(structmd,'petscoptions'), md.solver=structmd.petscoptions; end 238 222 if isfield(structmd,'g'), md.constants.g=structmd.g; end … … 328 312 if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end 329 313 if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end 314 if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end 315 if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end 316 if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end 317 if isfield(structmd,'lon'), md.mesh.lon=structmd.lon; end 318 if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end 319 if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end 330 320 331 321 %Field changes -
issm/trunk/src/m/kml/kml_mesh_elem.m
r8298 r9714 61 61 end 62 62 63 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...64 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))63 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 64 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 65 65 if ~exist('latsgn','var') 66 66 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 69 69 else 70 70 display('Converting x/y data to lat/long data.'); 71 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);71 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 72 72 end 73 73 end … … 160 160 161 161 for j=1:size(md.elements,2) 162 kring.coords(j,:)=[md. long(md.elements(i,j)) md.lat(md.elements(i,j)) alt];162 kring.coords(j,:)=[md.mesh.long(md.elements(i,j)) md.mesh.lat(md.elements(i,j)) alt]; 163 163 end 164 164 kring.coords(end,:)=kring.coords(1,:); -
issm/trunk/src/m/kml/kml_mesh_write.m
r9650 r9714 67 67 end 68 68 69 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...70 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))69 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 70 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 71 71 if ~exist('latsgn','var') 72 72 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 75 75 else 76 76 display('Converting x/y data to lat/long data.'); 77 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);77 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 78 78 end 79 79 end -
issm/trunk/src/m/kml/kml_part_edges.m
r9650 r9714 62 62 end 63 63 64 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...65 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))64 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 65 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 66 66 if ~exist('latsgn','var') 67 67 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 70 70 else 71 71 display('Converting x/y data to lat/long data.'); 72 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);72 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 73 73 end 74 74 end … … 234 234 % if first edge, write out first node 235 235 if ~elast 236 kline.coords(end+1,:)=[md. long(edgeper(j,1)) md.lat(edgeper(j,1)) alt];236 kline.coords(end+1,:)=[md.mesh.long(edgeper(j,1)) md.mesh.lat(edgeper(j,1)) alt]; 237 237 end 238 kline.coords(end+1,:)=[md. long(edgeper(j,2)) md.lat(edgeper(j,2)) alt];238 kline.coords(end+1,:)=[md.mesh.long(edgeper(j,2)) md.mesh.lat(edgeper(j,2)) alt]; 239 239 elast=elemper(j); 240 240 nlast=edgeper(j,2); … … 323 323 % write out first node of first side for half-edge to midpoint 324 324 % disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.']) 325 kline.coords(end+1,:)=[md. long(elemp(ielem,nlast)) ...326 md. lat(elemp(ielem,nlast)) alt];325 kline.coords(end+1,:)=[md.mesh.long(elemp(ielem,nlast)) ... 326 md.mesh.lat(elemp(ielem,nlast)) alt]; 327 327 end 328 328 nlast=0; 329 329 330 330 % write out midpoint of first side 331 kline.coords(end+1,:)=[(md. long(elemp(ielem,slast))...332 +md. long(elemp(ielem,mod(slast,3)+1)))/2. ...333 (md. lat(elemp(ielem,slast))...334 +md. lat(elemp(ielem,mod(slast,3)+1)))/2. alt];331 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,slast))... 332 +md.mesh.long(elemp(ielem,mod(slast,3)+1)))/2. ... 333 (md.mesh.lat(elemp(ielem,slast))... 334 +md.mesh.lat(elemp(ielem,mod(slast,3)+1)))/2. alt]; 335 335 end 336 336 … … 371 371 % write out half-edge from current node to midpoint of unshared side 372 372 % disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.']) 373 kline.coords(end+1,:)=[(md. long(elemp(ielem,nlast))...374 +md. long(elemp(ielem,nnext)))/2. ...375 (md. lat(elemp(ielem,nlast))...376 +md. lat(elemp(ielem,nnext)))/2. alt];373 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,nlast))... 374 +md.mesh.long(elemp(ielem,nnext)))/2. ... 375 (md.mesh.lat(elemp(ielem,nlast))... 376 +md.mesh.lat(elemp(ielem,nnext)))/2. alt]; 377 377 nlast=0; 378 378 … … 404 404 % all different, so cut through centroid 405 405 % disp(['element ielem=' int2str(ielem) ' centroid written.']) 406 kline.coords(end+1,:)=[sum(md. long(elemp(ielem,:)))/3. ...407 sum(md. lat(elemp(ielem,:)))/3. alt];406 kline.coords(end+1,:)=[sum(md.mesh.long(elemp(ielem,:)))/3. ... 407 sum(md.mesh.lat(elemp(ielem,:)))/3. alt]; 408 408 end 409 409 % one node is in current partition, so cut off other two nodes … … 419 419 % write out midpoint of opposite side 420 420 % disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.']) 421 kline.coords(end+1,:)=[(md. long(elemp(ielem,snext))...422 +md. long(elemp(ielem,mod(snext,3)+1)))/2. ...423 (md. lat(elemp(ielem,snext))...424 +md. lat(elemp(ielem,mod(snext,3)+1)))/2. alt];421 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,snext))... 422 +md.mesh.long(elemp(ielem,mod(snext,3)+1)))/2. ... 423 (md.mesh.lat(elemp(ielem,snext))... 424 +md.mesh.lat(elemp(ielem,mod(snext,3)+1)))/2. alt]; 425 425 elast=ielem; 426 426 nlast=0; … … 451 451 end 452 452 % disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.']) 453 kline.coords(end+1,:)=[md. long(elemp(elast,nnext)) ...454 md. lat(elemp(elast,nnext)) alt];453 kline.coords(end+1,:)=[md.mesh.long(elemp(elast,nnext)) ... 454 md.mesh.lat(elemp(elast,nnext)) alt]; 455 455 break 456 456 % if not unshared, advance perimeter list and watch for end -
issm/trunk/src/m/kml/kml_part_elems.m
r9650 r9714 62 62 end 63 63 64 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...65 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))64 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 65 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 66 66 if ~exist('latsgn','var') 67 67 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 70 70 else 71 71 display('Converting x/y data to lat/long data.'); 72 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);72 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 73 73 end 74 74 end … … 199 199 200 200 for j=1:size(elem,2) 201 kring.coords(j,:)=[md. long(elem(i,j)) md.lat(elem(i,j)) alt];201 kring.coords(j,:)=[md.mesh.long(elem(i,j)) md.mesh.lat(elem(i,j)) alt]; 202 202 end 203 203 kring.coords(end,:)=kring.coords(1,:); -
issm/trunk/src/m/kml/kml_part_flagedges.m
r9650 r9714 54 54 end 55 55 56 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...57 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))56 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 57 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 58 58 if ~exist('latsgn','var') 59 59 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 62 62 else 63 63 display('Converting x/y data to lat/long data.'); 64 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);64 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 65 65 end 66 66 end … … 74 74 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ... 75 75 md.qmu.numberofpartitions 76 [latseg,lonseg]=flagedges(md.elements,md. lat,md.long,md.qmu.partition);76 [latseg,lonseg]=flagedges(md.elements,md.mesh.lat,md.mesh.long,md.qmu.partition); 77 77 kfold=kml_folder(); 78 78 kfold.name ='Partition Segments'; -
issm/trunk/src/m/kml/kml_partitions.m
r9650 r9714 63 63 end 64 64 65 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...66 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))65 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 66 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 67 67 if ~exist('latsgn','var') 68 68 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 71 71 else 72 72 display('Converting x/y data to lat/long data.'); 73 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);73 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 74 74 end 75 75 end … … 236 236 % if first edge, write out first node 237 237 if ~elast 238 kring.coords(end+1,:)=[md. long(edgeper(j,1)) md.lat(edgeper(j,1)) alt];238 kring.coords(end+1,:)=[md.mesh.long(edgeper(j,1)) md.mesh.lat(edgeper(j,1)) alt]; 239 239 end 240 kring.coords(end+1,:)=[md. long(edgeper(j,2)) md.lat(edgeper(j,2)) alt];240 kring.coords(end+1,:)=[md.mesh.long(edgeper(j,2)) md.mesh.lat(edgeper(j,2)) alt]; 241 241 elast=elemper(j); 242 242 nlast=edgeper(j,2); … … 325 325 % write out first node of first side for half-edge to midpoint 326 326 % disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.']) 327 kring.coords(end+1,:)=[md. long(elemp(ielem,nlast)) ...328 md. lat(elemp(ielem,nlast)) alt];327 kring.coords(end+1,:)=[md.mesh.long(elemp(ielem,nlast)) ... 328 md.mesh.lat(elemp(ielem,nlast)) alt]; 329 329 end 330 330 nlast=0; 331 331 332 332 % write out midpoint of first side 333 kring.coords(end+1,:)=[(md. long(elemp(ielem,slast))...334 +md. long(elemp(ielem,mod(slast,3)+1)))/2. ...335 (md. lat(elemp(ielem,slast))...336 +md. lat(elemp(ielem,mod(slast,3)+1)))/2. alt];333 kring.coords(end+1,:)=[(md.mesh.long(elemp(ielem,slast))... 334 +md.mesh.long(elemp(ielem,mod(slast,3)+1)))/2. ... 335 (md.mesh.lat(elemp(ielem,slast))... 336 +md.mesh.lat(elemp(ielem,mod(slast,3)+1)))/2. alt]; 337 337 end 338 338 … … 373 373 % write out half-edge from current node to midpoint of unshared side 374 374 % disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.']) 375 kring.coords(end+1,:)=[(md. long(elemp(ielem,nlast))...376 +md. long(elemp(ielem,nnext)))/2. ...377 (md. lat(elemp(ielem,nlast))...378 +md. lat(elemp(ielem,nnext)))/2. alt];375 kring.coords(end+1,:)=[(md.mesh.long(elemp(ielem,nlast))... 376 +md.mesh.long(elemp(ielem,nnext)))/2. ... 377 (md.mesh.lat(elemp(ielem,nlast))... 378 +md.mesh.lat(elemp(ielem,nnext)))/2. alt]; 379 379 nlast=0; 380 380 … … 406 406 % all different, so cut through centroid 407 407 % disp(['element ielem=' int2str(ielem) ' centroid written.']) 408 kring.coords(end+1,:)=[sum(md. long(elemp(ielem,:)))/3. ...409 sum(md. lat(elemp(ielem,:)))/3. alt];408 kring.coords(end+1,:)=[sum(md.mesh.long(elemp(ielem,:)))/3. ... 409 sum(md.mesh.lat(elemp(ielem,:)))/3. alt]; 410 410 end 411 411 % one node is in current partition, so cut off other two nodes … … 421 421 % write out midpoint of opposite side 422 422 % disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.']) 423 kring.coords(end+1,:)=[(md. long(elemp(ielem,snext))...424 +md. long(elemp(ielem,mod(snext,3)+1)))/2. ...425 (md. lat(elemp(ielem,snext))...426 +md. lat(elemp(ielem,mod(snext,3)+1)))/2. alt];423 kring.coords(end+1,:)=[(md.mesh.long(elemp(ielem,snext))... 424 +md.mesh.long(elemp(ielem,mod(snext,3)+1)))/2. ... 425 (md.mesh.lat(elemp(ielem,snext))... 426 +md.mesh.lat(elemp(ielem,mod(snext,3)+1)))/2. alt]; 427 427 elast=ielem; 428 428 nlast=0; … … 453 453 end 454 454 % disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.']) 455 kring.coords(end+1,:)=[md. long(elemp(elast,nnext)) ...456 md. lat(elemp(elast,nnext)) alt];455 kring.coords(end+1,:)=[md.mesh.long(elemp(elast,nnext)) ... 456 md.mesh.lat(elemp(elast,nnext)) alt]; 457 457 break 458 458 % if not unshared, advance perimeter list and watch for end -
issm/trunk/src/m/kml/kml_unsh_edges.m
r9650 r9714 54 54 end 55 55 56 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...57 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))56 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 57 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 58 58 if ~exist('latsgn','var') 59 59 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']); … … 62 62 else 63 63 display('Converting x/y data to lat/long data.'); 64 [md. lat,md.long]=xy2ll(md.x,md.y,latsgn);64 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn); 65 65 end 66 66 end … … 104 104 105 105 for j=1:2 106 kline.coords(j,:)=[md. long(edgeuns(i,j)) md.lat(edgeuns(i,j)) alt];106 kline.coords(j,:)=[md.mesh.long(edgeuns(i,j)) md.mesh.lat(edgeuns(i,j)) alt]; 107 107 end 108 108 -
issm/trunk/src/m/model/bamg.m
r9691 r9714 321 321 md.elements=bamgmesh_out.Triangles(:,1:3); 322 322 md.edges=bamgmesh_out.IssmEdges; 323 md. segments=bamgmesh_out.IssmSegments(:,1:3);324 md. segmentmarkers=bamgmesh_out.IssmSegments(:,4);323 md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3); 324 md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4); 325 325 326 326 %Fill in rest of fields: … … 334 334 md.elementonbed=ones(md.numberofelements,1); 335 335 md.elementonsurface=ones(md.numberofelements,1); 336 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;336 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 337 337 md.elementconnectivity=md.private.bamg.mesh.ElementConnectivity; 338 338 md.elementconnectivity(find(isnan(md.elementconnectivity)))=0; -
issm/trunk/src/m/model/collapse.m
r9705 r9714 83 83 md.geometry.thickness=project2d(md,md.geometry.thickness,1); 84 84 md.geometry.bed=project2d(md,md.geometry.bed,1); 85 md. nodeonboundary=project2d(md,md.nodeonboundary,1);85 md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); 86 86 md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); 87 87 md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); -
issm/trunk/src/m/model/extrude.m
r9705 r9714 202 202 md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node'); 203 203 md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node'); 204 md. nodeonboundary=project3d(md,'vector',md.nodeonboundary,'type','node');204 md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); 205 205 md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element'); 206 206 md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node'); -
issm/trunk/src/m/model/ismodelselfconsistent.m
r9705 r9714 151 151 message(['model not consistent: md.rifts should be a structure!']); 152 152 end 153 if ~isempty(find(md. segmentmarkers>=2)),153 if ~isempty(find(md.mesh.segmentmarkers>=2)), 154 154 %We have segments with rift markers, but no rift structure! 155 155 message(['model not consistent: model ' md.miscellaneous.name ' should be processed for rifts (run meshprocessrifts)!']); -
issm/trunk/src/m/model/kmlimagesc.m
r9074 r9714 20 20 fontsize=getfieldvalue(options,'fontsize',12); 21 21 posting=getfieldvalue(options,'posting',.1); 22 minlong=getfieldvalue(options,'minlong',min(md. long));23 maxlong=getfieldvalue(options,'maxlong',max(md. long));24 minlat=getfieldvalue(options,'minlat',min(md. lat));25 maxlat=getfieldvalue(options,'maxlat',max(md. lat));22 minlong=getfieldvalue(options,'minlong',min(md.mesh.long)); 23 maxlong=getfieldvalue(options,'maxlong',max(md.mesh.long)); 24 minlat=getfieldvalue(options,'minlat',min(md.mesh.lat)); 25 maxlat=getfieldvalue(options,'maxlat',max(md.mesh.lat)); 26 26 minfield=getfieldvalue(options,'minfield',min(field)); 27 27 maxfield=getfieldvalue(options,'maxfield',max(field)); 28 28 29 29 %do we have hemisphere setup?: 30 if ~isstr(md. hemisphere),31 error('md. hemisphere should be ''s'' or ''n''');30 if ~isstr(md.mesh.hemisphere), 31 error('md.mesh.hemisphere should be ''s'' or ''n'''); 32 32 end 33 33 34 if strcmpi(md. hemisphere,'s'),34 if strcmpi(md.mesh.hemisphere,'s'), 35 35 hemisphere=1; 36 36 central_meridian=getfieldvalue(options,'central_meridian',45); 37 37 standard_parallel=getfieldvalue(options,'standard_parallel',70); 38 elseif strcmpi(md. hemisphere,'n'),38 elseif strcmpi(md.mesh.hemisphere,'n'), 39 39 hemisphere=-1; 40 40 central_meridian=getfieldvalue(options,'central_meridian',0); 41 41 standard_parallel=getfieldvalue(options,'standard_parallel',71); 42 42 else 43 error('md. hemisphere should be ''s'' or ''n''');43 error('md.mesh.hemisphere should be ''s'' or ''n'''); 44 44 end 45 45 … … 49 49 50 50 %regrid to lat,long grid 51 [x_m,y_m,field]=InterpFromMeshToGrid(md.elements,md. long,md.lat,field,minlong,maxlat,posting,posting,nlines,ncols,NaN);51 [x_m,y_m,field]=InterpFromMeshToGrid(md.elements,md.mesh.long,md.mesh.lat,field,minlong,maxlat,posting,posting,nlines,ncols,NaN); 52 52 field=flipud(field); 53 53 -
issm/trunk/src/m/model/mesh/meshconvert.m
r9644 r9714 32 32 md.elements=bamgmesh_out.Triangles(:,1:3); 33 33 md.edges=bamgmesh_out.IssmEdges; 34 md. segments=bamgmesh_out.IssmSegments(:,1:3);35 md. segmentmarkers=bamgmesh_out.IssmSegments(:,4);34 md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3); 35 md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4); 36 36 37 37 %Fill in rest of fields: … … 45 45 md.elementonbed=ones(md.numberofelements,1); 46 46 md.elementonsurface=ones(md.numberofelements,1); 47 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;47 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; -
issm/trunk/src/m/model/mesh/meshnodensity.m
r9703 r9714 25 25 %Mesh using TriMeshNoDensity 26 26 if strcmp(riftname,''), 27 [md.elements,md.x,md.y,md. segments,md.segmentmarkers]=TriMeshNoDensity(domainname);27 [md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshNoDensity(domainname); 28 28 else 29 29 [elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname); … … 50 50 md.y=y; 51 51 md.elements=elements; 52 md. segments=segments;53 md. segmentmarkers=segmentmarkers;52 md.mesh.segments=segments; 53 md.mesh.segmentmarkers=segmentmarkers; 54 54 end 55 55 … … 58 58 md.numberofnodes=length(md.x); 59 59 md.z=zeros(md.numberofnodes,1); 60 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;60 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 61 61 md.nodeonbed=ones(md.numberofnodes,1); 62 62 md.nodeonsurface=ones(md.numberofnodes,1); -
issm/trunk/src/m/model/mesh/meshrefine.m
r9619 r9714 16 16 17 17 %Refine using TriMeshRefine 18 [md.elements,md.x,md.y,md. segments,md.segmentmarkers]=TriMeshRefine(md.elements,md.x,md.y,md.segments,md.segmentmarkers,areas,'yes');18 [md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshRefine(md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers,areas,'yes'); 19 19 20 20 %Fill in rest of fields: … … 22 22 md.numberofnodes=length(md.x); 23 23 md.z=zeros(md.numberofnodes,1); 24 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;24 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 25 25 md.nodeonbed=ones(md.numberofnodes,1); 26 26 md.nodeonsurface=ones(md.numberofnodes,1); -
issm/trunk/src/m/model/mesh/meshyams.m
r9703 r9714 97 97 md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes); 98 98 md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity); 99 md. segments=findsegments(md);99 md.mesh.segments=findsegments(md); 100 100 md=meshyamsrecreateriftsegments(md); 101 101 end … … 110 110 111 111 %recreate segments 112 md. segments=findsegments(md);113 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;112 md.mesh.segments=findsegments(md); 113 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 114 114 115 115 %Fill in rest of fields: -
issm/trunk/src/m/model/mesh/reorder.m
r9451 r9714 21 21 %update all fields 22 22 md.elements=tnewnodes(md.elements(newelements,:)); 23 md. segments=[tnewnodes(md.segments(:,1)) tnewnodes(md.segments(:,2)) tnewelements(md.segments(:,3))];23 md.mesh.segments=[tnewnodes(md.mesh.segments(:,1)) tnewnodes(md.mesh.segments(:,2)) tnewelements(md.mesh.segments(:,3))]; 24 24 md.x=md.x(newnodes); 25 25 md.y=md.y(newnodes); 26 26 md.z=md.z(newnodes); 27 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;27 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; -
issm/trunk/src/m/model/mesh/rifts/meshaddrifts.m
r9619 r9714 74 74 end 75 75 76 md. segments(:,1:2)=nodeconv(md.segments(:,1:2));77 md. segments(:,3)=elconv(md.segments(:,3));76 md.mesh.segments(:,1:2)=nodeconv(md.mesh.segments(:,1:2)); 77 md.mesh.segments(:,3)=elconv(md.mesh.segments(:,3)); 78 78 79 79 end 80 80 81 81 %finish up "a la" mesh.h 82 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;82 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 83 83 md.nodeonbed=ones(md.numberofnodes,1); 84 84 md.nodeonsurface=ones(md.numberofnodes,1); -
issm/trunk/src/m/model/mesh/rifts/meshprocessoutsiderifts.m
r9661 r9714 68 68 69 69 %deal with segments 70 tipsegments=find((md. segments(:,1)==tip) | (md.segments(:,2)==tip));70 tipsegments=find((md.mesh.segments(:,1)==tip) | (md.mesh.segments(:,2)==tip)); 71 71 for k=1:length(tipsegments), 72 72 segment_index=tipsegments(k); 73 pos=find(md. segments(segment_index,1:2)~=tip);74 other_node=md. segments(segment_index,pos);73 pos=find(md.mesh.segments(segment_index,1:2)~=tip); 74 other_node=md.mesh.segments(segment_index,pos); 75 75 if ~isconnected(md.elements,other_node,tip), 76 pos=find(md. segments(segment_index,1:2)==tip);77 md. segments(segment_index,pos)=num;76 pos=find(md.mesh.segments(segment_index,1:2)==tip); 77 md.mesh.segments(segment_index,pos)=num; 78 78 end 79 79 end … … 86 86 md.numberofnodes=length(md.x); 87 87 md.z=zeros(md.numberofnodes,1); 88 md. nodeonboundary=zeros(length(md.x),1); md.nodeonboundary(md.segments(:,1:2))=1;88 md.mesh.vertexonboundary=zeros(length(md.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 89 89 md.rifts.numrifts=length(md.rifts.riftstruct); 90 90 md.flowequation.element_equation=3*ones(md.numberofelements,1); -
issm/trunk/src/m/model/mesh/rifts/meshprocessrifts.m
r9661 r9714 25 25 26 26 %Call MEX file 27 [md.elements,md.x,md.y,md. segments,md.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);27 [md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers); 28 28 if ~isstruct(md.rifts.riftstruct), 29 29 error('TriMeshProcessRifts did not find any rift'); … … 34 34 md.numberofnodes=length(md.x); 35 35 md.z=zeros(md.numberofnodes,1); 36 md. nodeonboundary=zeros(length(md.x),1); md.nodeonboundary(md.segments(:,1:2))=1;36 md.mesh.vertexonboundary=zeros(length(md.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 37 37 md.rifts.numrifts=length(md.rifts.riftstruct); 38 38 md.flowequation.element_equation=3*ones(md.numberofelements,1); -
issm/trunk/src/m/model/mesh/rifts/meshyamsrecreateriftsegments.m
r9619 r9714 11 11 12 12 %find tip1 and tip2 for this rift, in the new mesh created by yams. 13 pos=find_point(md.x(md. segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));14 tip1=md. segments(pos,1);15 pos=find_point(md.x(md. segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));16 tip2=md. segments(pos,1);13 pos=find_point(md.x(md.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 14 tip1=md.mesh.segments(pos,1); 15 pos=find_point(md.x(md.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2)); 16 tip2=md.mesh.segments(pos,1); 17 17 18 18 %start from tip1, and build segments of this rift. 19 pos=find_point(md.x(md. segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));19 pos=find_point(md.x(md.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 20 20 pos_record=[pos_record; pos]; 21 riftsegs=md. segments(pos,:);21 riftsegs=md.mesh.segments(pos,:); 22 22 while 1, 23 23 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3); 24 24 %find other segment that holds B. 25 pos=find(md. segments(:,1)==B);25 pos=find(md.mesh.segments(:,1)==B); 26 26 pos_record=[pos_record; pos]; 27 riftsegs=[riftsegs; md. segments(pos,:)];27 riftsegs=[riftsegs; md.mesh.segments(pos,:)]; 28 28 if riftsegs(end,2)==tip1, 29 29 break; … … 39 39 40 40 %find tip1 and tip2 for this rift, in the new mesh created by yams. 41 pos1=find_point(md.x(md. segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));42 tip1=md. segments(pos1,1);43 pos2=find_point(md.x(md. segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));44 tip2=md. segments(pos2,1);41 pos1=find_point(md.x(md.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 42 tip1=md.mesh.segments(pos1,1); 43 pos2=find_point(md.x(md.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2)); 44 tip2=md.mesh.segments(pos2,1); 45 45 if length(tip1)==2, 46 46 %swap. … … 53 53 54 54 pos_record=[pos_record; pos]; 55 riftsegs=md. segments(pos,:);55 riftsegs=md.mesh.segments(pos,:); 56 56 while 1, 57 57 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3); 58 58 %find other segment that holds B. 59 pos=find(md. segments(:,1)==B);59 pos=find(md.mesh.segments(:,1)==B); 60 60 pos_record=[pos_record; pos]; 61 riftsegs=[riftsegs; md. segments(pos,:)];61 riftsegs=[riftsegs; md.mesh.segments(pos,:)]; 62 62 if ((riftsegs(end,2)==tip2(1)) | (riftsegs(end,2)==tip2(2))), 63 63 %figure out which tip we reached … … 70 70 pos=pos2(index); 71 71 pos_record=[pos_record; pos]; 72 riftsegs=[riftsegs; md. segments(pos,:)];72 riftsegs=[riftsegs; md.mesh.segments(pos,:)]; 73 73 while 1, 74 74 A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3); 75 75 %find other segment that holds B. 76 pos=find(md. segments(:,1)==B);76 pos=find(md.mesh.segments(:,1)==B); 77 77 pos_record=[pos_record; pos]; 78 riftsegs=[riftsegs; md. segments(pos,:)];78 riftsegs=[riftsegs; md.mesh.segments(pos,:)]; 79 79 if riftsegs(end,2)==tip1, 80 80 break; … … 88 88 end 89 89 %take out rift segments from segments 90 md. segments(pos_record,:)=[];90 md.mesh.segments(pos_record,:)=[]; -
issm/trunk/src/m/model/mesh/setmesh.m
r9704 r9714 39 39 %Mesh using TriMesh 40 40 if strcmp(riftname,''), 41 [md.elements,md.x,md.y,md. segments,md.segmentmarkers]=TriMesh(domainname,area,'yes');41 [md.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,area,'yes'); 42 42 else 43 43 [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes'); … … 64 64 md.y=y; 65 65 md.elements=elements; 66 md. segments=segments;67 md. segmentmarkers=segmentmarkers;66 md.mesh.segments=segments; 67 md.mesh.segmentmarkers=segmentmarkers; 68 68 end 69 69 … … 72 72 md.numberofnodes=length(md.x); 73 73 md.z=zeros(md.numberofnodes,1); 74 md. nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;74 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 75 75 md.nodeonbed=ones(md.numberofnodes,1); 76 76 md.nodeonsurface=ones(md.numberofnodes,1); -
issm/trunk/src/m/model/modelextract.m
r9706 r9714 206 206 md2.elementconnectivity=ElementConnectivity(md2.elements,md2.nodeconnectivity); 207 207 md2.segments=contourenvelope(md2); 208 md2. nodeonboundary=zeros(numberofnodes2,1); md2.nodeonboundary(md2.segments(:,1:2))=1;208 md2.mesh.vertexonboundary=zeros(numberofnodes2,1); md2.nodeonboundary(md2.segments(:,1:2))=1; 209 209 end 210 210 -
issm/trunk/src/m/model/outflow.m
r9684 r9714 5 5 % flag=outflow(md); 6 6 7 A=md. segments(:,1);8 B=md. segments(:,2);7 A=md.mesh.segments(:,1); 8 B=md.mesh.segments(:,2); 9 9 Nx=-(md.y(A)-md.y(B)); 10 10 Ny= md.x(A)-md.x(B); -
issm/trunk/src/m/model/plot/applyoptions.m
r9650 r9714 9 9 10 10 %some defaults 11 if strcmpi(md. hemisphere,'n'), options=addfielddefault(options,'hemisphere','n');12 elseif strcmpi(md. hemisphere,'s'), options=addfielddefault(options,'hemisphere','s');11 if strcmpi(md.mesh.hemisphere,'n'), options=addfielddefault(options,'hemisphere','n'); 12 elseif strcmpi(md.mesh.hemisphere,'s'), options=addfielddefault(options,'hemisphere','s'); 13 13 end 14 14 … … 410 410 axis equal off 411 411 %box off 412 if strcmpi(md. hemisphere,'n') | strcmpi(md.hemisphere,'north'),412 if strcmpi(md.mesh.hemisphere,'n') | strcmpi(md.mesh.hemisphere,'north'), 413 413 A=expread([ issmdir() '/projects/Exp/GreenlandBoxFront.exp']); 414 414 [A.x A.y]=ll2xy(A.x,A.y,+1,45,70); 415 elseif strcmpi(md. hemisphere,'s') | strcmpi(md.hemisphere,'south'),415 elseif strcmpi(md.mesh.hemisphere,'s') | strcmpi(md.mesh.hemisphere,'south'), 416 416 A=expread([ issmdir() '/projects/Exp/Antarctica.exp']); 417 417 else -
issm/trunk/src/m/model/plot/kmlgroundoverlay.m
r7028 r9714 12 12 13 13 %first figure out if lat and long were computed! 14 if (isempty(md. lat) | isempty(md.long)),14 if (isempty(md.mesh.lat) | isempty(md.mesh.long)), 15 15 error('kmlgroundoverlay error message: project x,y onto lat,long fields of model!'); 16 16 end … … 28 28 29 29 %figure out min and max for lat and long of this image: 30 west=min(md. long);31 east=max(md. long);32 south=min(md. lat);33 north=max(md. lat);30 west=min(md.mesh.long); 31 east=max(md.mesh.long); 32 south=min(md.mesh.lat); 33 north=max(md.mesh.lat); 34 34 35 35 %print image at high resolution -
issm/trunk/src/m/model/plot/latlonoverlay.m
r8929 r9714 66 66 latitudes =lat*ones(size(longitudes)); 67 67 68 if strcmpi(md. hemisphere,'n'),68 if strcmpi(md.mesh.hemisphere,'n'), 69 69 if lat<0, continue; end 70 70 [x,y]=ll2xy(latitudes,longitudes,+1,45,70); 71 elseif strcmpi(md. hemisphere,'s'),71 elseif strcmpi(md.mesh.hemisphere,'s'), 72 72 if lat>0, continue; end 73 73 [x,y]=ll2xy(latitudes,longitudes,-1, 0,71); … … 105 105 for lon=-180:lonstep:180 106 106 107 if strcmpi(md. hemisphere,'n'),107 if strcmpi(md.mesh.hemisphere,'n'), 108 108 latitudes =0:resolution:90; 109 109 longitudes=lon*ones(size(latitudes)); 110 110 [x,y]=ll2xy(latitudes,longitudes,+1,45,70); 111 elseif strcmpi(md. hemisphere,'s'),111 elseif strcmpi(md.mesh.hemisphere,'s'), 112 112 latitudes =-90:resolution:0; 113 113 longitudes=lon*ones(size(latitudes)); -
issm/trunk/src/m/model/plot/plot_boundaries.m
r9619 r9714 15 15 [x y z elements is2d isplanet]=processmesh(md,[],options); 16 16 17 for i=1:size(md. segments,1),18 plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k.-');hold on;17 for i=1:size(md.mesh.segments,1), 18 plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k.-');hold on; 19 19 end 20 20 -
issm/trunk/src/m/model/plot/plot_riftfraction.m
r9619 r9714 19 19 20 20 %plot mesh boundaries 21 for i=1:size(md. segments,1),22 h1=plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k.-');21 for i=1:size(md.mesh.segments,1), 22 h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k.-'); 23 23 end 24 24 -
issm/trunk/src/m/model/plot/plot_riftnumbering.m
r9684 r9714 15 15 16 16 %plot mesh boundaries 17 for i=1:size(md. segments,1),18 plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k.-');17 for i=1:size(md.mesh.segments,1), 18 plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k.-'); 19 19 end 20 20 … … 24 24 if isstruct(md.rifts.riftstruct), 25 25 %plot mesh boundaries 26 for i=1:size(md. segments,1),27 h1=plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'b-');26 for i=1:size(md.mesh.segments,1), 27 h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'b-'); 28 28 end 29 29 for i=1:size(md.rifts.riftstruct,1), -
issm/trunk/src/m/model/plot/plot_riftpenetration.m
r9684 r9714 14 14 15 15 %plot mesh boundaries 16 for i=1:size(md. segments,1),17 plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k-');16 for i=1:size(md.mesh.segments,1), 17 plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k-'); 18 18 end 19 19 … … 23 23 if isstruct(md.rifts.riftstruct), 24 24 %plot mesh boundaries 25 for i=1:size(md. segments,1),26 h1=plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'b-');25 for i=1:size(md.mesh.segments,1), 26 h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'b-'); 27 27 end 28 28 for i=1:size(md.rifts.riftstruct,1), -
issm/trunk/src/m/model/plot/plot_riftrelvel.m
r9684 r9714 39 39 40 40 %plot mesh boundaries 41 for i=1:size(md. segments,1),42 plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k-');41 for i=1:size(md.mesh.segments,1), 42 plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k-'); 43 43 end 44 44 … … 59 59 60 60 %plot mesh boundaries 61 for i=1:size(md. segments,1),62 h1=plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'b-');61 for i=1:size(md.mesh.segments,1), 62 h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'b-'); 63 63 end 64 64 for i=1:md.rifts.numrifts, -
issm/trunk/src/m/model/plot/plot_riftvel.m
r9684 r9714 35 35 36 36 %plot mesh boundaries 37 for i=1:size(md. segments,1),38 plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'k.-');37 for i=1:size(md.mesh.segments,1), 38 plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'k.-'); 39 39 end 40 40 … … 53 53 54 54 %plot mesh boundaries 55 for i=1:size(md. segments,1),56 h1=plot(x(md. segments(i,1:2)),y(md.segments(i,1:2)),'b-');55 for i=1:size(md.mesh.segments,1), 56 h1=plot(x(md.mesh.segments(i,1:2)),y(md.mesh.segments(i,1:2)),'b-'); 57 57 end 58 58 -
issm/trunk/src/m/model/plot/plot_segments.m
r8472 r9714 12 12 %process mesh and data 13 13 [x y z elements is2d isplanet]=processmesh(md,[],options); 14 segments=md. segments;14 segments=md.mesh.segments; 15 15 16 16 if (md.dim==2), … … 27 27 if strcmpi(getfieldvalue(options,'segmentnumbering','off'),'on'), 28 28 text(sum(x(segments(:,1:2)),2)/2,sum(y(segments(:,1:2)),2)/2,sum(z(segments(:,1:2)),2)/2+1,... 29 num2str(md. segmentmarkers),...29 num2str(md.mesh.segmentmarkers),... 30 30 'backgroundcolor',[0.8 0.9 0.8],'HorizontalAlignment','center','VerticalAlignment','middle'); 31 31 end -
issm/trunk/src/m/model/plot/processmesh.m
r9532 r9714 21 21 y2d=md.y2d; 22 22 else 23 x=md. long;23 x=md.mesh.long; 24 24 x2d=md.x2d; 25 y=md. lat;25 y=md.mesh.lat; 26 26 y2d=md.y2d; 27 27 end -
issm/trunk/src/m/model/plot/showregion.m
r6860 r9714 15 15 axis equal off 16 16 %box off 17 if strcmpi(md. hemisphere,'n') | strcmpi(md.hemisphere,'north'),17 if strcmpi(md.mesh.hemisphere,'n') | strcmpi(md.mesh.hemisphere,'north'), 18 18 A=expread([issmdir() 'projects/Exp/Greenland.exp']); 19 elseif strcmpi(md. hemisphere,'s') | strcmpi(md.hemisphere,'south'),19 elseif strcmpi(md.mesh.hemisphere,'s') | strcmpi(md.mesh.hemisphere,'south'), 20 20 A=expread([issmdir() '/projects/Exp/Antarctica.exp']); 21 21 else -
issm/trunk/src/m/model/radarpower.m
r9620 r9714 29 29 %figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one. 30 30 if ~exist(options,'overlay_image'), 31 if strcmpi(md. hemisphere,'n'),31 if strcmpi(md.mesh.hemisphere,'n'), 32 32 if ~exist([issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpg']), 33 33 error(['radarpower error message: file ' issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']); … … 56 56 md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax); 57 57 58 elseif strcmpi(md. hemisphere,'s'),58 elseif strcmpi(md.mesh.hemisphere,'s'), 59 59 if highres, 60 60 if ~exist([issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']), -
issm/trunk/src/m/model/setmask2.m
r9641 r9714 120 120 121 121 %call findsegments to build segment using THIS conectivity 122 md. segments=findsegments(md,'elementconnectivity',elementconnectivity);122 md.mesh.segments=findsegments(md,'elementconnectivity',elementconnectivity); 123 123 124 124 %some final checks: … … 146 146 md.mask.elementongroundedice=elementongroundedice; 147 147 148 md. segmentmarkers(:)=1;148 md.mesh.segmentmarkers(:)=1; -
issm/trunk/src/m/solutions/flaim.m
r9650 r9714 50 50 % calculate latitude and longitude, if necessary 51 51 52 if isempty(md. lat) || ((numel(md.lat) == 1) && isnan(md.lat)) || ...53 isempty(md. long) || ((numel(md.long) == 1) && isnan(md.long))52 if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ... 53 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long)) 54 54 if ~isfield(options,'latsgn') 55 55 error(['Missing ''latsgn'' parameter to calculate missing lat/long values.']); … … 58 58 else 59 59 display('Converting x/y values to lat/long values.'); 60 [md. lat,md.long]=xy2ll(md.x,md.y,options.latsgn);60 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,options.latsgn); 61 61 end 62 62 end … … 80 80 81 81 display('Calling KMLMeshWrite.'); 82 KMLMeshWrite(md.miscellaneous.name,md.miscellaneous.notes,md.elements,md.nodeconnectivity,md. lat,md.long,md.qmu.partition,md.flaim.criterion,options.cmap,filekml);82 KMLMeshWrite(md.miscellaneous.name,md.miscellaneous.notes,md.elements,md.nodeconnectivity,md.mesh.lat,md.mesh.long,md.qmu.partition,md.flaim.criterion,options.cmap,filekml); 83 83 % for testing 84 84 %filekml='issm-split-geikie1-targets.kml'; -
issm/trunk/src/m/utils/BC/SetIceSheetBC.m
r9696 r9714 8 8 9 9 %node on Dirichlet 10 pos=find(md. nodeonboundary);10 pos=find(md.mesh.vertexonboundary); 11 11 md.diagnostic.spcvx=NaN*ones(md.numberofnodes,1); 12 12 md.diagnostic.spcvy=NaN*ones(md.numberofnodes,1); -
issm/trunk/src/m/utils/BC/SetIceShelfBC.m
r9684 r9714 20 20 if ~exist(icefrontfile), error(['SetIceShelfBC error message: ice front file ' icefrontfile ' not found']); end 21 21 nodeinsideicefront=ContourToMesh(md.elements,md.x,md.y,icefrontfile,'node',2); 22 nodeonicefront=double(md. nodeonboundary & nodeinsideicefront);22 nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront); 23 23 elseif nargin==1, 24 24 nodeonicefront=zeros(md.numberofnodes,1); … … 27 27 error('bad usage'); 28 28 end 29 pos=find(md. nodeonboundary & ~nodeonicefront);29 pos=find(md.mesh.vertexonboundary & ~nodeonicefront); 30 30 md.diagnostic.spcvx=NaN*ones(md.numberofnodes,1); 31 31 md.diagnostic.spcvy=NaN*ones(md.numberofnodes,1); … … 47 47 %segment on Ice Front 48 48 %segment on Neumann (Ice Front) 49 pos=find(nodeonicefront(md. segments(:,1)) | nodeonicefront(md.segments(:,2)));49 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); 50 50 if (md.dim==2) 51 pressureload=md. segments(pos,:);51 pressureload=md.mesh.segments(pos,:); 52 52 elseif md.dim==3 53 pressureload_layer1=[md. segments(pos,1:2) md.segments(pos,2)+md.numberofnodes2d md.segments(pos,1)+md.numberofnodes2d md.segments(pos,3)];53 pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.numberofnodes2d md.mesh.segments(pos,1)+md.numberofnodes2d md.mesh.segments(pos,3)]; 54 54 pressureload=[]; 55 55 for i=1:md.numlayers-1, -
issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m
r9696 r9714 24 24 end 25 25 nodeinsideicefront=ContourToMesh(md.elements,md.x,md.y,icefrontfile,'node',2); 26 vertexonicefront=double(md. nodeonboundary & nodeinsideicefront);26 vertexonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront); 27 27 else 28 28 %Guess where the ice front is 29 29 vertexonfloatingice=zeros(md.numberofnodes,1); 30 30 vertexonfloatingice(md.elements(find(md.mask.elementonfloatingice),:))=1; 31 vertexonicefront=double(md. nodeonboundary & vertexonfloatingice);31 vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice); 32 32 end 33 pos=find(md. nodeonboundary & ~vertexonicefront);33 pos=find(md.mesh.vertexonboundary & ~vertexonicefront); 34 34 if isempty(pos), 35 35 warning('SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually') … … 53 53 54 54 md.hydrology.spcwatercolumn=zeros(md.numberofnodes,2); 55 pos=find(md. nodeonboundary);55 pos=find(md.mesh.vertexonboundary); 56 56 md.hydrology.spcwatercolumn(pos,1)=1; 57 57 58 58 %segment on Neumann (Ice Front) 59 pos=find(vertexonicefront(md. segments(:,1)) | vertexonicefront(md.segments(:,2)));59 pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); 60 60 if (md.dim==2) 61 pressureload=md. segments(pos,:);61 pressureload=md.mesh.segments(pos,:); 62 62 elseif md.dim==3 63 pressureload_layer1=[md. segments(pos,1:2) md.segments(pos,2)+md.numberofnodes2d md.segments(pos,1)+md.numberofnodes2d md.segments(pos,3)];63 pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.numberofnodes2d md.mesh.segments(pos,1)+md.numberofnodes2d md.mesh.segments(pos,3)]; 64 64 pressureload=[]; 65 65 for i=1:md.numlayers-1, -
issm/trunk/src/m/utils/Mesh/argusmesh.m
r9451 r9714 91 91 92 92 %Add segments and nodes on boundary 93 md. segments=findsegments(md);94 md. nodeonboundary=zeros(md.numberofnodes,1);95 md. nodeonboundary(md.segments(:,1))=1;96 md. nodeonboundary(md.segments(:,2))=1;93 md.mesh.segments=findsegments(md); 94 md.mesh.vertexonboundary=zeros(md.numberofnodes,1); 95 md.mesh.vertexonboundary(md.mesh.segments(:,1))=1; 96 md.mesh.vertexonboundary(md.mesh.segments(:,2))=1; -
issm/trunk/src/m/utils/Mesh/squaremesh.m
r9604 r9714 60 60 md.z=zeros(nods,1); 61 61 md.numberofnodes=nods; 62 md. nodeonboundary=zeros(nods,1);md.nodeonboundary(segments(:,1:2))=1;62 md.mesh.vertexonboundary=zeros(nods,1);md.mesh.vertexonboundary(segments(:,1:2))=1; 63 63 md.nodeonbed=ones(nods,1); 64 64 md.nodeonsurface=ones(nods,1); … … 66 66 %plug elements 67 67 md.elements=index; 68 md. segments=segments;68 md.mesh.segments=segments; 69 69 md.numberofelements=nel; 70 70 md.elementonbed=ones(nel,1); -
issm/trunk/template
r9702 r9714 2 2 3 3 projection 4 lat5 long6 hemisphere7 8 4 x 9 5 y … … 34 30 nodeconnectivity 35 31 elementconnectivity 36 connectivity -> rename average_vertex_connectivity37 38 segments39 segmentmarkers40 nodeonboundary41 42 extractednodes43 extractedelements44 32 }}} 45 33 -
issm/trunk/test/NightlyRun/test131.m
r9703 r9714 3 3 md=parameterize(md,'../Par/SquareShelfConstrained.par'); 4 4 %Add boundary conditions on thickness on the border 5 pos=find(md. nodeonboundary);5 pos=find(md.mesh.vertexonboundary); 6 6 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 7 7 md=setflowequation(md,'macayeal','all'); -
issm/trunk/test/NightlyRun/test132.m
r9703 r9714 3 3 md=parameterize(md,'../Par/SquareShelfConstrained.par'); 4 4 %Add boundary conditions on thickness on the border 5 pos=find(md. nodeonboundary);5 pos=find(md.mesh.vertexonboundary); 6 6 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 7 7 md=setflowequation(md,'macayeal','all'); -
issm/trunk/test/NightlyRun/test133.m
r9703 r9714 4 4 md=extrude(md,5,1); 5 5 %Add boundary conditions on thickness on the border 6 pos=find(md. nodeonboundary);6 pos=find(md.mesh.vertexonboundary); 7 7 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 8 8 md=setflowequation(md,'macayeal','all'); -
issm/trunk/test/NightlyRun/test134.m
r9703 r9714 4 4 md=extrude(md,5,1); 5 5 %Add boundary conditions on thickness on the border 6 pos=find(md. nodeonboundary);6 pos=find(md.mesh.vertexonboundary); 7 7 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 8 8 md=setflowequation(md,'macayeal','all'); -
issm/trunk/test/NightlyRun/test233.m
r9703 r9714 54 54 nodeonicefront=zeros(md.numberofnodes,1); 55 55 pos=find(md.y==ymax); nodeonicefront(pos)=1; 56 pos=find(nodeonicefront(md. segments(:,1)) | nodeonicefront(md.segments(:,2))); diagnostic.icefront=md.segments(pos,:);56 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); diagnostic.icefront=md.mesh.segments(pos,:); 57 57 diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]; 58 58 md.diagnostic.icefront=diagnostic.icefront; -
issm/trunk/test/NightlyRun/test234.m
r9703 r9714 54 54 nodeonicefront=zeros(md.numberofnodes,1); 55 55 pos=find(md.y==ymax); nodeonicefront(pos)=1; 56 pos=find(nodeonicefront(md. segments(:,1)) | nodeonicefront(md.segments(:,2))); diagnostic.icefront=md.segments(pos,:);56 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); diagnostic.icefront=md.mesh.segments(pos,:); 57 57 diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]; 58 58 md.diagnostic.icefront=diagnostic.icefront; -
issm/trunk/test/NightlyRun/test235.m
r9691 r9714 51 51 nodeonicefront=zeros(md.numberofnodes,1); 52 52 pos=find(md.y==ymax); nodeonicefront(pos)=1; 53 pos=find(nodeonicefront(md. segments(:,1)) | nodeonicefront(md.segments(:,2))); diagnostic.icefront=md.segments(pos,:);53 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); diagnostic.icefront=md.mesh.segments(pos,:); 54 54 diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]; 55 55 md.diagnostic.icefront=diagnostic.icefront; -
issm/trunk/test/NightlyRun/test236.m
r9691 r9714 51 51 nodeonicefront=zeros(md.numberofnodes,1); 52 52 pos=find(md.y==ymax); nodeonicefront(pos)=1; 53 pos=find(nodeonicefront(md. segments(:,1)) | nodeonicefront(md.segments(:,2))); diagnostic.icefront=md.segments(pos,:);53 pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); diagnostic.icefront=md.mesh.segments(pos,:); 54 54 diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]; 55 55 md.diagnostic.icefront=diagnostic.icefront; -
issm/trunk/test/NightlyRun/test527.m
r9691 r9714 28 28 %refine existing mesh 3 29 29 hVertices=NaN*ones(md.numberofnodes,1); 30 hVertices(find(md. nodeonboundary))=500;30 hVertices(find(md.mesh.vertexonboundary))=500; 31 31 md2=bamg(md,'metric',md.miscellaneous.dummy,'hmin',1000,'hmax',20000,'gradation',3,'geometricalmetric',1,'anisomax',1,'hVertices',hVertices); 32 32 x5=md2.x; -
issm/trunk/test/NightlyRun/test625.m
r9703 r9714 7 7 %Ice sheet only 8 8 md=modelextract(md,md.mask.elementongroundedice); 9 pos=find(md. nodeonboundary);9 pos=find(md.mesh.vertexonboundary); 10 10 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 11 11 -
issm/trunk/test/NightlyRun/test626.m
r9703 r9714 7 7 %Ice sheet only 8 8 md=modelextract(md,md.mask.elementongroundedice); 9 pos=find(md. nodeonboundary);9 pos=find(md.mesh.vertexonboundary); 10 10 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 11 11 -
issm/trunk/test/Par/79North.par
r9702 r9714 44 44 %Boundary conditions: 45 45 md=SetMarineIceSheetBC(md); 46 pos=find(md. nodeonboundary);46 pos=find(md.mesh.vertexonboundary); 47 47 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); 48 48 md.prognostic.spcthickness(pos)=md.geometry.thickness(pos); -
issm/trunk/test/Par/RoundSheetShelf.par
r9702 r9714 83 83 md.diagnostic.spcvy(pos)=0; 84 84 85 pressureload=md. segments;85 pressureload=md.mesh.segments; 86 86 pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end)) + 0*md.mask.elementongroundedice(pressureload(:,end))]; 87 87 md.diagnostic.icefront=pressureload;
Note:
See TracChangeset
for help on using the changeset viewer.