Index: /issm/trunk-jpl/src/m/classes/modellist.m =================================================================== --- /issm/trunk-jpl/src/m/classes/modellist.m (revision 13011) +++ /issm/trunk-jpl/src/m/classes/modellist.m (revision 13012) @@ -292,2 +292,96 @@ end end + +function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)% {{{ +%BUILDMULTIPLEQUEUEINGSCRIPT - +% +% Usage: +% BuildMultipleQueueingScript(executionpath,codepath) + +disp('building queueing script'); + +%First try and figure out if there is a special script for this particular cluster +function_name=['BuildMultipleQueueingScript' cluster]; + +%some specific treatment of identical cluster, gemini, castor and pollux +if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'), + function_name='BuildMultipleQueueingScriptgemini'; +end + +if exist(function_name,'file'), + %Call this function: + eval([function_name '(name,executionpath,codepath);']); +else + %Call the generic BuildQueueingScript: + BuildMultipleQueueingScriptGeneric(name,executionpath,codepath); +end +end % }}} +function BuildQueueingScriptgemini(name,executionpath,codepath)% {{{ +%BUILDQUEUEINGSCRIPTGEMINI - ... +% +% Usage: +% BuildQueueingScriptgemini(md,executionpath,codepath) + +scriptname=[name '.queue']; + +fid=fopen(scriptname,'w'); +if fid==-1, + error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']); +end + +fprintf(fid,'#!/bin/sh\n'); +fprintf(fid,'cd %s\n',executionpath); +fprintf(fid,'mkdir %s\n',name); +fprintf(fid,'cd %s\n',name); +fprintf(fid,'mv ../ModelList.tar.gz ./\n'); +fprintf(fid,'tar -zxvf ModelList.tar.gz\n'); +fprintf(fid,'foreach i (%s-*vs*.queue)\n',name); +fprintf(fid,'qsub $i\n'); +fprintf(fid,'end\n'); +fclose(fid); +end% }}} +function LaunchMultipleQueueJob(cluster,name,executionpath)% {{{ +%LAUNCHMULTIPLEQUEUEJOB - ... +% +% Usage: +% LaunchMultipleQueueJob(executionpath) + +%First try and figure out if there is a special script for thie particular cluster +function_name=['LaunchMultipleQueueJob' cluster]; + +%some specific treatment of identical cluster, gemini, castor and pollux +if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'), + function_name='LaunchMultipleQueueJobgemini'; +end + +if exist(function_name,'file'), + %Call this function: + eval([function_name '(cluster,name,executionpath);']); +else + %Call the generic LaunchMultipleQueueJob: + LaunchMultipleQueueJobGeneric(cluster,name,executionpath); +end +end% }}} +function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)% {{{ +%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster +% +% Usage: +% LaunchMultipleQueueJobgemini(cluster,name,executionpath) + + +%first, check we have the binary file and the queueing script +if ~exist([ name '.queue'],'file'), + error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward'); +end + +if ~exist('ModelList.tar.gz','file'), + error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward'); +end + +%upload both files to cluster +disp('uploading input file, queueing script and variables script'); +eval(['!scp ModelList.tar.gz ' name '.queue ' cluster ':' executionpath]); + +disp('launching solution sequence on remote cluster'); +issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']); +end% }}} Index: /issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m (revision 13012) +++ /issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m (revision 13012) @@ -0,0 +1,84 @@ +function md=BamgCall(md,field,hmin,hmax,gradation,epsilon), +%BAMGCALL - call bam +% +% build a metric using the Hessian of the given field +% call Bamg and the output mesh is plugged onto the model +% -hmin = minimum edge length (m) +% -hmax = maximum edge length (m) +% -gradation = maximum edge length gradation between 2 elements +% -epsilon = average error on each element (m/yr) +% +% Usage: +% md=BamgCall(md,field,hmin,hmax,gradation,epsilon); +% +% Example: +% md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9); + +%2d geometric parameter (do not change) +scale=2/9; + +%Compute Hessian +t1=clock; fprintf('%s',' computing Hessian...'); +hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%Compute metric +t1=clock; fprintf('%s',' computing metric...'); +if length(md.nodeonwater)==md.mesh.numberofvertices, + pos=find(md.nodeonwater); +else + pos=[]; +end +metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%write files +t1=clock; fprintf('%s',' writing initial mesh files...'); +fid=fopen('carre0.met','w'); +fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3); +fprintf(fid,'%i %i %i\n',metric'); +fclose(fid); + +fid=fopen('carre0.mesh','w'); + +%initialiation +fprintf(fid,'%s %i\n','MeshVersionFormatted',0); + +%dimension +fprintf(fid,'\n%s\n%i\n','Dimension',2); + +%Vertices +fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); +fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); + +%Triangles +fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); +fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]'); +numberofelements1=md.mesh.numberofelements; + +%close +fclose(fid); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%call bamg +fprintf('%s\n',' call Bamg...'); +system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']); + +%plug new mesh +t1=clock; fprintf('\n%s',' reading final mesh files...'); +A=meshread('carre1.mesh'); +md.mesh.x=A.x; +md.mesh.y=A.y; +md.z=zeros(A.nods,1); +md.mesh.elements=A.index; +md.mesh.numberofvertices=A.nods; +md.mesh.numberofelements=A.nels; +numberofelements2=md.mesh.numberofelements; +t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%display number of elements +fprintf('\n%s %i',' inital number of elements:',numberofelements1); +fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); + +%clean up: +system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh'); Index: /issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m (revision 13012) +++ /issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m (revision 13012) @@ -0,0 +1,65 @@ +function md=BamgCallFromMetric(md,metric,gradation), +%BAMGCALL - call bam +% +% call Bamg and the output mesh is plugged onto the model +% -gradation = maximum edge length gradation between 2 elements +% +% Usage: +% md=BamgCallFromMetric(md,metric,gradation); +% +% Example: +% md=BamgCall(md,metric,1500,10^8,1.3,0.9); + +%2d geometric parameter (do not change) +scale=2/9; + +%write files +t1=clock; fprintf('%s',' writing initial mesh files...'); +fid=fopen('carre0.met','w'); +fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3); +fprintf(fid,'%i %i %i\n',metric'); +fclose(fid); + +fid=fopen('carre0.mesh','w'); + +%initialiation +fprintf(fid,'%s %i\n','MeshVersionFormatted',0); + +%dimension +fprintf(fid,'\n%s\n%i\n','Dimension',2); + +%Vertices +fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); +fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); + +%Triangles +fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); +fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]'); +numberofelements1=md.mesh.numberofelements; + +%close +fclose(fid); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%call bamg +fprintf('%s\n',' call Bamg...'); +system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']); + +%plug new mesh +t1=clock; fprintf('\n%s',' reading final mesh files...'); +A=meshread('carre1.mesh'); +md.mesh.x=A.x; +md.mesh.y=A.y; +md.z=zeros(A.nods,1); +md.mesh.elements=A.index; +md.mesh.numberofvertices=A.nods; +md.mesh.numberofelements=A.nels; +numberofelements2=md.mesh.numberofelements; +t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%display number of elements +fprintf('\n%s %i',' inital number of elements:',numberofelements1); +fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); + +%clean up: +system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh'); Index: /issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m (revision 13012) +++ /issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m (revision 13012) @@ -0,0 +1,104 @@ +function md=YamsCall(md,field,hmin,hmax,gradation,epsilon), +%YAMSCALL - call yams +% +% build a metric using the Hessian of the given field +% call Yams and the output mesh is plugged onto the model +% -hmin = minimum edge length (m) +% -hmax = maximum edge length (m) +% -gradation = maximum edge length gradation between 2 elements +% -epsilon = average error on each element (m/yr) +% +% Usage: +% md=YamsCall(md,field,hmin,hmax,gradation,epsilon); +% +% Example: +% md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9); + +%2d geometric parameter (do not change) +scale=2/9; + +%Compute Hessian +t1=clock; fprintf('%s',' computing Hessian...'); +hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%Compute metric +t1=clock; fprintf('%s',' computing metric...'); +if length(md.mask.vertexonwater)==md.mesh.numberofvertices, + pos=find(md.mask.vertexonwater); +else + pos=[]; +end +metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%write files +t1=clock; fprintf('%s',' writing initial mesh files...'); +save -ascii carre0.met metric + +fid=fopen('carre0.mesh','w'); + +%initialiation +fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1); + +%dimension +fprintf(fid,'\n%s\n%i\n','Dimension',2); + +%Vertices +fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); +fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]'); + +%Triangles +fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); +fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]'); +numberofelements1=md.mesh.numberofelements; + +%Deal with rifts +if ~isnan(md.rifts.riftstruct), + + %we have the list of triangles that make up the rift. keep those triangles around during refinement. + triangles=[]; + for i=1:size(md.rifts.riftstruct,1), + triangles=[triangles md.rifts(i).segments(:,3)']; + end + + fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles)); + fprintf(fid,'%i\n',triangles); +end + +%close +fclose(fid); +t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%call yams +fprintf('%s\n',' call Yams...'); +if ispc + %windows + system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); +elseif ismac + %Macosx + system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); +else + %Linux + system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); +end + +%plug new mesh +t1=clock; fprintf('\n%s',' reading final mesh files...'); +Tria=load('carre1.tria'); +Coor=load('carre1.coor'); +md.mesh.x=Coor(:,1); +md.mesh.y=Coor(:,2); +md.mesh.z=zeros(size(Coor,1),1); +md.mesh.elements=Tria; +md.mesh.numberofvertices=size(Coor,1); +md.mesh.numberofelements=size(Tria,1); +numberofelements2=md.mesh.numberofelements; +t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); + +%display number of elements +fprintf('\n%s %i',' inital number of elements:',numberofelements1); +fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); + +%clean up: +system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb'); Index: /issm/trunk-jpl/src/m/contrib/bamg/meshread.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/bamg/meshread.m (revision 13012) +++ /issm/trunk-jpl/src/m/contrib/bamg/meshread.m (revision 13012) @@ -0,0 +1,41 @@ +function Struct=meshread(filename); + +%some checks +if ~exist(filename), + error(['meshread error message: file ' filename ' not found!']); +end + +fid=fopen(filename,'r'); + +while (~feof(fid)), + + A=fscanf(fid,'%s',1); + + if strcmp(A,'MeshVersionFormatted'); + Struct.Version=fscanf(fid,'%s',1); + + elseif strcmp(A,'Dimension'), + Struct.Dimension=fscanf(fid,'%i',1); + + elseif strcmp(A,'Vertices'), + Struct.nods=fscanf(fid,'%i',1); + A=fscanf(fid,'%f %f %f',[3 Struct.nods]); + Struct.x=A(1,:)'; + Struct.y=A(2,:)'; + + elseif strcmp(A,'Triangles'), + Struct.nels=fscanf(fid,'%i',1); + A=fscanf(fid,'%i %i %i',[4 Struct.nels]); + Struct.index=A(1:3,:)'; + + elseif strcmp(A,'Quadrilaterals'), + Struct.nels=fscanf(fid,'%i',1); + A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]); + Struct.index=A(1:4,:)'; + else + %do nothing + + end +end + +fclose(fid); Index: /issm/trunk-jpl/src/m/contrib/bamg/yams.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/bamg/yams.m (revision 13012) +++ /issm/trunk-jpl/src/m/contrib/bamg/yams.m (revision 13012) @@ -0,0 +1,171 @@ +function md=yams(md,varargin); +%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator +% +% Usage: +% md=yams(md,varargin); +% where varargin is a lit of paired arguments. +% arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed +% arguments can be: 'velocities': matlab file containing the velocities [m/yr] +% optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice +% this option is used to minimize the metric on water (no refinement) +% optional arguments: 'resolution': initial mesh resolution [m] +% optional arguments: 'nsteps': number of steps of mesh adaptation +% optional arguments: 'epsilon': average interpolation error wished [m/yr] +% optional arguments: 'hmin': minimum edge length +% optional arguments: 'hmanx': maximum edge +% optional arguments: 'riftoutline': if rifts are present, specifies rift outline file. +% +% +% Examples: +% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat'); +% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp'); +% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000); + +%recover options +options=pairoptions(varargin{:}); +options=deleteduplicates(options,1); + +%recover some fields +disp('MeshYams Options:') +domainoutline=getfieldvalue(options,'domainoutline'); +disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline)); +riftoutline=getfieldvalue(options,'riftoutline','N/A'); +disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline)); +groundeddomain=getfieldvalue(options,'groundeddomain','N/A'); +disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain)); +velocities=getfieldvalue(options,'velocities'); +disp(sprintf(' %-15s: ''%s''','Velocities',velocities)); +resolution=getfieldvalue(options,'resolution',5000); +disp(sprintf(' %-15s: %f','Resolution',resolution)); +nsteps=getfieldvalue(options,'nsteps',6); +disp(sprintf(' %-15s: %i','nsteps',nsteps)); +gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1)); +disp(sprintf(' %-15s: %g','gradation',gradation(1))); +epsilon=getfieldvalue(options,'epsilon',3); +disp(sprintf(' %-15s: %f','epsilon',epsilon)); +hmin=getfieldvalue(options,'hmin',500); +disp(sprintf(' %-15s: %f','hmin',hmin)); +hmax=getfieldvalue(options,'hmax',150*10^3); +disp(sprintf(' %-15s: %f\n','hmax',hmax)); + +%mesh with initial resolution +disp('Initial mesh generation...'); +if strcmpi(riftoutline,'N/A'); + md=setmesh(md,domainoutline,resolution); +else + md=setmesh(md,domainoutline,riftoutline,resolution); + md=meshprocessrifts(md,domainoutline); +end +disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]); + +%load velocities +disp('loading velocities...'); +Names=VelFindVarNames(velocities); +Vel=load(velocities); + +%start mesh adaptation +for i=1:nsteps, + disp(['Iteration #' num2str(i) '/' num2str(nsteps)]); + + %interpolate velocities onto mesh + disp(' interpolating velocities...'); + if strcmpi(Names.interp,'node'), + vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); + vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); + else + vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); + vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); + end + field=sqrt(vx_obs.^2+vy_obs.^2); + + %set mask.vertexonwater field + if ~strcmp(groundeddomain,'N/A'), + nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); + md.mask.vertexonwater=ones(md.mesh.numberofvertices,1); + md.mask.vertexonwater(find(nodeground))=0; + else + md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); + end + + %adapt according to velocities + disp(' adapting...'); + md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon); + + %if we have rifts, we just messed them up, we need to recreate the segments that constitute those + %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement. + if md.rifts.numrifts, + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); + md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); + md.mesh.segments=findsegments(md); + md=yamsrecreateriftsegments(md); + end + +end + +disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]); + +%Now, build the connectivity tables for this mesh. +md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); +md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); + +%recreate segments +md.mesh.segments=findsegments(md); +md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; + +%Fill in rest of fields: +md.mesh.z=zeros(md.mesh.numberofvertices,1); +md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); +md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1); +md.mesh.elementonbed=ones(md.mesh.numberofelements,1); +md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); +if ~strcmp(groundeddomain,'N/A'), + nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); + md.mask.vertexonwater=ones(md.mesh.numberofvertices,1); + md.mask.vertexonwater(find(nodeground))=0; +else + md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); +end +if strcmpi(Names.interp,'node'), + md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); + md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); +else + md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); + md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); +end +md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2); + +%deal with rifts +if md.rifts.numrifts, + %first, recreate rift segments + md=meshyamsrecreateriftsegments(md); + + %using the segments, recreate the penaltypairs + for j=1:md.rifts.numrifts, + rift=md.rifts.riftstruct(j); + + %build normals and lengths of segments: + 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 ); + 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))))); + 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))))); + + %ok, build penaltypairs: + numpenaltypairs=length(rift.segments)/2-1; + rift.penaltypairs=zeros(numpenaltypairs,7); + + for i=1:numpenaltypairs, + rift.penaltypairs(i,1)=rift.segments(i,2); + rift.penaltypairs(i,2)=rift.segments(end-i,2); + rift.penaltypairs(i,3)=rift.segments(i,3); + rift.penaltypairs(i,4)=rift.segments(end-i,3); + rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1); + rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1); + rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2; + end + %renormalize norms: + norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2); + rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms; + rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms; + + md.rifts.riftstruct(j)=rift; + end +end Index: /issm/trunk-jpl/src/m/geometry/GetAreas.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/GetAreas.m (revision 13012) +++ /issm/trunk-jpl/src/m/geometry/GetAreas.m (revision 13012) @@ -0,0 +1,51 @@ +function areas=GetAreas(index,x,y,varargin) +%GETAREAS - compute areas or volumes of elements +% +% compute areas of triangular elements or volumes +% of pentahedrons +% +% Usage: +% areas =GetAreas(index,x,y); +% volumes=GetAreas(index,x,y,z); +% +% Examples: +% areas =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); +% volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z); + +%get number of elements and number of nodes +nels=size(index,1); +nods=length(x); +if nargin==4, z=varargin{1}; end + +%some checks +if nargout~=1 | (nargin~=3 & nargin~=4), + help GetAreas + error('GetAreas error message: bad usage') +end +if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)), + error('GetAreas error message: x,y and z do not have the same length') +end +if max(index(:))>nods, + error(['GetAreas error message: index should not have values above ' num2str(nods) ]) +end +if (nargin==3 & size(index,2)~=3), + error('GetAreas error message: index should have 3 columns for 2d meshes.') +end +if (nargin==4 & size(index,2)~=6), + error('GetAreas error message: index should have 6 columns for 3d meshes.') +end + +%initialization +areas=zeros(nels,1); +x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); +y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3)); + +%compute the volume of each element +if nargin==3, + %compute the surface of the triangle + areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))); +else + %V=area(triangle)*1/3(z1+z2+z3) + thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2); + areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness; +end Index: sm/trunk-jpl/src/m/mesh/BamgCall.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/BamgCall.m (revision 13011) +++ (revision ) @@ -1,84 +1,0 @@ -function md=BamgCall(md,field,hmin,hmax,gradation,epsilon), -%BAMGCALL - call bam -% -% build a metric using the Hessian of the given field -% call Bamg and the output mesh is plugged onto the model -% -hmin = minimum edge length (m) -% -hmax = maximum edge length (m) -% -gradation = maximum edge length gradation between 2 elements -% -epsilon = average error on each element (m/yr) -% -% Usage: -% md=BamgCall(md,field,hmin,hmax,gradation,epsilon); -% -% Example: -% md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9); - -%2d geometric parameter (do not change) -scale=2/9; - -%Compute Hessian -t1=clock; fprintf('%s',' computing Hessian...'); -hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%Compute metric -t1=clock; fprintf('%s',' computing metric...'); -if length(md.nodeonwater)==md.mesh.numberofvertices, - pos=find(md.nodeonwater); -else - pos=[]; -end -metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%write files -t1=clock; fprintf('%s',' writing initial mesh files...'); -fid=fopen('carre0.met','w'); -fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3); -fprintf(fid,'%i %i %i\n',metric'); -fclose(fid); - -fid=fopen('carre0.mesh','w'); - -%initialiation -fprintf(fid,'%s %i\n','MeshVersionFormatted',0); - -%dimension -fprintf(fid,'\n%s\n%i\n','Dimension',2); - -%Vertices -fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); -fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); - -%Triangles -fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); -fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]'); -numberofelements1=md.mesh.numberofelements; - -%close -fclose(fid); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%call bamg -fprintf('%s\n',' call Bamg...'); -system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']); - -%plug new mesh -t1=clock; fprintf('\n%s',' reading final mesh files...'); -A=meshread('carre1.mesh'); -md.mesh.x=A.x; -md.mesh.y=A.y; -md.z=zeros(A.nods,1); -md.mesh.elements=A.index; -md.mesh.numberofvertices=A.nods; -md.mesh.numberofelements=A.nels; -numberofelements2=md.mesh.numberofelements; -t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%display number of elements -fprintf('\n%s %i',' inital number of elements:',numberofelements1); -fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); - -%clean up: -system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh'); Index: sm/trunk-jpl/src/m/mesh/BamgCallFromMetric.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/BamgCallFromMetric.m (revision 13011) +++ (revision ) @@ -1,65 +1,0 @@ -function md=BamgCallFromMetric(md,metric,gradation), -%BAMGCALL - call bam -% -% call Bamg and the output mesh is plugged onto the model -% -gradation = maximum edge length gradation between 2 elements -% -% Usage: -% md=BamgCallFromMetric(md,metric,gradation); -% -% Example: -% md=BamgCall(md,metric,1500,10^8,1.3,0.9); - -%2d geometric parameter (do not change) -scale=2/9; - -%write files -t1=clock; fprintf('%s',' writing initial mesh files...'); -fid=fopen('carre0.met','w'); -fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3); -fprintf(fid,'%i %i %i\n',metric'); -fclose(fid); - -fid=fopen('carre0.mesh','w'); - -%initialiation -fprintf(fid,'%s %i\n','MeshVersionFormatted',0); - -%dimension -fprintf(fid,'\n%s\n%i\n','Dimension',2); - -%Vertices -fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); -fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]'); - -%Triangles -fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); -fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]'); -numberofelements1=md.mesh.numberofelements; - -%close -fclose(fid); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%call bamg -fprintf('%s\n',' call Bamg...'); -system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']); - -%plug new mesh -t1=clock; fprintf('\n%s',' reading final mesh files...'); -A=meshread('carre1.mesh'); -md.mesh.x=A.x; -md.mesh.y=A.y; -md.z=zeros(A.nods,1); -md.mesh.elements=A.index; -md.mesh.numberofvertices=A.nods; -md.mesh.numberofelements=A.nels; -numberofelements2=md.mesh.numberofelements; -t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%display number of elements -fprintf('\n%s %i',' inital number of elements:',numberofelements1); -fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); - -%clean up: -system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh'); Index: /issm/trunk-jpl/src/m/mesh/FixMesh.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/FixMesh.m (revision 13011) +++ /issm/trunk-jpl/src/m/mesh/FixMesh.m (revision 13012) @@ -1,10 +1,10 @@ function [index2 x2 y2 value2]=FixMesh(index,x,y,value) -%FixMesh fix mesh with broken triangles, orphan vertices, etc ... +% FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ... % -% Usage: -% [index2 x2 y2 value2]=FixMesh(index,x,y,value) -% where index,x,y is a delaunay triangulation, -% value is a field on the input triangulation, with values at the vertices -% index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices +% Usage: +% [index2 x2 y2 value2]=FixMesh(index,x,y,value) +% where index,x,y is a delaunay triangulation, +% value is a field on the input triangulation, with values at the vertices +% index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices % % Index: sm/trunk-jpl/src/m/mesh/GetAreas.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/GetAreas.m (revision 13011) +++ (revision ) @@ -1,51 +1,0 @@ -function areas=GetAreas(index,x,y,varargin) -%GETAREAS - compute areas or volumes of elements -% -% compute areas of triangular elements or volumes -% of pentahedrons -% -% Usage: -% areas =GetAreas(index,x,y); -% volumes=GetAreas(index,x,y,z); -% -% Examples: -% areas =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); -% volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z); - -%get number of elements and number of nodes -nels=size(index,1); -nods=length(x); -if nargin==4, z=varargin{1}; end - -%some checks -if nargout~=1 | (nargin~=3 & nargin~=4), - help GetAreas - error('GetAreas error message: bad usage') -end -if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)), - error('GetAreas error message: x,y and z do not have the same length') -end -if max(index(:))>nods, - error(['GetAreas error message: index should not have values above ' num2str(nods) ]) -end -if (nargin==3 & size(index,2)~=3), - error('GetAreas error message: index should have 3 columns for 2d meshes.') -end -if (nargin==4 & size(index,2)~=6), - error('GetAreas error message: index should have 6 columns for 3d meshes.') -end - -%initialization -areas=zeros(nels,1); -x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); -y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3)); - -%compute the volume of each element -if nargin==3, - %compute the surface of the triangle - areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))); -else - %V=area(triangle)*1/3(z1+z2+z3) - thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2); - areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness; -end Index: sm/trunk-jpl/src/m/mesh/GetCharacteristicLength.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/GetCharacteristicLength.m (revision 13011) +++ (revision ) @@ -1,45 +1,0 @@ -function length=GetCharacteristicLength(index,x,y,varargin) -%GETCHARACTERISTICLENGTH - compute characteristic length for a mesh -% -% compute characteristic lengths of every element of a mesh. -% -% Usage: -% length =GetCharacteristicLength(index,x,y); -% length =GetCharacteristicLength(index,x,y,z); -% -% Examples: -% length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y); -% length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y,md.z); - - -%get number of elements and number of nodes -nels=size(index,1); -nods=numel(x); - -%some checks -if nargout~=1 | (nargin~=3 & nargin~=4), - help GetCharacteristicLength - error('GetCharacteristicLength error message: bad usage') -end -if ((numel(y)~=nods) | (nargin==4 & numel(z)~=nods)), - error('GetCharacteristicLength error message: x,y and z do not have the same length') -end -if max(index(:))>nods, - error(['GetCharacteristicLength error message: index should not have values above ' num2str(nods) ]) -end -if (nargin==3 & size(index,2)~=3), - error('GetCharacteristicLength error message: index should have 3 columns for 2d meshes.') -end -if (nargin==4 & size(index,2)~=6), - error('GetCharacteristicLength error message: index should have 6 columns for 3d meshes.') -end - -%get areas or volumes. -areas=GetAreas(index,x,y,varargin{:}); - -%for a 2d mesh: -if nargin==3, - length=sqrt(2*areas); -else - error('not supported yet'); -end Index: /issm/trunk-jpl/src/m/mesh/NodeInElement.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/NodeInElement.m (revision 13011) +++ /issm/trunk-jpl/src/m/mesh/NodeInElement.m (revision 13012) @@ -1,6 +1,7 @@ function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity); -%NODEINELEMENT: find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to. +% NODEINELEMENT - find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to. % -% Usage: node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity); +% Usage: +% node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity); % % See also Nodeconnectivity Index: sm/trunk-jpl/src/m/mesh/YamsCall.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/YamsCall.m (revision 13011) +++ (revision ) @@ -1,104 +1,0 @@ -function md=YamsCall(md,field,hmin,hmax,gradation,epsilon), -%YAMSCALL - call yams -% -% build a metric using the Hessian of the given field -% call Yams and the output mesh is plugged onto the model -% -hmin = minimum edge length (m) -% -hmax = maximum edge length (m) -% -gradation = maximum edge length gradation between 2 elements -% -epsilon = average error on each element (m/yr) -% -% Usage: -% md=YamsCall(md,field,hmin,hmax,gradation,epsilon); -% -% Example: -% md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9); - -%2d geometric parameter (do not change) -scale=2/9; - -%Compute Hessian -t1=clock; fprintf('%s',' computing Hessian...'); -hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node'); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%Compute metric -t1=clock; fprintf('%s',' computing metric...'); -if length(md.mask.vertexonwater)==md.mesh.numberofvertices, - pos=find(md.mask.vertexonwater); -else - pos=[]; -end -metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%write files -t1=clock; fprintf('%s',' writing initial mesh files...'); -save -ascii carre0.met metric - -fid=fopen('carre0.mesh','w'); - -%initialiation -fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1); - -%dimension -fprintf(fid,'\n%s\n%i\n','Dimension',2); - -%Vertices -fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices); -fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]'); - -%Triangles -fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements); -fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]'); -numberofelements1=md.mesh.numberofelements; - -%Deal with rifts -if ~isnan(md.rifts.riftstruct), - - %we have the list of triangles that make up the rift. keep those triangles around during refinement. - triangles=[]; - for i=1:size(md.rifts.riftstruct,1), - triangles=[triangles md.rifts(i).segments(:,3)']; - end - - fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles)); - fprintf(fid,'%i\n',triangles); -end - -%close -fclose(fid); -t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%call yams -fprintf('%s\n',' call Yams...'); -if ispc - %windows - system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); -elseif ismac - %Macosx - system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); -else - %Linux - system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']); -end - -%plug new mesh -t1=clock; fprintf('\n%s',' reading final mesh files...'); -Tria=load('carre1.tria'); -Coor=load('carre1.coor'); -md.mesh.x=Coor(:,1); -md.mesh.y=Coor(:,2); -md.mesh.z=zeros(size(Coor,1),1); -md.mesh.elements=Tria; -md.mesh.numberofvertices=size(Coor,1); -md.mesh.numberofelements=size(Tria,1); -numberofelements2=md.mesh.numberofelements; -t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']); - -%display number of elements -fprintf('\n%s %i',' inital number of elements:',numberofelements1); -fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2); - -%clean up: -system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb'); Index: sm/trunk-jpl/src/m/mesh/isconnected.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/isconnected.m (revision 13011) +++ (revision ) @@ -1,13 +1,0 @@ -function flag=isconnected(elements,A,B) -%ISCONNECTED: are two nodes connected by a triangulation? -% -% Usage: flag=isconnected(elements,A,B) -% -% - -elements=ElementsFromEdge(elements,A,B); -if isempty(elements), - flag=0; -else - flag=1; -end Index: sm/trunk-jpl/src/m/mesh/meshread.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshread.m (revision 13011) +++ (revision ) @@ -1,41 +1,0 @@ -function Struct=meshread(filename); - -%some checks -if ~exist(filename), - error(['meshread error message: file ' filename ' not found!']); -end - -fid=fopen(filename,'r'); - -while (~feof(fid)), - - A=fscanf(fid,'%s',1); - - if strcmp(A,'MeshVersionFormatted'); - Struct.Version=fscanf(fid,'%s',1); - - elseif strcmp(A,'Dimension'), - Struct.Dimension=fscanf(fid,'%i',1); - - elseif strcmp(A,'Vertices'), - Struct.nods=fscanf(fid,'%i',1); - A=fscanf(fid,'%f %f %f',[3 Struct.nods]); - Struct.x=A(1,:)'; - Struct.y=A(2,:)'; - - elseif strcmp(A,'Triangles'), - Struct.nels=fscanf(fid,'%i',1); - A=fscanf(fid,'%i %i %i',[4 Struct.nels]); - Struct.index=A(1:3,:)'; - - elseif strcmp(A,'Quadrilaterals'), - Struct.nels=fscanf(fid,'%i',1); - A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]); - Struct.index=A(1:4,:)'; - else - %do nothing - - end -end - -fclose(fid); Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13011) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13012) @@ -93,2 +93,18 @@ md.mesh.elementonbed=ones(md.mesh.numberofelements,1); md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); +end + +function flag=isconnected(elements,A,B)% {{{ + %ISCONNECTED: are two nodes connected by a triangulation? + % + % Usage: flag=isconnected(elements,A,B) + % + % + + elements=ElementsFromEdge(elements,A,B); + if isempty(elements), + flag=0; + else + flag=1; + end +end % }}} Index: /issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m (revision 13012) +++ /issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m (revision 13012) @@ -0,0 +1,26 @@ +function tips=rifttipsonmesh(md,riftoutline) +%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of +% rifts. + +%read rifts from outline file +rifts=expread(riftoutline,1); + +tips=[]; + +for i=1:length(rifts), + rift=rifts(i); + + x_tip=rift.x(1); + y_tip=rift.y(1); + + index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); + tips(end+1)=index; + + x_tip=rift.x(end); + y_tip=rift.y(end); + + index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); + tips(end+1)=index; + +end + Index: sm/trunk-jpl/src/m/mesh/rifttipsonmesh.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifttipsonmesh.m (revision 13011) +++ (revision ) @@ -1,26 +1,0 @@ -function tips=rifttipsonmesh(md,riftoutline) -%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of -% rifts. - -%read rifts from outline file -rifts=expread(riftoutline,1); - -tips=[]; - -for i=1:length(rifts), - rift=rifts(i); - - x_tip=rift.x(1); - y_tip=rift.y(1); - - index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); - tips(end+1)=index; - - x_tip=rift.x(end); - y_tip=rift.y(end); - - index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip); - tips(end+1)=index; - -end - Index: sm/trunk-jpl/src/m/mesh/yams.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/yams.m (revision 13011) +++ (revision ) @@ -1,171 +1,0 @@ -function md=yams(md,varargin); -%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator -% -% Usage: -% md=yams(md,varargin); -% where varargin is a lit of paired arguments. -% arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed -% arguments can be: 'velocities': matlab file containing the velocities [m/yr] -% optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice -% this option is used to minimize the metric on water (no refinement) -% optional arguments: 'resolution': initial mesh resolution [m] -% optional arguments: 'nsteps': number of steps of mesh adaptation -% optional arguments: 'epsilon': average interpolation error wished [m/yr] -% optional arguments: 'hmin': minimum edge length -% optional arguments: 'hmanx': maximum edge -% optional arguments: 'riftoutline': if rifts are present, specifies rift outline file. -% -% -% Examples: -% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat'); -% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp'); -% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000); - -%recover options -options=pairoptions(varargin{:}); -options=deleteduplicates(options,1); - -%recover some fields -disp('MeshYams Options:') -domainoutline=getfieldvalue(options,'domainoutline'); -disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline)); -riftoutline=getfieldvalue(options,'riftoutline','N/A'); -disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline)); -groundeddomain=getfieldvalue(options,'groundeddomain','N/A'); -disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain)); -velocities=getfieldvalue(options,'velocities'); -disp(sprintf(' %-15s: ''%s''','Velocities',velocities)); -resolution=getfieldvalue(options,'resolution',5000); -disp(sprintf(' %-15s: %f','Resolution',resolution)); -nsteps=getfieldvalue(options,'nsteps',6); -disp(sprintf(' %-15s: %i','nsteps',nsteps)); -gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1)); -disp(sprintf(' %-15s: %g','gradation',gradation(1))); -epsilon=getfieldvalue(options,'epsilon',3); -disp(sprintf(' %-15s: %f','epsilon',epsilon)); -hmin=getfieldvalue(options,'hmin',500); -disp(sprintf(' %-15s: %f','hmin',hmin)); -hmax=getfieldvalue(options,'hmax',150*10^3); -disp(sprintf(' %-15s: %f\n','hmax',hmax)); - -%mesh with initial resolution -disp('Initial mesh generation...'); -if strcmpi(riftoutline,'N/A'); - md=setmesh(md,domainoutline,resolution); -else - md=setmesh(md,domainoutline,riftoutline,resolution); - md=meshprocessrifts(md,domainoutline); -end -disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]); - -%load velocities -disp('loading velocities...'); -Names=VelFindVarNames(velocities); -Vel=load(velocities); - -%start mesh adaptation -for i=1:nsteps, - disp(['Iteration #' num2str(i) '/' num2str(nsteps)]); - - %interpolate velocities onto mesh - disp(' interpolating velocities...'); - if strcmpi(Names.interp,'node'), - vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); - vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); - else - vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); - vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); - end - field=sqrt(vx_obs.^2+vy_obs.^2); - - %set mask.vertexonwater field - if ~strcmp(groundeddomain,'N/A'), - nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); - md.mask.vertexonwater=ones(md.mesh.numberofvertices,1); - md.mask.vertexonwater(find(nodeground))=0; - else - md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); - end - - %adapt according to velocities - disp(' adapting...'); - md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon); - - %if we have rifts, we just messed them up, we need to recreate the segments that constitute those - %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement. - if md.rifts.numrifts, - md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); - md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); - md.mesh.segments=findsegments(md); - md=yamsrecreateriftsegments(md); - end - -end - -disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]); - -%Now, build the connectivity tables for this mesh. -md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); -md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); - -%recreate segments -md.mesh.segments=findsegments(md); -md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; - -%Fill in rest of fields: -md.mesh.z=zeros(md.mesh.numberofvertices,1); -md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1); -md.mesh.elementonbed=ones(md.mesh.numberofelements,1); -md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); -if ~strcmp(groundeddomain,'N/A'), - nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2); - md.mask.vertexonwater=ones(md.mesh.numberofvertices,1); - md.mask.vertexonwater(find(nodeground))=0; -else - md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1); -end -if strcmpi(Names.interp,'node'), - md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); - md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); -else - md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0); - md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0); -end -md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2); - -%deal with rifts -if md.rifts.numrifts, - %first, recreate rift segments - md=meshyamsrecreateriftsegments(md); - - %using the segments, recreate the penaltypairs - for j=1:md.rifts.numrifts, - rift=md.rifts.riftstruct(j); - - %build normals and lengths of segments: - 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 ); - 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))))); - 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))))); - - %ok, build penaltypairs: - numpenaltypairs=length(rift.segments)/2-1; - rift.penaltypairs=zeros(numpenaltypairs,7); - - for i=1:numpenaltypairs, - rift.penaltypairs(i,1)=rift.segments(i,2); - rift.penaltypairs(i,2)=rift.segments(end-i,2); - rift.penaltypairs(i,3)=rift.segments(i,3); - rift.penaltypairs(i,4)=rift.segments(end-i,3); - rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1); - rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1); - rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2; - end - %renormalize norms: - norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2); - rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms; - rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms; - - md.rifts.riftstruct(j)=rift; - end -end Index: sm/trunk-jpl/src/m/miscellaneous/create_region.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/create_region.m (revision 13011) +++ (revision ) @@ -1,15 +1,0 @@ -function create_region(name) -%CREATE_REGION - create region ???? -% -% very temporary function. -% -% Usage: -% create_region(name) - -eval(['mkdir ' name]); -eval(['cd ' name ]); -!mkdir Delivery Exp_Par Results -cd Exp_Par -!cp ../../RonneShelf/Exp_Par/* ./ -!rm -rf Hole* -eval(['!mv Ronne.par ' name '.par']); Index: sm/trunk-jpl/src/m/miscellaneous/findarg.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/findarg.m (revision 13011) +++ (revision ) @@ -1,51 +1,0 @@ -function vals=findarg(arglist,field) -%FINDARG - find argument associated to a field in a list -% -% This function parses through an argument list (typically varargin in a routine) -% looking for a character array equal to field. Once this is found, we return the -% next value in the varargin (if possible). -% Because field might appear several times in the argument list, we return a structure -% holding all these values. -% Note that all comparisons to field value are case independent. -% -% Usage: -% vals=findarg(arglist,field) -% -% Example: -% routine foobar calls vals=findarg('Data',varargin) -% with varargin='Data',1,'Data','foo','Plot','velocity','Arrow',4 -% findarg would return the following structure: vals(1).value=1, vals(2).value='foo'; - -%some argument checking: -if ((nargin==0) | (nargout==0)), - help findarg; - error('findarg error message'); -end - -if ~ischar(field), - error('findarg error message: field should be a string'); -end - -if ~iscell(arglist), - error('findarg error message: argument list should be a cell array.'); -end - -%Recover data to plot -founddata=0; - -for i=1:(length(arglist)-1), %data in arglist comes in pairs, hence the -1. - if ischar(arglist{i}), - if (strcmpi(arglist{i},field)), - founddata=founddata+1; - if founddata==1, - vals.value=arglist{i+1}; - else - vals(end+1).value=arglist{i+1}; - end - end - end -end - -if founddata==0, - vals=[]; -end Index: sm/trunk-jpl/src/m/miscellaneous/netcdf.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/netcdf.m (revision 13011) +++ (revision ) @@ -1,136 +1,0 @@ -function S = netcdf(File,varargin) -% Function to read NetCDF files -% S = netcdf(File) -% Input Arguments -% File = NetCDF file to read -% Optional Input Arguments: -% 'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)] -% 'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs] -% Output Arguments: -% S = Structure of NetCDF data organised as per NetCDF definition -% Notes: -% Only version 1, classic 32bit, NetCDF files are supported. By default -% data are extracted into the S.VarArray().Data field for all variables. -% To read the header only call S = netcdf(File,'Var',[]); -% -% SEE ALSO -% --------------------------------------------------------------------------- -S = []; - -try - if exist(File,'file') fp = fopen(File,'r','b'); - else fp = []; error('File not found'); end - if fp == -1 error('Unable to open file'); end - -% Read header - Magic = fread(fp,4,'uint8=>char'); - if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end - if uint8(Magic(4))~=1 error('Version not supported'); end - S.NumRecs = fread(fp,1,'uint32=>uint32'); - S.DimArray = DimArray(fp); - S.AttArray = AttArray(fp); - S.VarArray = VarArray(fp); - -% Setup indexing to arrays and records - Var = ones(1,length(S.VarArray)); - Rec = ones(1,S.NumRecs); - for i = 1:2:length(varargin) - if strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1; - elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1; - else error('Optional input argument not recognised'); end - end - if sum(Var)==0 fclose(fp); return; end - -% Read non-record variables - Dim = double(cat(2,S.DimArray.Dim)); - ID = double(cat(2,S.VarArray.Type)); - - for i = 1:length(S.VarArray) - D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0); - if isempty(RecID{i}) - if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end - if Var(i) - S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D); - fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); - else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end - else S.VarArray(i).Data = []; end - end - -% Read record variables - for k = 1:S.NumRecs - for i = 1:length(S.VarArray) - if ~isempty(RecID{i}) - D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D); - if length(D)==1 D=[D,1]; end - if Var(i) & Rec(k) - S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,... - ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D)); - if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end - else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end - end - end - end - - fclose(fp); -catch - Err = lasterror; fprintf('%s\n',Err.message); - if ~isempty(fp) && fp ~= -1 fclose(fp); end -end - -% --------------------------------------------------------------------------------------- -% Utility functions - -function S = Size(ID) -% Size of NetCDF data type, ID, in bytes - S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}})); - -function T = Type(ID) -% Matlab string for CDF data type, ID - T = subsref({'int8','char','int16','int32','single','double'},... - struct('type','{}','subs',{{ID}})); - -function N = Pad(Num,ID) -% Number of elements to read after padding to 4 bytes for type ID - N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0); - -function S = String(fp) -% Read a CDF string; Size,[String,[Padding]] - S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').'; - -function A = ReOrder(A,S) -% Rearrange CDF array A to size S with matlab ordering - A = permute(reshape(A,fliplr(S)),fliplr(1:length(S))); - -function S = DimArray(fp) -% Read DimArray into structure - if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION - for i = 1:fread(fp,1,'uint32=>uint32') - S(i).Str = String(fp); - S(i).Dim = fread(fp,1,'uint32=>uint32'); - end - else fread(fp,1,'uint32=>uint32'); S = []; end - -function S = AttArray(fp) -% Read AttArray into structure - if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE - for i = 1:fread(fp,1,'uint32=>uint32') - S(i).Str = String(fp); - ID = fread(fp,1,'uint32=>uint32'); - Num = fread(fp,1,'uint32=>uint32'); - S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).'; - end - else fread(fp,1,'uint32=>uint32'); S = []; end - -function S = VarArray(fp) -% Read VarArray into structure - if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE - for i = 1:fread(fp,1,'uint32=>uint32') - S(i).Str = String(fp); - Num = double(fread(fp,1,'uint32=>uint32')); - S(i).DimID = double(fread(fp,Num,'uint32=>uint32')); - S(i).AttArray = AttArray(fp); - S(i).Type = fread(fp,1,'uint32=>uint32'); - S(i).VSize = fread(fp,1,'uint32=>uint32'); - S(i).Begin = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only - end - else fread(fp,1,'uint32=>uint32'); S = []; end Index: /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m (revision 13011) +++ /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m (revision 13012) @@ -6,5 +6,5 @@ %Read netcdf file -data=netcdf(File); +data=readnetcdf(File); %initialize output @@ -26,2 +26,141 @@ S.(fieldname)=fieldvalue; end +end + +function S = readnetcdf(File,varargin) +% Function to read NetCDF files +% S = netcdf(File) +% Input Arguments +% File = NetCDF file to read +% Optional Input Arguments: +% 'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)] +% 'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs] +% Output Arguments: +% S = Structure of NetCDF data organised as per NetCDF definition +% Notes: +% Only version 1, classic 32bit, NetCDF files are supported. By default +% data are extracted into the S.VarArray().Data field for all variables. +% To read the header only call S = netcdf(File,'Var',[]); +% +% SEE ALSO +% --------------------------------------------------------------------------- +S = []; + +try + if exist(File,'file') fp = fopen(File,'r','b'); + else fp = []; error('File not found'); end + if fp == -1 error('Unable to open file'); end + +% Read header + Magic = fread(fp,4,'uint8=>char'); + if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end + if uint8(Magic(4))~=1 error('Version not supported'); end + S.NumRecs = fread(fp,1,'uint32=>uint32'); + S.DimArray = DimArray(fp); + S.AttArray = AttArray(fp); + S.VarArray = VarArray(fp); + +% Setup indexing to arrays and records + Var = ones(1,length(S.VarArray)); + Rec = ones(1,S.NumRecs); + for i = 1:2:length(varargin) + if strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1; + elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1; + else error('Optional input argument not recognised'); end + end + if sum(Var)==0 fclose(fp); return; end + +% Read non-record variables + Dim = double(cat(2,S.DimArray.Dim)); + ID = double(cat(2,S.VarArray.Type)); + + for i = 1:length(S.VarArray) + D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0); + if isempty(RecID{i}) + if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end + if Var(i) + S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D); + fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); + else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end + else S.VarArray(i).Data = []; end + end + +% Read record variables + for k = 1:S.NumRecs + for i = 1:length(S.VarArray) + if ~isempty(RecID{i}) + D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D); + if length(D)==1 D=[D,1]; end + if Var(i) & Rec(k) + S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,... + ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D)); + if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end + else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end + end + end + end + + fclose(fp); +catch + Err = lasterror; fprintf('%s\n',Err.message); + if ~isempty(fp) && fp ~= -1 fclose(fp); end +end + +% --------------------------------------------------------------------------------------- +% Utility functions + +function S = Size(ID) +% Size of NetCDF data type, ID, in bytes + S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}})); + +function T = Type(ID) +% Matlab string for CDF data type, ID + T = subsref({'int8','char','int16','int32','single','double'},... + struct('type','{}','subs',{{ID}})); + +function N = Pad(Num,ID) +% Number of elements to read after padding to 4 bytes for type ID + N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0); + +function S = String(fp) +% Read a CDF string; Size,[String,[Padding]] + S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').'; + +function A = ReOrder(A,S) +% Rearrange CDF array A to size S with matlab ordering + A = permute(reshape(A,fliplr(S)),fliplr(1:length(S))); + +function S = DimArray(fp) +% Read DimArray into structure + if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION + for i = 1:fread(fp,1,'uint32=>uint32') + S(i).Str = String(fp); + S(i).Dim = fread(fp,1,'uint32=>uint32'); + end + else fread(fp,1,'uint32=>uint32'); S = []; end + +function S = AttArray(fp) +% Read AttArray into structure + if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE + for i = 1:fread(fp,1,'uint32=>uint32') + S(i).Str = String(fp); + ID = fread(fp,1,'uint32=>uint32'); + Num = fread(fp,1,'uint32=>uint32'); + S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).'; + end + else fread(fp,1,'uint32=>uint32'); S = []; end + +function S = VarArray(fp) +% Read VarArray into structure + if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE + for i = 1:fread(fp,1,'uint32=>uint32') + S(i).Str = String(fp); + Num = double(fread(fp,1,'uint32=>uint32')); + S(i).DimID = double(fread(fp,Num,'uint32=>uint32')); + S(i).AttArray = AttArray(fp); + S(i).Type = fread(fp,1,'uint32=>uint32'); + S(i).VSize = fread(fp,1,'uint32=>uint32'); + S(i).Begin = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only + end + else fread(fp,1,'uint32=>uint32'); S = []; end +end