Merge lp:~alexandros-avdis/meshing/boundaries_from_bathymetry into lp:meshing

[Alexandros Avdis]

Some code clean-up:Added doc string and a first attempt to divide into modules.

[lp:~alexandros-avdis/meshing/boundaries_from_bathymetry]

9. By Alexandros Avdis on 2012-05-18

Modified the example-usage message to reflect some changes.

[Alexandros Avdis (alexandros-avdis)]

Tried running the example-usage commands and the output geo file can be opened with Gmsh. However did not tried meshing.

[Adam Candy (asc)]

Can we create a (simple) test framework before these changes go in please? (since they directly modify the original script rtopo_mask_to_stereographic.py, and possibly the original functionality).

Alternatively, we could create a new main script, called for example, 'boundaries_from_bathymetry.py', which could exist alongside the original script until we can test both properly.

[Simon Mouradian (mouradian)]

I agree with this. This is a fairly large commit, so we need to have some 'testing' in place.

Alex, is it ok if we leave this merge for later?

[Alexandros Avdis (alexandros-avdis)]

No the merge has to take place, I would rather the testing framework be put in place now. The reason I see this as important is that the commit is basically a clean-up, and will be useful for the UROP students as we must have a rough idea of the module layout before we start asking people to commit work

[lp:~alexandros-avdis/meshing/boundaries_from_bathymetry]

10. By Alexandros Avdis on 2012-05-18

Commit following merge from trunk.

11. By Alexandros Avdis on 2012-05-18

Commit following merge from branch.

12. By Alexandros Avdis on 2012-05-19

Some changes to make the test cases pass.

[Alexandros Avdis (alexandros-avdis)]

The modifications appear not to change the output, the branch passes the test cases.

review: Approve

=== added file 'cartographic_projections.py'
--- cartographic_projections.py	1970-01-01 00:00:00 +0000
+++ cartographic_projections.py	2012-05-19 20:09:18 +0000
@@ -0,0 +1,13 @@
+def lotlat_deg_stereographic(longitude_deg, latitude_deg, central_longitude_deg=0.0, central_latitude_deg=-90.0, geoid_radious=6.37101e+06, central_point_buffer_deg=5.0):
+   '''Function converting lognitude and latitude into x & y coordinates on a stereographic
+      projection plane tangent to the south pole.'''
+   from scipy import pi, cos, sin
+   #Convert longitude, latitude, central longitude and central latitude to radians.
+   longitude_rad = longitude_deg*pi/180
+   latitude_rad = latitude_deg*pi/180
+   central_longitude_rad = central_longitude_deg*pi/180
+   central_latitude_rad = central_latitude_deg*pi/180
+   #Calculate coondinates in new projection.
+   stereographic_x = ((2*geoid_radious)/(1 + sin(central_latitude_rad)*sin(latitude_rad) + cos(central_latitude_rad)*cos(latitude_rad)*cos(longitude_rad-central_longitude_rad)))*cos(latitude_rad)*sin(longitude_rad-central_longitude_rad)
+   stereographic_y = ((2*geoid_radious)/(1 + sin(central_latitude_rad)*sin(latitude_rad) + cos(central_latitude_rad)*cos(latitude_rad)*cos(longitude_rad-central_longitude_rad)))*(cos(central_latitude_rad)*sin(latitude_rad) - sin(central_latitude_rad)*cos(latitude_rad)*cos(longitude_rad-central_longitude_rad))
+   return stereographic_x, stereographic_y

