Changeset 9734
- Timestamp:
- 09/09/11 15:10:56 (14 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 added
- 3 deleted
- 134 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Inputs.cpp
r9320 r9734 412 412 /*}}}*/ 413 413 /*FUNCTION Inputs::AXPY{{{1*/ 414 void Inputs::AXPY(int YEnum, double scalar, intXEnum){414 void Inputs::AXPY(int MeshYEnum, double scalar, int MeshXEnum){ 415 415 416 416 Input* xinput=NULL; … … 418 418 419 419 /*Find x and y inputs: */ 420 xinput=(Input*)this->GetInput( XEnum);421 yinput=(Input*)this->GetInput( YEnum);420 xinput=(Input*)this->GetInput(MeshXEnum); 421 yinput=(Input*)this->GetInput(MeshYEnum); 422 422 423 423 /*some checks: */ 424 if(!xinput) _error_(" input %s could not be found!",EnumToStringx( XEnum));425 if(!yinput) _error_(" input %s could not be found!",EnumToStringx( YEnum));424 if(!xinput) _error_(" input %s could not be found!",EnumToStringx(MeshXEnum)); 425 if(!yinput) _error_(" input %s could not be found!",EnumToStringx(MeshYEnum)); 426 426 427 427 /*Apply AXPY: */ -
issm/trunk/src/c/Container/Inputs.h
r8363 r9734 37 37 Input* GetInput(int enum_name); 38 38 Inputs* SpawnTriaInputs(int* indices); 39 void AXPY(int YEnum, double scalar, intXEnum);39 void AXPY(int MeshYEnum, double scalar, int MeshXEnum); 40 40 double InfinityNorm(int enumtype); 41 41 double Max(int enumtype); -
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
r9733 r9734 165 165 MeshElementsEnum, 166 166 MeshEdgesEnum, 167 MeshYEnum, 168 MeshZEnum, 169 MeshXEnum, 167 170 /*}}}*/ 168 171 /*Datasets {{{1*/ … … 473 476 TransientInputEnum, 474 477 /*Temporary*/ 475 YEnum,476 ZEnum,477 XEnum,478 478 OutputfilenameEnum, 479 479 WaterfractionEnum, -
issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
r9733 r9734 169 169 case MeshElementsEnum : return "MeshElements"; 170 170 case MeshEdgesEnum : return "MeshEdges"; 171 case MeshYEnum : return "MeshY"; 172 case MeshZEnum : return "MeshZ"; 173 case MeshXEnum : return "MeshX"; 171 174 case ConstraintsEnum : return "Constraints"; 172 175 case LoadsEnum : return "Loads"; … … 416 419 case DragCoefficientAbsGradientEnum : return "DragCoefficientAbsGradient"; 417 420 case TransientInputEnum : return "TransientInput"; 418 case YEnum : return "Y";419 case ZEnum : return "Z";420 case XEnum : return "X";421 421 case OutputfilenameEnum : return "Outputfilename"; 422 422 case WaterfractionEnum : return "Waterfraction"; -
issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
r9733 r9734 72 72 73 73 /*Fetch data:*/ 74 iomodel->FetchData(5, XEnum,YEnum,ZEnum,BedEnum,ThicknessEnum);74 iomodel->FetchData(5,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum); 75 75 76 76 for (i=0;i<numberofvertices;i++){ … … 85 85 86 86 /*Free data: */ 87 iomodel->DeleteData(5, XEnum,YEnum,ZEnum,BedEnum,ThicknessEnum);87 iomodel->DeleteData(5,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum); 88 88 89 89 /*Assign output pointer: */ -
issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
r9725 r9734 69 69 iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum); 70 70 iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum); 71 iomodel->FetchData(&z,NULL,NULL, ZEnum);71 iomodel->FetchData(&z,NULL,NULL,MeshZEnum); 72 72 73 73 /*Initialize counter: */ -
issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
r9733 r9734 167 167 else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; 168 168 else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum; 169 else if (strcmp(name,"MeshY")==0) return MeshYEnum; 170 else if (strcmp(name,"MeshZ")==0) return MeshZEnum; 171 else if (strcmp(name,"MeshX")==0) return MeshXEnum; 169 172 else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; 170 173 else if (strcmp(name,"Loads")==0) return LoadsEnum; … … 414 417 else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; 415 418 else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; 416 else if (strcmp(name,"Y")==0) return YEnum;417 else if (strcmp(name,"Z")==0) return ZEnum;418 else if (strcmp(name,"X")==0) return XEnum;419 419 else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum; 420 420 else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; -
issm/trunk/src/c/objects/Vertex.cpp
r9405 r9734 31 31 Vertex::Vertex(int vertex_id, int vertex_sid,int i, IoModel* iomodel){ 32 32 33 this->Init(vertex_id, vertex_sid, iomodel->Data( XEnum)[i],iomodel->Data(YEnum)[i],iomodel->Data(ZEnum)[i],(iomodel->Data(ZEnum)[i]-iomodel->Data(BedEnum)[i])/(iomodel->Data(ThicknessEnum)[i]));33 this->Init(vertex_id, vertex_sid, iomodel->Data(MeshXEnum)[i],iomodel->Data(MeshYEnum)[i],iomodel->Data(MeshZEnum)[i],(iomodel->Data(MeshZEnum)[i]-iomodel->Data(BedEnum)[i])/(iomodel->Data(ThicknessEnum)[i])); 34 34 35 35 } -
issm/trunk/src/m/classes/mesh.m
r9733 r9734 37 37 elements = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',2); 38 38 edges = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3); 39 x = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1); 40 y = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1); 41 z = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1); 39 42 end 40 43 methods -
issm/trunk/src/m/classes/model/model.m
r9733 r9734 41 41 autodiff = modelfield('default',0,'marshall',true); 42 42 mesh = modelfield('default',0,'marshall',true); 43 44 %FIXME: all other fields should belong to other classes45 x = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);46 y = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);47 z = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',1);48 43 49 44 %}}} … … 265 260 if isfield(structmd,'elements'), md.mesh.elements=structmd.mesh.elements; end 266 261 if isfield(structmd,'edges'), md.mesh.edges=structmd.mesh.edges; end 262 if isfield(structmd,'y'), md.mesh.y=structmd.y; end 263 if isfield(structmd,'x'), md.mesh.x=structmd.x; end 264 if isfield(structmd,'z'), md.mesh.z=structmd.z; end 267 265 268 266 %Field changes -
issm/trunk/src/m/model/BasinConstrain.m
r9733 r9734 33 33 end 34 34 %ok, flag elements and nodes 35 [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md. x,md.y,domain,'element and node',2);35 [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); 36 36 end 37 37 if invert, -
issm/trunk/src/m/model/BasinConstrainShelf.m
r9733 r9734 33 33 end 34 34 %ok, flag elements and nodes 35 [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md. x,md.y,domain,'element and node',2);35 [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); 36 36 end 37 37 if invert, -
issm/trunk/src/m/model/DepthAverage.m
r9725 r9734 17 17 vector_average=zeros(md.mesh.numberofvertices2d,1); 18 18 for i=1:md.mesh.numberoflayers-1, 19 vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md. z,i+1)-project2d(md,md.z,i));19 vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i)); 20 20 end 21 21 vector_average=vector_average./project2d(md,md.geometry.thickness,1); … … 25 25 vector_average=zeros(md.mesh.numberofelements2d,1); 26 26 for i=1:md.mesh.numberoflayers-1, 27 vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md. z,i+1)-project2d(md,md.z,i));27 vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i)); 28 28 end 29 29 vector_average=vector_average./project2d(md,md.geometry.thickness,1); -
issm/trunk/src/m/model/MeltingGroundingLines.m
r9733 r9734 17 17 18 18 %search the node on ice sheet the closest to i 19 [d posd]=min(sqrt((md. x(pos_shelf(i))-md.x(pos_GL)).^2+(md.y(pos_shelf(i))-md.y(pos_GL)).^2));19 [d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2)); 20 20 21 21 if d<distance, -
issm/trunk/src/m/model/PropagateFlagsUntilDistance.m
r9733 r9734 13 13 %make 3d work in 2d: 14 14 if md.mesh.dimension==3, 15 md. x=md.mesh.x2d;16 md. y=md.mesh.y2d;15 md.mesh.x=md.mesh.x2d; 16 md.mesh.y=md.mesh.y2d; 17 17 md.mesh.elements=md.mesh.elements2d; 18 18 end … … 27 27 28 28 %average x and y over elements: 29 x_elem=md. x(md.mesh.elements)*[1;1;1]/3;30 y_elem=md. y(md.mesh.elements)*[1;1;1]/3;29 x_elem=md.mesh.x(md.mesh.elements)*[1;1;1]/3; 30 y_elem=md.mesh.y(md.mesh.elements)*[1;1;1]/3; 31 31 32 32 while 1, -
issm/trunk/src/m/model/SectionValues.m
r9733 r9734 82 82 83 83 %Interpolation of data on specified points 84 data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md. x,md.y,data,X,Y);85 %data_interp=griddata(md. x,md.y,data,X,Y);84 data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y); 85 %data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y); 86 86 87 87 %Compute index … … 120 120 121 121 %Interpolation of data on specified points 122 data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md. x,md.y,md.z,data,X3,Y3,Z3,NaN);122 data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,NaN); 123 123 124 124 %build outputs -
issm/trunk/src/m/model/ThicknessCorrection.m
r9733 r9734 48 48 49 49 %search the node on ice sheet the closest to i 50 [d posd]=min(sqrt((md. x(pos_shelf(i))-md.x(pos_GL)).^2+(md.y(pos_shelf(i))-md.y(pos_GL)).^2));50 [d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2)); 51 51 52 52 if d>distance, -
issm/trunk/src/m/model/averaging.m
r9733 r9734 33 33 if md.mesh.dimension==3 34 34 rep=6; 35 areas=GetAreas(index,md. x,md.y,md.z);35 areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z); 36 36 else 37 37 rep=3; 38 areas=GetAreas(index,md. x,md.y);38 areas=GetAreas(index,md.mesh.x,md.mesh.y); 39 39 end 40 40 summation=1/rep*ones(rep,1); -
issm/trunk/src/m/model/bamg.m
r9733 r9734 270 270 bamg_mesh=bamgmesh(md.private.bamg.mesh); 271 271 else 272 bamg_mesh.Vertices=[md. x md.y ones(md.mesh.numberofvertices,1)];272 bamg_mesh.Vertices=[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]; 273 273 bamg_mesh.Triangles=[md.mesh.elements ones(md.mesh.numberofelements,1)]; 274 274 end … … 317 317 md.private.bamg.mesh=bamgmesh(bamgmesh_out); 318 318 md.private.bamg.geometry=bamggeom(bamggeom_out); 319 md. x=bamgmesh_out.Vertices(:,1);320 md. y=bamgmesh_out.Vertices(:,2);319 md.mesh.x=bamgmesh_out.Vertices(:,1); 320 md.mesh.y=bamgmesh_out.Vertices(:,2); 321 321 md.mesh.elements=bamgmesh_out.Triangles(:,1:3); 322 322 md.mesh.edges=bamgmesh_out.IssmEdges; … … 327 327 md.mesh.dimension=2; 328 328 md.mesh.numberofelements=size(md.mesh.elements,1); 329 md.mesh.numberofvertices=length(md. x);329 md.mesh.numberofvertices=length(md.mesh.x); 330 330 md.mesh.numberofedges=size(md.mesh.edges,1); 331 md. z=zeros(md.mesh.numberofvertices,1);331 md.mesh.z=zeros(md.mesh.numberofvertices,1); 332 332 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); 333 333 md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); -
issm/trunk/src/m/model/basevert.m
r9733 r9734 14 14 15 15 for n=1:md.mesh.numberofelements 16 X=inv([md. x(md.mesh.elements(n,:)) md.y(md.mesh.elements(n,:)) ones(3,1)]);16 X=inv([md.mesh.x(md.mesh.elements(n,:)) md.mesh.y(md.mesh.elements(n,:)) ones(3,1)]); 17 17 alpha(n,:)=X(1,:); 18 18 beta(n,:)=X(2,:); -
issm/trunk/src/m/model/bedslope.m
r9733 r9734 10 10 numberofnodes=md.mesh.numberofvertices; 11 11 index=md.mesh.elements; 12 x=md. x; y=md.y;12 x=md.mesh.x; y=md.mesh.y; 13 13 else 14 14 numberofelements=md.mesh.numberofelements2d; -
issm/trunk/src/m/model/collapse.m
r9733 r9734 90 90 91 91 %Initialize with the 2d mesh 92 md. x=md.mesh.x2d;93 md. y=md.mesh.y2d;94 md. z=zeros(size(md.mesh.x2d));92 md.mesh.x=md.mesh.x2d; 93 md.mesh.y=md.mesh.y2d; 94 md.mesh.z=zeros(size(md.mesh.x2d)); 95 95 md.mesh.numberofvertices=md.mesh.numberofvertices2d; 96 96 md.mesh.numberofelements=md.mesh.numberofelements2d; -
issm/trunk/src/m/model/contourenvelope.m
r9733 r9734 47 47 if isfile, 48 48 %get flag list of elements and nodes inside the contour 49 nodein=ContourToMesh(md.mesh.elements,md. x,md.y,file,'node',1);49 nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1); 50 50 elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2)); 51 51 %modify element connectivity -
issm/trunk/src/m/model/contourmassbalance.m
r9733 r9734 22 22 23 23 %get flag list of elements and nodes inside the contour 24 nodein=ContourToMesh(md.mesh.elements,md. x,md.y,file,'node',1);24 nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1); 25 25 elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2)); 26 26 27 27 %conputing Mass flux 28 x=md. x;29 y=md. y;28 x=md.mesh.x; 29 y=md.mesh.y; 30 30 vx=mean(md.initialization.vx(segments(:,1:end-1)),2); 31 31 vy=mean(md.initialization.vy(segments(:,1:end-1)),2); … … 36 36 flux = - md.materials.rho_ice*sum(L.*H.*(vx.*nx+vy.*ny)); %outflux is negative! 37 37 disp(['mass outflux on ' file ' = ' num2str(-flux/10^9) ' Gt/yr']); 38 areas=GetAreas(md.mesh.elements,md. x,md.y);38 areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); 39 39 dhdt=flux/(sum(areas(find(elemin)))*md.materials.rho_ice); 40 40 disp(['dhdt on ' file ' (Flux method) = ' num2str(dhdt) ' m/yr']); -
issm/trunk/src/m/model/divergence.m
r9733 r9734 9 9 numberofnodes=md.mesh.numberofvertices; 10 10 index=md.mesh.elements; 11 x=md. x; y=md.y; z=md.z;11 x=md.mesh.x; y=md.mesh.y; z=md.mesh.z; 12 12 else 13 13 numberofelements=md.mesh.numberofelements2d; -
issm/trunk/src/m/model/extrude.m
r9733 r9734 72 72 %Create the new layers 73 73 for i=1:numlayers, 74 x3d=[x3d; md. x];75 y3d=[y3d; md. y];74 x3d=[x3d; md.mesh.x]; 75 y3d=[y3d; md.mesh.y]; 76 76 %nodes are distributed between bed and surface accordingly to the given exponent 77 77 z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; … … 103 103 104 104 %Save old mesh 105 md.mesh.x2d=md. x;106 md.mesh.y2d=md. y;105 md.mesh.x2d=md.mesh.x; 106 md.mesh.y2d=md.mesh.y; 107 107 md.mesh.elements2d=md.mesh.elements; 108 108 md.mesh.numberofelements2d=md.mesh.numberofelements; … … 114 114 %Build global 3d mesh 115 115 md.mesh.elements=elements3d; 116 md. x=x3d;117 md. y=y3d;118 md. z=z3d;116 md.mesh.x=x3d; 117 md.mesh.y=y3d; 118 md.mesh.z=z3d; 119 119 md.mesh.numberofelements=number_el3d; 120 120 md.mesh.numberofvertices=number_nodes3d; … … 213 213 %Put lithostatic pressure is there is an existing pressure 214 214 if ~isnan(md.initialization.pressure), 215 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md. z);215 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); 216 216 end 217 217 -
issm/trunk/src/m/model/mechanicalproperties.m
r9733 r9734 34 34 35 35 %compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma 36 [alpha beta]=GetNodalFunctionsCoeff(index,md. x,md.y);36 [alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y); 37 37 38 38 %compute shear -
issm/trunk/src/m/model/mesh/meshadaptation.m
r9733 r9734 41 41 42 42 %get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma 43 [alpha beta]=GetNodalFunctionsCoeff(index,md. x,md.y);44 areas=GetAreas(index,md. x,md.y);43 [alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y); 44 areas=GetAreas(index,md.mesh.x,md.mesh.y); 45 45 46 46 %Compute gradient for each element -
issm/trunk/src/m/model/mesh/meshbamg.m
r9733 r9734 77 77 disp(' interpolating velocities...'); 78 78 if strcmpi(NamesV.interp,'node'), 79 vx_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md. x,md.y,0);80 vy_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md. x,md.y,0);79 vx_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.mesh.x,md.mesh.y,0); 80 vy_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.mesh.x,md.mesh.y,0); 81 81 else 82 vx_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md. x,md.y,0);83 vy_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md. x,md.y,0);82 vx_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.mesh.x,md.mesh.y,0); 83 vy_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.mesh.x,md.mesh.y,0); 84 84 end 85 85 field=sqrt(vx_obs.^2+vy_obs.^2); … … 88 88 disp(' interpolating thickness...'); 89 89 if strcmpi(NamesH.interp,'node'), 90 h=InterpFromGridToMesh(Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md. x,md.y,0);90 h=InterpFromGridToMesh(Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.mesh.x,md.mesh.y,0); 91 91 else 92 h=InterpFromMeshToMesh2d(Thi.(NamesH.indexname),Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md. x,md.y,0);92 h=InterpFromMeshToMesh2d(Thi.(NamesH.indexname),Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.mesh.x,md.mesh.y,0); 93 93 end 94 94 field=[field h]; … … 97 97 %set nodeonwater field 98 98 if ~strcmp(groundeddomain,'N/A'), 99 nodeground=ContourToMesh(md.mesh.elements,md. x,md.y,groundeddomain,'node',2);99 nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); 100 100 md.nodeonwater=ones(md.mesh.numberofvertices,1); 101 101 md.nodeonwater(find(nodeground))=0; -
issm/trunk/src/m/model/mesh/meshconvert.m
r9733 r9734 12 12 13 13 if nargin==1, 14 x=md. x;15 y=md. y;14 x=md.mesh.x; 15 y=md.mesh.y; 16 16 index=md.mesh.elements; 17 17 else … … 28 28 md.private.bamg.mesh=bamgmesh(bamgmesh_out); 29 29 md.private.bamg.geometry=bamggeom(bamggeom_out); 30 md. x=bamgmesh_out.Vertices(:,1);31 md. y=bamgmesh_out.Vertices(:,2);30 md.mesh.x=bamgmesh_out.Vertices(:,1); 31 md.mesh.y=bamgmesh_out.Vertices(:,2); 32 32 md.mesh.elements=bamgmesh_out.Triangles(:,1:3); 33 33 md.mesh.edges=bamgmesh_out.IssmEdges; … … 38 38 md.mesh.dimension=2; 39 39 md.mesh.numberofelements=size(md.mesh.elements,1); 40 md.mesh.numberofvertices=length(md. x);40 md.mesh.numberofvertices=length(md.mesh.x); 41 41 md.mesh.numberofedges=size(md.mesh.edges,1); 42 md. z=zeros(md.mesh.numberofvertices,1);42 md.mesh.z=zeros(md.mesh.numberofvertices,1); 43 43 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); 44 44 md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); -
issm/trunk/src/m/model/mesh/meshexprefine.m
r9733 r9734 35 35 36 36 %Read domainame file into a matlab array (x,y): 37 refinearea=ContourToMesh(md.mesh.elements,md. x,md.y,domainname,'element',1);38 aires=GetAreas(md.mesh.elements,md. x,md.y);37 refinearea=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,domainname,'element',1); 38 aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); 39 39 40 40 %flags areas within the domain -
issm/trunk/src/m/model/mesh/meshnodensity.m
r9733 r9734 25 25 %Mesh using TriMeshNoDensity 26 26 if strcmp(riftname,''), 27 [md.mesh.elements,md. x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshNoDensity(domainname);27 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshNoDensity(domainname); 28 28 else 29 29 [elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname); … … 47 47 48 48 %plug into md 49 md. x=x;50 md. y=y;49 md.mesh.x=x; 50 md.mesh.y=y; 51 51 md.mesh.elements=elements; 52 52 md.mesh.segments=segments; … … 56 56 %Fill in rest of fields: 57 57 md.mesh.numberofelements=length(md.mesh.elements); 58 md.mesh.numberofvertices=length(md. x);59 md. z=zeros(md.mesh.numberofvertices,1);58 md.mesh.numberofvertices=length(md.mesh.x); 59 md.mesh.z=zeros(md.mesh.numberofvertices,1); 60 60 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 61 61 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -
issm/trunk/src/m/model/mesh/meshrefine.m
r9733 r9734 16 16 17 17 %Refine using TriMeshRefine 18 [md.mesh.elements,md. x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshRefine(md.mesh.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers,areas,'yes');18 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshRefine(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,areas,'yes'); 19 19 20 20 %Fill in rest of fields: 21 21 md.mesh.numberofelements=length(md.mesh.elements); 22 md.mesh.numberofvertices=length(md. x);23 md. z=zeros(md.mesh.numberofvertices,1);22 md.mesh.numberofvertices=length(md.mesh.x); 23 md.mesh.z=zeros(md.mesh.numberofvertices,1); 24 24 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 25 25 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -
issm/trunk/src/m/model/mesh/meshyams.m
r9733 r9734 71 71 disp(' interpolating velocities...'); 72 72 if strcmpi(Names.interp,'node'), 73 vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,0);74 vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,0);73 vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); 74 vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); 75 75 else 76 vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,0);77 vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,0);76 vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); 77 vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); 78 78 end 79 79 field=sqrt(vx_obs.^2+vy_obs.^2); … … 81 81 %set nodeonwater field 82 82 if ~strcmp(groundeddomain,'N/A'), 83 nodeground=ContourToMesh(md.mesh.elements,md. x,md.y,groundeddomain,'node',2);83 nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); 84 84 md.nodeonwater=ones(md.mesh.numberofvertices,1); 85 85 md.nodeonwater(find(nodeground))=0; … … 114 114 115 115 %Fill in rest of fields: 116 md. z=zeros(md.mesh.numberofvertices,1);116 md.mesh.z=zeros(md.mesh.numberofvertices,1); 117 117 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); 118 118 md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1); … … 120 120 md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); 121 121 if ~strcmp(groundeddomain,'N/A'), 122 nodeground=ContourToMesh(md.mesh.elements,md. x,md.y,groundeddomain,'node',2);122 nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); 123 123 md.nodeonwater=ones(md.mesh.numberofvertices,1); 124 124 md.nodeonwater(find(nodeground))=0; … … 127 127 end 128 128 if strcmpi(Names.interp,'node'), 129 md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,0);130 md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,0);129 md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); 130 md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); 131 131 else 132 md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,0);133 md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,0);132 md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); 133 md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); 134 134 end 135 135 md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2); … … 145 145 146 146 %build normals and lengths of segments: 147 lengths=sqrt((md. x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 );148 normalsx=cos(atan2((md. x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));149 normalsy=sin(atan2((md. x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));147 lengths=sqrt((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))).^2 + (md.mesh.y(rift.segments(:,1))-md.mesh.y(rift.segments(:,2))).^2 ); 148 normalsx=cos(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1))))); 149 normalsy=sin(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1))))); 150 150 151 151 %ok, build penaltypairs: -
issm/trunk/src/m/model/mesh/reorder.m
r9733 r9734 22 22 md.mesh.elements=tnewnodes(md.mesh.elements(newelements,:)); 23 23 md.mesh.segments=[tnewnodes(md.mesh.segments(:,1)) tnewnodes(md.mesh.segments(:,2)) tnewelements(md.mesh.segments(:,3))]; 24 md. x=md.x(newnodes);25 md. y=md.y(newnodes);26 md. z=md.z(newnodes);24 md.mesh.x=md.mesh.x(newnodes); 25 md.mesh.y=md.mesh.y(newnodes); 26 md.mesh.z=md.mesh.z(newnodes); 27 27 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; -
issm/trunk/src/m/model/mesh/rifts/meshaddrifts.m
r9733 r9734 56 56 57 57 %plug md2 mesh into md mesh: 58 [md.mesh.elements,md. x,md.y,md.z,md.mesh.numberofelements,md.mesh.numberofvertices,elconv,nodeconv,elconv2,nodeconv2]=meshplug(md.mesh.elements,md.x,md.y,md.z,...58 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,md.mesh.numberofelements,md.mesh.numberofvertices,elconv,nodeconv,elconv2,nodeconv2]=meshplug(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,... 59 59 md2.mesh.elements,md2.x,md2.y,md2.z,md2.extractednodes,md2.extractedelements,domain_index); 60 60 -
issm/trunk/src/m/model/mesh/rifts/meshprocessoutsiderifts.m
r9733 r9734 12 12 13 13 %first, flag nodes that belong to the domain outline 14 flags=ContourToMesh(md.mesh.elements,md. x,md.y,domainoutline,'node',0);14 flags=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,domainoutline,'node',0); 15 15 16 16 rift=md.rifts.riftstruct(i); … … 55 55 %take the list of elements on one side of the rift that connect to the tip, 56 56 %and duplicate the tip on them, so as to open the rift to the outside. 57 num=length(md. x)+1;58 md. x=[md.x;md.x(tip)];59 md. y=[md.y;md.y(tip)];57 num=length(md.mesh.x)+1; 58 md.mesh.x=[md.mesh.x;md.mesh.x(tip)]; 59 md.mesh.y=[md.mesh.y;md.mesh.y(tip)]; 60 60 md.mesh.numberofvertices=num; 61 61 … … 84 84 %Fill in rest of fields: 85 85 md.mesh.numberofelements=length(md.mesh.elements); 86 md.mesh.numberofvertices=length(md. x);87 md. z=zeros(md.mesh.numberofvertices,1);88 md.mesh.vertexonboundary=zeros(length(md. x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;86 md.mesh.numberofvertices=length(md.mesh.x); 87 md.mesh.z=zeros(md.mesh.numberofvertices,1); 88 md.mesh.vertexonboundary=zeros(length(md.mesh.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.mesh.numberofelements,1); -
issm/trunk/src/m/model/mesh/rifts/meshprocessrifts.m
r9733 r9734 25 25 26 26 %Call MEX file 27 [md.mesh.elements,md. x,md.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers);27 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers); 28 28 md.rifts 29 29 md.rifts.riftstruct … … 34 34 %Fill in rest of fields: 35 35 md.mesh.numberofelements=length(md.mesh.elements); 36 md.mesh.numberofvertices=length(md. x);37 md. z=zeros(md.mesh.numberofvertices,1);38 md.mesh.vertexonboundary=zeros(length(md. x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;36 md.mesh.numberofvertices=length(md.mesh.x); 37 md.mesh.z=zeros(md.mesh.numberofvertices,1); 38 md.mesh.vertexonboundary=zeros(length(md.mesh.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 39 39 md.rifts.numrifts=length(md.rifts.riftstruct); 40 40 md.flowequation.element_equation=3*ones(md.mesh.numberofelements,1); … … 46 46 %get coordinates of rift tips 47 47 for i=1:md.rifts.numrifts, 48 md.rifts.riftstruct(i).tip1coordinates=[md. x(md.rifts.riftstruct(i).tips(1)) md.y(md.rifts.riftstruct(i).tips(1))];49 md.rifts.riftstruct(i).tip2coordinates=[md. x(md.rifts.riftstruct(i).tips(2)) md.y(md.rifts.riftstruct(i).tips(2))];48 md.rifts.riftstruct(i).tip1coordinates=[md.mesh.x(md.rifts.riftstruct(i).tips(1)) md.mesh.y(md.rifts.riftstruct(i).tips(1))]; 49 md.rifts.riftstruct(i).tip2coordinates=[md.mesh.x(md.rifts.riftstruct(i).tips(2)) md.mesh.y(md.rifts.riftstruct(i).tips(2))]; 50 50 end 51 51 52 52 %In case we have rifts that open up the domain outline, we need to open them: 53 flags=ContourToMesh(md.mesh.elements,md. x,md.y,domainoutline,'node',0);53 flags=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,domainoutline,'node',0); 54 54 found=0; 55 55 for i=1:md.rifts.numrifts, … … 68 68 69 69 %get elements that are not correctly oriented in the correct direction: 70 aires=GetAreas(md.mesh.elements,md. x,md.y);70 aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); 71 71 pos=find(aires<0); 72 72 md.mesh.elements(pos,:)=[md.mesh.elements(pos,2) md.mesh.elements(pos,1) md.mesh.elements(pos,3)]; -
issm/trunk/src/m/model/mesh/rifts/meshyamsrecreateriftsegments.m
r9714 r9734 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.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));13 pos=find_point(md.mesh.x(md.mesh.segments(:,1)),md.mesh.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 14 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));15 pos=find_point(md.mesh.x(md.mesh.segments(:,1)),md.mesh.y(md.mesh.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2)); 16 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.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));19 pos=find_point(md.mesh.x(md.mesh.segments(:,1)),md.mesh.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 20 20 pos_record=[pos_record; pos]; 21 21 riftsegs=md.mesh.segments(pos,:); … … 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.mesh.segments(:,1)),md.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));41 pos1=find_point(md.mesh.x(md.mesh.segments(:,1)),md.mesh.y(md.mesh.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2)); 42 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));43 pos2=find_point(md.mesh.x(md.mesh.segments(:,1)),md.mesh.y(md.mesh.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2)); 44 44 tip2=md.mesh.segments(pos2,1); 45 45 if length(tip1)==2, -
issm/trunk/src/m/model/mesh/setmesh.m
r9733 r9734 39 39 %Mesh using TriMesh 40 40 if strcmp(riftname,''), 41 [md.mesh.elements,md. x,md.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,area,'yes');41 [md.mesh.elements,md.mesh.x,md.mesh.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'); … … 61 61 62 62 %plug into md 63 md. x=x;64 md. y=y;63 md.mesh.x=x; 64 md.mesh.y=y; 65 65 md.mesh.elements=elements; 66 66 md.mesh.segments=segments; … … 70 70 %Fill in rest of fields: 71 71 md.mesh.numberofelements=length(md.mesh.elements); 72 md.mesh.numberofvertices=length(md. x);73 md. z=zeros(md.mesh.numberofvertices,1);72 md.mesh.numberofvertices=length(md.mesh.x); 73 md.mesh.z=zeros(md.mesh.numberofvertices,1); 74 74 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; 75 75 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -
issm/trunk/src/m/model/misfit.m
r9733 r9734 11 11 if md.mesh.dimension==2, 12 12 elements=md.mesh.elements; 13 x=md. x;14 y=md. y;13 x=md.mesh.x; 14 y=md.mesh.y; 15 15 vx=md.initialization.vx; 16 16 vy=md.initialization.vy; -
issm/trunk/src/m/model/outflow.m
r9725 r9734 7 7 A=md.mesh.segments(:,1); 8 8 B=md.mesh.segments(:,2); 9 Nx=-(md. y(A)-md.y(B));10 Ny= md. x(A)-md.x(B);9 Nx=-(md.mesh.y(A)-md.mesh.y(B)); 10 Ny= md.mesh.x(A)-md.mesh.x(B); 11 11 Vx=(md.initialization.vx(A)+md.initialization.vx(B))/2; 12 12 Vy=(md.initialization.vy(A)+md.initialization.vy(B))/2; -
issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m
r9733 r9734 17 17 18 18 md.mesh.elements=md.mesh.elements2d; 19 md. x=md.mesh.x2d;20 md. y=md.mesh.y2d;19 md.mesh.x=md.mesh.x2d; 20 md.mesh.y=md.mesh.y2d; 21 21 md.mesh.numberofvertices=md.mesh.numberofvertices2d; 22 22 md.mesh.numberofelements=md.mesh.numberofelements2d; -
issm/trunk/src/m/model/partition/adjacency.m
r9733 r9734 16 16 17 17 %now, build vwgt: 18 areas=GetAreas(md.mesh.elements,md. x,md.y);18 areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); 19 19 20 20 %get node connectivity -
issm/trunk/src/m/model/partition/partitioner.m
r9733 r9734 35 35 md3d=md; %save for later 36 36 md.mesh.elements=md.mesh.elements2d; 37 md. x=md.mesh.x2d;38 md. y=md.mesh.y2d;37 md.mesh.x=md.mesh.x2d; 38 md.mesh.y=md.mesh.y2d; 39 39 md.mesh.numberofvertices=md.mesh.numberofvertices2d; 40 40 md.mesh.numberofelements=md.mesh.numberofelements2d; … … 70 70 71 71 % partition into nparts 72 part=Chaco(md.qmu.adjacency,weights,[],md. x, md.y ,md.z,method,npart,[])'+1; %index partitions from 1 up. like metis.72 part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y ,md.mesh.z,method,npart,[])'+1; %index partitions from 1 up. like metis. 73 73 74 74 elseif strcmpi(package,'scotch'), -
issm/trunk/src/m/model/plot/applyoptions.m
r9733 r9734 424 424 [mdx mdy]=basinzoom(options); 425 425 else 426 mdx=[min(md. x) max(md.x)];427 mdy=[min(md. y) max(md.y)];426 mdx=[min(md.mesh.x) max(md.mesh.x)]; 427 mdy=[min(md.mesh.y) max(md.mesh.y)]; 428 428 end 429 429 line(A.x,A.y,ones(size(A.x)),'color','b'); … … 437 437 %flag edges of a partition 438 438 if exist(options,'partitionedges') 439 [xsegments ysegments]=flagedges(md.mesh.elements,md. x,md.y,md.qmu.partition);439 [xsegments ysegments]=flagedges(md.mesh.elements,md.mesh.x,md.mesh.y,md.qmu.partition); 440 440 xsegments=xsegments*getfieldvalue(options,'unit',1); 441 441 ysegments=ysegments*getfieldvalue(options,'unit',1); -
issm/trunk/src/m/model/plot/checkplotoptions.m
r9641 r9734 180 180 if strcmpi(getfieldvalue(options,'northarrow'),'on') 181 181 %default values 182 Lx=max(md. y)-min(md.y);183 Ly=max(md. y)-min(md.y);184 %default values 185 options=changefieldvalue(options,'northarrow',[min(md. x)+1/6*Lx min(md.y)+5/6*Ly 1/15*Ly 0.25 1/250*Ly]);182 Lx=max(md.mesh.y)-min(md.mesh.y); 183 Ly=max(md.mesh.y)-min(md.mesh.y); 184 %default values 185 options=changefieldvalue(options,'northarrow',[min(md.mesh.x)+1/6*Lx min(md.mesh.y)+5/6*Ly 1/15*Ly 0.25 1/250*Ly]); 186 186 end 187 187 end … … 191 191 if strcmpi(getfieldvalue(options,'scaleruler'),'on') 192 192 %default values 193 Lx=max(md. y)-min(md.y);194 Ly=max(md. y)-min(md.y);195 %default values 196 options=changefieldvalue(options,'scaleruler',[min(md. x)+6/8*Lx min(md.y)+1/10*Ly 10^(ceil(log10(Lx)))/5 floor(Lx/100) 5]);193 Lx=max(md.mesh.y)-min(md.mesh.y); 194 Ly=max(md.mesh.y)-min(md.mesh.y); 195 %default values 196 options=changefieldvalue(options,'scaleruler',[min(md.mesh.x)+6/8*Lx min(md.mesh.y)+1/10*Ly 10^(ceil(log10(Lx)))/5 floor(Lx/100) 5]); 197 197 end 198 198 end -
issm/trunk/src/m/model/plot/plot_BC.m
r9679 r9734 7 7 8 8 %plot dirichlets 9 h1=plot(md. x(find(~isnan(md.diagnostic.spcvx))),md.y(find(~isnan(md.diagnostic.spcvx))),'ro','MarkerSize',14,'MarkerFaceColor','r');10 h2=plot(md. x(find(~isnan(md.diagnostic.spcvy))),md.y(find(~isnan(md.diagnostic.spcvy))),'bo','MarkerSize',10,'MarkerFaceColor','b');11 h3=plot(md. x(find(~isnan(md.diagnostic.spcvz))),md.y(find(~isnan(md.diagnostic.spcvz))),'yo','MarkerSize',6 ,'MarkerFaceColor','y');9 h1=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvx))),md.mesh.y(find(~isnan(md.diagnostic.spcvx))),'ro','MarkerSize',14,'MarkerFaceColor','r'); 10 h2=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvy))),md.mesh.y(find(~isnan(md.diagnostic.spcvy))),'bo','MarkerSize',10,'MarkerFaceColor','b'); 11 h3=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvz))),md.mesh.y(find(~isnan(md.diagnostic.spcvz))),'yo','MarkerSize',6 ,'MarkerFaceColor','y'); 12 12 13 13 %update legend -
issm/trunk/src/m/model/plot/plot_quiver.m
r9684 r9734 8 8 % 9 9 % Example: 10 % plot_quiver(md. x,md.y,md.initialization.vx,md.initialization.vy,options);10 % plot_quiver(md.mesh.x,md.mesh.y,md.initialization.vx,md.initialization.vy,options); 11 11 12 12 %process fields -
issm/trunk/src/m/model/plot/plot_quiver3.m
r9684 r9734 8 8 % 9 9 % Example: 10 % plot_quiver(md. x,md.y,md.z,md.initialization.vx,md.initialization.vy,md.initialization.vz,options);10 % plot_quiver(md.mesh.x,md.mesh.y,md.mesh.z,md.initialization.vx,md.initialization.vy,md.initialization.vz,options); 11 11 12 12 %keep only non NaN elements -
issm/trunk/src/m/model/plot/plot_sarpwr.m
r9620 r9734 13 13 if exist(options,'unit'), 14 14 unit=getfieldvalue(options,'unit'); 15 md. x=md.x*unit;16 md. y=md.y*unit;17 md. z=md.z*unit;15 md.mesh.x=md.mesh.x*unit; 16 md.mesh.y=md.mesh.y*unit; 17 md.mesh.z=md.mesh.z*unit; 18 18 end 19 19 -
issm/trunk/src/m/model/plot/plot_section.m
r9733 r9734 32 32 md3d=md; 33 33 if exist(options,'layer') 34 md. x=md.mesh.x2d; md.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.mesh.dimension=2;34 md.mesh.x=md.mesh.x2d; md.mesh.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.mesh.dimension=2; 35 35 end 36 36 … … 96 96 text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) 97 97 plot(x,y,'-r') 98 axis([min(md. x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])98 axis([min(md.mesh.x)-1 max(md.mesh.x)+1 min(md.mesh.y)-1 max(md.mesh.y)+1]) 99 99 view(2) 100 100 end … … 138 138 text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) 139 139 plot(x,y,'-r') 140 axis([min(md. x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])140 axis([min(md.mesh.x)-1 max(md.mesh.x)+1 min(md.mesh.y)-1 max(md.mesh.y)+1]) 141 141 view(2) 142 142 end … … 173 173 text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8]) 174 174 plot(x,y,'-r') 175 axis([min(md. x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])175 axis([min(md.mesh.x)-1 max(md.mesh.x)+1 min(md.mesh.y)-1 max(md.mesh.y)+1]) 176 176 view(2) 177 177 end -
issm/trunk/src/m/model/plot/plot_tensor_principalaxis.m
r9733 r9734 29 29 %take the center of each element if ~isonnode 30 30 if datatype==1, 31 x=mean(md. x(md.mesh.elements'))'; y=mean(md.y(md.mesh.elements'))'; z=mean(md.z(md.mesh.elements'))';31 x=mean(md.mesh.x(md.mesh.elements'))'; y=mean(md.mesh.y(md.mesh.elements'))'; z=mean(md.mesh.z(md.mesh.elements'))'; 32 32 end 33 33 -
issm/trunk/src/m/model/plot/processmesh.m
r9733 r9734 16 16 17 17 if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon'), 18 x=md. x;18 x=md.mesh.x; 19 19 x2d=md.mesh.x2d; 20 y=md. y;20 y=md.mesh.y; 21 21 y2d=md.mesh.y2d; 22 22 else … … 27 27 end 28 28 29 z_field=getfieldvalue(options,'z',md. z);29 z_field=getfieldvalue(options,'z',md.mesh.z); 30 30 if ischar(z_field), 31 31 z=md.(z_field); … … 33 33 z=z_field; 34 34 else 35 z=md. z;35 z=md.mesh.z; 36 36 end 37 37 … … 68 68 if (md.mesh.dimension==2), 69 69 elements=transpose(reshape(1:3*md.mesh.numberofelements,3,md.mesh.numberofelements)); 70 x=transpose(reshape(md. x(data.index)',1,3*md.mesh.numberofelements));71 y=transpose(reshape(md. y(data.index)',1,3*md.mesh.numberofelements));70 x=transpose(reshape(md.mesh.x(data.index)',1,3*md.mesh.numberofelements)); 71 y=transpose(reshape(md.mesh.y(data.index)',1,3*md.mesh.numberofelements)); 72 72 z=zeros(3*md.mesh.numberofelements,1); 73 73 is2d=1; 74 74 else 75 75 elements=transpose(reshape(1:6*md.mesh.numberofelements,6,md.mesh.numberofelements)); 76 x=transpose(reshape(md. x(data.index)',1,6*md.mesh.numberofelements));77 y=transpose(reshape(md. y(data.index)',1,6*md.mesh.numberofelements));78 z=transpose(reshape(md. z(data.index)',1,6*md.mesh.numberofelements));76 x=transpose(reshape(md.mesh.x(data.index)',1,6*md.mesh.numberofelements)); 77 y=transpose(reshape(md.mesh.y(data.index)',1,6*md.mesh.numberofelements)); 78 z=transpose(reshape(md.mesh.z(data.index)',1,6*md.mesh.numberofelements)); 79 79 is2d=0; 80 80 end -
issm/trunk/src/m/model/plot/showregion.m
r9714 r9734 26 26 Ay=[min(A.y) max(A.y)]; 27 27 28 mdx=[min(md. x) max(md.x)];29 mdy=[min(md. y) max(md.y)];28 mdx=[min(md.mesh.x) max(md.mesh.x)]; 29 mdy=[min(md.mesh.y) max(md.mesh.y)]; 30 30 31 31 line(A.x,A.y,'color','b'); -
issm/trunk/src/m/model/radarpower.m
r9714 r9734 19 19 20 20 highres=getfieldvalue(options,'highres',0); 21 xlim=getfieldvalue(options,'xlim',[min(md. x) max(md.x)]);22 ylim=getfieldvalue(options,'ylim',[min(md. y) max(md.y)]);21 xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]); 22 ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]); 23 23 posting=getfieldvalue(options,'posting',0); % 0 -> image posting default 24 24 -
issm/trunk/src/m/model/setmask.m
r9733 r9734 21 21 22 22 %Get assigned fields 23 x=md. x;24 y=md. y;23 x=md.mesh.x; 24 y=md.mesh.y; 25 25 elements=md.mesh.elements; 26 26 -
issm/trunk/src/m/model/setmask2.m
r9733 r9734 9 9 10 10 %Get assigned fields 11 x=md. x;12 y=md. y;11 x=md.mesh.x; 12 y=md.mesh.y; 13 13 elements=md.mesh.elements; 14 14 -
issm/trunk/src/m/model/shear2d.m
r9733 r9734 8 8 % s=shear2d(md); 9 9 10 [alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md. x,md.y);10 [alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 11 11 12 12 summation=[1;1;1]; -
issm/trunk/src/m/model/slope.m
r9733 r9734 10 10 numberofnodes=md.mesh.numberofvertices; 11 11 index=md.mesh.elements; 12 x=md. x; y=md.y;12 x=md.mesh.x; y=md.mesh.y; 13 13 else 14 14 numberofelements=md.mesh.numberofelements2d; -
issm/trunk/src/m/model/thicknessevolution.m
r9733 r9734 20 20 21 21 %compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma 22 [alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md. x,md.y);22 [alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 23 23 24 24 %compute dhdt=div(Hu) -
issm/trunk/src/m/utils/BC/SetIceShelfBC.m
r9733 r9734 19 19 icefrontfile=varargin{1}; 20 20 if ~exist(icefrontfile), error(['SetIceShelfBC error message: ice front file ' icefrontfile ' not found']); end 21 nodeinsideicefront=ContourToMesh(md.mesh.elements,md. x,md.y,icefrontfile,'node',2);21 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2); 22 22 nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront); 23 23 elseif nargin==1, -
issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m
r9733 r9734 23 23 error(['SetMarineIceSheetBC error message: ice front file ' icefrontfile ' not found']); 24 24 end 25 nodeinsideicefront=ContourToMesh(md.mesh.elements,md. x,md.y,icefrontfile,'node',2);25 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2); 26 26 vertexonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront); 27 27 else -
issm/trunk/src/m/utils/Exp/clicktoflowline.m
r9733 r9734 4 4 % Usage: clicktoflowline(index,x,y,u,v,x0,y0,filename) 5 5 % 6 % Ex: clicktoflowline(md.mesh.elements,md. x,md.y,md.inversion.vx_obs,md.inversion.vy_obs,'flowline1.exp')6 % Ex: clicktoflowline(md.mesh.elements,md.mesh.x,md.mesh.y,md.inversion.vx_obs,md.inversion.vy_obs,'flowline1.exp') 7 7 8 8 -
issm/trunk/src/m/utils/Exp/downstreamflowlines.m
r9733 r9734 10 10 % 11 11 % Example: 12 % flowpath=downstreamflowlines(md.mesh.elements,md. x,md.y,md.vx,md.initialization.vy,x0,y0)12 % flowpath=downstreamflowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.vy,x0,y0) 13 13 14 14 %check input size -
issm/trunk/src/m/utils/Exp/flowlines.m
r9733 r9734 10 10 % 11 11 % Example: 12 % flowpath=flowlines(md.mesh.elements,md. x,md.y,md.vx,md.initialization.vy,x0,y0)12 % flowpath=flowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.vy,x0,y0) 13 13 14 14 %check input size -
issm/trunk/src/m/utils/Geometry/FlagElements.m
r9733 r9734 33 33 if ~exist(region,'file'), 34 34 [xlim,ylim]=basinzoom('basin',region); 35 flag_nodes=double(md. x<xlim(2) & md.x>xlim(1) & md.y<ylim(2) & md.y>ylim(1));35 flag_nodes=double(md.mesh.x<xlim(2) & md.mesh.x>xlim(1) & md.mesh.y<ylim(2) & md.mesh.y>ylim(1)); 36 36 flag=prod(flag_nodes(md.mesh.elements),2); 37 37 else 38 38 %ok, flag elements 39 flag=ContourToMesh(md.mesh.elements(:,1:3),md. x,md.y,region,'element',1);39 flag=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,region,'element',1); 40 40 end 41 41 end -
issm/trunk/src/m/utils/Interp/FillHole.m
r9733 r9734 6 6 % 7 7 % Example: 8 % md.geometry.surface=FillHole(md.mesh.elements,x,md. y,md.geometry.surface)8 % md.geometry.surface=FillHole(md.mesh.elements,x,md.mesh.y,md.geometry.surface) 9 9 % 10 10 -
issm/trunk/src/m/utils/Interp/InterpFromFile.m
r9691 r9734 19 19 % 20 20 % Example: 21 % md.geometry.surface=InterpFromFile(md. x,md.y,'surfacefile.mat',0);21 % md.geometry.surface=InterpFromFile(md.mesh.x,md.mesh.y,'surfacefile.mat',0); 22 22 % 23 23 % See also: PLUGVELOCITIES, INTERPFROMGRID, INTERPFROMMESH2D, INTERPFROMMESH3D -
issm/trunk/src/m/utils/Interp/plugvelocities.m
r9684 r9734 29 29 %Interpolation 30 30 if strcmpi(Names.interp,'node'), 31 md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,default_value);32 md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,default_value);31 md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,default_value); 32 md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,default_value); 33 33 else 34 md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md. x,md.y,default_value);35 md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md. x,md.y,default_value);34 md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,default_value); 35 md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,default_value); 36 36 end 37 37 -
issm/trunk/src/m/utils/Mesh/BamgCall.m
r9733 r9734 20 20 %Compute Hessian 21 21 t1=clock; fprintf('%s',' computing Hessian...'); 22 hessian=ComputeHessian(md.mesh.elements,md. x,md.y,field,'node');22 hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); 23 23 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); 24 24 … … 50 50 %Vertices 51 51 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); 52 fprintf(fid,'%8g %8g %i\n',[md. x md.y ones(md.mesh.numberofvertices,1)]');52 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); 53 53 54 54 %Triangles … … 68 68 t1=clock; fprintf('\n%s',' reading final mesh files...'); 69 69 A=meshread('carre1.mesh'); 70 md. x=A.x;71 md. y=A.y;70 md.mesh.x=A.x; 71 md.mesh.y=A.y; 72 72 md.z=zeros(A.nods,1); 73 73 md.mesh.elements=A.index; -
issm/trunk/src/m/utils/Mesh/BamgCallFromMetric.m
r9733 r9734 31 31 %Vertices 32 32 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); 33 fprintf(fid,'%8g %8g %i\n',[md. x md.y ones(md.mesh.numberofvertices,1)]');33 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); 34 34 35 35 %Triangles … … 49 49 t1=clock; fprintf('\n%s',' reading final mesh files...'); 50 50 A=meshread('carre1.mesh'); 51 md. x=A.x;52 md. y=A.y;51 md.mesh.x=A.x; 52 md.mesh.y=A.y; 53 53 md.z=zeros(A.nods,1); 54 54 md.mesh.elements=A.index; -
issm/trunk/src/m/utils/Mesh/ComputeHessian.m
r9733 r9734 10 10 % 11 11 % Example: 12 % hessian=ComputeHessian(md.mesh.elements,md. x,md.y,md.inversion.vel_obs,'node')12 % hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,md.inversion.vel_obs,'node') 13 13 14 14 %some variables -
issm/trunk/src/m/utils/Mesh/GetAreas.m
r9733 r9734 10 10 % 11 11 % Examples: 12 % areas =GetAreas(md.mesh.elements,md. x,md.y);13 % volumes=GetAreas(md.mesh.elements,md. x,md.y,md.z);12 % areas =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); 13 % volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z); 14 14 15 15 %get number of elements and number of nodes -
issm/trunk/src/m/utils/Mesh/GetCharacteristicLength.m
r9733 r9734 9 9 % 10 10 % Examples: 11 % length =GetCharacteristicLength(md.mesh.elements,md. x,md.y);12 % length =GetCharacteristicLength(md.mesh.elements,md. x,md.y,md.z);11 % length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y); 12 % length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y,md.z); 13 13 14 14 -
issm/trunk/src/m/utils/Mesh/GetNodalFunctionsCoeff.m
r9733 r9734 12 12 % 13 13 % Example: 14 % [alpha beta gamma]=GetNodalFunctionsCoeff(md.mesh.elements,md. x,md.y);14 % [alpha beta gamma]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 15 15 16 16 %make columns out of x and y -
issm/trunk/src/m/utils/Mesh/MeshQuality.m
r9733 r9734 7 7 %Get some variables from the model 8 8 index=md.mesh.elements; 9 x=md. x;10 y=md. y;9 x=md.mesh.x; 10 y=md.mesh.y; 11 11 12 12 %2d geometric parameter (do not change) -
issm/trunk/src/m/utils/Mesh/ProfileProjectOntoMesh.m
r9733 r9734 9 9 %make a curve out of the mesh, to use the intersections routine. 10 10 rows=[md.mesh.elements md.mesh.elements(:,1)]'; rows=rows(:); 11 x=md. x(rows);12 y=md. y(rows);11 x=md.mesh.x(rows); 12 y=md.mesh.y(rows); 13 13 14 14 %[x0,y0] = intersections(profile.x,profile.y,x,y,1); … … 48 48 49 49 %now, for each node, figure out which element it belongs to. 50 node_in_element=NodeInElement(newx,newy,md.mesh.elements,md. x,md.y,md.mesh.vertexconnectivity);50 node_in_element=NodeInElement(newx,newy,md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.vertexconnectivity); 51 51 52 52 % eliminate nodes that don't fall in any element -
issm/trunk/src/m/utils/Mesh/YamsCall.m
r9733 r9734 20 20 %Compute Hessian 21 21 t1=clock; fprintf('%s',' computing Hessian...'); 22 hessian=ComputeHessian(md.mesh.elements,md. x,md.y,field,'node');22 hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); 23 23 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); 24 24 … … 47 47 %Vertices 48 48 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); 49 fprintf(fid,'%8g %8g %i\n',[md. x md.y zeros(md.mesh.numberofvertices,1)]');49 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]'); 50 50 51 51 %Triangles … … 88 88 Tria=load('carre1.tria'); 89 89 Coor=load('carre1.coor'); 90 md. x=Coor(:,1);91 md. y=Coor(:,2);90 md.mesh.x=Coor(:,1); 91 md.mesh.y=Coor(:,2); 92 92 md.z=zeros(size(Coor,1),1); 93 93 md.mesh.elements=Tria; -
issm/trunk/src/m/utils/Mesh/argusmesh.m
r9733 r9734 78 78 %Finally, use model constructor to build a complete model: 79 79 md.mesh.elements=elements; 80 md. x=x;81 md. y=y;80 md.mesh.x=x; 81 md.mesh.y=y; 82 82 md.z=z; 83 md.mesh.numberofvertices=size(md. x,1);83 md.mesh.numberofvertices=size(md.mesh.x,1); 84 84 md.mesh.numberofelements=size(md.mesh.elements,1); 85 85 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -
issm/trunk/src/m/utils/Mesh/rifttipsonmesh.m
r5819 r9734 14 14 y_tip=rift.y(1); 15 15 16 index=find_point(md. x,md.y,x_tip,y_tip);16 index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); 17 17 tips(end+1)=index; 18 18 … … 20 20 y_tip=rift.y(end); 21 21 22 index=find_point(md. x,md.y,x_tip,y_tip);22 index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); 23 23 tips(end+1)=index; 24 24 -
issm/trunk/src/m/utils/Mesh/roundmesh.m
r9703 r9734 27 27 28 28 %move the closest node to the center 29 [mini pos]=min(md. x.^2+md.y.^2);30 md. x(pos)=0;31 md. y(pos)=0;29 [mini pos]=min(md.mesh.x.^2+md.mesh.y.^2); 30 md.mesh.x(pos)=0; 31 md.mesh.y(pos)=0; 32 32 33 33 %delete domain -
issm/trunk/src/m/utils/Mesh/squaremesh.m
r9733 r9734 56 56 57 57 %plug coordinates and nodes 58 md. x=x;59 md. y=y;58 md.mesh.x=x; 59 md.mesh.y=y; 60 60 md.z=zeros(nods,1); 61 61 md.mesh.numberofvertices=nods; -
issm/trunk/src/m/utils/Numerics/cfl_step.m
r9733 r9734 16 16 17 17 index=md.mesh.elements; 18 edgex=max(md. x(index),[],2)-min(md.x(index),[],2);19 edgey=max(md. y(index),[],2)-min(md.y(index),[],2);18 edgex=max(md.mesh.x(index),[],2)-min(md.mesh.x(index),[],2); 19 edgey=max(md.mesh.y(index),[],2)-min(md.mesh.y(index),[],2); 20 20 vx=max(abs(vx(index)),[],2); 21 21 vy=max(abs(vy(index)),[],2); -
issm/trunk/template
r9729 r9734 1 mesh {{{2 3 projection4 x5 y6 z7 elements8 edges9 10 elements2d11 x2d12 y2d13 }}}14 15 1 %To be completed 16 2 cluster{{{ -
issm/trunk/test/Miscellaneous/Bump/Bump.par
r9733 r9734 5 5 hmin=1000; 6 6 hmax=1000; 7 ymin=min(md. y);8 ymax=max(md. y);9 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);7 ymin=min(md.mesh.y); 8 ymax=max(md.mesh.y); 9 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 10 10 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+10; 11 11 md.geometry.surface=md.geometry.bed+md.geometry.thickness; 12 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+height*exp(-((md. x-50000).^2+(md.y-50000).^2)/(4000)^2)+10;12 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+height*exp(-((md.mesh.x-50000).^2+(md.mesh.y-50000).^2)/(4000)^2)+10; 13 13 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 14 14 -
issm/trunk/test/Miscellaneous/GJM_test1/SquareShelf.par
r9733 r9734 4 4 hmin=300; 5 5 hmax=1000; 6 ymin=min(md. y);7 ymax=max(md. y);8 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);6 ymin=min(md.mesh.y); 7 ymax=max(md.mesh.y); 8 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 9 9 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; 10 10 md.geometry.surface=md.geometry.bed+md.geometry.thickness; -
issm/trunk/test/Miscellaneous/connectivity/Square.par
r9733 r9734 9 9 hmin=300; 10 10 hmax=1000; 11 ymin=min(md. y);12 ymax=max(md. y);13 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);11 ymin=min(md.mesh.y); 12 ymax=max(md.mesh.y); 13 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 14 14 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; 15 15 md.geometry.surface=md.geometry.bed+md.geometry.thickness; -
issm/trunk/test/NightlyRun/test1101.m
r9729 r9734 30 30 31 31 %Create MPCs to have periodic boundary conditions 32 posx=find(md. x==0);33 posx2=find(md. x==max(md.x));32 posx=find(md.mesh.x==0); 33 posx2=find(md.mesh.x==max(md.mesh.x)); 34 34 35 posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same nodes two times36 posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x));35 posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times 36 posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); 37 37 38 38 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1102.m
r9733 r9734 17 17 18 18 % %Find elements at the corner and extract model 19 % posnodes=find((md. x==0 | md.x==max(md.x)) & (md.y==0 | md.y==max(md.y)));19 % posnodes=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y))); 20 20 % [a,b]=find(ismember(md.mesh.elements,posnodes)); 21 21 % elements=ones(md.mesh.numberofelements,1); … … 34 34 35 35 %Create MPCs to have periodic boundary conditions 36 %posx=find(md. x==0);37 %posx2=find(md. x==max(md.x));38 %posx=find(md. x==0 & md.y~=0 & md.y~=max(md.y) & ~md.mesh.vertexonbed);39 %posx2=find(md. x==max(md.x) & md.y~=0 & md.y~=max(md.y) & ~md.mesh.vertexonbed);36 %posx=find(md.mesh.x==0); 37 %posx2=find(md.mesh.x==max(md.mesh.x)); 38 %posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed); 39 %posx2=find(md.mesh.x==max(md.mesh.x) & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed); 40 40 41 %posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times42 %posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x) & ~md.mesh.vertexonbed);41 %posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times 42 %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); 43 43 44 44 %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1103.m
r9729 r9734 27 27 28 28 %Create MPCs to have periodic boundary conditions 29 posx=find(md. x==0);30 posx2=find(md. x==max(md.x));29 posx=find(md.mesh.x==0); 30 posx2=find(md.mesh.x==max(md.mesh.x)); 31 31 32 posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same nodes two times33 posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x));32 posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times 33 posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); 34 34 35 35 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1104.m
r9729 r9734 27 27 28 28 %Create MPCs to have periodic boundary conditions 29 posx=find(md. x==0);30 posx2=find(md. x==max(md.x));29 posx=find(md.mesh.x==0); 30 posx2=find(md.mesh.x==max(md.mesh.x)); 31 31 32 posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same nodes two times33 posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x));32 posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times 33 posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); 34 34 35 35 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1105.m
r9729 r9734 25 25 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); 26 26 27 posx=find(md. x==0 & md.y~=0 & md.y~=L);28 posx2=find(md. x==L & md.y~=0 & md.y~=L);27 posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=L); 28 posx2=find(md.mesh.x==L & md.mesh.y~=0 & md.mesh.y~=L); 29 29 30 posy=find(md. y==0 & md.x~=0 & md.x~=L); %Don't take the same nodes two times31 posy2=find(md. y==L & md.x~=0 & md.x~=L);30 posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=L); %Don't take the same nodes two times 31 posy2=find(md.mesh.y==L & md.mesh.x~=0 & md.mesh.x~=L); 32 32 33 33 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; 34 34 35 35 %Add spc on the corners 36 pos=find((md. x==0 | md.x==L) & (md.y==0 | md.y==L) & md.mesh.vertexonbed);36 pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed); 37 37 md.diagnostic.spcvx(pos)=0; 38 38 md.diagnostic.spcvy(pos)=0; -
issm/trunk/test/NightlyRun/test1106.m
r9725 r9734 19 19 %Create MPCs to have periodic boundary conditions 20 20 21 %posx=find(md. x==0);22 %posx2=find(md. x==L);21 %posx=find(md.mesh.x==0); 22 %posx2=find(md.mesh.x==L); 23 23 24 %posy=find(md. y==0 & md.x~=0 & md.x~=L); %Don't take the same nodes two times25 %posy2=find(md. y==L & md.x~=0 & md.x~=L);24 %posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=L); %Don't take the same nodes two times 25 %posy2=find(md.mesh.y==L & md.mesh.x~=0 & md.mesh.x~=L); 26 26 27 27 %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1107.m
r9729 r9734 26 26 27 27 %Create MPCs to have periodic boundary conditions 28 posx=find(md. x==0 & ~(md.y==0 & md.mesh.vertexonbed) & ~(md.y==L & md.mesh.vertexonbed));29 posx2=find(md. x==max(md.x) & ~(md.y==0 & md.mesh.vertexonbed) & ~(md.y==L & md.mesh.vertexonbed));28 posx=find(md.mesh.x==0 & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed)); 29 posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed)); 30 30 31 posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same nodes two times32 posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x));31 posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times 32 posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); 33 33 34 34 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; 35 35 36 36 %Add spc on the corners 37 pos=find((md. x==0 | md.x==L) & (md.y==0 | md.y==L) & md.mesh.vertexonbed);37 pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed); 38 38 md.diagnostic.spcvy(:)=0; 39 39 md.diagnostic.spcvx(pos)=0; -
issm/trunk/test/NightlyRun/test1108.m
r9725 r9734 22 22 %md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); 23 23 24 %pos=find((md. x==0 | md.x==max(md.x)) & (md.y==0 | md.y==max(md.y)));24 %pos=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y))); 25 25 %md.diagnostic.spcvx(pos)=0; 26 26 %md.diagnostic.spcvy(pos)=0; … … 28 28 29 29 %%Create MPCs to have periodic boundary conditions 30 %posx=find(md. x==0);31 %posx2=find(md. x==max(md.x));30 %posx=find(md.mesh.x==0); 31 %posx2=find(md.mesh.x==max(md.mesh.x)); 32 32 33 %posy=find(md. y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same nodes two times34 %posy2=find(md. y==max(md.y) & md.x~=0 & md.x~=max(md.x));33 %posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times 34 %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); 35 35 36 36 %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1109.m
r9729 r9734 26 26 27 27 %Create MPCs to have periodic boundary conditions 28 posx=find(md. x==0);29 posx2=find(md. x==max(md.x));28 posx=find(md.mesh.x==0); 29 posx2=find(md.mesh.x==max(md.mesh.x)); 30 30 md.diagnostic.vertex_pairing=[posx,posx2]; 31 31 … … 41 41 %Remove the spc where there is some sliding (case 3 and 4): 42 42 if i==3 | i==4, 43 pos=find(md. y/max(md.y)>=0.44 & md.y/max(md.y)<=0.5);43 pos=find(md.mesh.y/max(md.mesh.y)>=0.44 & md.mesh.y/max(md.mesh.y)<=0.5); 44 44 md.diagnostic.spcvx(pos)=NaN; 45 45 md.diagnostic.spcvy(pos)=NaN; -
issm/trunk/test/NightlyRun/test1110.m
r9729 r9734 30 30 md.diagnostic.spcvz(pos)=0; 31 31 else 32 pos=find(md.mesh.vertexonbed & (md. x==0 | md.x==max(md.x)) & (md.y==0 | md.y==max(md.y)));32 pos=find(md.mesh.vertexonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y))); 33 33 md.diagnostic.spcvx(pos)=100; %because we need a dirichlet somewhere 34 34 md.diagnostic.spcvy(pos)=0; … … 39 39 40 40 %Create MPCs to have periodic boundary conditions 41 posx=find(md. x==0);42 posx2=find(md. x==max(md.x));41 posx=find(md.mesh.x==0); 42 posx2=find(md.mesh.x==max(md.mesh.x)); 43 43 44 posy=find(md. y==0);45 posy2=find(md. y==max(md.y));44 posy=find(md.mesh.y==0); 45 posy2=find(md.mesh.y==max(md.mesh.y)); 46 46 47 47 md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; -
issm/trunk/test/NightlyRun/test1201.m
r9725 r9734 27 27 28 28 %spc thickness 29 pos=find(md. y>199999.9);29 pos=find(md.mesh.y>199999.9); 30 30 times=0:1:500; 31 31 md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices+1,length(times)); … … 48 48 [elements,x,y,z,s,h2]=SectionValues(md,results{2},'../Exp/CrossLineEISMINT.exp',100); 49 49 [elements,x,y,z,s,h3]=SectionValues(md,results{3},'../Exp/CrossLineEISMINT.exp',100); 50 [elements,x,y,z,s,hth]=SectionValues(md, 500+100*sin(2*pi/200*(500-md. y/400)),'../Exp/CrossLineEISMINT.exp',100);50 [elements,x,y,z,s,hth]=SectionValues(md, 500+100*sin(2*pi/200*(500-md.mesh.y/400)),'../Exp/CrossLineEISMINT.exp',100); 51 51 plot(s,h1,'r',s,h2,'b',s,h3,'g',s,hth,'k') 52 52 legend('Art. diff.','No Art. diff.','D.G.','Theoretical') -
issm/trunk/test/NightlyRun/test1203.m
r9703 r9734 10 10 11 11 %Impose a non zero velocity on the upper boundary condition (y=max(y)) 12 pos=find(md. y==max(md.y));13 md.diagnostic.spcvy(pos)=400*(((md. x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);12 pos=find(md.mesh.y==max(md.mesh.y)); 13 md.diagnostic.spcvy(pos)=400*(((md.mesh.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.mesh.x(pos)-100000)/25000).^2); 14 14 15 15 %Compute solution for MacAyeal's model -
issm/trunk/test/NightlyRun/test1204.m
r9703 r9734 10 10 11 11 %Impose a non zero velocity on the upper boundary condition (y=max(y)) 12 pos=find(md. y==max(md.y));13 md.diagnostic.spcvy(pos)=400*(((md. x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);12 pos=find(md.mesh.y==max(md.mesh.y)); 13 md.diagnostic.spcvy(pos)=400*(((md.mesh.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.mesh.x(pos)-100000)/25000).^2); 14 14 15 15 %Compute solution for MacAyeal's model -
issm/trunk/test/NightlyRun/test1205.m
r9733 r9734 13 13 %Calculation of the analytical 2d velocity field 14 14 constant=0.3; 15 vx_obs=constant/2*md. x.*(md.geometry.thickness).^-1;16 vy_obs=constant/2*md. y.*(md.geometry.thickness).^-1;15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; 16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; 17 17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); 18 18 … … 51 51 set(gcf,'Position',[1 1 1580 1150]) 52 52 subplot(2,2,1) 53 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...53 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 54 54 vel,'FaceColor','interp','EdgeColor','none'); 55 55 title('Modelled velocity','FontSize',14,'FontWeight','bold') … … 58 58 59 59 subplot(2,2,2) 60 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...60 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 61 61 vel_obs,'FaceColor','interp','EdgeColor','none'); 62 62 title('Analytical velocity','FontSize',14,'FontWeight','bold') … … 66 66 subplot(2,2,3) 67 67 hold on; 68 plot(sqrt((md. x(1:md.mesh.numberofvertices2d)).^2+(md.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.');69 plot(sqrt((md. x2d).^2+(md.y2d).^2),vel_obs,'b.');68 plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.'); 69 plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.'); 70 70 title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold'); 71 71 xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold'); … … 76 76 77 77 subplot(2,2,4) 78 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...78 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 79 79 abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none'); 80 80 title('Relative misfit [%]','FontSize',14,'FontWeight','bold') -
issm/trunk/test/NightlyRun/test1206.m
r9733 r9734 13 13 %Calculation of the analytical 2d velocity field 14 14 constant=0.3; 15 vx_obs=constant/2*md. x.*(md.geometry.thickness).^-1;16 vy_obs=constant/2*md. y.*(md.geometry.thickness).^-1;15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; 16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; 17 17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); 18 18 … … 50 50 figure(1) 51 51 subplot(2,2,1) 52 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...52 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 53 53 vel,'FaceColor','interp','EdgeColor','none'); 54 54 title('Modelled velocity','FontSize',14,'FontWeight','bold') … … 57 57 58 58 subplot(2,2,2) 59 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...59 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 60 60 vel_obs,'FaceColor','interp','EdgeColor','none'); 61 61 title('Analytical velocity','FontSize',14,'FontWeight','bold') … … 65 65 subplot(2,2,3) 66 66 hold on; 67 plot(sqrt((md. x(1:md.mesh.numberofvertices2d)).^2+(md.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.');68 plot(sqrt((md. x2d).^2+(md.y2d).^2),vel_obs,'b.');67 plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.'); 68 plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.'); 69 69 title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold'); 70 70 xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold'); … … 75 75 76 76 subplot(2,2,4) 77 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...77 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 78 78 abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none'); 79 79 title('Relative misfit [%]','FontSize',14,'FontWeight','bold') -
issm/trunk/test/NightlyRun/test1207.m
r9733 r9734 13 13 %Calculation of the analytical 2d velocity field 14 14 constant=0.3; 15 vx_obs=constant/2*md. x.*(md.geometry.thickness).^-1;16 vy_obs=constant/2*md. y.*(md.geometry.thickness).^-1;15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; 16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; 17 17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); 18 18 … … 50 50 figure(1) 51 51 subplot(2,2,1) 52 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...52 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 53 53 vel,'FaceColor','interp','EdgeColor','none'); 54 54 title('Modelled velocity','FontSize',14,'FontWeight','bold') … … 57 57 58 58 subplot(2,2,2) 59 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...59 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 60 60 vel_obs,'FaceColor','interp','EdgeColor','none'); 61 61 title('Analytical velocity','FontSize',14,'FontWeight','bold') … … 65 65 subplot(2,2,3) 66 66 hold on; 67 plot(sqrt((md. x(1:md.mesh.numberofvertices2d)).^2+(md.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.');68 plot(sqrt((md. x2d).^2+(md.y2d).^2),vel_obs,'b.');67 plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.'); 68 plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.'); 69 69 title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold'); 70 70 xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold'); … … 75 75 76 76 subplot(2,2,4) 77 p=patch('Faces',md.mesh.elements2d,'Vertices',[md. x2d md.y2d],'FaceVertexCData',...77 p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',... 78 78 abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none'); 79 79 title('Relative misfit [%]','FontSize',14,'FontWeight','bold') -
issm/trunk/test/NightlyRun/test1302.m
r9733 r9734 23 23 A=10/(exp(alpha*(-1000))-1); %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10 24 24 B=-A; 25 md.initialization.temperature=A*exp(alpha*md. z)+B;25 md.initialization.temperature=A*exp(alpha*md.mesh.z)+B; 26 26 27 27 %modeled results -
issm/trunk/test/NightlyRun/test1303.m
r9733 r9734 18 18 %d2T/dz2=0 T(bed)=10 T(surface)=0 => T=0*(z-bed)/thickness+10*(surface-z)/thickness 19 19 %each layer of the 3d mesh must have a constant value 20 md.initialization.temperature=10*(md.geometry.surface-md. z)./md.geometry.thickness;20 md.initialization.temperature=10*(md.geometry.surface-md.mesh.z)./md.geometry.thickness; 21 21 22 22 %modeled results -
issm/trunk/test/NightlyRun/test1304.m
r9733 r9734 19 19 %the result is linear with depth and is equal to 0 on the upper surface (See BC) 20 20 %d2T/dz2=0 -k*dT/dz(bed)=G T(surface)=0 => T=-G/k*(z-surface) 21 md.initialization.temperature=-0.1/md.thermalconductivity*(md. z-md.geometry.surface); %G=0.1 W/m221 md.initialization.temperature=-0.1/md.thermalconductivity*(md.mesh.z-md.geometry.surface); %G=0.1 W/m2 22 22 23 23 %modeled results -
issm/trunk/test/NightlyRun/test137.m
r9725 r9734 2 2 %Simple mesh 3 3 md=bamg(model,'domain','../Exp/Square.exp','hmax',100000); 4 x1=md. x;5 y1=md. y;4 x1=md.mesh.x; 5 y1=md.mesh.y; 6 6 7 7 %hVertices 8 8 md=bamg(model,'domain','../Exp/Square.exp','hmax',300000,'hvertices',[10000 100000 400000 100000]'); 9 x2=md. x;10 y2=md. y;9 x2=md.mesh.x; 10 y2=md.mesh.y; 11 11 12 12 %big mesh -
issm/trunk/test/NightlyRun/test1401.m
r9681 r9734 1 1 %test the anisotropic mesh adaptation 2 %function to capture = exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;2 %function to capture = exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 3 3 printingflag=false; 4 4 … … 11 11 %mesh adaptation loop YAMS 12 12 md=squaremesh(md,L,L,nx,ny); 13 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;13 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 14 14 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 15 15 if printingflag, … … 20 20 21 21 md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,1.3,10^-4); 22 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;22 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 23 23 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 24 24 if printingflag, … … 29 29 30 30 md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,2.5,0.008); 31 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;31 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 32 32 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 33 33 if printingflag, … … 36 36 system(['mv mesh1_yams3.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Mesh/ ']); 37 37 end 38 x1=md. x;39 y1=md. y;38 x1=md.mesh.x; 39 y1=md.mesh.y; 40 40 41 41 %mesh adaptation loop BAMG 42 42 md=squaremesh(md,L,L,nx,ny); 43 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;43 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 44 44 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 45 45 if printingflag, … … 51 51 md.bamg=NaN; 52 52 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10^-4); 53 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;53 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 54 54 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 55 55 if printingflag, … … 61 61 md.bamg=NaN; 62 62 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2.5,'err',0.008); 63 md.inversion.vel_obs=exp(-(sqrt((md. x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2;63 md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2; 64 64 plotmodel(md,'data','vel_obs','data','vel_obs','nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5); 65 65 if printingflag, … … 68 68 system(['mv mesh1_bamg3.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Mesh/ ']); 69 69 end 70 x2=md. x;71 y2=md. y;70 x2=md.mesh.x; 71 y2=md.mesh.y; 72 72 73 73 %Fields and tolerances to track changes -
issm/trunk/test/NightlyRun/test1402.m
r9681 r9734 10 10 %mesh adaptation loop YAMS 11 11 md=squaremesh(md,L,L,nx,ny); 12 u=4*md. x-2; v=4*md.y-2;12 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 13 13 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 14 14 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 22 22 23 23 md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,2.3,10^-2); 24 u=4*md. x-2; v=4*md.y-2;24 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 25 25 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 26 26 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 34 34 35 35 md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,3,0.005); 36 u=4*md. x-2; v=4*md.y-2;36 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 37 37 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 38 38 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 44 44 system(['mv mesh2_yams3.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Mesh/ ']); 45 45 end 46 x1=md. x;47 y1=md. y;46 x1=md.mesh.x; 47 y1=md.mesh.y; 48 48 49 49 %mesh adaptation loop BAMG 50 50 md=squaremesh(md,L,L,nx,ny); 51 u=4*md. x-2; v=4*md.y-2;51 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 52 52 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 53 53 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 62 62 md.bamg=NaN; 63 63 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2.3,'err',10^-2); 64 u=4*md. x-2; v=4*md.y-2;64 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 65 65 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 66 66 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 75 75 md.bamg=NaN; 76 76 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',3,'err',0.005); 77 u=4*md. x-2; v=4*md.y-2;77 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 78 78 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 79 79 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 88 88 md.bamg=NaN; 89 89 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',1.5,'err',0.003,'anisomax',1); 90 u=4*md. x-2; v=4*md.y-2;90 u=4*md.mesh.x-2; v=4*md.mesh.y-2; 91 91 md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... 92 92 +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... … … 98 98 system(['mv mesh2_bamgiso.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Mesh/ ']); 99 99 end 100 x2=md. x;101 y2=md. y;100 x2=md.mesh.x; 101 y2=md.mesh.y; 102 102 103 103 %Fields and tolerances to track changes -
issm/trunk/test/NightlyRun/test1501.m
r9733 r9734 97 97 98 98 index = md.mesh.elements; 99 x1=md. x(index(:,1)); x2=md.x(index(:,2)); x3=md.x(index(:,3));100 y1=md. y(index(:,1)); y2=md.y(index(:,2)); y3=md.y(index(:,3));99 x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3)); 100 y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3)); 101 101 areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))); 102 102 -
issm/trunk/test/NightlyRun/test1502.m
r9733 r9734 103 103 104 104 index = md.mesh.elements; 105 x1=md. x(index(:,1)); x2=md.x(index(:,2)); x3=md.x(index(:,3));106 y1=md. y(index(:,1)); y2=md.y(index(:,2)); y3=md.y(index(:,3));105 x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3)); 106 y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3)); 107 107 areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))); 108 108 -
issm/trunk/test/NightlyRun/test233.m
r9725 r9734 13 13 14 14 %needed later 15 ymin=min(md. y);16 ymax=max(md. y);17 xmin=min(md. x);18 xmax=max(md. x);15 ymin=min(md.mesh.y); 16 ymax=max(md.mesh.y); 17 xmin=min(md.mesh.x); 18 xmax=max(md.mesh.x); 19 19 20 20 di=md.materials.rho_ice/md.materials.rho_water; … … 42 42 43 43 %constrain flanks to 0 normal velocity 44 pos=find(md. x==xmin | md.x==xmax);44 pos=find(md.mesh.x==xmin | md.mesh.x==xmax); 45 45 md.diagnostic.spcvx(pos)=0; 46 46 md.diagnostic.spcvz(pos)=NaN; 47 47 48 48 %constrain grounding line to 0 velocity 49 pos=find(md. y==ymin);49 pos=find(md.mesh.y==ymin); 50 50 md.diagnostic.spcvx(pos)=0; 51 51 md.diagnostic.spcvy(pos)=0; … … 53 53 %icefront 54 54 nodeonicefront=zeros(md.mesh.numberofvertices,1); 55 pos=find(md. y==ymax); nodeonicefront(pos)=1;55 pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; 56 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))]; … … 62 62 %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292. 63 63 %ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3; 64 %vy_c=ey_c.*md. y*md.constants.yts;64 %vy_c=ey_c.*md.mesh.y*md.constants.yts; 65 65 66 66 %Fields and tolerances to track changes -
issm/trunk/test/NightlyRun/test234.m
r9725 r9734 13 13 14 14 %needed later 15 ymin=min(md. y);16 ymax=max(md. y);17 xmin=min(md. x);18 xmax=max(md. x);15 ymin=min(md.mesh.y); 16 ymax=max(md.mesh.y); 17 xmin=min(md.mesh.x); 18 xmax=max(md.mesh.x); 19 19 20 20 di=md.materials.rho_ice/md.materials.rho_water; … … 42 42 43 43 %constrain flanks to 0 normal velocity 44 pos=find(md. x==xmin | md.x==xmax);44 pos=find(md.mesh.x==xmin | md.mesh.x==xmax); 45 45 md.diagnostic.spcvx(pos)=0; 46 46 md.diagnostic.spcvz(pos)=NaN; 47 47 48 48 %constrain grounding line to 0 velocity 49 pos=find(md. y==ymin);49 pos=find(md.mesh.y==ymin); 50 50 md.diagnostic.spcvx(pos)=0; 51 51 md.diagnostic.spcvy(pos)=0; … … 53 53 %icefront 54 54 nodeonicefront=zeros(md.mesh.numberofvertices,1); 55 pos=find(md. y==ymax); nodeonicefront(pos)=1;55 pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; 56 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))]; … … 62 62 %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292. 63 63 %ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).^3; 64 %vy_c=ey_c.*md. y*md.constants.yts;64 %vy_c=ey_c.*md.mesh.y*md.constants.yts; 65 65 66 66 %Fields and tolerances to track changes -
issm/trunk/test/NightlyRun/test235.m
r9725 r9734 10 10 11 11 %needed later 12 ymin=min(md. y);13 ymax=max(md. y);14 xmin=min(md. x);15 xmax=max(md. x);12 ymin=min(md.mesh.y); 13 ymax=max(md.mesh.y); 14 xmin=min(md.mesh.x); 15 xmax=max(md.mesh.x); 16 16 17 17 di=md.materials.rho_ice/md.materials.rho_water; … … 39 39 40 40 %constrain flanks to 0 normal velocity 41 pos=find(md. x==xmin | md.x==xmax);41 pos=find(md.mesh.x==xmin | md.mesh.x==xmax); 42 42 md.diagnostic.spcvx(pos)=0; 43 43 md.diagnostic.spcvz(pos)=NaN; 44 44 45 45 %constrain grounding line to 0 velocity 46 pos=find(md. y==ymin);46 pos=find(md.mesh.y==ymin); 47 47 md.diagnostic.spcvx(pos)=0; 48 48 md.diagnostic.spcvy(pos)=0; … … 50 50 %icefront 51 51 nodeonicefront=zeros(md.mesh.numberofvertices,1); 52 pos=find(md. y==ymax); nodeonicefront(pos)=1;52 pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; 53 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))]; -
issm/trunk/test/NightlyRun/test236.m
r9725 r9734 10 10 11 11 %needed later 12 ymin=min(md. y);13 ymax=max(md. y);14 xmin=min(md. x);15 xmax=max(md. x);12 ymin=min(md.mesh.y); 13 ymax=max(md.mesh.y); 14 xmin=min(md.mesh.x); 15 xmax=max(md.mesh.x); 16 16 17 17 di=md.materials.rho_ice/md.materials.rho_water; … … 39 39 40 40 %constrain flanks to 0 normal velocity 41 pos=find(md. x==xmin | md.x==xmax);41 pos=find(md.mesh.x==xmin | md.mesh.x==xmax); 42 42 md.diagnostic.spcvx(pos)=0; 43 43 md.diagnostic.spcvz(pos)=NaN; 44 44 45 45 %constrain grounding line to 0 velocity 46 pos=find(md. y==ymin);46 pos=find(md.mesh.y==ymin); 47 47 md.diagnostic.spcvx(pos)=0; 48 48 md.diagnostic.spcvy(pos)=0; … … 50 50 %icefront 51 51 nodeonicefront=zeros(md.mesh.numberofvertices,1); 52 pos=find(md. y==ymax); nodeonicefront(pos)=1;52 pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; 53 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))]; -
issm/trunk/test/NightlyRun/test446.m
r9733 r9734 5 5 md=roundmesh(model,radius,resolution); 6 6 %fix center node to 0,0 7 rad=sqrt((md. x).*md.x+(md.y).*md.y);7 rad=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); 8 8 pos=find(rad==min(rad)); 9 md. x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center9 md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center 10 10 %}}} 11 11 %put our grounding line 'shelfextent' meters from the icefront {{{1 12 xelem=md. x(md.mesh.elements)*[1;1;1]/3;13 yelem=md. y(md.mesh.elements)*[1;1;1]/3;12 xelem=md.mesh.x(md.mesh.elements)*[1;1;1]/3; 13 yelem=md.mesh.y(md.mesh.elements)*[1;1;1]/3; 14 14 rad=sqrt((xelem).*xelem+(yelem).*yelem); 15 15 flags=zeros(md.mesh.numberofelements,1); … … 24 24 %plug holes into the ice sheet, to test for grounding line migration. {{{1 25 25 di=md.materials.rho_ice/md.materials.rho_water; 26 rad=sqrt((md. x).*md.x+(md.y).*md.y);26 rad=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); 27 27 pos=find(rad<200000); 28 28 md.geometry.thickness(pos)=100; … … 30 30 md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); 31 31 32 pos=find(md. x<.2*1e6 & md.x>-.2*1e6 & md.y>0);32 pos=find(md.mesh.x<.2*1e6 & md.mesh.x>-.2*1e6 & md.mesh.y>0); 33 33 md.geometry.thickness(pos)=100; 34 34 md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; 35 35 md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); 36 36 37 pos=find(md. x<.1*1e6 & md.x>-.1*1e6 & md.y<-.5*1e6 & md.y>-.6*1e6);37 pos=find(md.mesh.x<.1*1e6 & md.mesh.x>-.1*1e6 & md.mesh.y<-.5*1e6 & md.mesh.y>-.6*1e6); 38 38 md.geometry.thickness(pos)=100; 39 39 md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; -
issm/trunk/test/NightlyRun/test527.m
r9733 r9734 3 3 hVertices(1:5)=1000; 4 4 md=bamg(model,'domain','../Exp/Pig.exp','hmax',20000,'hVertices',hVertices,'gradation',3,'geometricalmetric',1); 5 x1=md. x;6 y1=md. y;5 x1=md.mesh.x; 6 y1=md.mesh.y; 7 7 8 8 %Simple mesh 2 … … 10 10 md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); 11 11 md=parameterize(md,'../Par/Pig.par'); 12 x2=md. x;13 y2=md. y;12 x2=md.mesh.x; 13 y2=md.mesh.y; 14 14 15 15 %refine existing mesh 1 16 hessian=ComputeHessian(md.mesh.elements,md. x,md.y,md.inversion.vy_obs,'node');16 hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,md.inversion.vy_obs,'node'); 17 17 metric=ComputeMetric(hessian,2/9,1,1000,25*10^3,[]); 18 18 md.miscellaneous.dummy=metric; -
issm/trunk/test/Par/79North.par
r9733 r9734 3 3 %Geometry and observation 4 4 load('../Data/79North.data','-mat'); 5 md.initialization.vx =InterpFromMeshToMesh2d(index,x,y,vx,md. x,md.y);6 md.initialization.vy =InterpFromMeshToMesh2d(index,x,y,vy,md. x,md.y);7 md.geometry.surface =InterpFromMeshToMesh2d(index,x,y,surface,md. x,md.y);8 md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md. x,md.y);5 md.initialization.vx =InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); 6 md.initialization.vy =InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); 7 md.geometry.surface =InterpFromMeshToMesh2d(index,x,y,surface,md.mesh.x,md.mesh.y); 8 md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y); 9 9 md.geometry.bed=md.geometry.surface-md.geometry.thickness; 10 10 clear surface thickness vx vy x y index; -
issm/trunk/test/Par/ISMIPA.par
r9733 r9734 2 2 3 3 disp(' creating thickness'); 4 md.geometry.surface=-md. x*tan(0.5*pi/180);5 md.geometry.bed=md.geometry.surface-1000+500*sin(md. x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x));4 md.geometry.surface=-md.mesh.x*tan(0.5*pi/180); 5 md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x)); 6 6 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 7 7 -
issm/trunk/test/Par/ISMIPB.par
r9733 r9734 2 2 3 3 disp(' creating thickness'); 4 md.geometry.surface=-md. x*tan(0.5*pi/180);5 md.geometry.bed=md.geometry.surface-1000+500*sin(md. x*2*pi/max(md.x));4 md.geometry.surface=-md.mesh.x*tan(0.5*pi/180); 5 md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)); 6 6 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 7 7 -
issm/trunk/test/Par/ISMIPC.par
r9733 r9734 2 2 3 3 disp(' creating thickness'); 4 md.geometry.surface=2000-md. x*tan(0.1*pi/180); %to have z>04 md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0 5 5 md.geometry.bed=md.geometry.surface-1000; 6 6 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 7 7 8 8 disp(' creating drag'); 9 %md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md. x*2*pi/max(md.x/2)).*sin(md.y*2*pi/max(md.x/2)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));10 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md. x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x))));9 %md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x/2)).*sin(md.mesh.y*2*pi/max(md.mesh.x/2)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed))); 10 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x)))); 11 11 %Take care of iceshelves: no basal drag 12 12 pos=find(md.mask.elementonfloatingice); -
issm/trunk/test/Par/ISMIPD.par
r9733 r9734 2 2 3 3 disp(' creating thickness'); 4 md.geometry.surface=2000-md. x*tan(0.1*pi/180); %to have z>04 md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0 5 5 md.geometry.bed=md.geometry.surface-1000; 6 6 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 7 7 8 8 disp(' creating drag'); 9 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md. x*2*pi/max(md.x))));9 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x)))); 10 10 %Take care of iceshelves: no basal drag 11 11 pos=find(md.mask.elementonfloatingice); -
issm/trunk/test/Par/ISMIPE.par
r9725 r9734 7 7 md.geometry.bed=zeros(md.mesh.numberofvertices,1); 8 8 for i=1:md.mesh.numberofvertices 9 y=md. y(i);9 y=md.mesh.y(i); 10 10 point1=floor(y/100)+1; 11 11 point2=min(point1+1,51); -
issm/trunk/test/Par/ISMIPF.par
r9725 r9734 3 3 4 4 disp(' creating thickness'); 5 md.geometry.surface=-md. x*tan(3*pi/180);5 md.geometry.surface=-md.mesh.x*tan(3*pi/180); 6 6 %md.geometry.bed=md.geometry.surface-1000; 7 md.geometry.bed=md.geometry.surface-1000+100*exp(-((md. x-max(md.x)/2).^2+(md.y-max(md.y)/2).^2)/(10000^2));7 md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-max(md.mesh.x)/2).^2+(md.mesh.y-max(md.mesh.y)/2).^2)/(10000^2)); 8 8 md.geometry.thickness=md.geometry.surface-md.geometry.bed; 9 9 … … 28 28 md.initialization.pressure=zeros(md.mesh.numberofvertices,1); 29 29 md.initialization.temperature=255*ones(md.mesh.numberofvertices,1); 30 pos=find(md. x==min(md.x) | md.x==max(md.x) | md.y==min(md.y) | md.y==max(md.y));30 pos=find(md.mesh.x==min(md.mesh.x) | md.mesh.x==max(md.mesh.x) | md.mesh.y==min(md.mesh.y) | md.mesh.y==max(md.mesh.y)); 31 31 md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 32 32 md.balancethickness.spcthickness(pos)=md.geometry.thickness(pos); -
issm/trunk/test/Par/Pig.par
r9733 r9734 3 3 %Geometry and observation 4 4 load('../Data/Pig.data','-mat'); 5 md.inversion.vx_obs =InterpFromMeshToMesh2d(index,x,y,vx_obs,md. x,md.y);6 md.inversion.vy_obs =InterpFromMeshToMesh2d(index,x,y,vy_obs,md. x,md.y);7 md.geometry.surface =InterpFromMeshToMesh2d(index,x,y,surface,md. x,md.y);8 md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md. x,md.y);5 md.inversion.vx_obs =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y); 6 md.inversion.vy_obs =InterpFromMeshToMesh2d(index,x,y,vy_obs,md.mesh.x,md.mesh.y); 7 md.geometry.surface =InterpFromMeshToMesh2d(index,x,y,surface,md.mesh.x,md.mesh.y); 8 md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y); 9 9 md.geometry.bed=md.geometry.surface-md.geometry.thickness; 10 10 clear surface thickness vx_obs vy_obs x y index; -
issm/trunk/test/Par/RoundSheetEISMINT.par
r9725 r9734 13 13 tmin=238.15; %K 14 14 st=1.67*10^-2/1000; %k/m; 15 radius=sqrt((md. x).^2+(md.y).^2);15 radius=sqrt((md.mesh.x).^2+(md.mesh.y).^2); 16 16 md.initialization.temperature=(tmin+st*radius); 17 17 md.basalforcings.geothermalflux=4.2*10^-2*ones(md.mesh.numberofvertices,1); … … 29 29 disp(' creating velocities'); 30 30 constant=0.3; 31 md.inversion.vx_obs=constant/2*md. x.*(md.geometry.thickness).^-1;32 md.inversion.vy_obs=constant/2*md. y.*(md.geometry.thickness).^-1;31 md.inversion.vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; 32 md.inversion.vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; 33 33 md.inversion.vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); 34 34 md.initialization.vx=zeros(md.mesh.numberofvertices,1); … … 41 41 md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp'); 42 42 43 radius=sqrt((md. x).*md.x+(md.y).*md.y);43 radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); 44 44 pos=find(radius==min(radius)); 45 md. x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center45 md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center 46 46 47 47 md.diagnostic.spcvx(pos)=0; -
issm/trunk/test/Par/RoundSheetShelf.par
r9733 r9734 7 7 hmin=300; 8 8 hmax=1000; 9 radius=sqrt((md. x).*md.x+(md.y).*md.y);9 radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); 10 10 ymin=min(radius); 11 11 ymax=max(radius); … … 23 23 24 24 25 pos=find(md. x<.2*1e6 & md.x>-.2*1e6 & md.y>0);25 pos=find(md.mesh.x<.2*1e6 & md.mesh.x>-.2*1e6 & md.mesh.y>0); 26 26 md.geometry.thickness(pos)=100; 27 27 md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; 28 28 md.geometry.surface(pos)=md.geometry.bed(pos)+md.geometry.thickness(pos); 29 29 30 pos=find(md. x<.1*1e6 & md.x>-.1*1e6 & md.y<-.5*1e6 & md.y>-.6*1e6);30 pos=find(md.mesh.x<.1*1e6 & md.mesh.x>-.1*1e6 & md.mesh.y<-.5*1e6 & md.mesh.y>-.6*1e6); 31 31 md.geometry.thickness(pos)=100; 32 32 md.geometry.bed(pos)=-di*md.geometry.thickness(pos)-20; … … 79 79 md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); 80 80 81 pos=find(md. x==0 & md.y==0);81 pos=find(md.mesh.x==0 & md.mesh.y==0); 82 82 md.diagnostic.spcvx(pos)=0; 83 83 md.diagnostic.spcvy(pos)=0; -
issm/trunk/test/Par/RoundSheetStaticEISMINT.par
r9733 r9734 2 2 hmin=0.01; 3 3 hmax=2756.7; 4 radius=(sqrt((md. x).^2+(md.y).^2));4 radius=(sqrt((md.mesh.x).^2+(md.mesh.y).^2)); 5 5 radiusmax=max(radius); 6 md.geometry.thickness=hmin*ones(size(md. x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8);6 md.geometry.thickness=hmin*ones(size(md.mesh.x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.mesh.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8); 7 7 md.geometry.bed=0*md.geometry.thickness; 8 8 md.geometry.surface=md.geometry.bed+md.geometry.thickness; … … 34 34 disp(' creating velocities'); 35 35 constant=0.3; 36 md.inversion.vx_obs=constant/2*md. x.*(md.geometry.thickness).^-1;37 md.inversion.vy_obs=constant/2*md. y.*(md.geometry.thickness).^-1;36 md.inversion.vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1; 37 md.inversion.vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1; 38 38 md.inversion.vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2)); 39 39 md.initialization.vx=zeros(md.mesh.numberofvertices,1); … … 46 46 md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp'); 47 47 48 radius=sqrt((md. x).*md.x+(md.y).*md.y);48 radius=sqrt((md.mesh.x).*md.mesh.x+(md.mesh.y).*md.mesh.y); 49 49 pos=find(radius==min(radius)); 50 md. x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center50 md.mesh.x(pos)=0; md.mesh.y(pos)=0; %the closest node to the center is changed to be exactly at the center 51 51 52 52 md.diagnostic.spcvx(pos)=0; -
issm/trunk/test/Par/SquareEISMINT.par
r9733 r9734 2 2 3 3 disp(' creating thickness'); 4 ymin=min(md. y);5 ymax=max(md. y);4 ymin=min(md.mesh.y); 5 ymax=max(md.mesh.y); 6 6 md.geometry.thickness=500*ones(md.mesh.numberofvertices,1); 7 7 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; … … 31 31 32 32 %Evolution of the ice shelf 33 pos=find(md. y==200000); %nodes on the upper boundary condition33 pos=find(md.mesh.y==200000); %nodes on the upper boundary condition 34 34 md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 35 35 md.balancethickness.spcthickness(pos)=500; -
issm/trunk/test/Par/SquareSheetConstrained.par
r9733 r9734 4 4 hmin=300; 5 5 hmax=1000; 6 ymin=min(md. y);7 ymax=max(md. y);8 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);6 ymin=min(md.mesh.y); 7 ymax=max(md.mesh.y); 8 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 9 9 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20; 10 10 md.geometry.surface=md.geometry.bed+md.geometry.thickness; … … 12 12 %Initial velocity 13 13 load('../Data/SquareSheetConstrained.data','-mat'); 14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md. x,md.y);15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md. x,md.y);14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); 15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); 16 16 clear vx vy x y index; 17 17 md.initialization.vz=zeros(md.mesh.numberofvertices,1); -
issm/trunk/test/Par/SquareSheetShelf.par
r9733 r9734 4 4 hmin=300; 5 5 hmax=1000; 6 ymin=min(md. y);7 ymax=max(md. y);8 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);6 ymin=min(md.mesh.y); 7 ymax=max(md.mesh.y); 8 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 9 9 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; 10 10 bed_sheet=-md.materials.rho_ice/md.materials.rho_water*(hmax+(hmin-hmax)*(ymax/2-ymin)/(ymax-ymin)); 11 pos=find(md. y<=ymax/2);11 pos=find(md.mesh.y<=ymax/2); 12 12 md.geometry.bed(pos)=bed_sheet; 13 13 md.geometry.surface=md.geometry.bed+md.geometry.thickness; … … 15 15 %Initial velocity 16 16 load('../Data/SquareSheetShelf.data','-mat'); 17 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md. x,md.y);18 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md. x,md.y);17 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); 18 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); 19 19 clear vx vy x y index; 20 20 md.initialization.vz=zeros(md.mesh.numberofvertices,1); -
issm/trunk/test/Par/SquareShelf.par
r9733 r9734 4 4 hmin=300; 5 5 hmax=1000; 6 ymin=min(md. y);7 ymax=max(md. y);8 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);6 ymin=min(md.mesh.y); 7 ymax=max(md.mesh.y); 8 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 9 9 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; 10 10 md.geometry.surface=md.geometry.bed+md.geometry.thickness; … … 12 12 %Initial velocity and pressure 13 13 load('../Data/SquareShelf.data','-mat'); 14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md. x,md.y);15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md. x,md.y);14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); 15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); 16 16 clear vx vy x y index; 17 17 md.initialization.vz=zeros(md.mesh.numberofvertices,1); -
issm/trunk/test/Par/SquareShelfConstrained.par
r9733 r9734 4 4 hmin=300; 5 5 hmax=1000; 6 ymin=min(md. y);7 ymax=max(md. y);8 md.geometry.thickness=hmax+(hmin-hmax)*(md. y-ymin)/(ymax-ymin);6 ymin=min(md.mesh.y); 7 ymax=max(md.mesh.y); 8 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin); 9 9 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness; 10 10 md.geometry.surface=md.geometry.bed+md.geometry.thickness; … … 12 12 %Initial velocity 13 13 load('../Data/SquareShelfConstrained.data','-mat'); 14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md. x,md.y);15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md. x,md.y);14 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y); 15 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y); 16 16 clear vx vy x y index; 17 17 md.initialization.vz=zeros(md.mesh.numberofvertices,1);
Note:
See TracChangeset
for help on using the changeset viewer.