Index: /issm/trunk-jpl/src/m/contrib/bamg/gmsh3d.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/bamg/gmsh3d.m	(revision 18153)
+++ /issm/trunk-jpl/src/m/contrib/bamg/gmsh3d.m	(revision 18153)
@@ -0,0 +1,103 @@
+function md=gmsh3d(md,domainfile,h),
+
+%Read domain
+domain=expread(domainfile);
+x   = domain.x(1:end-1);
+y   = domain.y(1:end-1);
+nbv = numel(x);
+
+%write files
+t1=clock; fprintf('%s','      writing initial mesh files...');
+fid=fopen('model.geo','w');
+fprintf(fid,['// Gmsh input file, created by ISSM on ' date '\n']);
+
+%Lower points
+for i=1:nbv,
+	fprintf(fid,'Point(%i) = {%8g, %8g, %8g};\n',i,x(i),y(i),0.);
+end
+%Upper points
+for i=1:nbv,
+	fprintf(fid,'Point(%i) = {%8g, %8g, %8g};\n',nbv+i,x(i),y(i),1.);
+end
+
+%Lower lines
+for i=1:nbv-1
+	fprintf(fid,'Line(%i) = {%i, %i};\n',i,i,i+1);
+end
+fprintf(fid,'Line(%i) = {%i, %i};\n',nbv,nbv,1);
+%Upper lines
+for i=nbv+1:2*nbv-1
+	fprintf(fid,'Line(%i) = {%i, %i};\n',i,i,i+1);
+end
+fprintf(fid,'Line(%i) = {%i, %i};\n',2*nbv,2*nbv,nbv+1);
+%Side lines
+for i=1:nbv
+	fprintf(fid,'Line(%i) = {%i, %i};\n',2*nbv+i,i,nbv+i);
+end
+
+counter = 3*nbv;
+ps = zeros(nbv+2,1);
+%Lower surface Loop and surface
+counter = counter+1;
+fprintf(fid,['Line Loop(' num2str(counter) ') = {']);
+for i=1:nbv-1
+	fprintf(fid,'%i,',i);
+end
+fprintf(fid,'%i};\n',nbv);
+fprintf(fid,['Plane Surface(' num2str(counter+1) ') = {' num2str(counter) '};\n']);
+ps(1)=counter+1;
+counter = counter+1;
+%Upper surface Loop and surface
+counter = counter+1;
+fprintf(fid,['Line Loop(' num2str(counter) ') = {']);
+for i=nbv+1:2*nbv-1
+	fprintf(fid,'%i,',i);
+end
+fprintf(fid,'%i};\n',2*nbv);
+ps(2)=counter+1;
+fprintf(fid,['Plane Surface(' num2str(counter+1) ') = {' num2str(counter) '};\n']);
+counter = counter+2;
+%Sides surfaces
+for i=1:nbv-1,
+	fprintf(fid,['Line Loop(' num2str(counter) ') = {' num2str(i) ',' num2str(2*nbv+i+1) ',-' num2str(nbv+i) ',-' num2str(2*nbv+i) '};\n']);
+	fprintf(fid,['Plane Surface(' num2str(counter+1) ') = {' num2str(counter) '};\n']);
+	ps(2+i)=counter+1;
+	counter=counter+2;
+end
+fprintf(fid,['Line Loop(' num2str(counter) ') = {' num2str(nbv) ',' num2str(2*nbv+1) ',-' num2str(2*nbv) ',-' num2str(3*nbv) '};\n']);
+fprintf(fid,['Plane Surface(' num2str(counter+1) ') = {' num2str(counter) '};\n']);
+ps(2+nbv)=counter+1;
+counter=counter+2;
+
+%Physical surfaces
+counter = counter+1;
+fprintf(fid,['Surface Loop(' num2str(counter) ') = {']);
+for i=1:numel(ps)-1
+	fprintf(fid,'%i,',ps(i));
+end
+fprintf(fid,'%i};\n',ps(end));
+fprintf(fid,['Physical Surface(1) = {' num2str(ps(2)) '};\n']);
+fprintf(fid,['Physical Surface(2) = {' num2str(ps(1)) '};\n']);
+fprintf(fid,['Physical Surface(3) = {']);
+for i=3:numel(ps)-1
+	fprintf(fid,'%i,',ps(i));
+end
+fprintf(fid,'%i};\n',ps(end));
+
+%Volume
+fprintf(fid,['Volume(1) = {' num2str(counter) '};\n']);
+fprintf(fid,['Physical Volume(2) = {1};\n']);
+
+%resolution
+fprintf(fid,'Mesh.CharacteristicLengthMax = %g;',h);
+fclose(fid);
+
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%call gmsh
+fprintf('%s\n','      call gmsh...');
+system([issmdir() '/externalpackages/gmsh/install/gmsh -3 model.geo']);
+
+%plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
+md=importgmsh('model.msh',3);