=== added file 'gmsh.py'
--- gmsh.py	1970-01-01 00:00:00 +0000
+++ gmsh.py	2012-05-19 20:09:18 +0000
@@ -0,0 +1,236 @@
+def gmsh_geo_header(output):
+  '''Function writing the header to the Gmsh .geo file. The header consits of creation
+     of new IDs for point, line, line-loops, surface and field entities. Then the first
+     point is drawn on the centre of the planet and the second point is drawn on the North
+     pole. The last statement in the header instructs Gmsh that the following geometry
+     commands are on the strereographic projection plane. Note that the planet is
+     modelled as a sphere.'''
+  earth_radius = 6.37101e+06
+  output.write( '''
+IP = newp;
+IL = newl;
+ILL = newll;
+IS = news;
+IFI = newf;
+Point ( IP + 0 ) = { 0, 0, 0 };
+Point ( IP + 1 ) = { 0, 0, %(earth_radius)g };
+PolarSphere ( IS + 0 ) = { IP, IP + 1 };
+
+''' % { 'earth_radius': earth_radius } )
+
+
+def gmsh_geo_footer(output, loopstart, loopend):
+  output.write( '''
+Field [ IFI + 0 ]  = Attractor;
+Field [ IFI + 0 ].NodesList  = { IP + %(loopstart)i : IP + %(loopend)i };
+''' % { 'loopstart':loopstart, 'loopend':loopend } )
+
+def gmsh_geo_remove_projection_points(output):
+  '''Function deleting the points placed at the centre of the geoid and at the north pole.'''
+  output.write( '''
+Delete { Point{1}; }
+Delete { Point{2}; }
+''' )
+
+
+def gmsh_geo_comment(output, comment):
+  '''Function witing a single-line comment to Gmsh geo script file.'''
+  output.write( '// ' + comment + '\n')
+
+
+def gmsh_geo_draw_point(output, index, loc, z):
+  '''Function writing (drawing) a point to the Gmsh geo script file.'''
+  accuracy = '.8'
+  format = 'Point ( IP + %%i ) = { %%%(dp)sf, %%%(dp)sf, %%%(dp)sf };\n' % { 'dp': accuracy }
+  output.write(format % (index, loc[0], loc[1], z))
+
+
+def gmsh_geo_draw_loop(output, index, loopstartpoint, last, open):
+  '''Function writing (drawing) a line-loop to the Gmsh geo script file.'''
+  if (index.point <= index.start):
+    return index
+  #pointstart = indexstart
+  #pointend   = index.point
+  #loopnumber = index.loop
+  if (last):
+    closure = ', IP + %(pointstart)i' % { 'pointstart':loopstartpoint }
+  else:
+    closure = ''
+  if (open):
+    index.open.append(index.path)
+    type = 'open'
+  else:
+    index.contour.append(index.path)
+    type = 'contour'
+
+  index.pathsinloop.append(index.path)
+
+#//Line Loop( ILL + %(loopnumber)i ) = { IL + %(loopnumber)i };
+#// Identified as a %(type)s path
+
+  output.write( '''LoopStart%(loopnumber)i = IP + %(pointstart)i;
+LoopEnd%(loopnumber)i = IP + %(pointend)i;
+BSpline ( IL + %(loopnumber)i ) = { IP + %(pointstart)i : IP + %(pointend)i%(loopstartpoint)s };
+Physical Line( IL + %(loopnumber)i ) = { IL + %(loopnumber)i };
+
+''' % { 'pointstart':index.start, 'pointend':index.point, 'loopnumber':index.path, 'loopstartpoint':closure, 'type':type } )
+
+  if (last):
+    output.write( '''Line Loop( ILL + %(loop)i ) = { %(loopnumbers)s };
+''' % { 'loop':index.loop , 'loopnumbers':list_to_comma_separated(index.pathsinloop, prefix = 'IL + ') } )
+    index.loops.append(index.loop)
+    index.loop += 1
+    index.pathsinloop = []
+
+  index.path +=1
+  index.start = index.point
+  return index
+
+
+def gmsh_geo_define_point(output, name, location):
+  '''Function writing a point definition in longitude-latitude as wall as in cartesian
+     coordinates in the stereographic projection plane. No points are drawn in Gmsh,
+     variables are assigned'''
+  # location [long, lat]
+  output.write('''
+//Point %(name)s is located at, %(longitude).2f deg, %(latitude).2f deg.
+Point_%(name)s_longitude_rad = (%(longitude)f + (00/60))*(Pi/180);
+Point_%(name)s_latitude_rad  = (%(latitude)f + (00/60))*(Pi/180);
+Point_%(name)s_stereographic_y = Cos(Point_%(name)s_longitude_rad)*Cos(Point_%(name)s_latitude_rad)  / ( 1 + Sin(Point_%(name)s_latitude_rad) );
+Point_%(name)s_stereographic_x = Cos(Point_%(name)s_latitude_rad) *Sin(Point_%(name)s_longitude_rad) / ( 1 + Sin(Point_%(name)s_latitude_rad) );
+''' % { 'name':name, 'longitude':location[0], 'latitude':location[1] } )
+
+
+def gmsh_geo_draw_parallel(output, startn, endn, start, end, points=200):
+  '''Function for drawing a parallel-segment through specified points.'''
+  startp = project(start)
+  endp = project(end)
+
+  output.write('''
+pointsOnParallel = %(points)i;
+parallelSectionStartingX = %(start_x)g;
+parallelSectionStartingY = %(start_y)g;
+firstPointOnParallel = IP + %(start_n)i;
+parallelSectionEndingX = %(end_x)g;
+parallelSectionEndingY = %(end_y)g;
+lastPointOnParallel = IP + %(end_n)i;
+newParallelID = IL + 10100;
+Call DrawParallel;
+''' % { 'start_x':startp[0], 'start_y':startp[1], 'end_x':endp[0], 'end_y':endp[1], 'start_n':startn, 'end_n':endn, 'points':points })
+
+
+def gmsh_geo_define_surfaces(output, index, boundary):
+  '''Function declaring plane surfaces in gmsh geo script file.'''
+  boundary_list = list_to_comma_separated(index.contour + index.open)
+#//Line Loop( ILL + %(loopnumber)i ) = { %(boundary_list)s };
+#//Plane Surface( %(surface)i ) = { ILL + %(loopnumber)i };
+  output.write('''
+Plane Surface( %(surface)i ) = { %(boundary_list)s };
+Physical Surface( %(surface)i ) = { %(surface)i };
+''' % { 'loopnumber':index.path, 'surface':boundary.surface + 1, 'boundary_list':list_to_comma_separated(index.loops, prefix = 'ILL + ') } )
+
+
+def gmsh_geo_output_fields(output, index,boundary):
+  '''Function writing fields controlling mesh size to gmsh geo file. This function also writes
+     some other options to the geo file, an operation that should probably be performed somewhere else.'''
+  if (index.contour is not None):
+    output.write('''
+Printf("Assigning characteristic mesh sizes...");
+// Field[ IFI + 1] = Attractor;
+// Field[ IFI + 1].EdgesList = { 999999, %(boundary_list)s };
+// Field [ IFI + 1 ].NNodesByEdge = 5e4;
+// 
+// Field[ IFI + 2] = Threshold;
+// Field[ IFI + 2].DistMax = 2e6;
+// Field[ IFI + 2].DistMin = 3e4;
+// Field[ IFI + 2].IField = IFI + 1;
+// Field[ IFI + 2].LcMin = 5e4;
+// Field[ IFI + 2].LcMax = 2e5;
+//
+// Background Field = IFI + 2;
+// Dont extent the elements sizes from the boundary inside the domain
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
+
+Field[ IFI + 1] = MathEval;
+Field[ IFI + 1].F = "1.0E5";
+Background Field = IFI + 1;
+
+//Set some options for better png output
+General.Color.Background = {255,255,255};
+General.Color.BackgroundGradient = {255,255,255};
+General.Color.Foreground = Black;
+Mesh.Color.Lines = {0,0,0};
+
+General.Trackball = 0 ;
+General.RotationX = 180;
+General.RotationY = 0;
+General.RotationZ = 270;
+''' % { 'boundary_list':list_to_comma_separated(index.contour, prefix = 'IL + ') } )
+
+
+def gmsh_geo_define_loop(output, index, loopstartpoint, last, open):
+  '''Function writing line-loops to gmsh geo file.'''
+  if (index.point <= index.start):
+    return index
+  #pointstart = indexstart
+  #pointend   = index.point
+  #loopnumber = index.loop
+  if (last):
+    closure = ', IP + %(pointstart)i' % { 'pointstart':loopstartpoint }
+  else:
+    closure = ''
+  if (open):
+    index.open.append(index.path)
+    type = 'open'
+  else:
+    index.contour.append(index.path)
+    type = 'contour'
+
+  index.pathsinloop.append(index.path)
+
+#//Line Loop( ILL + %(loopnumber)i ) = { IL + %(loopnumber)i };
+#// Identified as a %(type)s path
+
+  output.write( '''LoopStart%(loopnumber)i = IP + %(pointstart)i;
+LoopEnd%(loopnumber)i = IP + %(pointend)i;
+BSpline ( IL + %(loopnumber)i ) = { IP + %(pointstart)i : IP + %(pointend)i%(loopstartpoint)s };
+Physical Line( IL + %(loopnumber)i ) = { IL + %(loopnumber)i };
+
+''' % { 'pointstart':index.start, 'pointend':index.point, 'loopnumber':index.path, 'loopstartpoint':closure, 'type':type } )
+
+  if (last):
+    output.write( '''Line Loop( ILL + %(loop)i ) = { %(loopnumbers)s };
+''' % { 'loop':index.loop , 'loopnumbers':list_to_comma_separated(index.pathsinloop, prefix = 'IL + ') } )
+    index.loops.append(index.loop)
+    index.loop += 1
+    index.pathsinloop = []
+
+  index.path +=1
+  index.start = index.point
+  return index
+
+def list_to_comma_separated(numbers, prefix='', add=0):
+  '''Function converting a python list to a comma-seperated string.'''
+  requirecomma = False
+  string = ''
+  for number in numbers:
+    if (requirecomma):
+      string += ', '
+    else:
+      requirecomma = True
+    string += prefix
+    string += str(number + add)
+  return string
+
+def list_to_space_separated(numbers, prefix='', add=0):
+  '''Function converting a python list to a space-separated string.'''
+  requirespace = False
+  string = ''
+  for number in numbers:
+    if (requirespace):
+      string += ' '
+    else:
+      requirespace = True
+    string += prefix
+    string += str(number + add)
+  return string

