source: issm/oecreview/Archive/12678-13393/ISSM-13011-13012.diff

Last change on this file was 13394, checked in by Mathieu Morlighem, 13 years ago

Added 12678-13393
File size: 71.7 KB

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJob.m

     
    1 function LaunchMultipleQueueJob(cluster,name,executionpath)
    2 %LAUNCHMULTIPLEQUEUEJOB - ...
    3 %
    4 %   Usage:
    5 %      LaunchMultipleQueueJob(executionpath)
    6 
    7 %First try and figure out if there is a special script for thie particular cluster
    8 function_name=['LaunchMultipleQueueJob' cluster];
    9 
    10 %some specific treatment of identical cluster, gemini, castor and pollux
    11 if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
    12         function_name='LaunchMultipleQueueJobgemini';
    13 end
    14 
    15 if exist(function_name,'file'),
    16         %Call this function:
    17         eval([function_name '(cluster,name,executionpath);']);
    18 else
    19         %Call the generic LaunchMultipleQueueJob:
    20         LaunchMultipleQueueJobGeneric(cluster,name,executionpath);
    21 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptgemini.m

     
    1 function BuildQueueingScriptgemini(name,executionpath,codepath)
    2 %BUILDQUEUEINGSCRIPTGEMINI - ...
    3 %
    4 %   Usage:
    5 %      BuildQueueingScriptgemini(md,executionpath,codepath)
    6 
    7 scriptname=[name '.queue'];
    8 
    9 fid=fopen(scriptname,'w');
    10 if fid==-1,
    11         error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']);
    12 end
    13 
    14 fprintf(fid,'#!/bin/sh\n');
    15 fprintf(fid,'cd %s\n',executionpath);
    16 fprintf(fid,'mkdir %s\n',name);
    17 fprintf(fid,'cd %s\n',name);
    18 fprintf(fid,'mv ../ModelList.tar.gz ./\n');
    19 fprintf(fid,'tar -zxvf ModelList.tar.gz\n');
    20 fprintf(fid,'foreach i (%s-*vs*.queue)\n',name);
    21 fprintf(fid,'qsub $i\n');
    22 fprintf(fid,'end\n');
    23 fclose(fid);

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobgemini.m

     
    1 function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)
    2 %LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster
    3 %
    4 %   Usage:
    5 %      LaunchMultipleQueueJobgemini(cluster,name,executionpath)
    6 
    7 
    8 %first, check we have the binary file and the queueing script
    9 if ~exist([ name '.queue'],'file'),
    10         error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward');
    11 end
    12 
    13 if ~exist('ModelList.tar.gz','file'),
    14         error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward');
    15 end
    16 
    17 %upload both files to cluster
    18 disp('uploading input file,  queueing script and variables script');
    19 eval(['!scp ModelList.tar.gz ' name '.queue '  cluster ':' executionpath]);
    20 
    21 disp('launching solution sequence on remote cluster');
    22 issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']);

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptGeneric.m

     
    1 function BuildMultipleQueueingScriptGeneric(name,executionpath,codepath)
    2 %BUILDMULTIPLEQUEUEINGSCRIPTGENERIC - ...
    3 %
    4 %   Usage:
    5 %      BuildMultipleQueueingScriptGeneric(executionpath,codepath)
    6 
    7 %not done yet
    8 error('BuildMultipleQueueingScriptGenericerror message: not supported yet!');
    9 

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScript.m

     
    1 function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)
    2 %BUILDMULTIPLEQUEUEINGSCRIPT -
    3 %
    4 %   Usage:
    5 %      BuildMultipleQueueingScript(executionpath,codepath)
    6 
    7 disp('building queueing script');
    8 
    9 %First try and figure out if there is a special script for this particular cluster
    10 function_name=['BuildMultipleQueueingScript' cluster];
    11 
    12 %some specific treatment of identical cluster, gemini, castor and pollux
    13 if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
    14         function_name='BuildMultipleQueueingScriptgemini';
    15 end
    16 
    17 if exist(function_name,'file'),
    18         %Call this function:
    19         eval([function_name '(name,executionpath,codepath);']);
    20 else
    21         %Call the generic BuildQueueingScript:
    22         BuildMultipleQueueingScriptGeneric(name,executionpath,codepath);
    23 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobGeneric.m

     
    1 function LaunchMultipleQueueJobGeneric(cluster,name,executionpath)
    2 %LAUNCHMULTIPLEQUEUEJOBGENERIC - Generic routine to launch multiple queueing job
    3 %
    4 %   Usage:
    5 %      LaunchMultipleQueueJobGeneric(cluster,name,executionpath)
    6 
    7 error('LaunchMultipleQueueJobGeneric error message: not supported yet!');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/geometry/GetAreas.m

     
     1function areas=GetAreas(index,x,y,varargin)
     2%GETAREAS - compute areas or volumes of elements
     3%
     4%   compute areas of triangular elements or volumes
     5%   of pentahedrons
     6%
     7%   Usage:
     8%      areas  =GetAreas(index,x,y);
     9%      volumes=GetAreas(index,x,y,z);
     10%
     11%   Examples:
     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
     15%get number of elements and number of nodes
     16nels=size(index,1);
     17nods=length(x);
     18if nargin==4, z=varargin{1}; end
     19
     20%some checks
     21if nargout~=1 | (nargin~=3 & nargin~=4),
     22        help GetAreas
     23        error('GetAreas error message: bad usage')
     24end
     25if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
     26        error('GetAreas error message: x,y and z do not have the same length')
     27end
     28if max(index(:))>nods,
     29        error(['GetAreas error message: index should not have values above ' num2str(nods) ])
     30end
     31if (nargin==3 & size(index,2)~=3),
     32        error('GetAreas error message: index should have 3 columns for 2d meshes.')
     33end
     34if (nargin==4 & size(index,2)~=6),
     35        error('GetAreas error message: index should have 6 columns for 3d meshes.')
     36end
     37
     38%initialization
     39areas=zeros(nels,1);
     40x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
     41y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
     42
     43%compute the volume of each element
     44if nargin==3,
     45        %compute the surface of the triangle
     46        areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
     47else
     48        %V=area(triangle)*1/3(z1+z2+z3)
     49        thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
     50        areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
     51end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/meshread.m

     
    1 function Struct=meshread(filename);
    2 
    3 %some checks
    4 if ~exist(filename),
    5         error(['meshread error message: file ' filename ' not found!']);
    6 end
    7 
    8 fid=fopen(filename,'r');
    9 
    10 while (~feof(fid)),
    11 
    12         A=fscanf(fid,'%s',1);
    13 
    14         if strcmp(A,'MeshVersionFormatted');
    15                 Struct.Version=fscanf(fid,'%s',1);
    16 
    17         elseif strcmp(A,'Dimension'),
    18                 Struct.Dimension=fscanf(fid,'%i',1);
    19 
    20         elseif strcmp(A,'Vertices'),
    21                 Struct.nods=fscanf(fid,'%i',1);
    22                 A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
    23                 Struct.x=A(1,:)';
    24                 Struct.y=A(2,:)';
    25 
    26         elseif strcmp(A,'Triangles'),
    27                 Struct.nels=fscanf(fid,'%i',1);
    28                 A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
    29                 Struct.index=A(1:3,:)';
    30 
    31         elseif strcmp(A,'Quadrilaterals'),
    32                 Struct.nels=fscanf(fid,'%i',1);
    33                 A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
    34                 Struct.index=A(1:4,:)';
    35         else
    36                 %do nothing
    37 
    38         end
    39 end
    40 
    41 fclose(fid);

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCallFromMetric.m

     
    1 function md=BamgCallFromMetric(md,metric,gradation),
    2 %BAMGCALL - call bam
    3 %
    4 %   call Bamg and the output mesh is plugged onto the model
    5 %   -gradation = maximum edge length gradation between 2 elements
    6 %
    7 %   Usage:
    8 %      md=BamgCallFromMetric(md,metric,gradation);
    9 %
    10 %   Example:
    11 %      md=BamgCall(md,metric,1500,10^8,1.3,0.9);
    12 
    13 %2d geometric parameter (do not change)
    14 scale=2/9;
    15 
    16 %write files
    17 t1=clock; fprintf('%s','      writing initial mesh files...');
    18 fid=fopen('carre0.met','w');
    19 fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
    20 fprintf(fid,'%i %i %i\n',metric');
    21 fclose(fid);
    22 
    23 fid=fopen('carre0.mesh','w');
    24 
    25 %initialiation
    26 fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
    27 
    28 %dimension
    29 fprintf(fid,'\n%s\n%i\n','Dimension',2);
    30 
    31 %Vertices
    32 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
    33 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
    34 
    35 %Triangles
    36 fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
    37 fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
    38 numberofelements1=md.mesh.numberofelements;
    39 
    40 %close
    41 fclose(fid);
    42 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    43 
    44 %call bamg
    45 fprintf('%s\n','      call Bamg...');
    46 system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
    47 
    48 %plug new mesh
    49 t1=clock; fprintf('\n%s','      reading final mesh files...');
    50 A=meshread('carre1.mesh');
    51 md.mesh.x=A.x;
    52 md.mesh.y=A.y;
    53 md.z=zeros(A.nods,1);
    54 md.mesh.elements=A.index;
    55 md.mesh.numberofvertices=A.nods;
    56 md.mesh.numberofelements=A.nels;
    57 numberofelements2=md.mesh.numberofelements;
    58 t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    59 
    60 %display number of elements
    61 fprintf('\n%s %i','      inital number of elements:',numberofelements1);
    62 fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
    63 
    64 %clean up:
    65 system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetCharacteristicLength.m

     
    1 function length=GetCharacteristicLength(index,x,y,varargin)
    2 %GETCHARACTERISTICLENGTH - compute characteristic length for a mesh
    3 %
    4 %   compute characteristic lengths of every element of a mesh.
    5 %
    6 %   Usage:
    7 %      length  =GetCharacteristicLength(index,x,y);
    8 %      length  =GetCharacteristicLength(index,x,y,z);
    9 %
    10 %   Examples:
    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 
    14 
    15 %get number of elements and number of nodes
    16 nels=size(index,1);
    17 nods=numel(x);
    18 
    19 %some checks
    20 if nargout~=1 | (nargin~=3 & nargin~=4),
    21         help GetCharacteristicLength
    22         error('GetCharacteristicLength error message: bad usage')
    23 end
    24 if ((numel(y)~=nods) | (nargin==4 & numel(z)~=nods)),
    25         error('GetCharacteristicLength error message: x,y and z do not have the same length')
    26 end
    27 if max(index(:))>nods,
    28         error(['GetCharacteristicLength error message: index should not have values above ' num2str(nods) ])
    29 end
    30 if (nargin==3 & size(index,2)~=3),
    31         error('GetCharacteristicLength error message: index should have 3 columns for 2d meshes.')
    32 end
    33 if (nargin==4 & size(index,2)~=6),
    34         error('GetCharacteristicLength error message: index should have 6 columns for 3d meshes.')
    35 end
    36 
    37 %get areas or volumes.
    38 areas=GetAreas(index,x,y,varargin{:});
    39 
    40 %for a 2d mesh:
    41 if nargin==3,
    42         length=sqrt(2*areas);
    43 else
    44         error('not supported yet');
    45 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetAreas.m

     
    1 function areas=GetAreas(index,x,y,varargin)
    2 %GETAREAS - compute areas or volumes of elements
    3 %
    4 %   compute areas of triangular elements or volumes
    5 %   of pentahedrons
    6 %
    7 %   Usage:
    8 %      areas  =GetAreas(index,x,y);
    9 %      volumes=GetAreas(index,x,y,z);
    10 %
    11 %   Examples:
    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 
    15 %get number of elements and number of nodes
    16 nels=size(index,1);
    17 nods=length(x);
    18 if nargin==4, z=varargin{1}; end
    19 
    20 %some checks
    21 if nargout~=1 | (nargin~=3 & nargin~=4),
    22         help GetAreas
    23         error('GetAreas error message: bad usage')
    24 end
    25 if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
    26         error('GetAreas error message: x,y and z do not have the same length')
    27 end
    28 if max(index(:))>nods,
    29         error(['GetAreas error message: index should not have values above ' num2str(nods) ])
    30 end
    31 if (nargin==3 & size(index,2)~=3),
    32         error('GetAreas error message: index should have 3 columns for 2d meshes.')
    33 end
    34 if (nargin==4 & size(index,2)~=6),
    35         error('GetAreas error message: index should have 6 columns for 3d meshes.')
    36 end
    37 
    38 %initialization
    39 areas=zeros(nels,1);
    40 x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
    41 y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
    42 
    43 %compute the volume of each element
    44 if nargin==3,
    45         %compute the surface of the triangle
    46         areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
    47 else
    48         %V=area(triangle)*1/3(z1+z2+z3)
    49         thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
    50         areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
    51 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/isconnected.m

     
    1 function flag=isconnected(elements,A,B)
    2 %ISCONNECTED: are two nodes connected by a triangulation?
    3 %
    4 %   Usage: flag=isconnected(elements,A,B)
    5 %
    6 %
    7 
    8 elements=ElementsFromEdge(elements,A,B);
    9 if isempty(elements),
    10         flag=0;
    11 else
    12         flag=1;
    13 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCall.m

     
    1 function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
    2 %BAMGCALL - call bam
    3 %
    4 %   build a metric using the Hessian of the given field
    5 %   call Bamg and the output mesh is plugged onto the model
    6 %   -hmin = minimum edge length (m)
    7 %   -hmax = maximum edge length (m)
    8 %   -gradation = maximum edge length gradation between 2 elements
    9 %   -epsilon = average error on each element (m/yr)
    10 %
    11 %   Usage:
    12 %      md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
    13 %
    14 %   Example:
    15 %      md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
    16 
    17 %2d geometric parameter (do not change)
    18 scale=2/9;
    19 
    20 %Compute Hessian
    21 t1=clock; fprintf('%s','      computing Hessian...');
    22 hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
    23 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    24 
    25 %Compute metric
    26 t1=clock; fprintf('%s','      computing metric...');
    27 if length(md.nodeonwater)==md.mesh.numberofvertices,
    28         pos=find(md.nodeonwater);
    29 else
    30         pos=[];
    31 end
    32 metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
    33 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    34 
    35 %write files
    36 t1=clock; fprintf('%s','      writing initial mesh files...');
    37 fid=fopen('carre0.met','w');
    38 fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
    39 fprintf(fid,'%i %i %i\n',metric');
    40 fclose(fid);
    41 
    42 fid=fopen('carre0.mesh','w');
    43 
    44 %initialiation
    45 fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
    46 
    47 %dimension
    48 fprintf(fid,'\n%s\n%i\n','Dimension',2);
    49 
    50 %Vertices
    51 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
    52 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
    53 
    54 %Triangles
    55 fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
    56 fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
    57 numberofelements1=md.mesh.numberofelements;
    58 
    59 %close
    60 fclose(fid);
    61 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    62 
    63 %call bamg
    64 fprintf('%s\n','      call Bamg...');
    65 system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
    66 
    67 %plug new mesh
    68 t1=clock; fprintf('\n%s','      reading final mesh files...');
    69 A=meshread('carre1.mesh');
    70 md.mesh.x=A.x;
    71 md.mesh.y=A.y;
    72 md.z=zeros(A.nods,1);
    73 md.mesh.elements=A.index;
    74 md.mesh.numberofvertices=A.nods;
    75 md.mesh.numberofelements=A.nels;
    76 numberofelements2=md.mesh.numberofelements;
    77 t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    78 
    79 %display number of elements
    80 fprintf('\n%s %i','      inital number of elements:',numberofelements1);
    81 fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
    82 
    83 %clean up:
    84 system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/YamsCall.m

     
    1 function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
    2 %YAMSCALL - call yams
    3 %
    4 %   build a metric using the Hessian of the given field
    5 %   call Yams and the output mesh is plugged onto the model
    6 %   -hmin = minimum edge length (m)
    7 %   -hmax = maximum edge length (m)
    8 %   -gradation = maximum edge length gradation between 2 elements
    9 %   -epsilon = average error on each element (m/yr)
    10 %
    11 %   Usage:
    12 %      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
    13 %
    14 %   Example:
    15 %      md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
    16 
    17 %2d geometric parameter (do not change)
    18 scale=2/9;
    19 
    20 %Compute Hessian
    21 t1=clock; fprintf('%s','      computing Hessian...');
    22 hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
    23 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    24 
    25 %Compute metric
    26 t1=clock; fprintf('%s','      computing metric...');
    27 if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
    28         pos=find(md.mask.vertexonwater);
    29 else
    30         pos=[];
    31 end
    32 metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
    33 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    34 
    35 %write files
    36 t1=clock; fprintf('%s','      writing initial mesh files...');
    37 save -ascii carre0.met  metric
    38 
    39 fid=fopen('carre0.mesh','w');
    40 
    41 %initialiation
    42 fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
    43 
    44 %dimension
    45 fprintf(fid,'\n%s\n%i\n','Dimension',2);
    46 
    47 %Vertices
    48 fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
    49 fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
    50 
    51 %Triangles
    52 fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
    53 fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
    54 numberofelements1=md.mesh.numberofelements;
    55        
    56 %Deal with rifts
    57 if ~isnan(md.rifts.riftstruct),
    58        
    59         %we have the list of triangles that make up the rift. keep those triangles around during refinement.
    60         triangles=[];
    61         for i=1:size(md.rifts.riftstruct,1),
    62                 triangles=[triangles md.rifts(i).segments(:,3)'];
    63         end
    64 
    65         fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
    66         fprintf(fid,'%i\n',triangles);
    67 end
    68 
    69 %close
    70 fclose(fid);
    71 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    72 
    73 %call yams
    74 fprintf('%s\n','      call Yams...');
    75 if ispc
    76         %windows
    77         system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
    78 elseif ismac
    79         %Macosx
    80         system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
    81 else
    82         %Linux
    83         system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
    84 end
    85 
    86 %plug new mesh
    87 t1=clock; fprintf('\n%s','      reading final mesh files...');
    88 Tria=load('carre1.tria');
    89 Coor=load('carre1.coor');
    90 md.mesh.x=Coor(:,1);
    91 md.mesh.y=Coor(:,2);
    92 md.mesh.z=zeros(size(Coor,1),1);
    93 md.mesh.elements=Tria;
    94 md.mesh.numberofvertices=size(Coor,1);
    95 md.mesh.numberofelements=size(Tria,1);
    96 numberofelements2=md.mesh.numberofelements;
    97 t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
    98 
    99 %display number of elements
    100 fprintf('\n%s %i','      inital number of elements:',numberofelements1);
    101 fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
    102 
    103 %clean up:
    104 system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/yams.m

     
    1 function md=yams(md,varargin);
    2 %MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
    3 %
    4 %   Usage:
    5 %      md=yams(md,varargin);
    6 %      where varargin is a lit of paired arguments.
    7 %      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
    8 %      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
    9 %      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
    10 %                          this option is used to minimize the metric on water (no refinement)
    11 %      optional arguments: 'resolution': initial mesh resolution [m]
    12 %      optional arguments: 'nsteps': number of steps of mesh adaptation
    13 %      optional arguments: 'epsilon': average interpolation error wished [m/yr]
    14 %      optional arguments: 'hmin': minimum edge length
    15 %      optional arguments: 'hmanx': maximum edge
    16 %      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
    17 %     
    18 %
    19 %   Examples:
    20 %      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
    21 %      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
    22 %      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
    23 
    24 %recover options
    25 options=pairoptions(varargin{:});
    26 options=deleteduplicates(options,1);
    27 
    28 %recover some fields
    29 disp('MeshYams Options:')
    30 domainoutline=getfieldvalue(options,'domainoutline');
    31 disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
    32 riftoutline=getfieldvalue(options,'riftoutline','N/A');
    33 disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
    34 groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
    35 disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
    36 velocities=getfieldvalue(options,'velocities');
    37 disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
    38 resolution=getfieldvalue(options,'resolution',5000);
    39 disp(sprintf('   %-15s: %f','Resolution',resolution));
    40 nsteps=getfieldvalue(options,'nsteps',6);
    41 disp(sprintf('   %-15s: %i','nsteps',nsteps));
    42 gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
    43 disp(sprintf('   %-15s: %g','gradation',gradation(1)));
    44 epsilon=getfieldvalue(options,'epsilon',3);
    45 disp(sprintf('   %-15s: %f','epsilon',epsilon));
    46 hmin=getfieldvalue(options,'hmin',500);
    47 disp(sprintf('   %-15s: %f','hmin',hmin));
    48 hmax=getfieldvalue(options,'hmax',150*10^3);
    49 disp(sprintf('   %-15s: %f\n','hmax',hmax));
    50 
    51 %mesh with initial resolution
    52 disp('Initial mesh generation...');
    53 if strcmpi(riftoutline,'N/A');
    54         md=setmesh(md,domainoutline,resolution);
    55 else
    56         md=setmesh(md,domainoutline,riftoutline,resolution);
    57         md=meshprocessrifts(md,domainoutline);
    58 end
    59 disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
    60 
    61 %load velocities
    62 disp('loading velocities...');
    63 Names=VelFindVarNames(velocities);
    64 Vel=load(velocities);
    65 
    66 %start mesh adaptation
    67 for i=1:nsteps,
    68         disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
    69 
    70         %interpolate velocities onto mesh
    71         disp('   interpolating velocities...');
    72         if strcmpi(Names.interp,'node'),
    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         else
    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         end
    79         field=sqrt(vx_obs.^2+vy_obs.^2);
    80 
    81         %set mask.vertexonwater  field
    82         if ~strcmp(groundeddomain,'N/A'),
    83                 nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
    84                 md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
    85                 md.mask.vertexonwater(find(nodeground))=0;
    86         else
    87                 md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
    88         end
    89 
    90         %adapt according to velocities
    91         disp('   adapting...');
    92         md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
    93 
    94         %if we have rifts, we just messed them up, we need to recreate the segments that constitute those
    95         %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
    96         if md.rifts.numrifts,
    97                 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
    98                 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
    99                 md.mesh.segments=findsegments(md);
    100                 md=yamsrecreateriftsegments(md);
    101         end
    102 
    103 end
    104        
    105 disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
    106 
    107 %Now, build the connectivity tables for this mesh.
    108 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
    109 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
    110 
    111 %recreate segments
    112 md.mesh.segments=findsegments(md);
    113 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
    114 
    115 %Fill in rest of fields:
    116 md.mesh.z=zeros(md.mesh.numberofvertices,1);
    117 md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
    118 md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
    119 md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
    120 md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
    121 if ~strcmp(groundeddomain,'N/A'),
    122         nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
    123         md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
    124         md.mask.vertexonwater(find(nodeground))=0;
    125 else
    126         md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
    127 end
    128 if strcmpi(Names.interp,'node'),
    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 else
    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 end
    135 md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
    136 
    137 %deal with rifts
    138 if md.rifts.numrifts,
    139         %first, recreate rift segments
    140         md=meshyamsrecreateriftsegments(md);
    141 
    142         %using the segments, recreate the penaltypairs
    143         for j=1:md.rifts.numrifts,
    144                 rift=md.rifts.riftstruct(j);
    145 
    146                 %build normals and lengths of segments:
    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 
    151                 %ok, build penaltypairs:
    152                 numpenaltypairs=length(rift.segments)/2-1;
    153                 rift.penaltypairs=zeros(numpenaltypairs,7);
    154 
    155                 for i=1:numpenaltypairs,
    156                         rift.penaltypairs(i,1)=rift.segments(i,2);
    157                         rift.penaltypairs(i,2)=rift.segments(end-i,2);
    158                         rift.penaltypairs(i,3)=rift.segments(i,3);
    159                         rift.penaltypairs(i,4)=rift.segments(end-i,3);
    160                         rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
    161                         rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
    162                         rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
    163                 end
    164                 %renormalize norms:
    165                 norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
    166                 rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
    167                 rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
    168 
    169                 md.rifts.riftstruct(j)=rift;
    170         end
    171 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifttipsonmesh.m

     
    1 function tips=rifttipsonmesh(md,riftoutline)
    2 %RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of
    3 %                rifts.
    4 
    5 %read rifts from outline file
    6 rifts=expread(riftoutline,1);
    7 
    8 tips=[];
    9 
    10 for i=1:length(rifts),
    11         rift=rifts(i);
    12 
    13         x_tip=rift.x(1);
    14         y_tip=rift.y(1);
    15        
    16         index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
    17         tips(end+1)=index;
    18 
    19         x_tip=rift.x(end);
    20         y_tip=rift.y(end);
    21        
    22         index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
    23         tips(end+1)=index;
    24 
    25 end
    26 

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/NodeInElement.m

     
    11function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity);
    2 %NODEINELEMENT: find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
     2% NODEINELEMENT - find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
    33%
    4 %  Usage: node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
     4%    Usage:
     5%      node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
    56%
    67%  See also Nodeconnectivity
    78%

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m

     
     1function tips=rifttipsonmesh(md,riftoutline)
     2%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of
     3%                rifts.
     4
     5%read rifts from outline file
     6rifts=expread(riftoutline,1);
     7
     8tips=[];
     9
     10for i=1:length(rifts),
     11        rift=rifts(i);
     12
     13        x_tip=rift.x(1);
     14        y_tip=rift.y(1);
     15       
     16        index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
     17        tips(end+1)=index;
     18
     19        x_tip=rift.x(end);
     20        y_tip=rift.y(end);
     21       
     22        index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
     23        tips(end+1)=index;
     24
     25end
     26

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m

     
    9292md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
    9393md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
    9494md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
     95end
     96
     97function flag=isconnected(elements,A,B)% {{{
     98        %ISCONNECTED: are two nodes connected by a triangulation?
     99        %
     100        %   Usage: flag=isconnected(elements,A,B)
     101        %
     102        %
     103
     104        elements=ElementsFromEdge(elements,A,B);
     105        if isempty(elements),
     106                flag=0;
     107        else
     108                flag=1;
     109        end
     110end % }}}

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/FixMesh.m

     
    11function  [index2 x2 y2 value2]=FixMesh(index,x,y,value)
    2 %FixMesh fix mesh with broken triangles, orphan vertices, etc ...
     2% FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ...
    33%
    4 % Usage:
    5 %            [index2 x2 y2 value2]=FixMesh(index,x,y,value)
    6 %            where index,x,y is a delaunay triangulation,
    7 %                  value is a field on the input triangulation, with values at the vertices
    8 %                  index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
     4%   Usage:
     5%      [index2 x2 y2 value2]=FixMesh(index,x,y,value)
     6%      where index,x,y is a delaunay triangulation,
     7%      value is a field on the input triangulation, with values at the vertices
     8%      index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
    99%
    1010%
    1111

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/modellist.m

     
    291291                end % }}}
    292292        end
    293293end
     294
     295function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)% {{{
     296%BUILDMULTIPLEQUEUEINGSCRIPT -
     297%
     298%   Usage:
     299%      BuildMultipleQueueingScript(executionpath,codepath)
     300
     301disp('building queueing script');
     302
     303%First try and figure out if there is a special script for this particular cluster
     304function_name=['BuildMultipleQueueingScript' cluster];
     305
     306%some specific treatment of identical cluster, gemini, castor and pollux
     307if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
     308        function_name='BuildMultipleQueueingScriptgemini';
     309end
     310
     311if exist(function_name,'file'),
     312        %Call this function:
     313        eval([function_name '(name,executionpath,codepath);']);
     314else
     315        %Call the generic BuildQueueingScript:
     316        BuildMultipleQueueingScriptGeneric(name,executionpath,codepath);
     317end
     318end % }}}
     319function BuildQueueingScriptgemini(name,executionpath,codepath)% {{{
     320%BUILDQUEUEINGSCRIPTGEMINI - ...
     321%
     322%   Usage:
     323%      BuildQueueingScriptgemini(md,executionpath,codepath)
     324
     325scriptname=[name '.queue'];
     326
     327fid=fopen(scriptname,'w');
     328if fid==-1,
     329        error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']);
     330end
     331
     332fprintf(fid,'#!/bin/sh\n');
     333fprintf(fid,'cd %s\n',executionpath);
     334fprintf(fid,'mkdir %s\n',name);
     335fprintf(fid,'cd %s\n',name);
     336fprintf(fid,'mv ../ModelList.tar.gz ./\n');
     337fprintf(fid,'tar -zxvf ModelList.tar.gz\n');
     338fprintf(fid,'foreach i (%s-*vs*.queue)\n',name);
     339fprintf(fid,'qsub $i\n');
     340fprintf(fid,'end\n');
     341fclose(fid);
     342end% }}}
     343function LaunchMultipleQueueJob(cluster,name,executionpath)% {{{
     344%LAUNCHMULTIPLEQUEUEJOB - ...
     345%
     346%   Usage:
     347%      LaunchMultipleQueueJob(executionpath)
     348
     349%First try and figure out if there is a special script for thie particular cluster
     350function_name=['LaunchMultipleQueueJob' cluster];
     351
     352%some specific treatment of identical cluster, gemini, castor and pollux
     353if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
     354        function_name='LaunchMultipleQueueJobgemini';
     355end
     356
     357if exist(function_name,'file'),
     358        %Call this function:
     359        eval([function_name '(cluster,name,executionpath);']);
     360else
     361        %Call the generic LaunchMultipleQueueJob:
     362        LaunchMultipleQueueJobGeneric(cluster,name,executionpath);
     363end
     364end% }}}
     365function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)% {{{
     366%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster
     367%
     368%   Usage:
     369%      LaunchMultipleQueueJobgemini(cluster,name,executionpath)
     370
     371
     372%first, check we have the binary file and the queueing script
     373if ~exist([ name '.queue'],'file'),
     374        error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward');
     375end
     376
     377if ~exist('ModelList.tar.gz','file'),
     378        error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward');
     379end
     380
     381%upload both files to cluster
     382disp('uploading input file,  queueing script and variables script');
     383eval(['!scp ModelList.tar.gz ' name '.queue '  cluster ':' executionpath]);
     384
     385disp('launching solution sequence on remote cluster');
     386issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']);
     387end% }}}

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/YamsCall.m

     
     1function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
     2%YAMSCALL - call yams
     3%
     4%   build a metric using the Hessian of the given field
     5%   call Yams and the output mesh is plugged onto the model
     6%   -hmin = minimum edge length (m)
     7%   -hmax = maximum edge length (m)
     8%   -gradation = maximum edge length gradation between 2 elements
     9%   -epsilon = average error on each element (m/yr)
     10%
     11%   Usage:
     12%      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
     13%
     14%   Example:
     15%      md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
     16
     17%2d geometric parameter (do not change)
     18scale=2/9;
     19
     20%Compute Hessian
     21t1=clock; fprintf('%s','      computing Hessian...');
     22hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
     23t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     24
     25%Compute metric
     26t1=clock; fprintf('%s','      computing metric...');
     27if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
     28        pos=find(md.mask.vertexonwater);
     29else
     30        pos=[];
     31end
     32metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
     33t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     34
     35%write files
     36t1=clock; fprintf('%s','      writing initial mesh files...');
     37save -ascii carre0.met  metric
     38
     39fid=fopen('carre0.mesh','w');
     40
     41%initialiation
     42fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
     43
     44%dimension
     45fprintf(fid,'\n%s\n%i\n','Dimension',2);
     46
     47%Vertices
     48fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
     49fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
     50
     51%Triangles
     52fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
     53fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
     54numberofelements1=md.mesh.numberofelements;
     55       
     56%Deal with rifts
     57if ~isnan(md.rifts.riftstruct),
     58       
     59        %we have the list of triangles that make up the rift. keep those triangles around during refinement.
     60        triangles=[];
     61        for i=1:size(md.rifts.riftstruct,1),
     62                triangles=[triangles md.rifts(i).segments(:,3)'];
     63        end
     64
     65        fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
     66        fprintf(fid,'%i\n',triangles);
     67end
     68
     69%close
     70fclose(fid);
     71t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     72
     73%call yams
     74fprintf('%s\n','      call Yams...');
     75if ispc
     76        %windows
     77        system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
     78elseif ismac
     79        %Macosx
     80        system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
     81else
     82        %Linux
     83        system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
     84end
     85
     86%plug new mesh
     87t1=clock; fprintf('\n%s','      reading final mesh files...');
     88Tria=load('carre1.tria');
     89Coor=load('carre1.coor');
     90md.mesh.x=Coor(:,1);
     91md.mesh.y=Coor(:,2);
     92md.mesh.z=zeros(size(Coor,1),1);
     93md.mesh.elements=Tria;
     94md.mesh.numberofvertices=size(Coor,1);
     95md.mesh.numberofelements=size(Tria,1);
     96numberofelements2=md.mesh.numberofelements;
     97t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     98
     99%display number of elements
     100fprintf('\n%s %i','      inital number of elements:',numberofelements1);
     101fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
     102
     103%clean up:
     104system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/meshread.m

     
     1function Struct=meshread(filename);
     2
     3%some checks
     4if ~exist(filename),
     5        error(['meshread error message: file ' filename ' not found!']);
     6end
     7
     8fid=fopen(filename,'r');
     9
     10while (~feof(fid)),
     11
     12        A=fscanf(fid,'%s',1);
     13
     14        if strcmp(A,'MeshVersionFormatted');
     15                Struct.Version=fscanf(fid,'%s',1);
     16
     17        elseif strcmp(A,'Dimension'),
     18                Struct.Dimension=fscanf(fid,'%i',1);
     19
     20        elseif strcmp(A,'Vertices'),
     21                Struct.nods=fscanf(fid,'%i',1);
     22                A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
     23                Struct.x=A(1,:)';
     24                Struct.y=A(2,:)';
     25
     26        elseif strcmp(A,'Triangles'),
     27                Struct.nels=fscanf(fid,'%i',1);
     28                A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
     29                Struct.index=A(1:3,:)';
     30
     31        elseif strcmp(A,'Quadrilaterals'),
     32                Struct.nels=fscanf(fid,'%i',1);
     33                A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
     34                Struct.index=A(1:4,:)';
     35        else
     36                %do nothing
     37
     38        end
     39end
     40
     41fclose(fid);

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/yams.m

     
     1function md=yams(md,varargin);
     2%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
     3%
     4%   Usage:
     5%      md=yams(md,varargin);
     6%      where varargin is a lit of paired arguments.
     7%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
     8%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
     9%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
     10%                          this option is used to minimize the metric on water (no refinement)
     11%      optional arguments: 'resolution': initial mesh resolution [m]
     12%      optional arguments: 'nsteps': number of steps of mesh adaptation
     13%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
     14%      optional arguments: 'hmin': minimum edge length
     15%      optional arguments: 'hmanx': maximum edge
     16%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
     17%     
     18%
     19%   Examples:
     20%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
     21%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
     22%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
     23
     24%recover options
     25options=pairoptions(varargin{:});
     26options=deleteduplicates(options,1);
     27
     28%recover some fields
     29disp('MeshYams Options:')
     30domainoutline=getfieldvalue(options,'domainoutline');
     31disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
     32riftoutline=getfieldvalue(options,'riftoutline','N/A');
     33disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
     34groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
     35disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
     36velocities=getfieldvalue(options,'velocities');
     37disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
     38resolution=getfieldvalue(options,'resolution',5000);
     39disp(sprintf('   %-15s: %f','Resolution',resolution));
     40nsteps=getfieldvalue(options,'nsteps',6);
     41disp(sprintf('   %-15s: %i','nsteps',nsteps));
     42gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
     43disp(sprintf('   %-15s: %g','gradation',gradation(1)));
     44epsilon=getfieldvalue(options,'epsilon',3);
     45disp(sprintf('   %-15s: %f','epsilon',epsilon));
     46hmin=getfieldvalue(options,'hmin',500);
     47disp(sprintf('   %-15s: %f','hmin',hmin));
     48hmax=getfieldvalue(options,'hmax',150*10^3);
     49disp(sprintf('   %-15s: %f\n','hmax',hmax));
     50
     51%mesh with initial resolution
     52disp('Initial mesh generation...');
     53if strcmpi(riftoutline,'N/A');
     54        md=setmesh(md,domainoutline,resolution);
     55else
     56        md=setmesh(md,domainoutline,riftoutline,resolution);
     57        md=meshprocessrifts(md,domainoutline);
     58end
     59disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
     60
     61%load velocities
     62disp('loading velocities...');
     63Names=VelFindVarNames(velocities);
     64Vel=load(velocities);
     65
     66%start mesh adaptation
     67for i=1:nsteps,
     68        disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
     69
     70        %interpolate velocities onto mesh
     71        disp('   interpolating velocities...');
     72        if strcmpi(Names.interp,'node'),
     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        else
     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        end
     79        field=sqrt(vx_obs.^2+vy_obs.^2);
     80
     81        %set mask.vertexonwater  field
     82        if ~strcmp(groundeddomain,'N/A'),
     83                nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
     84                md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
     85                md.mask.vertexonwater(find(nodeground))=0;
     86        else
     87                md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
     88        end
     89
     90        %adapt according to velocities
     91        disp('   adapting...');
     92        md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
     93
     94        %if we have rifts, we just messed them up, we need to recreate the segments that constitute those
     95        %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
     96        if md.rifts.numrifts,
     97                md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
     98                md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
     99                md.mesh.segments=findsegments(md);
     100                md=yamsrecreateriftsegments(md);
     101        end
     102
     103end
     104       
     105disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
     106
     107%Now, build the connectivity tables for this mesh.
     108md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
     109md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
     110
     111%recreate segments
     112md.mesh.segments=findsegments(md);
     113md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
     114
     115%Fill in rest of fields:
     116md.mesh.z=zeros(md.mesh.numberofvertices,1);
     117md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
     118md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
     119md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
     120md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
     121if ~strcmp(groundeddomain,'N/A'),
     122        nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
     123        md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
     124        md.mask.vertexonwater(find(nodeground))=0;
     125else
     126        md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
     127end
     128if strcmpi(Names.interp,'node'),
     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);
     131else
     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);
     134end
     135md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
     136
     137%deal with rifts
     138if md.rifts.numrifts,
     139        %first, recreate rift segments
     140        md=meshyamsrecreateriftsegments(md);
     141
     142        %using the segments, recreate the penaltypairs
     143        for j=1:md.rifts.numrifts,
     144                rift=md.rifts.riftstruct(j);
     145
     146                %build normals and lengths of segments:
     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
     151                %ok, build penaltypairs:
     152                numpenaltypairs=length(rift.segments)/2-1;
     153                rift.penaltypairs=zeros(numpenaltypairs,7);
     154
     155                for i=1:numpenaltypairs,
     156                        rift.penaltypairs(i,1)=rift.segments(i,2);
     157                        rift.penaltypairs(i,2)=rift.segments(end-i,2);
     158                        rift.penaltypairs(i,3)=rift.segments(i,3);
     159                        rift.penaltypairs(i,4)=rift.segments(end-i,3);
     160                        rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
     161                        rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
     162                        rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
     163                end
     164                %renormalize norms:
     165                norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
     166                rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
     167                rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
     168
     169                md.rifts.riftstruct(j)=rift;
     170        end
     171end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m

     
     1function md=BamgCallFromMetric(md,metric,gradation),
     2%BAMGCALL - call bam
     3%
     4%   call Bamg and the output mesh is plugged onto the model
     5%   -gradation = maximum edge length gradation between 2 elements
     6%
     7%   Usage:
     8%      md=BamgCallFromMetric(md,metric,gradation);
     9%
     10%   Example:
     11%      md=BamgCall(md,metric,1500,10^8,1.3,0.9);
     12
     13%2d geometric parameter (do not change)
     14scale=2/9;
     15
     16%write files
     17t1=clock; fprintf('%s','      writing initial mesh files...');
     18fid=fopen('carre0.met','w');
     19fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
     20fprintf(fid,'%i %i %i\n',metric');
     21fclose(fid);
     22
     23fid=fopen('carre0.mesh','w');
     24
     25%initialiation
     26fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
     27
     28%dimension
     29fprintf(fid,'\n%s\n%i\n','Dimension',2);
     30
     31%Vertices
     32fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
     33fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
     34
     35%Triangles
     36fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
     37fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
     38numberofelements1=md.mesh.numberofelements;
     39
     40%close
     41fclose(fid);
     42t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     43
     44%call bamg
     45fprintf('%s\n','      call Bamg...');
     46system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
     47
     48%plug new mesh
     49t1=clock; fprintf('\n%s','      reading final mesh files...');
     50A=meshread('carre1.mesh');
     51md.mesh.x=A.x;
     52md.mesh.y=A.y;
     53md.z=zeros(A.nods,1);
     54md.mesh.elements=A.index;
     55md.mesh.numberofvertices=A.nods;
     56md.mesh.numberofelements=A.nels;
     57numberofelements2=md.mesh.numberofelements;
     58t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     59
     60%display number of elements
     61fprintf('\n%s %i','      inital number of elements:',numberofelements1);
     62fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
     63
     64%clean up:
     65system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCall.m

     
     1function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
     2%BAMGCALL - call bam
     3%
     4%   build a metric using the Hessian of the given field
     5%   call Bamg and the output mesh is plugged onto the model
     6%   -hmin = minimum edge length (m)
     7%   -hmax = maximum edge length (m)
     8%   -gradation = maximum edge length gradation between 2 elements
     9%   -epsilon = average error on each element (m/yr)
     10%
     11%   Usage:
     12%      md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
     13%
     14%   Example:
     15%      md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
     16
     17%2d geometric parameter (do not change)
     18scale=2/9;
     19
     20%Compute Hessian
     21t1=clock; fprintf('%s','      computing Hessian...');
     22hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
     23t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     24
     25%Compute metric
     26t1=clock; fprintf('%s','      computing metric...');
     27if length(md.nodeonwater)==md.mesh.numberofvertices,
     28        pos=find(md.nodeonwater);
     29else
     30        pos=[];
     31end
     32metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
     33t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     34
     35%write files
     36t1=clock; fprintf('%s','      writing initial mesh files...');
     37fid=fopen('carre0.met','w');
     38fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
     39fprintf(fid,'%i %i %i\n',metric');
     40fclose(fid);
     41
     42fid=fopen('carre0.mesh','w');
     43
     44%initialiation
     45fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
     46
     47%dimension
     48fprintf(fid,'\n%s\n%i\n','Dimension',2);
     49
     50%Vertices
     51fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
     52fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
     53
     54%Triangles
     55fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
     56fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
     57numberofelements1=md.mesh.numberofelements;
     58
     59%close
     60fclose(fid);
     61t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     62
     63%call bamg
     64fprintf('%s\n','      call Bamg...');
     65system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
     66
     67%plug new mesh
     68t1=clock; fprintf('\n%s','      reading final mesh files...');
     69A=meshread('carre1.mesh');
     70md.mesh.x=A.x;
     71md.mesh.y=A.y;
     72md.z=zeros(A.nods,1);
     73md.mesh.elements=A.index;
     74md.mesh.numberofvertices=A.nods;
     75md.mesh.numberofelements=A.nels;
     76numberofelements2=md.mesh.numberofelements;
     77t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
     78
     79%display number of elements
     80fprintf('\n%s %i','      inital number of elements:',numberofelements1);
     81fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
     82
     83%clean up:
     84system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/create_region.m

     
    1 function create_region(name)
    2 %CREATE_REGION - create region  ????
    3 %
    4 %   very temporary function.
    5 %   
    6 %   Usage:
    7 %      create_region(name)
    8 
    9 eval(['mkdir ' name]);
    10 eval(['cd ' name ]);
    11 !mkdir Delivery Exp_Par Results
    12 cd Exp_Par
    13 !cp ../../RonneShelf/Exp_Par/* ./
    14 !rm -rf Hole*
    15 eval(['!mv Ronne.par ' name '.par']);

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/findarg.m

     
    1 function  vals=findarg(arglist,field)
    2 %FINDARG - find argument associated to a field in a list
    3 %
    4 %   This function parses through an argument list (typically varargin in a routine)
    5 %   looking for a character array equal to field. Once this is found, we return the
    6 %   next value in the varargin (if possible).
    7 %   Because field might appear several times in the argument list, we return a structure
    8 %   holding all these values.
    9 %   Note that all comparisons to field value are case independent.
    10 %
    11 %   Usage:
    12 %      vals=findarg(arglist,field)
    13 %
    14 %   Example:
    15 %      routine foobar calls vals=findarg('Data',varargin)
    16 %      with varargin='Data',1,'Data','foo','Plot','velocity','Arrow',4
    17 %      findarg would return the following structure: vals(1).value=1, vals(2).value='foo';
    18 
    19 %some argument checking:
    20 if ((nargin==0) | (nargout==0)),
    21         help findarg;
    22         error('findarg error message');
    23 end
    24 
    25 if ~ischar(field),
    26         error('findarg error message: field should be a string');
    27 end
    28 
    29 if ~iscell(arglist),
    30         error('findarg error message: argument list should be a cell array.');
    31 end
    32 
    33 %Recover data to plot
    34 founddata=0;
    35 
    36 for i=1:(length(arglist)-1), %data in arglist comes in pairs, hence the -1.
    37         if ischar(arglist{i}),
    38                 if (strcmpi(arglist{i},field)),
    39                         founddata=founddata+1;
    40                         if founddata==1,
    41                                 vals.value=arglist{i+1};
    42                         else
    43                                 vals(end+1).value=arglist{i+1};
    44                         end
    45                 end
    46         end
    47 end
    48 
    49 if founddata==0,
    50         vals=[];
    51 end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf.m

     
    1 function S = netcdf(File,varargin)
    2 % Function to read NetCDF files
    3 %   S = netcdf(File)
    4 % Input Arguments
    5 %   File = NetCDF file to read
    6 % Optional Input Arguments:
    7 %   'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
    8 %   'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
    9 % Output Arguments:
    10 %   S    = Structure of NetCDF data organised as per NetCDF definition
    11 % Notes:
    12 %   Only version 1, classic 32bit, NetCDF files are supported. By default
    13 % data are extracted into the S.VarArray().Data field for all variables.
    14 % To read the header only call S = netcdf(File,'Var',[]);
    15 %
    16 % SEE ALSO
    17 % ---------------------------------------------------------------------------
    18 S = [];
    19 
    20 try
    21    if exist(File,'file') fp = fopen(File,'r','b');
    22    else fp = []; error('File not found'); end
    23    if fp == -1   error('Unable to open file'); end
    24 
    25 % Read header
    26    Magic = fread(fp,4,'uint8=>char');
    27    if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
    28    if uint8(Magic(4))~=1       error('Version not supported'); end
    29    S.NumRecs  = fread(fp,1,'uint32=>uint32');
    30    S.DimArray = DimArray(fp);
    31    S.AttArray = AttArray(fp);
    32    S.VarArray = VarArray(fp);
    33 
    34 % Setup indexing to arrays and records
    35    Var = ones(1,length(S.VarArray));
    36    Rec = ones(1,S.NumRecs);
    37    for i = 1:2:length(varargin)
    38       if     strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
    39       elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
    40       else error('Optional input argument not recognised'); end
    41    end
    42    if sum(Var)==0 fclose(fp); return; end
    43 
    44 % Read non-record variables
    45    Dim = double(cat(2,S.DimArray.Dim));
    46    ID  = double(cat(2,S.VarArray.Type));
    47 
    48    for i = 1:length(S.VarArray)
    49       D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
    50       if isempty(RecID{i})
    51          if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
    52          if Var(i)
    53             S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
    54             fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
    55          else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
    56       else S.VarArray(i).Data = []; end
    57    end
    58 
    59 % Read record variables
    60    for k = 1:S.NumRecs
    61       for i = 1:length(S.VarArray)
    62          if ~isempty(RecID{i})
    63             D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
    64             if length(D)==1 D=[D,1]; end
    65             if Var(i) & Rec(k)
    66                S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
    67                   ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
    68                if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
    69             else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
    70          end
    71       end
    72    end
    73 
    74    fclose(fp);
    75 catch
    76    Err = lasterror; fprintf('%s\n',Err.message);
    77    if ~isempty(fp) && fp ~= -1 fclose(fp); end
    78 end
    79 
    80 % ---------------------------------------------------------------------------------------
    81 % Utility functions
    82 
    83 function S = Size(ID)
    84 % Size of NetCDF data type, ID, in bytes
    85    S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
    86 
    87 function T = Type(ID)
    88 % Matlab string for CDF data type, ID
    89    T = subsref({'int8','char','int16','int32','single','double'},...
    90                struct('type','{}','subs',{{ID}}));
    91 
    92 function N = Pad(Num,ID)
    93 % Number of elements to read after padding to 4 bytes for type ID
    94    N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
    95 
    96 function S = String(fp)
    97 % Read a CDF string; Size,[String,[Padding]]
    98    S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
    99 
    100 function A = ReOrder(A,S)
    101 % Rearrange CDF array A to size S with matlab ordering
    102    A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
    103 
    104 function S = DimArray(fp)
    105 % Read DimArray into structure
    106    if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
    107       for i = 1:fread(fp,1,'uint32=>uint32')
    108          S(i).Str = String(fp);
    109          S(i).Dim = fread(fp,1,'uint32=>uint32');
    110       end
    111    else fread(fp,1,'uint32=>uint32'); S = []; end
    112 
    113 function S = AttArray(fp)
    114 % Read AttArray into structure
    115    if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
    116       for i = 1:fread(fp,1,'uint32=>uint32')
    117          S(i).Str = String(fp);
    118          ID       = fread(fp,1,'uint32=>uint32');
    119          Num      = fread(fp,1,'uint32=>uint32');
    120          S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
    121       end
    122    else fread(fp,1,'uint32=>uint32'); S = []; end
    123 
    124 function S = VarArray(fp)
    125 % Read VarArray into structure
    126    if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
    127       for i = 1:fread(fp,1,'uint32=>uint32')
    128          S(i).Str      = String(fp);
    129          Num           = double(fread(fp,1,'uint32=>uint32'));
    130          S(i).DimID    = double(fread(fp,Num,'uint32=>uint32'));
    131          S(i).AttArray = AttArray(fp);
    132          S(i).Type     = fread(fp,1,'uint32=>uint32');
    133          S(i).VSize    = fread(fp,1,'uint32=>uint32');
    134          S(i).Begin    = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
    135       end
    136    else fread(fp,1,'uint32=>uint32'); S = []; end

  • u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf2struct.m

     
    55%      S=netcdf2struct(File);
    66
    77%Read netcdf file
    8 data=netcdf(File);
     8data=readnetcdf(File);
    99
    1010%initialize output
    1111S=struct();
    … …  
    2525        fieldvalue=double(variables(i).Val);
    2626        S.(fieldname)=fieldvalue;
    2727end
     28end
     29
     30function S = readnetcdf(File,varargin)
     31% Function to read NetCDF files
     32%   S = netcdf(File)
     33% Input Arguments
     34%   File = NetCDF file to read
     35% Optional Input Arguments:
     36%   'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
     37%   'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
     38% Output Arguments:
     39%   S    = Structure of NetCDF data organised as per NetCDF definition
     40% Notes:
     41%   Only version 1, classic 32bit, NetCDF files are supported. By default
     42% data are extracted into the S.VarArray().Data field for all variables.
     43% To read the header only call S = netcdf(File,'Var',[]);
     44%
     45% SEE ALSO
     46% ---------------------------------------------------------------------------
     47S = [];
     48
     49try
     50   if exist(File,'file') fp = fopen(File,'r','b');
     51   else fp = []; error('File not found'); end
     52   if fp == -1   error('Unable to open file'); end
     53
     54% Read header
     55   Magic = fread(fp,4,'uint8=>char');
     56   if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
     57   if uint8(Magic(4))~=1       error('Version not supported'); end
     58   S.NumRecs  = fread(fp,1,'uint32=>uint32');
     59   S.DimArray = DimArray(fp);
     60   S.AttArray = AttArray(fp);
     61   S.VarArray = VarArray(fp);
     62
     63% Setup indexing to arrays and records
     64   Var = ones(1,length(S.VarArray));
     65   Rec = ones(1,S.NumRecs);
     66   for i = 1:2:length(varargin)
     67      if     strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
     68      elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
     69      else error('Optional input argument not recognised'); end
     70   end
     71   if sum(Var)==0 fclose(fp); return; end
     72
     73% Read non-record variables
     74   Dim = double(cat(2,S.DimArray.Dim));
     75   ID  = double(cat(2,S.VarArray.Type));
     76
     77   for i = 1:length(S.VarArray)
     78      D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
     79      if isempty(RecID{i})
     80         if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
     81         if Var(i)
     82            S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
     83            fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
     84         else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
     85      else S.VarArray(i).Data = []; end
     86   end
     87
     88% Read record variables
     89   for k = 1:S.NumRecs
     90      for i = 1:length(S.VarArray)
     91         if ~isempty(RecID{i})
     92            D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
     93            if length(D)==1 D=[D,1]; end
     94            if Var(i) & Rec(k)
     95               S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
     96                  ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
     97               if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
     98            else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
     99         end
     100      end
     101   end
     102
     103   fclose(fp);
     104catch
     105   Err = lasterror; fprintf('%s\n',Err.message);
     106   if ~isempty(fp) && fp ~= -1 fclose(fp); end
     107end
     108
     109% ---------------------------------------------------------------------------------------
     110% Utility functions
     111
     112function S = Size(ID)
     113% Size of NetCDF data type, ID, in bytes
     114   S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
     115
     116function T = Type(ID)
     117% Matlab string for CDF data type, ID
     118   T = subsref({'int8','char','int16','int32','single','double'},...
     119               struct('type','{}','subs',{{ID}}));
     120
     121function N = Pad(Num,ID)
     122% Number of elements to read after padding to 4 bytes for type ID
     123   N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
     124
     125function S = String(fp)
     126% Read a CDF string; Size,[String,[Padding]]
     127   S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
     128
     129function A = ReOrder(A,S)
     130% Rearrange CDF array A to size S with matlab ordering
     131   A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
     132
     133function S = DimArray(fp)
     134% Read DimArray into structure
     135   if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
     136      for i = 1:fread(fp,1,'uint32=>uint32')
     137         S(i).Str = String(fp);
     138         S(i).Dim = fread(fp,1,'uint32=>uint32');
     139      end
     140   else fread(fp,1,'uint32=>uint32'); S = []; end
     141
     142function S = AttArray(fp)
     143% Read AttArray into structure
     144   if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
     145      for i = 1:fread(fp,1,'uint32=>uint32')
     146         S(i).Str = String(fp);
     147         ID       = fread(fp,1,'uint32=>uint32');
     148         Num      = fread(fp,1,'uint32=>uint32');
     149         S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
     150      end
     151   else fread(fp,1,'uint32=>uint32'); S = []; end
     152
     153function S = VarArray(fp)
     154% Read VarArray into structure
     155   if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
     156      for i = 1:fread(fp,1,'uint32=>uint32')
     157         S(i).Str      = String(fp);
     158         Num           = double(fread(fp,1,'uint32=>uint32'));
     159         S(i).DimID    = double(fread(fp,Num,'uint32=>uint32'));
     160         S(i).AttArray = AttArray(fp);
     161         S(i).Type     = fread(fp,1,'uint32=>uint32');
     162         S(i).VSize    = fread(fp,1,'uint32=>uint32');
     163         S(i).Begin    = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
     164      end
     165   else fread(fp,1,'uint32=>uint32'); S = []; end
     166end
Note: See TracBrowser for help on using the repository browser.