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
RevLine 
[13394]1Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJob.m
2===================================================================
3--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJob.m (revision 13011)
4+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJob.m (revision 13012)
5@@ -1,21 +0,0 @@
6-function LaunchMultipleQueueJob(cluster,name,executionpath)
7-%LAUNCHMULTIPLEQUEUEJOB - ...
8-%
9-% Usage:
10-% LaunchMultipleQueueJob(executionpath)
11-
12-%First try and figure out if there is a special script for thie particular cluster
13-function_name=['LaunchMultipleQueueJob' cluster];
14-
15-%some specific treatment of identical cluster, gemini, castor and pollux
16-if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
17- function_name='LaunchMultipleQueueJobgemini';
18-end
19-
20-if exist(function_name,'file'),
21- %Call this function:
22- eval([function_name '(cluster,name,executionpath);']);
23-else
24- %Call the generic LaunchMultipleQueueJob:
25- LaunchMultipleQueueJobGeneric(cluster,name,executionpath);
26-end
27Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptgemini.m
28===================================================================
29--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptgemini.m (revision 13011)
30+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptgemini.m (revision 13012)
31@@ -1,23 +0,0 @@
32-function BuildQueueingScriptgemini(name,executionpath,codepath)
33-%BUILDQUEUEINGSCRIPTGEMINI - ...
34-%
35-% Usage:
36-% BuildQueueingScriptgemini(md,executionpath,codepath)
37-
38-scriptname=[name '.queue'];
39-
40-fid=fopen(scriptname,'w');
41-if fid==-1,
42- error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']);
43-end
44-
45-fprintf(fid,'#!/bin/sh\n');
46-fprintf(fid,'cd %s\n',executionpath);
47-fprintf(fid,'mkdir %s\n',name);
48-fprintf(fid,'cd %s\n',name);
49-fprintf(fid,'mv ../ModelList.tar.gz ./\n');
50-fprintf(fid,'tar -zxvf ModelList.tar.gz\n');
51-fprintf(fid,'foreach i (%s-*vs*.queue)\n',name);
52-fprintf(fid,'qsub $i\n');
53-fprintf(fid,'end\n');
54-fclose(fid);
55Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobgemini.m
56===================================================================
57--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobgemini.m (revision 13011)
58+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobgemini.m (revision 13012)
59@@ -1,22 +0,0 @@
60-function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)
61-%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster
62-%
63-% Usage:
64-% LaunchMultipleQueueJobgemini(cluster,name,executionpath)
65-
66-
67-%first, check we have the binary file and the queueing script
68-if ~exist([ name '.queue'],'file'),
69- error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward');
70-end
71-
72-if ~exist('ModelList.tar.gz','file'),
73- error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward');
74-end
75-
76-%upload both files to cluster
77-disp('uploading input file, queueing script and variables script');
78-eval(['!scp ModelList.tar.gz ' name '.queue ' cluster ':' executionpath]);
79-
80-disp('launching solution sequence on remote cluster');
81-issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']);
82Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptGeneric.m
83===================================================================
84--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptGeneric.m (revision 13011)
85+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScriptGeneric.m (revision 13012)
86@@ -1,9 +0,0 @@
87-function BuildMultipleQueueingScriptGeneric(name,executionpath,codepath)
88-%BUILDMULTIPLEQUEUEINGSCRIPTGENERIC - ...
89-%
90-% Usage:
91-% BuildMultipleQueueingScriptGeneric(executionpath,codepath)
92-
93-%not done yet
94-error('BuildMultipleQueueingScriptGenericerror message: not supported yet!');
95-
96Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScript.m
97===================================================================
98--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScript.m (revision 13011)
99+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/BuildMultipleQueueingScript.m (revision 13012)
100@@ -1,23 +0,0 @@
101-function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)
102-%BUILDMULTIPLEQUEUEINGSCRIPT -
103-%
104-% Usage:
105-% BuildMultipleQueueingScript(executionpath,codepath)
106-
107-disp('building queueing script');
108-
109-%First try and figure out if there is a special script for this particular cluster
110-function_name=['BuildMultipleQueueingScript' cluster];
111-
112-%some specific treatment of identical cluster, gemini, castor and pollux
113-if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
114- function_name='BuildMultipleQueueingScriptgemini';
115-end
116-
117-if exist(function_name,'file'),
118- %Call this function:
119- eval([function_name '(name,executionpath,codepath);']);
120-else
121- %Call the generic BuildQueueingScript:
122- BuildMultipleQueueingScriptGeneric(name,executionpath,codepath);
123-end
124Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobGeneric.m
125===================================================================
126--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobGeneric.m (revision 13011)
127+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/multiplequeue/LaunchMultipleQueueJobGeneric.m (revision 13012)
128@@ -1,7 +0,0 @@
129-function LaunchMultipleQueueJobGeneric(cluster,name,executionpath)
130-%LAUNCHMULTIPLEQUEUEJOBGENERIC - Generic routine to launch multiple queueing job
131-%
132-% Usage:
133-% LaunchMultipleQueueJobGeneric(cluster,name,executionpath)
134-
135-error('LaunchMultipleQueueJobGeneric error message: not supported yet!');
136Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/geometry/GetAreas.m
137===================================================================
138--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/geometry/GetAreas.m (revision 0)
139+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/geometry/GetAreas.m (revision 13012)
140@@ -0,0 +1,51 @@
141+function areas=GetAreas(index,x,y,varargin)
142+%GETAREAS - compute areas or volumes of elements
143+%
144+% compute areas of triangular elements or volumes
145+% of pentahedrons
146+%
147+% Usage:
148+% areas =GetAreas(index,x,y);
149+% volumes=GetAreas(index,x,y,z);
150+%
151+% Examples:
152+% areas =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
153+% volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
154+
155+%get number of elements and number of nodes
156+nels=size(index,1);
157+nods=length(x);
158+if nargin==4, z=varargin{1}; end
159+
160+%some checks
161+if nargout~=1 | (nargin~=3 & nargin~=4),
162+ help GetAreas
163+ error('GetAreas error message: bad usage')
164+end
165+if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
166+ error('GetAreas error message: x,y and z do not have the same length')
167+end
168+if max(index(:))>nods,
169+ error(['GetAreas error message: index should not have values above ' num2str(nods) ])
170+end
171+if (nargin==3 & size(index,2)~=3),
172+ error('GetAreas error message: index should have 3 columns for 2d meshes.')
173+end
174+if (nargin==4 & size(index,2)~=6),
175+ error('GetAreas error message: index should have 6 columns for 3d meshes.')
176+end
177+
178+%initialization
179+areas=zeros(nels,1);
180+x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
181+y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
182+
183+%compute the volume of each element
184+if nargin==3,
185+ %compute the surface of the triangle
186+ areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
187+else
188+ %V=area(triangle)*1/3(z1+z2+z3)
189+ thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
190+ areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
191+end
192Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/meshread.m
193===================================================================
194--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/meshread.m (revision 13011)
195+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/meshread.m (revision 13012)
196@@ -1,41 +0,0 @@
197-function Struct=meshread(filename);
198-
199-%some checks
200-if ~exist(filename),
201- error(['meshread error message: file ' filename ' not found!']);
202-end
203-
204-fid=fopen(filename,'r');
205-
206-while (~feof(fid)),
207-
208- A=fscanf(fid,'%s',1);
209-
210- if strcmp(A,'MeshVersionFormatted');
211- Struct.Version=fscanf(fid,'%s',1);
212-
213- elseif strcmp(A,'Dimension'),
214- Struct.Dimension=fscanf(fid,'%i',1);
215-
216- elseif strcmp(A,'Vertices'),
217- Struct.nods=fscanf(fid,'%i',1);
218- A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
219- Struct.x=A(1,:)';
220- Struct.y=A(2,:)';
221-
222- elseif strcmp(A,'Triangles'),
223- Struct.nels=fscanf(fid,'%i',1);
224- A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
225- Struct.index=A(1:3,:)';
226-
227- elseif strcmp(A,'Quadrilaterals'),
228- Struct.nels=fscanf(fid,'%i',1);
229- A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
230- Struct.index=A(1:4,:)';
231- else
232- %do nothing
233-
234- end
235-end
236-
237-fclose(fid);
238Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCallFromMetric.m
239===================================================================
240--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCallFromMetric.m (revision 13011)
241+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCallFromMetric.m (revision 13012)
242@@ -1,65 +0,0 @@
243-function md=BamgCallFromMetric(md,metric,gradation),
244-%BAMGCALL - call bam
245-%
246-% call Bamg and the output mesh is plugged onto the model
247-% -gradation = maximum edge length gradation between 2 elements
248-%
249-% Usage:
250-% md=BamgCallFromMetric(md,metric,gradation);
251-%
252-% Example:
253-% md=BamgCall(md,metric,1500,10^8,1.3,0.9);
254-
255-%2d geometric parameter (do not change)
256-scale=2/9;
257-
258-%write files
259-t1=clock; fprintf('%s',' writing initial mesh files...');
260-fid=fopen('carre0.met','w');
261-fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
262-fprintf(fid,'%i %i %i\n',metric');
263-fclose(fid);
264-
265-fid=fopen('carre0.mesh','w');
266-
267-%initialiation
268-fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
269-
270-%dimension
271-fprintf(fid,'\n%s\n%i\n','Dimension',2);
272-
273-%Vertices
274-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
275-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
276-
277-%Triangles
278-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
279-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
280-numberofelements1=md.mesh.numberofelements;
281-
282-%close
283-fclose(fid);
284-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
285-
286-%call bamg
287-fprintf('%s\n',' call Bamg...');
288-system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
289-
290-%plug new mesh
291-t1=clock; fprintf('\n%s',' reading final mesh files...');
292-A=meshread('carre1.mesh');
293-md.mesh.x=A.x;
294-md.mesh.y=A.y;
295-md.z=zeros(A.nods,1);
296-md.mesh.elements=A.index;
297-md.mesh.numberofvertices=A.nods;
298-md.mesh.numberofelements=A.nels;
299-numberofelements2=md.mesh.numberofelements;
300-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
301-
302-%display number of elements
303-fprintf('\n%s %i',' inital number of elements:',numberofelements1);
304-fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
305-
306-%clean up:
307-system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
308Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetCharacteristicLength.m
309===================================================================
310--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetCharacteristicLength.m (revision 13011)
311+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetCharacteristicLength.m (revision 13012)
312@@ -1,45 +0,0 @@
313-function length=GetCharacteristicLength(index,x,y,varargin)
314-%GETCHARACTERISTICLENGTH - compute characteristic length for a mesh
315-%
316-% compute characteristic lengths of every element of a mesh.
317-%
318-% Usage:
319-% length =GetCharacteristicLength(index,x,y);
320-% length =GetCharacteristicLength(index,x,y,z);
321-%
322-% Examples:
323-% length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y);
324-% length =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
325-
326-
327-%get number of elements and number of nodes
328-nels=size(index,1);
329-nods=numel(x);
330-
331-%some checks
332-if nargout~=1 | (nargin~=3 & nargin~=4),
333- help GetCharacteristicLength
334- error('GetCharacteristicLength error message: bad usage')
335-end
336-if ((numel(y)~=nods) | (nargin==4 & numel(z)~=nods)),
337- error('GetCharacteristicLength error message: x,y and z do not have the same length')
338-end
339-if max(index(:))>nods,
340- error(['GetCharacteristicLength error message: index should not have values above ' num2str(nods) ])
341-end
342-if (nargin==3 & size(index,2)~=3),
343- error('GetCharacteristicLength error message: index should have 3 columns for 2d meshes.')
344-end
345-if (nargin==4 & size(index,2)~=6),
346- error('GetCharacteristicLength error message: index should have 6 columns for 3d meshes.')
347-end
348-
349-%get areas or volumes.
350-areas=GetAreas(index,x,y,varargin{:});
351-
352-%for a 2d mesh:
353-if nargin==3,
354- length=sqrt(2*areas);
355-else
356- error('not supported yet');
357-end
358Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetAreas.m
359===================================================================
360--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetAreas.m (revision 13011)
361+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/GetAreas.m (revision 13012)
362@@ -1,51 +0,0 @@
363-function areas=GetAreas(index,x,y,varargin)
364-%GETAREAS - compute areas or volumes of elements
365-%
366-% compute areas of triangular elements or volumes
367-% of pentahedrons
368-%
369-% Usage:
370-% areas =GetAreas(index,x,y);
371-% volumes=GetAreas(index,x,y,z);
372-%
373-% Examples:
374-% areas =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
375-% volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
376-
377-%get number of elements and number of nodes
378-nels=size(index,1);
379-nods=length(x);
380-if nargin==4, z=varargin{1}; end
381-
382-%some checks
383-if nargout~=1 | (nargin~=3 & nargin~=4),
384- help GetAreas
385- error('GetAreas error message: bad usage')
386-end
387-if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
388- error('GetAreas error message: x,y and z do not have the same length')
389-end
390-if max(index(:))>nods,
391- error(['GetAreas error message: index should not have values above ' num2str(nods) ])
392-end
393-if (nargin==3 & size(index,2)~=3),
394- error('GetAreas error message: index should have 3 columns for 2d meshes.')
395-end
396-if (nargin==4 & size(index,2)~=6),
397- error('GetAreas error message: index should have 6 columns for 3d meshes.')
398-end
399-
400-%initialization
401-areas=zeros(nels,1);
402-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
403-y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
404-
405-%compute the volume of each element
406-if nargin==3,
407- %compute the surface of the triangle
408- areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
409-else
410- %V=area(triangle)*1/3(z1+z2+z3)
411- thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
412- areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
413-end
414Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/isconnected.m
415===================================================================
416--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/isconnected.m (revision 13011)
417+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/isconnected.m (revision 13012)
418@@ -1,13 +0,0 @@
419-function flag=isconnected(elements,A,B)
420-%ISCONNECTED: are two nodes connected by a triangulation?
421-%
422-% Usage: flag=isconnected(elements,A,B)
423-%
424-%
425-
426-elements=ElementsFromEdge(elements,A,B);
427-if isempty(elements),
428- flag=0;
429-else
430- flag=1;
431-end
432Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCall.m
433===================================================================
434--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCall.m (revision 13011)
435+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/BamgCall.m (revision 13012)
436@@ -1,84 +0,0 @@
437-function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
438-%BAMGCALL - call bam
439-%
440-% build a metric using the Hessian of the given field
441-% call Bamg and the output mesh is plugged onto the model
442-% -hmin = minimum edge length (m)
443-% -hmax = maximum edge length (m)
444-% -gradation = maximum edge length gradation between 2 elements
445-% -epsilon = average error on each element (m/yr)
446-%
447-% Usage:
448-% md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
449-%
450-% Example:
451-% md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
452-
453-%2d geometric parameter (do not change)
454-scale=2/9;
455-
456-%Compute Hessian
457-t1=clock; fprintf('%s',' computing Hessian...');
458-hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
459-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
460-
461-%Compute metric
462-t1=clock; fprintf('%s',' computing metric...');
463-if length(md.nodeonwater)==md.mesh.numberofvertices,
464- pos=find(md.nodeonwater);
465-else
466- pos=[];
467-end
468-metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
469-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
470-
471-%write files
472-t1=clock; fprintf('%s',' writing initial mesh files...');
473-fid=fopen('carre0.met','w');
474-fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
475-fprintf(fid,'%i %i %i\n',metric');
476-fclose(fid);
477-
478-fid=fopen('carre0.mesh','w');
479-
480-%initialiation
481-fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
482-
483-%dimension
484-fprintf(fid,'\n%s\n%i\n','Dimension',2);
485-
486-%Vertices
487-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
488-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
489-
490-%Triangles
491-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
492-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
493-numberofelements1=md.mesh.numberofelements;
494-
495-%close
496-fclose(fid);
497-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
498-
499-%call bamg
500-fprintf('%s\n',' call Bamg...');
501-system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
502-
503-%plug new mesh
504-t1=clock; fprintf('\n%s',' reading final mesh files...');
505-A=meshread('carre1.mesh');
506-md.mesh.x=A.x;
507-md.mesh.y=A.y;
508-md.z=zeros(A.nods,1);
509-md.mesh.elements=A.index;
510-md.mesh.numberofvertices=A.nods;
511-md.mesh.numberofelements=A.nels;
512-numberofelements2=md.mesh.numberofelements;
513-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
514-
515-%display number of elements
516-fprintf('\n%s %i',' inital number of elements:',numberofelements1);
517-fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
518-
519-%clean up:
520-system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
521Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/YamsCall.m
522===================================================================
523--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/YamsCall.m (revision 13011)
524+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/YamsCall.m (revision 13012)
525@@ -1,104 +0,0 @@
526-function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
527-%YAMSCALL - call yams
528-%
529-% build a metric using the Hessian of the given field
530-% call Yams and the output mesh is plugged onto the model
531-% -hmin = minimum edge length (m)
532-% -hmax = maximum edge length (m)
533-% -gradation = maximum edge length gradation between 2 elements
534-% -epsilon = average error on each element (m/yr)
535-%
536-% Usage:
537-% md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
538-%
539-% Example:
540-% md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
541-
542-%2d geometric parameter (do not change)
543-scale=2/9;
544-
545-%Compute Hessian
546-t1=clock; fprintf('%s',' computing Hessian...');
547-hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
548-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
549-
550-%Compute metric
551-t1=clock; fprintf('%s',' computing metric...');
552-if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
553- pos=find(md.mask.vertexonwater);
554-else
555- pos=[];
556-end
557-metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
558-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
559-
560-%write files
561-t1=clock; fprintf('%s',' writing initial mesh files...');
562-save -ascii carre0.met metric
563-
564-fid=fopen('carre0.mesh','w');
565-
566-%initialiation
567-fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
568-
569-%dimension
570-fprintf(fid,'\n%s\n%i\n','Dimension',2);
571-
572-%Vertices
573-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
574-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
575-
576-%Triangles
577-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
578-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
579-numberofelements1=md.mesh.numberofelements;
580-
581-%Deal with rifts
582-if ~isnan(md.rifts.riftstruct),
583-
584- %we have the list of triangles that make up the rift. keep those triangles around during refinement.
585- triangles=[];
586- for i=1:size(md.rifts.riftstruct,1),
587- triangles=[triangles md.rifts(i).segments(:,3)'];
588- end
589-
590- fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
591- fprintf(fid,'%i\n',triangles);
592-end
593-
594-%close
595-fclose(fid);
596-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
597-
598-%call yams
599-fprintf('%s\n',' call Yams...');
600-if ispc
601- %windows
602- system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
603-elseif ismac
604- %Macosx
605- system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
606-else
607- %Linux
608- system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
609-end
610-
611-%plug new mesh
612-t1=clock; fprintf('\n%s',' reading final mesh files...');
613-Tria=load('carre1.tria');
614-Coor=load('carre1.coor');
615-md.mesh.x=Coor(:,1);
616-md.mesh.y=Coor(:,2);
617-md.mesh.z=zeros(size(Coor,1),1);
618-md.mesh.elements=Tria;
619-md.mesh.numberofvertices=size(Coor,1);
620-md.mesh.numberofelements=size(Tria,1);
621-numberofelements2=md.mesh.numberofelements;
622-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
623-
624-%display number of elements
625-fprintf('\n%s %i',' inital number of elements:',numberofelements1);
626-fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
627-
628-%clean up:
629-system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');
630Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/yams.m
631===================================================================
632--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/yams.m (revision 13011)
633+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/yams.m (revision 13012)
634@@ -1,171 +0,0 @@
635-function md=yams(md,varargin);
636-%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
637-%
638-% Usage:
639-% md=yams(md,varargin);
640-% where varargin is a lit of paired arguments.
641-% arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
642-% arguments can be: 'velocities': matlab file containing the velocities [m/yr]
643-% optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
644-% this option is used to minimize the metric on water (no refinement)
645-% optional arguments: 'resolution': initial mesh resolution [m]
646-% optional arguments: 'nsteps': number of steps of mesh adaptation
647-% optional arguments: 'epsilon': average interpolation error wished [m/yr]
648-% optional arguments: 'hmin': minimum edge length
649-% optional arguments: 'hmanx': maximum edge
650-% optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
651-%
652-%
653-% Examples:
654-% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
655-% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
656-% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
657-
658-%recover options
659-options=pairoptions(varargin{:});
660-options=deleteduplicates(options,1);
661-
662-%recover some fields
663-disp('MeshYams Options:')
664-domainoutline=getfieldvalue(options,'domainoutline');
665-disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline));
666-riftoutline=getfieldvalue(options,'riftoutline','N/A');
667-disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline));
668-groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
669-disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain));
670-velocities=getfieldvalue(options,'velocities');
671-disp(sprintf(' %-15s: ''%s''','Velocities',velocities));
672-resolution=getfieldvalue(options,'resolution',5000);
673-disp(sprintf(' %-15s: %f','Resolution',resolution));
674-nsteps=getfieldvalue(options,'nsteps',6);
675-disp(sprintf(' %-15s: %i','nsteps',nsteps));
676-gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
677-disp(sprintf(' %-15s: %g','gradation',gradation(1)));
678-epsilon=getfieldvalue(options,'epsilon',3);
679-disp(sprintf(' %-15s: %f','epsilon',epsilon));
680-hmin=getfieldvalue(options,'hmin',500);
681-disp(sprintf(' %-15s: %f','hmin',hmin));
682-hmax=getfieldvalue(options,'hmax',150*10^3);
683-disp(sprintf(' %-15s: %f\n','hmax',hmax));
684-
685-%mesh with initial resolution
686-disp('Initial mesh generation...');
687-if strcmpi(riftoutline,'N/A');
688- md=setmesh(md,domainoutline,resolution);
689-else
690- md=setmesh(md,domainoutline,riftoutline,resolution);
691- md=meshprocessrifts(md,domainoutline);
692-end
693-disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
694-
695-%load velocities
696-disp('loading velocities...');
697-Names=VelFindVarNames(velocities);
698-Vel=load(velocities);
699-
700-%start mesh adaptation
701-for i=1:nsteps,
702- disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
703-
704- %interpolate velocities onto mesh
705- disp(' interpolating velocities...');
706- if strcmpi(Names.interp,'node'),
707- vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
708- vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
709- else
710- vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
711- vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
712- end
713- field=sqrt(vx_obs.^2+vy_obs.^2);
714-
715- %set mask.vertexonwater field
716- if ~strcmp(groundeddomain,'N/A'),
717- nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
718- md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
719- md.mask.vertexonwater(find(nodeground))=0;
720- else
721- md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
722- end
723-
724- %adapt according to velocities
725- disp(' adapting...');
726- md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
727-
728- %if we have rifts, we just messed them up, we need to recreate the segments that constitute those
729- %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
730- if md.rifts.numrifts,
731- md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
732- md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
733- md.mesh.segments=findsegments(md);
734- md=yamsrecreateriftsegments(md);
735- end
736-
737-end
738-
739-disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
740-
741-%Now, build the connectivity tables for this mesh.
742-md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
743-md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
744-
745-%recreate segments
746-md.mesh.segments=findsegments(md);
747-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
748-
749-%Fill in rest of fields:
750-md.mesh.z=zeros(md.mesh.numberofvertices,1);
751-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
752-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
753-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
754-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
755-if ~strcmp(groundeddomain,'N/A'),
756- nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
757- md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
758- md.mask.vertexonwater(find(nodeground))=0;
759-else
760- md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
761-end
762-if strcmpi(Names.interp,'node'),
763- md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
764- md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
765-else
766- md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
767- md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
768-end
769-md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
770-
771-%deal with rifts
772-if md.rifts.numrifts,
773- %first, recreate rift segments
774- md=meshyamsrecreateriftsegments(md);
775-
776- %using the segments, recreate the penaltypairs
777- for j=1:md.rifts.numrifts,
778- rift=md.rifts.riftstruct(j);
779-
780- %build normals and lengths of segments:
781- 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 );
782- 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)))));
783- 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)))));
784-
785- %ok, build penaltypairs:
786- numpenaltypairs=length(rift.segments)/2-1;
787- rift.penaltypairs=zeros(numpenaltypairs,7);
788-
789- for i=1:numpenaltypairs,
790- rift.penaltypairs(i,1)=rift.segments(i,2);
791- rift.penaltypairs(i,2)=rift.segments(end-i,2);
792- rift.penaltypairs(i,3)=rift.segments(i,3);
793- rift.penaltypairs(i,4)=rift.segments(end-i,3);
794- rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
795- rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
796- rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
797- end
798- %renormalize norms:
799- norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
800- rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
801- rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
802-
803- md.rifts.riftstruct(j)=rift;
804- end
805-end
806Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifttipsonmesh.m
807===================================================================
808--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifttipsonmesh.m (revision 13011)
809+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifttipsonmesh.m (revision 13012)
810@@ -1,26 +0,0 @@
811-function tips=rifttipsonmesh(md,riftoutline)
812-%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of
813-% rifts.
814-
815-%read rifts from outline file
816-rifts=expread(riftoutline,1);
817-
818-tips=[];
819-
820-for i=1:length(rifts),
821- rift=rifts(i);
822-
823- x_tip=rift.x(1);
824- y_tip=rift.y(1);
825-
826- index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
827- tips(end+1)=index;
828-
829- x_tip=rift.x(end);
830- y_tip=rift.y(end);
831-
832- index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
833- tips(end+1)=index;
834-
835-end
836-
837Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/NodeInElement.m
838===================================================================
839--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/NodeInElement.m (revision 13011)
840+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/NodeInElement.m (revision 13012)
841@@ -1,7 +1,8 @@
842 function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity);
843-%NODEINELEMENT: find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
844+% NODEINELEMENT - find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
845 %
846-% Usage: node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
847+% Usage:
848+% node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
849 %
850 % See also Nodeconnectivity
851 %
852Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m
853===================================================================
854--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m (revision 0)
855+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m (revision 13012)
856@@ -0,0 +1,26 @@
857+function tips=rifttipsonmesh(md,riftoutline)
858+%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of
859+% rifts.
860+
861+%read rifts from outline file
862+rifts=expread(riftoutline,1);
863+
864+tips=[];
865+
866+for i=1:length(rifts),
867+ rift=rifts(i);
868+
869+ x_tip=rift.x(1);
870+ y_tip=rift.y(1);
871+
872+ index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
873+ tips(end+1)=index;
874+
875+ x_tip=rift.x(end);
876+ y_tip=rift.y(end);
877+
878+ index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
879+ tips(end+1)=index;
880+
881+end
882+
883Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m
884===================================================================
885--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13011)
886+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13012)
887@@ -92,3 +92,19 @@
888 md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
889 md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
890 md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
891+end
892+
893+function flag=isconnected(elements,A,B)% {{{
894+ %ISCONNECTED: are two nodes connected by a triangulation?
895+ %
896+ % Usage: flag=isconnected(elements,A,B)
897+ %
898+ %
899+
900+ elements=ElementsFromEdge(elements,A,B);
901+ if isempty(elements),
902+ flag=0;
903+ else
904+ flag=1;
905+ end
906+end % }}}
907Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/FixMesh.m
908===================================================================
909--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/FixMesh.m (revision 13011)
910+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/mesh/FixMesh.m (revision 13012)
911@@ -1,11 +1,11 @@
912 function [index2 x2 y2 value2]=FixMesh(index,x,y,value)
913-%FixMesh fix mesh with broken triangles, orphan vertices, etc ...
914+% FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ...
915 %
916-% Usage:
917-% [index2 x2 y2 value2]=FixMesh(index,x,y,value)
918-% where index,x,y is a delaunay triangulation,
919-% value is a field on the input triangulation, with values at the vertices
920-% index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
921+% Usage:
922+% [index2 x2 y2 value2]=FixMesh(index,x,y,value)
923+% where index,x,y is a delaunay triangulation,
924+% value is a field on the input triangulation, with values at the vertices
925+% index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
926 %
927 %
928
929Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/modellist.m
930===================================================================
931--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/modellist.m (revision 13011)
932+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/modellist.m (revision 13012)
933@@ -291,3 +291,97 @@
934 end % }}}
935 end
936 end
937+
938+function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)% {{{
939+%BUILDMULTIPLEQUEUEINGSCRIPT -
940+%
941+% Usage:
942+% BuildMultipleQueueingScript(executionpath,codepath)
943+
944+disp('building queueing script');
945+
946+%First try and figure out if there is a special script for this particular cluster
947+function_name=['BuildMultipleQueueingScript' cluster];
948+
949+%some specific treatment of identical cluster, gemini, castor and pollux
950+if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
951+ function_name='BuildMultipleQueueingScriptgemini';
952+end
953+
954+if exist(function_name,'file'),
955+ %Call this function:
956+ eval([function_name '(name,executionpath,codepath);']);
957+else
958+ %Call the generic BuildQueueingScript:
959+ BuildMultipleQueueingScriptGeneric(name,executionpath,codepath);
960+end
961+end % }}}
962+function BuildQueueingScriptgemini(name,executionpath,codepath)% {{{
963+%BUILDQUEUEINGSCRIPTGEMINI - ...
964+%
965+% Usage:
966+% BuildQueueingScriptgemini(md,executionpath,codepath)
967+
968+scriptname=[name '.queue'];
969+
970+fid=fopen(scriptname,'w');
971+if fid==-1,
972+ error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']);
973+end
974+
975+fprintf(fid,'#!/bin/sh\n');
976+fprintf(fid,'cd %s\n',executionpath);
977+fprintf(fid,'mkdir %s\n',name);
978+fprintf(fid,'cd %s\n',name);
979+fprintf(fid,'mv ../ModelList.tar.gz ./\n');
980+fprintf(fid,'tar -zxvf ModelList.tar.gz\n');
981+fprintf(fid,'foreach i (%s-*vs*.queue)\n',name);
982+fprintf(fid,'qsub $i\n');
983+fprintf(fid,'end\n');
984+fclose(fid);
985+end% }}}
986+function LaunchMultipleQueueJob(cluster,name,executionpath)% {{{
987+%LAUNCHMULTIPLEQUEUEJOB - ...
988+%
989+% Usage:
990+% LaunchMultipleQueueJob(executionpath)
991+
992+%First try and figure out if there is a special script for thie particular cluster
993+function_name=['LaunchMultipleQueueJob' cluster];
994+
995+%some specific treatment of identical cluster, gemini, castor and pollux
996+if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
997+ function_name='LaunchMultipleQueueJobgemini';
998+end
999+
1000+if exist(function_name,'file'),
1001+ %Call this function:
1002+ eval([function_name '(cluster,name,executionpath);']);
1003+else
1004+ %Call the generic LaunchMultipleQueueJob:
1005+ LaunchMultipleQueueJobGeneric(cluster,name,executionpath);
1006+end
1007+end% }}}
1008+function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)% {{{
1009+%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster
1010+%
1011+% Usage:
1012+% LaunchMultipleQueueJobgemini(cluster,name,executionpath)
1013+
1014+
1015+%first, check we have the binary file and the queueing script
1016+if ~exist([ name '.queue'],'file'),
1017+ error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward');
1018+end
1019+
1020+if ~exist('ModelList.tar.gz','file'),
1021+ error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward');
1022+end
1023+
1024+%upload both files to cluster
1025+disp('uploading input file, queueing script and variables script');
1026+eval(['!scp ModelList.tar.gz ' name '.queue ' cluster ':' executionpath]);
1027+
1028+disp('launching solution sequence on remote cluster');
1029+issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']);
1030+end% }}}
1031Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/YamsCall.m
1032===================================================================
1033--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/YamsCall.m (revision 0)
1034+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/YamsCall.m (revision 13012)
1035@@ -0,0 +1,104 @@
1036+function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
1037+%YAMSCALL - call yams
1038+%
1039+% build a metric using the Hessian of the given field
1040+% call Yams and the output mesh is plugged onto the model
1041+% -hmin = minimum edge length (m)
1042+% -hmax = maximum edge length (m)
1043+% -gradation = maximum edge length gradation between 2 elements
1044+% -epsilon = average error on each element (m/yr)
1045+%
1046+% Usage:
1047+% md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
1048+%
1049+% Example:
1050+% md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
1051+
1052+%2d geometric parameter (do not change)
1053+scale=2/9;
1054+
1055+%Compute Hessian
1056+t1=clock; fprintf('%s',' computing Hessian...');
1057+hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
1058+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1059+
1060+%Compute metric
1061+t1=clock; fprintf('%s',' computing metric...');
1062+if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
1063+ pos=find(md.mask.vertexonwater);
1064+else
1065+ pos=[];
1066+end
1067+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
1068+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1069+
1070+%write files
1071+t1=clock; fprintf('%s',' writing initial mesh files...');
1072+save -ascii carre0.met metric
1073+
1074+fid=fopen('carre0.mesh','w');
1075+
1076+%initialiation
1077+fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
1078+
1079+%dimension
1080+fprintf(fid,'\n%s\n%i\n','Dimension',2);
1081+
1082+%Vertices
1083+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
1084+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
1085+
1086+%Triangles
1087+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
1088+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
1089+numberofelements1=md.mesh.numberofelements;
1090+
1091+%Deal with rifts
1092+if ~isnan(md.rifts.riftstruct),
1093+
1094+ %we have the list of triangles that make up the rift. keep those triangles around during refinement.
1095+ triangles=[];
1096+ for i=1:size(md.rifts.riftstruct,1),
1097+ triangles=[triangles md.rifts(i).segments(:,3)'];
1098+ end
1099+
1100+ fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
1101+ fprintf(fid,'%i\n',triangles);
1102+end
1103+
1104+%close
1105+fclose(fid);
1106+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1107+
1108+%call yams
1109+fprintf('%s\n',' call Yams...');
1110+if ispc
1111+ %windows
1112+ system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
1113+elseif ismac
1114+ %Macosx
1115+ system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
1116+else
1117+ %Linux
1118+ system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation) ' carre0 carre1']);
1119+end
1120+
1121+%plug new mesh
1122+t1=clock; fprintf('\n%s',' reading final mesh files...');
1123+Tria=load('carre1.tria');
1124+Coor=load('carre1.coor');
1125+md.mesh.x=Coor(:,1);
1126+md.mesh.y=Coor(:,2);
1127+md.mesh.z=zeros(size(Coor,1),1);
1128+md.mesh.elements=Tria;
1129+md.mesh.numberofvertices=size(Coor,1);
1130+md.mesh.numberofelements=size(Tria,1);
1131+numberofelements2=md.mesh.numberofelements;
1132+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1133+
1134+%display number of elements
1135+fprintf('\n%s %i',' inital number of elements:',numberofelements1);
1136+fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
1137+
1138+%clean up:
1139+system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');
1140Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/meshread.m
1141===================================================================
1142--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/meshread.m (revision 0)
1143+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/meshread.m (revision 13012)
1144@@ -0,0 +1,41 @@
1145+function Struct=meshread(filename);
1146+
1147+%some checks
1148+if ~exist(filename),
1149+ error(['meshread error message: file ' filename ' not found!']);
1150+end
1151+
1152+fid=fopen(filename,'r');
1153+
1154+while (~feof(fid)),
1155+
1156+ A=fscanf(fid,'%s',1);
1157+
1158+ if strcmp(A,'MeshVersionFormatted');
1159+ Struct.Version=fscanf(fid,'%s',1);
1160+
1161+ elseif strcmp(A,'Dimension'),
1162+ Struct.Dimension=fscanf(fid,'%i',1);
1163+
1164+ elseif strcmp(A,'Vertices'),
1165+ Struct.nods=fscanf(fid,'%i',1);
1166+ A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
1167+ Struct.x=A(1,:)';
1168+ Struct.y=A(2,:)';
1169+
1170+ elseif strcmp(A,'Triangles'),
1171+ Struct.nels=fscanf(fid,'%i',1);
1172+ A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
1173+ Struct.index=A(1:3,:)';
1174+
1175+ elseif strcmp(A,'Quadrilaterals'),
1176+ Struct.nels=fscanf(fid,'%i',1);
1177+ A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
1178+ Struct.index=A(1:4,:)';
1179+ else
1180+ %do nothing
1181+
1182+ end
1183+end
1184+
1185+fclose(fid);
1186Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/yams.m
1187===================================================================
1188--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/yams.m (revision 0)
1189+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/yams.m (revision 13012)
1190@@ -0,0 +1,171 @@
1191+function md=yams(md,varargin);
1192+%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
1193+%
1194+% Usage:
1195+% md=yams(md,varargin);
1196+% where varargin is a lit of paired arguments.
1197+% arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
1198+% arguments can be: 'velocities': matlab file containing the velocities [m/yr]
1199+% optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
1200+% this option is used to minimize the metric on water (no refinement)
1201+% optional arguments: 'resolution': initial mesh resolution [m]
1202+% optional arguments: 'nsteps': number of steps of mesh adaptation
1203+% optional arguments: 'epsilon': average interpolation error wished [m/yr]
1204+% optional arguments: 'hmin': minimum edge length
1205+% optional arguments: 'hmanx': maximum edge
1206+% optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
1207+%
1208+%
1209+% Examples:
1210+% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
1211+% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
1212+% md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
1213+
1214+%recover options
1215+options=pairoptions(varargin{:});
1216+options=deleteduplicates(options,1);
1217+
1218+%recover some fields
1219+disp('MeshYams Options:')
1220+domainoutline=getfieldvalue(options,'domainoutline');
1221+disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline));
1222+riftoutline=getfieldvalue(options,'riftoutline','N/A');
1223+disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline));
1224+groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
1225+disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain));
1226+velocities=getfieldvalue(options,'velocities');
1227+disp(sprintf(' %-15s: ''%s''','Velocities',velocities));
1228+resolution=getfieldvalue(options,'resolution',5000);
1229+disp(sprintf(' %-15s: %f','Resolution',resolution));
1230+nsteps=getfieldvalue(options,'nsteps',6);
1231+disp(sprintf(' %-15s: %i','nsteps',nsteps));
1232+gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
1233+disp(sprintf(' %-15s: %g','gradation',gradation(1)));
1234+epsilon=getfieldvalue(options,'epsilon',3);
1235+disp(sprintf(' %-15s: %f','epsilon',epsilon));
1236+hmin=getfieldvalue(options,'hmin',500);
1237+disp(sprintf(' %-15s: %f','hmin',hmin));
1238+hmax=getfieldvalue(options,'hmax',150*10^3);
1239+disp(sprintf(' %-15s: %f\n','hmax',hmax));
1240+
1241+%mesh with initial resolution
1242+disp('Initial mesh generation...');
1243+if strcmpi(riftoutline,'N/A');
1244+ md=setmesh(md,domainoutline,resolution);
1245+else
1246+ md=setmesh(md,domainoutline,riftoutline,resolution);
1247+ md=meshprocessrifts(md,domainoutline);
1248+end
1249+disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
1250+
1251+%load velocities
1252+disp('loading velocities...');
1253+Names=VelFindVarNames(velocities);
1254+Vel=load(velocities);
1255+
1256+%start mesh adaptation
1257+for i=1:nsteps,
1258+ disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
1259+
1260+ %interpolate velocities onto mesh
1261+ disp(' interpolating velocities...');
1262+ if strcmpi(Names.interp,'node'),
1263+ vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
1264+ vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
1265+ else
1266+ vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
1267+ vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
1268+ end
1269+ field=sqrt(vx_obs.^2+vy_obs.^2);
1270+
1271+ %set mask.vertexonwater field
1272+ if ~strcmp(groundeddomain,'N/A'),
1273+ nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
1274+ md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
1275+ md.mask.vertexonwater(find(nodeground))=0;
1276+ else
1277+ md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
1278+ end
1279+
1280+ %adapt according to velocities
1281+ disp(' adapting...');
1282+ md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
1283+
1284+ %if we have rifts, we just messed them up, we need to recreate the segments that constitute those
1285+ %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
1286+ if md.rifts.numrifts,
1287+ md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
1288+ md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
1289+ md.mesh.segments=findsegments(md);
1290+ md=yamsrecreateriftsegments(md);
1291+ end
1292+
1293+end
1294+
1295+disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
1296+
1297+%Now, build the connectivity tables for this mesh.
1298+md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
1299+md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
1300+
1301+%recreate segments
1302+md.mesh.segments=findsegments(md);
1303+md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
1304+
1305+%Fill in rest of fields:
1306+md.mesh.z=zeros(md.mesh.numberofvertices,1);
1307+md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
1308+md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
1309+md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
1310+md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
1311+if ~strcmp(groundeddomain,'N/A'),
1312+ nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
1313+ md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
1314+ md.mask.vertexonwater(find(nodeground))=0;
1315+else
1316+ md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
1317+end
1318+if strcmpi(Names.interp,'node'),
1319+ md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
1320+ md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
1321+else
1322+ md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
1323+ md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
1324+end
1325+md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
1326+
1327+%deal with rifts
1328+if md.rifts.numrifts,
1329+ %first, recreate rift segments
1330+ md=meshyamsrecreateriftsegments(md);
1331+
1332+ %using the segments, recreate the penaltypairs
1333+ for j=1:md.rifts.numrifts,
1334+ rift=md.rifts.riftstruct(j);
1335+
1336+ %build normals and lengths of segments:
1337+ 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 );
1338+ 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)))));
1339+ 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)))));
1340+
1341+ %ok, build penaltypairs:
1342+ numpenaltypairs=length(rift.segments)/2-1;
1343+ rift.penaltypairs=zeros(numpenaltypairs,7);
1344+
1345+ for i=1:numpenaltypairs,
1346+ rift.penaltypairs(i,1)=rift.segments(i,2);
1347+ rift.penaltypairs(i,2)=rift.segments(end-i,2);
1348+ rift.penaltypairs(i,3)=rift.segments(i,3);
1349+ rift.penaltypairs(i,4)=rift.segments(end-i,3);
1350+ rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
1351+ rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
1352+ rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
1353+ end
1354+ %renormalize norms:
1355+ norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
1356+ rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
1357+ rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
1358+
1359+ md.rifts.riftstruct(j)=rift;
1360+ end
1361+end
1362Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m
1363===================================================================
1364--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m (revision 0)
1365+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m (revision 13012)
1366@@ -0,0 +1,65 @@
1367+function md=BamgCallFromMetric(md,metric,gradation),
1368+%BAMGCALL - call bam
1369+%
1370+% call Bamg and the output mesh is plugged onto the model
1371+% -gradation = maximum edge length gradation between 2 elements
1372+%
1373+% Usage:
1374+% md=BamgCallFromMetric(md,metric,gradation);
1375+%
1376+% Example:
1377+% md=BamgCall(md,metric,1500,10^8,1.3,0.9);
1378+
1379+%2d geometric parameter (do not change)
1380+scale=2/9;
1381+
1382+%write files
1383+t1=clock; fprintf('%s',' writing initial mesh files...');
1384+fid=fopen('carre0.met','w');
1385+fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
1386+fprintf(fid,'%i %i %i\n',metric');
1387+fclose(fid);
1388+
1389+fid=fopen('carre0.mesh','w');
1390+
1391+%initialiation
1392+fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
1393+
1394+%dimension
1395+fprintf(fid,'\n%s\n%i\n','Dimension',2);
1396+
1397+%Vertices
1398+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
1399+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
1400+
1401+%Triangles
1402+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
1403+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
1404+numberofelements1=md.mesh.numberofelements;
1405+
1406+%close
1407+fclose(fid);
1408+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1409+
1410+%call bamg
1411+fprintf('%s\n',' call Bamg...');
1412+system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
1413+
1414+%plug new mesh
1415+t1=clock; fprintf('\n%s',' reading final mesh files...');
1416+A=meshread('carre1.mesh');
1417+md.mesh.x=A.x;
1418+md.mesh.y=A.y;
1419+md.z=zeros(A.nods,1);
1420+md.mesh.elements=A.index;
1421+md.mesh.numberofvertices=A.nods;
1422+md.mesh.numberofelements=A.nels;
1423+numberofelements2=md.mesh.numberofelements;
1424+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1425+
1426+%display number of elements
1427+fprintf('\n%s %i',' inital number of elements:',numberofelements1);
1428+fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
1429+
1430+%clean up:
1431+system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
1432Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCall.m
1433===================================================================
1434--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCall.m (revision 0)
1435+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/contrib/bamg/BamgCall.m (revision 13012)
1436@@ -0,0 +1,84 @@
1437+function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
1438+%BAMGCALL - call bam
1439+%
1440+% build a metric using the Hessian of the given field
1441+% call Bamg and the output mesh is plugged onto the model
1442+% -hmin = minimum edge length (m)
1443+% -hmax = maximum edge length (m)
1444+% -gradation = maximum edge length gradation between 2 elements
1445+% -epsilon = average error on each element (m/yr)
1446+%
1447+% Usage:
1448+% md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
1449+%
1450+% Example:
1451+% md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
1452+
1453+%2d geometric parameter (do not change)
1454+scale=2/9;
1455+
1456+%Compute Hessian
1457+t1=clock; fprintf('%s',' computing Hessian...');
1458+hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
1459+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1460+
1461+%Compute metric
1462+t1=clock; fprintf('%s',' computing metric...');
1463+if length(md.nodeonwater)==md.mesh.numberofvertices,
1464+ pos=find(md.nodeonwater);
1465+else
1466+ pos=[];
1467+end
1468+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
1469+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1470+
1471+%write files
1472+t1=clock; fprintf('%s',' writing initial mesh files...');
1473+fid=fopen('carre0.met','w');
1474+fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
1475+fprintf(fid,'%i %i %i\n',metric');
1476+fclose(fid);
1477+
1478+fid=fopen('carre0.mesh','w');
1479+
1480+%initialiation
1481+fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
1482+
1483+%dimension
1484+fprintf(fid,'\n%s\n%i\n','Dimension',2);
1485+
1486+%Vertices
1487+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
1488+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
1489+
1490+%Triangles
1491+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
1492+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
1493+numberofelements1=md.mesh.numberofelements;
1494+
1495+%close
1496+fclose(fid);
1497+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1498+
1499+%call bamg
1500+fprintf('%s\n',' call Bamg...');
1501+system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
1502+
1503+%plug new mesh
1504+t1=clock; fprintf('\n%s',' reading final mesh files...');
1505+A=meshread('carre1.mesh');
1506+md.mesh.x=A.x;
1507+md.mesh.y=A.y;
1508+md.z=zeros(A.nods,1);
1509+md.mesh.elements=A.index;
1510+md.mesh.numberofvertices=A.nods;
1511+md.mesh.numberofelements=A.nels;
1512+numberofelements2=md.mesh.numberofelements;
1513+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
1514+
1515+%display number of elements
1516+fprintf('\n%s %i',' inital number of elements:',numberofelements1);
1517+fprintf('\n%s %i\n\n',' new number of elements:',numberofelements2);
1518+
1519+%clean up:
1520+system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
1521Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/create_region.m
1522===================================================================
1523--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/create_region.m (revision 13011)
1524+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/create_region.m (revision 13012)
1525@@ -1,15 +0,0 @@
1526-function create_region(name)
1527-%CREATE_REGION - create region ????
1528-%
1529-% very temporary function.
1530-%
1531-% Usage:
1532-% create_region(name)
1533-
1534-eval(['mkdir ' name]);
1535-eval(['cd ' name ]);
1536-!mkdir Delivery Exp_Par Results
1537-cd Exp_Par
1538-!cp ../../RonneShelf/Exp_Par/* ./
1539-!rm -rf Hole*
1540-eval(['!mv Ronne.par ' name '.par']);
1541Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/findarg.m
1542===================================================================
1543--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/findarg.m (revision 13011)
1544+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/findarg.m (revision 13012)
1545@@ -1,51 +0,0 @@
1546-function vals=findarg(arglist,field)
1547-%FINDARG - find argument associated to a field in a list
1548-%
1549-% This function parses through an argument list (typically varargin in a routine)
1550-% looking for a character array equal to field. Once this is found, we return the
1551-% next value in the varargin (if possible).
1552-% Because field might appear several times in the argument list, we return a structure
1553-% holding all these values.
1554-% Note that all comparisons to field value are case independent.
1555-%
1556-% Usage:
1557-% vals=findarg(arglist,field)
1558-%
1559-% Example:
1560-% routine foobar calls vals=findarg('Data',varargin)
1561-% with varargin='Data',1,'Data','foo','Plot','velocity','Arrow',4
1562-% findarg would return the following structure: vals(1).value=1, vals(2).value='foo';
1563-
1564-%some argument checking:
1565-if ((nargin==0) | (nargout==0)),
1566- help findarg;
1567- error('findarg error message');
1568-end
1569-
1570-if ~ischar(field),
1571- error('findarg error message: field should be a string');
1572-end
1573-
1574-if ~iscell(arglist),
1575- error('findarg error message: argument list should be a cell array.');
1576-end
1577-
1578-%Recover data to plot
1579-founddata=0;
1580-
1581-for i=1:(length(arglist)-1), %data in arglist comes in pairs, hence the -1.
1582- if ischar(arglist{i}),
1583- if (strcmpi(arglist{i},field)),
1584- founddata=founddata+1;
1585- if founddata==1,
1586- vals.value=arglist{i+1};
1587- else
1588- vals(end+1).value=arglist{i+1};
1589- end
1590- end
1591- end
1592-end
1593-
1594-if founddata==0,
1595- vals=[];
1596-end
1597Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf.m
1598===================================================================
1599--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf.m (revision 13011)
1600+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf.m (revision 13012)
1601@@ -1,136 +0,0 @@
1602-function S = netcdf(File,varargin)
1603-% Function to read NetCDF files
1604-% S = netcdf(File)
1605-% Input Arguments
1606-% File = NetCDF file to read
1607-% Optional Input Arguments:
1608-% 'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
1609-% 'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
1610-% Output Arguments:
1611-% S = Structure of NetCDF data organised as per NetCDF definition
1612-% Notes:
1613-% Only version 1, classic 32bit, NetCDF files are supported. By default
1614-% data are extracted into the S.VarArray().Data field for all variables.
1615-% To read the header only call S = netcdf(File,'Var',[]);
1616-%
1617-% SEE ALSO
1618-% ---------------------------------------------------------------------------
1619-S = [];
1620-
1621-try
1622- if exist(File,'file') fp = fopen(File,'r','b');
1623- else fp = []; error('File not found'); end
1624- if fp == -1 error('Unable to open file'); end
1625-
1626-% Read header
1627- Magic = fread(fp,4,'uint8=>char');
1628- if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
1629- if uint8(Magic(4))~=1 error('Version not supported'); end
1630- S.NumRecs = fread(fp,1,'uint32=>uint32');
1631- S.DimArray = DimArray(fp);
1632- S.AttArray = AttArray(fp);
1633- S.VarArray = VarArray(fp);
1634-
1635-% Setup indexing to arrays and records
1636- Var = ones(1,length(S.VarArray));
1637- Rec = ones(1,S.NumRecs);
1638- for i = 1:2:length(varargin)
1639- if strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
1640- elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
1641- else error('Optional input argument not recognised'); end
1642- end
1643- if sum(Var)==0 fclose(fp); return; end
1644-
1645-% Read non-record variables
1646- Dim = double(cat(2,S.DimArray.Dim));
1647- ID = double(cat(2,S.VarArray.Type));
1648-
1649- for i = 1:length(S.VarArray)
1650- D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
1651- if isempty(RecID{i})
1652- if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
1653- if Var(i)
1654- S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
1655- fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
1656- else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
1657- else S.VarArray(i).Data = []; end
1658- end
1659-
1660-% Read record variables
1661- for k = 1:S.NumRecs
1662- for i = 1:length(S.VarArray)
1663- if ~isempty(RecID{i})
1664- D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
1665- if length(D)==1 D=[D,1]; end
1666- if Var(i) & Rec(k)
1667- S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
1668- ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
1669- if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
1670- else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
1671- end
1672- end
1673- end
1674-
1675- fclose(fp);
1676-catch
1677- Err = lasterror; fprintf('%s\n',Err.message);
1678- if ~isempty(fp) && fp ~= -1 fclose(fp); end
1679-end
1680-
1681-% ---------------------------------------------------------------------------------------
1682-% Utility functions
1683-
1684-function S = Size(ID)
1685-% Size of NetCDF data type, ID, in bytes
1686- S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
1687-
1688-function T = Type(ID)
1689-% Matlab string for CDF data type, ID
1690- T = subsref({'int8','char','int16','int32','single','double'},...
1691- struct('type','{}','subs',{{ID}}));
1692-
1693-function N = Pad(Num,ID)
1694-% Number of elements to read after padding to 4 bytes for type ID
1695- N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
1696-
1697-function S = String(fp)
1698-% Read a CDF string; Size,[String,[Padding]]
1699- S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
1700-
1701-function A = ReOrder(A,S)
1702-% Rearrange CDF array A to size S with matlab ordering
1703- A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
1704-
1705-function S = DimArray(fp)
1706-% Read DimArray into structure
1707- if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
1708- for i = 1:fread(fp,1,'uint32=>uint32')
1709- S(i).Str = String(fp);
1710- S(i).Dim = fread(fp,1,'uint32=>uint32');
1711- end
1712- else fread(fp,1,'uint32=>uint32'); S = []; end
1713-
1714-function S = AttArray(fp)
1715-% Read AttArray into structure
1716- if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
1717- for i = 1:fread(fp,1,'uint32=>uint32')
1718- S(i).Str = String(fp);
1719- ID = fread(fp,1,'uint32=>uint32');
1720- Num = fread(fp,1,'uint32=>uint32');
1721- S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
1722- end
1723- else fread(fp,1,'uint32=>uint32'); S = []; end
1724-
1725-function S = VarArray(fp)
1726-% Read VarArray into structure
1727- if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
1728- for i = 1:fread(fp,1,'uint32=>uint32')
1729- S(i).Str = String(fp);
1730- Num = double(fread(fp,1,'uint32=>uint32'));
1731- S(i).DimID = double(fread(fp,Num,'uint32=>uint32'));
1732- S(i).AttArray = AttArray(fp);
1733- S(i).Type = fread(fp,1,'uint32=>uint32');
1734- S(i).VSize = fread(fp,1,'uint32=>uint32');
1735- S(i).Begin = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
1736- end
1737- else fread(fp,1,'uint32=>uint32'); S = []; end
1738Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf2struct.m
1739===================================================================
1740--- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf2struct.m (revision 13011)
1741+++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/miscellaneous/netcdf2struct.m (revision 13012)
1742@@ -5,7 +5,7 @@
1743 % S=netcdf2struct(File);
1744
1745 %Read netcdf file
1746-data=netcdf(File);
1747+data=readnetcdf(File);
1748
1749 %initialize output
1750 S=struct();
1751@@ -25,3 +25,142 @@
1752 fieldvalue=double(variables(i).Val);
1753 S.(fieldname)=fieldvalue;
1754 end
1755+end
1756+
1757+function S = readnetcdf(File,varargin)
1758+% Function to read NetCDF files
1759+% S = netcdf(File)
1760+% Input Arguments
1761+% File = NetCDF file to read
1762+% Optional Input Arguments:
1763+% 'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
1764+% 'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
1765+% Output Arguments:
1766+% S = Structure of NetCDF data organised as per NetCDF definition
1767+% Notes:
1768+% Only version 1, classic 32bit, NetCDF files are supported. By default
1769+% data are extracted into the S.VarArray().Data field for all variables.
1770+% To read the header only call S = netcdf(File,'Var',[]);
1771+%
1772+% SEE ALSO
1773+% ---------------------------------------------------------------------------
1774+S = [];
1775+
1776+try
1777+ if exist(File,'file') fp = fopen(File,'r','b');
1778+ else fp = []; error('File not found'); end
1779+ if fp == -1 error('Unable to open file'); end
1780+
1781+% Read header
1782+ Magic = fread(fp,4,'uint8=>char');
1783+ if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
1784+ if uint8(Magic(4))~=1 error('Version not supported'); end
1785+ S.NumRecs = fread(fp,1,'uint32=>uint32');
1786+ S.DimArray = DimArray(fp);
1787+ S.AttArray = AttArray(fp);
1788+ S.VarArray = VarArray(fp);
1789+
1790+% Setup indexing to arrays and records
1791+ Var = ones(1,length(S.VarArray));
1792+ Rec = ones(1,S.NumRecs);
1793+ for i = 1:2:length(varargin)
1794+ if strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
1795+ elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
1796+ else error('Optional input argument not recognised'); end
1797+ end
1798+ if sum(Var)==0 fclose(fp); return; end
1799+
1800+% Read non-record variables
1801+ Dim = double(cat(2,S.DimArray.Dim));
1802+ ID = double(cat(2,S.VarArray.Type));
1803+
1804+ for i = 1:length(S.VarArray)
1805+ D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
1806+ if isempty(RecID{i})
1807+ if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
1808+ if Var(i)
1809+ S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
1810+ fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
1811+ else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
1812+ else S.VarArray(i).Data = []; end
1813+ end
1814+
1815+% Read record variables
1816+ for k = 1:S.NumRecs
1817+ for i = 1:length(S.VarArray)
1818+ if ~isempty(RecID{i})
1819+ D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
1820+ if length(D)==1 D=[D,1]; end
1821+ if Var(i) & Rec(k)
1822+ S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
1823+ ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
1824+ if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
1825+ else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
1826+ end
1827+ end
1828+ end
1829+
1830+ fclose(fp);
1831+catch
1832+ Err = lasterror; fprintf('%s\n',Err.message);
1833+ if ~isempty(fp) && fp ~= -1 fclose(fp); end
1834+end
1835+
1836+% ---------------------------------------------------------------------------------------
1837+% Utility functions
1838+
1839+function S = Size(ID)
1840+% Size of NetCDF data type, ID, in bytes
1841+ S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
1842+
1843+function T = Type(ID)
1844+% Matlab string for CDF data type, ID
1845+ T = subsref({'int8','char','int16','int32','single','double'},...
1846+ struct('type','{}','subs',{{ID}}));
1847+
1848+function N = Pad(Num,ID)
1849+% Number of elements to read after padding to 4 bytes for type ID
1850+ N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
1851+
1852+function S = String(fp)
1853+% Read a CDF string; Size,[String,[Padding]]
1854+ S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
1855+
1856+function A = ReOrder(A,S)
1857+% Rearrange CDF array A to size S with matlab ordering
1858+ A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
1859+
1860+function S = DimArray(fp)
1861+% Read DimArray into structure
1862+ if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
1863+ for i = 1:fread(fp,1,'uint32=>uint32')
1864+ S(i).Str = String(fp);
1865+ S(i).Dim = fread(fp,1,'uint32=>uint32');
1866+ end
1867+ else fread(fp,1,'uint32=>uint32'); S = []; end
1868+
1869+function S = AttArray(fp)
1870+% Read AttArray into structure
1871+ if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
1872+ for i = 1:fread(fp,1,'uint32=>uint32')
1873+ S(i).Str = String(fp);
1874+ ID = fread(fp,1,'uint32=>uint32');
1875+ Num = fread(fp,1,'uint32=>uint32');
1876+ S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
1877+ end
1878+ else fread(fp,1,'uint32=>uint32'); S = []; end
1879+
1880+function S = VarArray(fp)
1881+% Read VarArray into structure
1882+ if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
1883+ for i = 1:fread(fp,1,'uint32=>uint32')
1884+ S(i).Str = String(fp);
1885+ Num = double(fread(fp,1,'uint32=>uint32'));
1886+ S(i).DimID = double(fread(fp,Num,'uint32=>uint32'));
1887+ S(i).AttArray = AttArray(fp);
1888+ S(i).Type = fread(fp,1,'uint32=>uint32');
1889+ S(i).VSize = fread(fp,1,'uint32=>uint32');
1890+ S(i).Begin = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
1891+ end
1892+ else fread(fp,1,'uint32=>uint32'); S = []; end
1893+end
Note: See TracBrowser for help on using the repository browser.