=== modified file 'rtopo_mask_to_stereographic.py'
--- rtopo_mask_to_stereographic.py	2012-05-18 21:28:13 +0000
+++ rtopo_mask_to_stereographic.py	2012-05-19 20:09:18 +0000
@@ -7,12 +7,14 @@
 from Scientific.IO import NetCDF
 import matplotlib
 matplotlib.use('Agg')
+import matplotlib.pyplot
 from pylab import contour
 #import matplotlib
 #matplotlib._cntr.Cntr
 #from matplotlib import contour
 #matplotlib.use('Agg')
 from numpy import zeros, array, append, exp
+import gmsh
 
 #contour = matplotlib.pyplot.contour
 
@@ -27,15 +29,12 @@
 def printv(text):
   if (arguments.verbose):
     print text
-  gmsh_comment(text)
+  gmsh.gmsh_geo_comment(output, text)
 
 def printvv(text):
   if (arguments.debug):
     print text
 
-def gmsh_comment(comment):
-  output.write( '// ' + comment + '\n')
-
 def expand_boxes(region, boxes):
   def error():
     print 'Error in argument for -b.'
@@ -104,15 +103,15 @@
 ------------------------------------------------------------
 Example usage:
 Include only the main Antarctic mass (path 1), and only parts which lie below 60S
-  rtopo_mask_to_stereographic.py -r 'latitude <= -60.0' -p 1
+  ./rtopo_mask_to_stereographic.py RTopo105b_50S.nc -r 'latitude <= -60.0' -p 1
 Filchner-Ronne extended out to the 65S parallel
-  rtopo_mask_to_stereographic.py -no -b -85.0:-20.0,-89.0:-75.0 -64.0:-30.0,-89.0:-70.0 -30.0:-20.0,-89.0:-75.0 -lat '-65.0'
+  ./rtopo_mask_to_stereographic.py RTopo105b_50S.nc -no -b -85.0:-20.0,-89.0:-75.0 -64.0:-30.0,-89.0:-70.0 -30.0:-20.0,-89.0:-75.0 -lat '-65.0'
 Antarctica, everything below the 60S parallel, coarse approximation to open boundary
-  rtopo_mask_to_stereographic.py -dx 2 -r 'latitude <= -60'
+  ./rtopo_mask_to_stereographic.py RTopo105b_50S.nc -dx 2 -r 'latitude <= -60'
 Small region close to the Filcher-Ronne ice shelf
-  rtopo_mask_to_stereographic.py -no -b -85.0:-20.0,-89.0:-75.0 -64.0:-30.0,-89.0:-70.0 -30.0:-20.0,-89.0:-75.0 -p 1 -r 'latitude <= -83'
+  ./rtopo_mask_to_stereographic.py RTopo105b_50S.nc -no -b -85.0:-20.0,-89.0:-75.0 -64.0:-30.0,-89.0:-70.0 -30.0:-20.0,-89.0:-75.0 -p 1 -r 'latitude <= -83'
 Amundsen Sea
-  rtopo_mask_to_stereographic.py -no -b -130.0:-85.0,-85.0:-60.0 -lat -64.0
+  ./rtopo_mask_to_stereographic.py RTopo105b_50S.nc -no -b -130.0:-85.0,-85.0:-60.0 -lat -64.0
 
 Small islands, single out, or group with -p
   312, 314
@@ -126,7 +125,7 @@
 argv = sys.argv[1:]
 dx_default = 0.1
 class arguments:
-  input  = '/d/dataset/rtopo/RTopo105b_50S.nc'
+  input  = './RTopo105b_50S.nc'
   #output = './stereographic_projection.geo'
   output = './shorelines.geo'
   boundaries = []
@@ -174,10 +173,10 @@
 arguments.region = expand_boxes(arguments.region, arguments.box)
  
 
-#source = file(arguments.input,'r')
+source = file(arguments.input,'r')
 output = file(arguments.output,'w')
 
-gmsh_comment('Arguments: ' + arguments.call)
+gmsh.gmsh_geo_comment(output, 'Arguments: ' + arguments.call)
 printv('Source netCDF located at ' + arguments.input)
 printv('Output to ' + arguments.output)
 if (len(arguments.boundaries) > 0):
@@ -189,21 +188,7 @@
 if (arguments.extendtolatitude is not None):
   printv('Extending region to meet parallel on latitude ' + str(arguments.extendtolatitude))
 
-gmsh_comment('')
-
-def gmsh_header():
-  earth_radius = 6.37101e+06
-  return '''
-IP = newp;
-IL = newl;
-ILL = newll;
-IS = news;
-IFI = newf;
-Point ( IP + 0 ) = { 0, 0, 0 };
-Point ( IP + 1 ) = { 0, 0, %(earth_radius)g };
-PolarSphere ( IS + 0 ) = { IP, IP + 1 };
-
-''' % { 'earth_radius': earth_radius }
+gmsh.gmsh_geo_comment(output, '')
 
 def smoothGaussian(list,degree,strippedXs=False):
   list = list.tolist()
@@ -221,25 +206,6 @@
     smoothed[i]=sum(array(list[i:i+window])*weight)/sum(weight)
   return array(smoothed)
 
-def gmsh_footer(loopstart, loopend):
-  output.write( '''
-Field [ IFI + 0 ]  = Attractor;
-Field [ IFI + 0 ].NodesList  = { IP + %(loopstart)i : IP + %(loopend)i };
-''' % { 'loopstart':loopstart, 'loopend':loopend } )
-
-def gmsh_remove_projection_points():
-  output.write( '''
-Delete { Point{1}; }
-Delete { Point{2}; }
-''' )
-
-
-def gmsh_format_point(index, loc, z):
-  accuracy = '.8'
-  format = 'Point ( IP + %%i ) = { %%%(dp)sf, %%%(dp)sf, %%%(dp)sf };\n' % { 'dp': accuracy }
-  output.write(format % (index, loc[0], loc[1], z))
-  #return "Point ( IP + %i ) = { %f, %f, %f }\n" % (index, x, y, z)
-
 def project(location):
   longitude = location[0]
   latitude  = location[1]
@@ -289,7 +255,7 @@
   return eval(region, globals)
 
 def array_to_gmsh_points(num, index, location, minarea, region, dx, latitude_max):
-  gmsh_comment('Ice-Land mass number %s' % (num))
+  gmsh.gmsh_geo_comment(output, 'Ice-Land mass number %s' % (num))
   count = 0 
   pointnumber = len(location[:,0])
   valid = [False]*pointnumber
@@ -314,7 +280,7 @@
     
   if (loopend is None):
     printvv('Path %i skipped (no points found in region)' % ( num ))
-    gmsh_comment('  Skipped (no points found in region)\n')
+    gmsh.gmsh_geo_comment(output, '  Skipped (no points found in region)\n')
     return index
   
   closelast=False
@@ -359,7 +325,7 @@
   area = area_enclosed(validlocation)
   if (area < minarea):
     printvv('Path %i skipped (area too small)' % ( num ))
-    gmsh_comment('  Skipped (area too small)\n')
+    gmsh.gmsh_geo_comment(output, '  Skipped (area too small)\n')
     return index
 
   printv('Path %i points %i/%i area %g%s' % ( num, validnumber, pointnumber, area_enclosed(validlocation), closingtext ))
@@ -372,16 +338,16 @@
   #     latitude  = validlocation[point,1]
   #     index += 1
   #     loc = project(longitude, latitude)
-  #     output.write( gmsh_format_point(index, loc, 0) )
+  #     gmsh.gmsh_geo_draw_point(output, index, loc, 0) )
   #     if (latitude_max is None):
   #       latitude_max = latitude
   #     else:
   #       latitude_max = max(latitude_max, latitude)
-  #   draw_parallel(index, index_start, [ validlocation[point,0], max(latitude_max, validlocation[point,1]) ], [ validlocation[0,0], max(latitude_max, validlocation[0,1]) ], points=200)
+  #   gmsh.gmsh_geo_draw_parallel(output, index, index_start, [ validlocation[point,0], max(latitude_max, validlocation[point,1]) ], [ validlocation[0,0], max(latitude_max, validlocation[0,1]) ], points=200)
   #   index += 200
   #   
   #   index += 1
-  #   output.write( gmsh_format_point(index, project(validlocation[0,0], validlocation[0,1]), 0) )
+  #   gmsh.gmsh_geo_draw_point(output, index, project(validlocation[0,0], validlocation[0,1]), 0) )
   #   
   # else:
   if (close[0]):
@@ -394,24 +360,24 @@
     #latitude  = validlocation[point,1]
     
     if ((close[point]) and (point == validnumber - 1) and (not (compare_points(validlocation[point], validlocation[0], dx)))):
-      gmsh_comment('**** END ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
-      index = gmsh_loop(index, loopstartpoint, False, False)
+      gmsh.gmsh_geo_comment(output, '**** END ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
+      index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, False, False)
       index = draw_parallel_explicit(validlocation[point], validlocation[0], index, latitude_max, dx)
-      index = gmsh_loop(index, loopstartpoint, True, True)
-      gmsh_comment('**** END end of loop ' + str(closelast) + str(point) + '/' + str(validnumber-1) + str(close[point]))
+      index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, True, True)
+      gmsh.gmsh_geo_comment(output, '**** END end of loop ' + str(closelast) + str(point) + '/' + str(validnumber-1) + str(close[point]))
     elif ((close[point]) and (point > 0) and (not (compare_points(validlocation[point], validlocation[0], dx)))):
-      gmsh_comment('**** NOT END ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
-      gmsh_comment(str(validlocation[point,:]) + str(validlocation[point,:]))
-      index = gmsh_loop(index, loopstartpoint, False, False)
+      gmsh.gmsh_geo_comment(output, '**** NOT END ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
+      gmsh.gmsh_geo_comment(output, str(validlocation[point,:]) + str(validlocation[point,:]))
+      index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, False, False)
       index = draw_parallel_explicit(validlocation[point - 1], validlocation[point], index, latitude_max, dx)
-      index = gmsh_loop(index, loopstartpoint, False, True)
-      gmsh_comment('**** NOT END end of loop ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
+      index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, False, True)
+      gmsh.gmsh_geo_comment(output, '**** NOT END end of loop ' + str(point) + '/' + str(validnumber-1) + str(close[point]))
     else:
       index.point += 1
-      gmsh_format_point(index.point, project(validlocation[point,:]), 0)
+      gmsh.gmsh_geo_draw_point(output, index.point, project(validlocation[point,:]), 0)
       index.contournodes.append(index.point)
 
-  index = gmsh_loop(index, loopstartpoint, (closelast and (point == validnumber - 1)), False)
+  index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, (closelast and (point == validnumber - 1)), False)
 
   return index
 
@@ -426,53 +392,10 @@
 #Line Loop( ILL + 20 ) = { IL + 2 };
 
 
-def gmsh_loop(index, loopstartpoint, last, open):
-  if (index.point <= index.start):
-    return index
-  #pointstart = indexstart
-  #pointend   = index.point
-  #loopnumber = index.loop
-  if (last):
-    closure = ', IP + %(pointstart)i' % { 'pointstart':loopstartpoint }
-  else:
-    closure = ''
-  if (open):
-    index.open.append(index.path)
-    type = 'open'
-  else:
-    index.contour.append(index.path)
-    type = 'contour'
-
-  index.pathsinloop.append(index.path)
-
-#//Line Loop( ILL + %(loopnumber)i ) = { IL + %(loopnumber)i };
-#// Identified as a %(type)s path
-
-  output.write( '''LoopStart%(loopnumber)i = IP + %(pointstart)i;
-LoopEnd%(loopnumber)i = IP + %(pointend)i;
-BSpline ( IL + %(loopnumber)i ) = { IP + %(pointstart)i : IP + %(pointend)i%(loopstartpoint)s };
-Physical Line( IL + %(loopnumber)i ) = { IL + %(loopnumber)i };
-
-''' % { 'pointstart':index.start, 'pointend':index.point, 'loopnumber':index.path, 'loopstartpoint':closure, 'type':type } )
-
-  if (last):
-    output.write( '''Line Loop( ILL + %(loop)i ) = { %(loopnumbers)s };
-''' % { 'loop':index.loop , 'loopnumbers':list_to_comma_separated(index.pathsinloop, prefix = 'IL + ') } )
-    index.loops.append(index.loop)
-    index.loop += 1
-    index.pathsinloop = []
-
-  index.path +=1
-  index.start = index.point
-  return index
-
-
-
-
-def output_boundaries(index, filename, paths=None, minarea=0, region='True', dx=dx_default, latitude_max=None):
+def output_boundaries(index, filename, paths=None, minarea=0, region='True', dx=0.1, latitude_max=None):
   pathall = read_rtopo(filename)
   printv('Paths found: ' + str(len(pathall)))
-  output.write( gmsh_header() )
+  gmsh.gmsh_geo_header(output)
   splinenumber = 0
   indexbase = 1
   index.point = indexbase
@@ -493,54 +416,41 @@
     index = array_to_gmsh_points(num, index, xy, minarea, region, dx, latitude_max)
   #for i in range(-85, 0, 5):
   #  indexend += 1
-  #  output.write( gmsh_format_point(indexend, project(0, i), 0) )
+  #  gmsh.gmsh_geo_draw_point(output, indexend, project(0, i), 0) )
   #for i in range(-85, 0, 5):
   #  indexend += 1
-  #  output.write( gmsh_format_point(indexend, project(45, i), 0) )
-  gmsh_remove_projection_points()
+  #  gmsh.gmsh_geo_draw_point(output, indexend, project(45, i), 0) )
+  gmsh.gmsh_geo_remove_projection_points(output)
   return index
 
-def define_point(name, location):
-  # location [long, lat]
-  output.write('''
-//Point %(name)s is located at, %(longitude).2f deg, %(latitude).2f deg.
-Point_%(name)s_longitude_rad = (%(longitude)f + (00/60))*(Pi/180);
-Point_%(name)s_latitude_rad  = (%(latitude)f + (00/60))*(Pi/180);
-Point_%(name)s_stereographic_y = Cos(Point_%(name)s_longitude_rad)*Cos(Point_%(name)s_latitude_rad)  / ( 1 + Sin(Point_%(name)s_latitude_rad) );
-Point_%(name)s_stereographic_x = Cos(Point_%(name)s_latitude_rad) *Sin(Point_%(name)s_longitude_rad) / ( 1 + Sin(Point_%(name)s_latitude_rad) );
-''' % { 'name':name, 'longitude':location[0], 'latitude':location[1] } )
-
-def draw_parallel(startn, endn, start, end, points=200):
-  startp = project(start)
-  endp = project(end)
-  
-  output.write('''
-pointsOnParallel = %(points)i;
-parallelSectionStartingX = %(start_x)g;
-parallelSectionStartingY = %(start_y)g;
-firstPointOnParallel = IP + %(start_n)i;
-parallelSectionEndingX = %(end_x)g;
-parallelSectionEndingY = %(end_y)g;
-lastPointOnParallel = IP + %(end_n)i;
-newParallelID = IL + 10100;
-Call DrawParallel;
-''' % { 'start_x':startp[0], 'start_y':startp[1], 'end_x':endp[0], 'end_y':endp[1], 'start_n':startn, 'end_n':endn, 'points':points })
-
 def compare_points(a, b, dx):
   tolerance = dx * 0.6
   if ( not (abs(a[1] - b[1]) < tolerance) ):
-    #gmsh_comment('lat differ')
+    #gmsh.gmsh_geo_comment(output, 'lat differ')
     return False
   elif (abs(a[0] - b[0]) < tolerance):
-    #gmsh_comment('long same')
+    #gmsh.gmsh_geo_comment(output, 'long same')
     return True
   elif ((abs(abs(a[0]) - 180) < tolerance) and (abs(abs(b[0]) - 180) < tolerance)):
-    #gmsh_comment('long +/-180')
+    #gmsh.gmsh_geo_comment(output, 'long +/-180')
     return True
   else:
-    #gmsh_comment('not same %g %g' % (abs(abs(a[0]) - 180), abs(abs(b[0]) - 180) ) )
+    #gmsh.gmsh_geo_comment(output, 'not same %g %g' % (abs(abs(a[0]) - 180), abs(abs(b[0]) - 180) ) )
     return False
 
+
+def output_open_boundaries(index, boundary, dx):
+  parallel = arguments.bounding_lat
+  index.start = index.point + 1
+  loopstartpoint = index.start
+  index = draw_parallel_explicit([   -1.0, parallel], [ 179.0, parallel], index, None, dx)
+  index = draw_parallel_explicit([-179.0,  parallel], [   1.0, parallel], index, None, dx)
+  
+  index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, True, True)
+
+  return index
+
+
 def draw_parallel_explicit(start, end, index, latitude_max, dx):
 
   #print start, end, index.point
@@ -549,16 +459,16 @@
     latitude_max = max(start[1], end[1])
   else:
     latitude_max = max(latitude_max, start[1], end[1])
-  current = start 
+  current = start
   tolerance = dx * 0.6
-  
-  gmsh_comment( 'Closing path with parallels and merdians, from (%.8f, %.8f) to  (%.8f, %.8f)' % ( start[0], start[1], end[0], end[1] ) )
-  
+
+  gmsh.gmsh_geo_comment(output, 'Closing path with parallels and merdians, from (%.8f, %.8f) to  (%.8f, %.8f)' % ( start[0], start[1], end[0], end[1] ) )
+
   if (compare_points(current, end, dx)):
-    gmsh_comment('Points already close enough, no need to draw parallels and meridians after all')
+    gmsh.gmsh_geo_comment(output, 'Points already close enough, no need to draw parallels and meridians after all')
     return index
-  
-  gmsh_comment('Drawing meridian to max latitude index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], latitude_max))
+
+  gmsh.gmsh_geo_comment(output, 'Drawing meridian to max latitude index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], latitude_max))
   while (current[1] != latitude_max):
     if (current[1] < latitude_max):
       current[1] = current[1] + dx
@@ -569,9 +479,9 @@
     index.point += 1
     printvv('Drawing meridian to max latitude index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], latitude_max))
     loc = project(current)
-    gmsh_format_point(index.point, loc, 0.0)
+    gmsh.gmsh_geo_draw_point(output, index.point, loc, 0.0)
 
-  gmsh_comment('Drawing parallel index %s at %f.2 (to match %f.2), %f.2' % (index.point, current[0], end[0], current[1]))
+  gmsh.gmsh_geo_comment(output, 'Drawing parallel index %s at %f.2 (to match %f.2), %f.2' % (index.point, current[0], end[0], current[1]))
   while (current[0] != end[0]):
     if (current[0] < end[0]):
       current[0] = current[0] + dx
@@ -582,9 +492,9 @@
     index.point += 1
     printvv('Drawing parallel index %s at %f.2 (to match %f.2), %f.2' % (index.point, current[0], end[0], current[1]))
     loc = project(current)
-    gmsh_format_point(index.point, loc, 0.0)
+    gmsh.gmsh_geo_draw_point(output, index.point, loc, 0.0)
 
-  gmsh_comment('Drawing meridian to end index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], end[1]))
+  gmsh.gmsh_geo_comment(output, 'Drawing meridian to end index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], end[1]))
   while (current[1] != end[1]):
     if (current[1] < end[1]):
       current[1] = current[1] + dx
@@ -595,57 +505,12 @@
     index.point += 1
     printvv('Drawing meridian to end index %s at %f.2, %f.2 (to match %f.2)' % (index.point, current[0], current[1], end[1]))
     loc = project(current)
-    gmsh_format_point(index.point, loc, 0.0)
-  
-  gmsh_comment( 'Closed path with parallels and merdians, from (%.8f, %.8f) to  (%.8f, %.8f)' % ( start[0], start[1], end[0], end[1] ) )
-
-  return index
-
-def list_to_comma_separated(numbers, prefix='', add=0):
-  requirecomma = False
-  string = ''
-  for number in numbers:
-    if (requirecomma):
-      string += ', '
-    else:
-      requirecomma = True
-    string += prefix
-    string += str(number + add)
-  return string
-
-def list_to_space_separated(numbers, prefix='', add=0):
-  requirespace = False
-  string = ''
-  for number in numbers:
-    if (requirespace):
-      string += ' '
-    else:
-      requirespace = True
-    string += prefix
-    string += str(number + add)
-  return string
-
-def output_open_boundaries(index, boundary, dx):
-  parallel = arguments.bounding_lat
-  index.start = index.point + 1
-  loopstartpoint = index.start
-  index = draw_parallel_explicit([   -1.0, parallel], [ 179.0, parallel], index, None, dx)
-  index = draw_parallel_explicit([-179.0,  parallel], [   1.0, parallel], index, None, dx)
-  
-  index = gmsh_loop(index, loopstartpoint, True, True)
-
-  return index
-
-def output_surfaces(index, boundary):
-  printv('Open boundaries   (id %i): %s' % (boundary.open, list_to_space_separated(index.open, add=1)))
-  printv('Closed boundaries (id %i): %s' % (boundary.contour, list_to_space_separated(index.contour, add=1)))
-  boundary_list = list_to_comma_separated(index.contour + index.open)
-#//Line Loop( ILL + %(loopnumber)i ) = { %(boundary_list)s };
-#//Plane Surface( %(surface)i ) = { ILL + %(loopnumber)i };
-  output.write('''
-Plane Surface( %(surface)i ) = { %(boundary_list)s };
-Physical Surface( %(surface)i ) = { %(surface)i };
-''' % { 'loopnumber':index.path, 'surface':boundary.surface + 1, 'boundary_list':list_to_comma_separated(index.loops, prefix = 'ILL + ') } )
+    gmsh.gmsh_geo_draw_point(output, index.point, loc, 0.0)
+
+  gmsh.gmsh_geo_comment(output, 'Closed path with parallels and merdians, from (%.8f, %.8f) to  (%.8f, %.8f)' % ( start[0], start[1], end[0], end[1] ) )
+
+  return index
+
 
 def acc_array():
   acc = array([[   1.0, -53.0 ],
@@ -691,13 +556,13 @@
 
 def draw_acc_old(index, boundary, dx):
   acc = acc_array()
-  gmsh_comment('ACC')
+  gmsh.gmsh_geo_comment(output, 'ACC')
   index.start = index.point + 1
   loopstartpoint = index.start
   for i in range(len(acc[:,0])):
     index.point += 1
     location = project(acc[i,:])
-    gmsh_format_point(index.point, location, 0.0)
+    gmsh.gmsh_geo_draw_point(output, index.point, location, 0.0)
 
   for i in range(len(acc[:,0])):
     a = index.start + i
@@ -715,62 +580,27 @@
   acc2 = acc[19:,:]
   print acc1
   print acc2
-  gmsh_comment('ACC')
+  gmsh.gmsh_geo_comment(output, 'ACC')
 
   index.start = index.point + 1
   loopstartpoint = index.start
   for i in range(len(acc1[:,0])):
     index.point += 1
     location = project(acc1[i,:])
-    gmsh_format_point(index.point, location, 0.0)
-  index = gmsh_loop(index, loopstartpoint, False, True)
+    gmsh.gmsh_geo_draw_point(output, index.point, location, 0.0)
+  index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, False, True)
 
   #index.start = index.point + 1
   #loopstartpoint = index.start
   for i in range(len(acc2[:,0])):
     index.point += 1
     location = project(acc2[i,:])
-    gmsh_format_point(index.point, location, 0.0)
-  index = gmsh_loop(index, loopstartpoint, True, True)
+    gmsh.gmsh_geo_draw_point(output, index.point, location, 0.0)
+  index = gmsh.gmsh_geo_draw_loop(output, index, loopstartpoint, True, True)
 
   return index
 
 
-def output_fields(index,boundary):
-  if (index.contour is not None):
-    output.write('''
-Printf("Assigning characteristic mesh sizes...");
-// Field[ IFI + 1] = Attractor;
-// Field[ IFI + 1].EdgesList = { 999999, %(boundary_list)s };
-// Field [ IFI + 1 ].NNodesByEdge = 5e4;
-// 
-// Field[ IFI + 2] = Threshold;
-// Field[ IFI + 2].DistMax = 2e6;
-// Field[ IFI + 2].DistMin = 3e4;
-// Field[ IFI + 2].IField = IFI + 1;
-// Field[ IFI + 2].LcMin = 5e4;
-// Field[ IFI + 2].LcMax = 2e5;
-//
-// Background Field = IFI + 2;
-// Dont extent the elements sizes from the boundary inside the domain
-//Mesh.CharacteristicLengthExtendFromBoundary = 0;
-
-Field[ IFI + 1] = MathEval;
-Field[ IFI + 1].F = "1.0E5";
-Background Field = IFI + 1;
-
-//Set some options for better png output
-General.Color.Background = {255,255,255};
-General.Color.BackgroundGradient = {255,255,255};
-General.Color.Foreground = Black;
-Mesh.Color.Lines = {0,0,0};
-
-General.Trackball = 0 ;
-General.RotationX = 180;
-General.RotationY = 0;
-General.RotationZ = 270;
-''' % { 'boundary_list':list_to_comma_separated(index.contour, prefix = 'IL + ') } )
-
 class index:
   point = 0
   path = 0
@@ -793,12 +623,14 @@
 
 
 if (arguments.open): index = output_open_boundaries(index, boundary, arguments.dx)
-output_surfaces(index, boundary)
+printv('Open boundaries   (id %i): %s' % (boundary.open, gmsh.list_to_space_separated(index.open, add=1)))
+printv('Closed boundaries (id %i): %s' % (boundary.contour, gmsh.list_to_space_separated(index.contour, add=1)))
+gmsh.gmsh_geo_define_surfaces(output, index, boundary)
 
 #index = draw_acc(index, boundary, arguments.dx)
 
 
-output_fields(index,boundary)
+gmsh.gmsh_geo_output_fields(output, index,boundary)
 
 
 if (len(index.skipped) > 0):