Index: /issm/trunk-jpl/src/m/classes/clusters/generic.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/clusters/generic.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/classes/clusters/generic.py	(revision 24115)
@@ -12,198 +12,198 @@
 
 class generic(object):
-	"""
-	GENERIC cluster class definition
- 
-	   Usage:
-	      cluster=generic('name','astrid','np',3);
-	      cluster=generic('name',gethostname(),'np',3,'login','username');
-	"""
-
-	def __init__(self,*args):    # {{{
-
-		self.name=''
-		self.login=''
-		self.np=1
-		self.port=0
-		self.interactive=1
-		self.codepath=IssmConfig('ISSM_PREFIX')[0]+'/bin'
-		self.executionpath=issmdir()+'/execution'
-		self.valgrind=issmdir()+'/externalpackages/valgrind/install/bin/valgrind'
-		self.valgrindlib=issmdir()+'/externalpackages/valgrind/install/lib/libmpidebug.so'
-		self.valgrindsup=issmdir()+'/externalpackages/valgrind/issm.supp'
-
-		#use provided options to change fields
-		options=pairoptions(*args)
-
-		#get name
-		self.name=socket.gethostname()
-
-		#initialize cluster using user settings if provided
-		if os.path.exists(self.name+'_settings.py'):
-			exec(compile(open(self.name+'_settings.py').read(), self.name+'_settings.py', 'exec'),globals())
-
-		#OK get other fields
-		self=options.AssignObjectFields(self)
-	# }}}
-	def __repr__(self):    # {{{
-		#  display the object
-		s ="class '%s' object '%s' = \n" % (type(self),'self')
-		s+="    name: %s\n" % self.name
-		s+="    login: %s\n" % self.login
-		s+="    np: %i\n" % self.np
-		s+="    port: %i\n" % self.port
-		s+="    codepath: %s\n" % self.codepath
-		s+="    executionpath: %s\n" % self.executionpath
-		s+="    valgrind: %s\n" % self.valgrind
-		s+="    valgrindlib: %s\n" % self.valgrindlib
-		s+="    valgrindsup: %s\n" % self.valgrindsup
-		return s
-	# }}}
-	def checkconsistency(self,md,solution,analyses):    # {{{
-		if self.np<1:
-			md = checkmessage(md,'number of processors should be at least 1')
-		if math.isnan(self.np):
-			md = checkmessage(md,'number of processors should not be NaN!')
-
-		return md
-	# }}}
-	def BuildQueueScript(self,dirname,modelname,solution,io_gather,isvalgrind,isgprof,isdakota,isoceancoupling):    # {{{
-
-		executable='issm.exe';
-		if isdakota:
-			version=IssmConfig('_DAKOTA_VERSION_')
-			version=float(version[0])
-			if version>=6:
-				executable='issm_dakota.exe'
-		if isoceancoupling:
-			executable='issm_ocean.exe'
-
-		#write queuing script 
-		if not m.ispc():
-
-			fid=open(modelname+'.queue','w')
-			fid.write('#!/bin/sh\n')
-			if not isvalgrind:
-				if self.interactive: 
-					if IssmConfig('_HAVE_MPI_')[0]:
-						fid.write('mpiexec -np %i %s/%s %s %s/%s %s ' % (self.np,self.codepath,executable,solution,self.executionpath,dirname,modelname))
-					else:
-						fid.write('%s/%s %s %s/%s %s ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname))
-				else:
-					if IssmConfig('_HAVE_MPI_')[0]:
-						fid.write('mpiexec -np %i %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % (self.np,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
-					else:
-						fid.write('%s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
-			elif isgprof:
-				fid.write('\n gprof %s/%s gmon.out > %s.performance' % (self.codepath,executable,modelname))
-			else:
-				#Add --gen-suppressions=all to get suppression lines
-				fid.write('LD_PRELOAD=%s \\\n' % self.valgrindlib)
-				if IssmConfig('_HAVE_MPI_')[0]:
-					fid.write('mpiexec -np %i %s --leak-check=full --suppressions=%s %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % \
-							(self.np,self.valgrind,self.valgrindsup,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
-				else:	
-					fid.write('%s --leak-check=full --suppressions=%s %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % \
-							(self.valgrind,self.valgrindsup,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
-
-			if not io_gather:    #concatenate the output files:
-				fid.write('\ncat %s.outbin.* > %s.outbin' % (modelname,modelname))
-			fid.close()
-
-		else:    # Windows
-
-			fid=open(modelname+'.bat','w')
-			fid.write('@echo off\n')
-			if self.interactive:
-				fid.write('"%s/%s" %s "%s/%s" %s ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname))
-			else:
-				fid.write('"%s/%s" %s "%s/%s" %s 2> %s.errlog >%s.outlog' % \
-					(self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
-			fid.close()
-
-		#in interactive mode, create a run file, and errlog and outlog file
-		if self.interactive:
-			fid=open(modelname+'.errlog','w')
-			fid.close()
-			fid=open(modelname+'.outlog','w')
-			fid.close()
-	# }}}
-	def BuildKrigingQueueScript(self,modelname,solution,io_gather,isvalgrind,isgprof):    # {{{
-
-		#write queuing script 
-		if not m.ispc():
-
-			fid=open(modelname+'.queue','w')
-			fid.write('#!/bin/sh\n')
-			if not isvalgrind:
-				if self.interactive:
-					fid.write('mpiexec -np %i %s/kriging.exe %s/%s %s ' % (self.np,self.codepath,self.executionpath,modelname,modelname))
-				else:
-					fid.write('mpiexec -np %i %s/kriging.exe %s/%s %s 2> %s.errlog >%s.outlog ' % (self.np,self.codepath,self.executionpath,modelname,modelname,modelname,modelname))
-			elif isgprof:
-				fid.write('\n gprof %s/kriging.exe gmon.out > %s.performance' & (self.codepath,modelname))
-			else:
-				#Add --gen-suppressions=all to get suppression lines
-				fid.write('LD_PRELOAD=%s \\\n' % self.valgrindlib)
-				fid.write('mpiexec -np %i %s --leak-check=full --suppressions=%s %s/kriging.exe %s/%s %s 2> %s.errlog >%s.outlog ' % \
-					(self.np,self.valgrind,self.valgrindsup,self.codepath,self.executionpath,modelname,modelname,modelname,modelname))
-			if not io_gather:    #concatenate the output files:
-				fid.write('\ncat %s.outbin.* > %s.outbin' % (modelname,modelname))
-			fid.close()
-
-		else:    # Windows
-
-			fid=open(modelname+'.bat','w')
-			fid.write('@echo off\n')
-			if self.interactive:
-				fid.write('"%s/issm.exe" %s "%s/%s" %s ' % (self.codepath,solution,self.executionpath,modelname,modelname))
-			else:
-				fid.write('"%s/issm.exe" %s "%s/%s" %s 2> %s.errlog >%s.outlog' % \
-					(self.codepath,solution,self.executionpath,modelname,modelname,modelname,modelname))
-			fid.close()
-
-		#in interactive mode, create a run file, and errlog and outlog file
-		if self.interactive:
-			fid=open(modelname+'.errlog','w')
-			fid.close()
-			fid=open(modelname+'.outlog','w')
-			fid.close()
-	# }}}
-	def UploadQueueJob(self,modelname,dirname,filelist):    # {{{
-
-		#compress the files into one zip.
-		compressstring='tar -zcf %s.tar.gz ' % dirname
-		for file in filelist:
-			compressstring += ' %s' % file
-		if self.interactive:
-			compressstring += ' %s.errlog %s.outlog ' % (modelname,modelname)
-		subprocess.call(compressstring,shell=True)
-
-		print('uploading input file and queueing script')
-		issmscpout(self.name,self.executionpath,self.login,self.port,[dirname+'.tar.gz'])
-
-	# }}}
-	def LaunchQueueJob(self,modelname,dirname,filelist,restart,batch):    # {{{
-
-		print('launching solution sequence on remote cluster')
-		if restart:
-			launchcommand='cd %s && cd %s chmod 777 %s.queue && ./%s.queue' % (self.executionpath,dirname,modelname,modelname)
-		else:
-			if batch:
-				launchcommand='cd %s && rm -rf ./%s && mkdir %s && cd %s && mv ../%s.tar.gz ./ && tar -zxf %s.tar.gz' % \
-						(self.executionpath,dirname,dirname,dirname,dirname,dirname)
-			else:
-				launchcommand='cd %s && rm -rf ./%s && mkdir %s && cd %s && mv ../%s.tar.gz ./ && tar -zxf %s.tar.gz  && chmod 777 %s.queue && ./%s.queue' % \
-					(self.executionpath,dirname,dirname,dirname,dirname,dirname,modelname,modelname)
-		issmssh(self.name,self.login,self.port,launchcommand)
-	# }}}
-	def Download(self,dirname,filelist):     # {{{
-
-		if m.ispc():
-			#do nothing
-			return
-
-		#copy files from cluster to current directory
-		directory='%s/%s/' % (self.executionpath,dirname)
-		issmscpin(self.name,self.login,self.port,directory,filelist)
-	# }}}
+    """
+    GENERIC cluster class definition
+
+       Usage:
+          cluster=generic('name','astrid','np',3);
+          cluster=generic('name',gethostname(),'np',3,'login','username');
+    """
+
+    def __init__(self,*args):    # {{{
+
+            self.name=''
+            self.login=''
+            self.np=1
+            self.port=0
+            self.interactive=1
+            self.codepath=IssmConfig('ISSM_PREFIX')[0]+'/bin'
+            self.executionpath=issmdir()+'/execution'
+            self.valgrind=issmdir()+'/externalpackages/valgrind/install/bin/valgrind'
+            self.valgrindlib=issmdir()+'/externalpackages/valgrind/install/lib/libmpidebug.so'
+            self.valgrindsup=issmdir()+'/externalpackages/valgrind/issm.supp'
+
+            #use provided options to change fields
+            options=pairoptions(*args)
+
+            #get name
+            self.name=socket.gethostname()
+
+            #initialize cluster using user settings if provided
+            if os.path.exists(self.name+'_settings.py'):
+                    exec(compile(open(self.name+'_settings.py').read(), self.name+'_settings.py', 'exec'),globals())
+
+            #OK get other fields
+            self=options.AssignObjectFields(self)
+    # }}}
+    def __repr__(self):    # {{{
+            #  display the object
+            s ="class '%s' object '%s' = \n" % (type(self),'self')
+            s+="    name: %s\n" % self.name
+            s+="    login: %s\n" % self.login
+            s+="    np: %i\n" % self.np
+            s+="    port: %i\n" % self.port
+            s+="    codepath: %s\n" % self.codepath
+            s+="    executionpath: %s\n" % self.executionpath
+            s+="    valgrind: %s\n" % self.valgrind
+            s+="    valgrindlib: %s\n" % self.valgrindlib
+            s+="    valgrindsup: %s\n" % self.valgrindsup
+            return s
+    # }}}
+    def checkconsistency(self,md,solution,analyses):    # {{{
+            if self.np<1:
+                    md = checkmessage(md,'number of processors should be at least 1')
+            if math.isnan(self.np):
+                    md = checkmessage(md,'number of processors should not be NaN!')
+
+            return md
+    # }}}
+    def BuildQueueScript(self,dirname,modelname,solution,io_gather,isvalgrind,isgprof,isdakota,isoceancoupling):    # {{{
+
+            executable='issm.exe';
+            if isdakota:
+                    version=IssmConfig('_DAKOTA_VERSION_')
+                    version=float(version[0])
+                    if version>=6:
+                            executable='issm_dakota.exe'
+            if isoceancoupling:
+                    executable='issm_ocean.exe'
+
+            #write queuing script
+            if not m.ispc():
+
+                    fid=open(modelname+'.queue','w')
+                    fid.write('#!/bin/sh\n')
+                    if not isvalgrind:
+                            if self.interactive:
+                                    if IssmConfig('_HAVE_MPI_')[0]:
+                                            fid.write('mpiexec -np %i %s/%s %s %s/%s %s ' % (self.np,self.codepath,executable,solution,self.executionpath,dirname,modelname))
+                                    else:
+                                            fid.write('%s/%s %s %s/%s %s ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname))
+                            else:
+                                    if IssmConfig('_HAVE_MPI_')[0]:
+                                            fid.write('mpiexec -np %i %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % (self.np,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
+                                    else:
+                                            fid.write('%s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
+                    elif isgprof:
+                            fid.write('\n gprof %s/%s gmon.out > %s.performance' % (self.codepath,executable,modelname))
+                    else:
+                            #Add --gen-suppressions=all to get suppression lines
+                            fid.write('LD_PRELOAD=%s \\\n' % self.valgrindlib)
+                            if IssmConfig('_HAVE_MPI_')[0]:
+                                    fid.write('mpiexec -np %i %s --leak-check=full --suppressions=%s %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % \
+                                                    (self.np,self.valgrind,self.valgrindsup,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
+                            else:
+                                    fid.write('%s --leak-check=full --suppressions=%s %s/%s %s %s/%s %s 2> %s.errlog >%s.outlog ' % \
+                                                    (self.valgrind,self.valgrindsup,self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
+
+                    if not io_gather:    #concatenate the output files:
+                            fid.write('\ncat %s.outbin.* > %s.outbin' % (modelname,modelname))
+                    fid.close()
+
+            else:    # Windows
+
+                    fid=open(modelname+'.bat','w')
+                    fid.write('@echo off\n')
+                    if self.interactive:
+                            fid.write('"%s/%s" %s "%s/%s" %s ' % (self.codepath,executable,solution,self.executionpath,dirname,modelname))
+                    else:
+                            fid.write('"%s/%s" %s "%s/%s" %s 2> %s.errlog >%s.outlog' % \
+                                    (self.codepath,executable,solution,self.executionpath,dirname,modelname,modelname,modelname))
+                    fid.close()
+
+            #in interactive mode, create a run file, and errlog and outlog file
+            if self.interactive:
+                    fid=open(modelname+'.errlog','w')
+                    fid.close()
+                    fid=open(modelname+'.outlog','w')
+                    fid.close()
+    # }}}
+    def BuildKrigingQueueScript(self,modelname,solution,io_gather,isvalgrind,isgprof):    # {{{
+
+            #write queuing script
+            if not m.ispc():
+
+                    fid=open(modelname+'.queue','w')
+                    fid.write('#!/bin/sh\n')
+                    if not isvalgrind:
+                            if self.interactive:
+                                    fid.write('mpiexec -np %i %s/kriging.exe %s/%s %s ' % (self.np,self.codepath,self.executionpath,modelname,modelname))
+                            else:
+                                    fid.write('mpiexec -np %i %s/kriging.exe %s/%s %s 2> %s.errlog >%s.outlog ' % (self.np,self.codepath,self.executionpath,modelname,modelname,modelname,modelname))
+                    elif isgprof:
+                            fid.write('\n gprof %s/kriging.exe gmon.out > %s.performance' & (self.codepath,modelname))
+                    else:
+                            #Add --gen-suppressions=all to get suppression lines
+                            fid.write('LD_PRELOAD=%s \\\n' % self.valgrindlib)
+                            fid.write('mpiexec -np %i %s --leak-check=full --suppressions=%s %s/kriging.exe %s/%s %s 2> %s.errlog >%s.outlog ' % \
+                                    (self.np,self.valgrind,self.valgrindsup,self.codepath,self.executionpath,modelname,modelname,modelname,modelname))
+                    if not io_gather:    #concatenate the output files:
+                            fid.write('\ncat %s.outbin.* > %s.outbin' % (modelname,modelname))
+                    fid.close()
+
+            else:    # Windows
+
+                    fid=open(modelname+'.bat','w')
+                    fid.write('@echo off\n')
+                    if self.interactive:
+                            fid.write('"%s/issm.exe" %s "%s/%s" %s ' % (self.codepath,solution,self.executionpath,modelname,modelname))
+                    else:
+                            fid.write('"%s/issm.exe" %s "%s/%s" %s 2> %s.errlog >%s.outlog' % \
+                                    (self.codepath,solution,self.executionpath,modelname,modelname,modelname,modelname))
+                    fid.close()
+
+            #in interactive mode, create a run file, and errlog and outlog file
+            if self.interactive:
+                    fid=open(modelname+'.errlog','w')
+                    fid.close()
+                    fid=open(modelname+'.outlog','w')
+                    fid.close()
+    # }}}
+    def UploadQueueJob(self,modelname,dirname,filelist):    # {{{
+
+        #compress the files into one zip.
+        compressstring='tar -zcf %s.tar.gz ' % dirname
+        for file in filelist:
+            compressstring += ' %s' % file
+        if self.interactive:
+            compressstring += ' %s.errlog %s.outlog ' % (modelname,modelname)
+        subprocess.call(compressstring,shell=True)
+
+        print('uploading input file and queueing script')
+        issmscpout(self.name,self.executionpath,self.login,self.port,[dirname+'.tar.gz'])
+
+    # }}}
+    def LaunchQueueJob(self,modelname,dirname,filelist,restart,batch):    # {{{
+
+        print('launching solution sequence on remote cluster')
+        if restart:
+            launchcommand='cd %s && cd %s chmod 777 %s.queue && ./%s.queue' % (self.executionpath,dirname,modelname,modelname)
+        else:
+            if batch:
+                launchcommand='cd %s && rm -rf ./%s && mkdir %s && cd %s && mv ../%s.tar.gz ./ && tar -zxf %s.tar.gz' % \
+                               (self.executionpath,dirname,dirname,dirname,dirname,dirname)
+            else:
+                launchcommand='cd %s && rm -rf ./%s && mkdir %s && cd %s && mv ../%s.tar.gz ./ && tar -zxf %s.tar.gz  && chmod 777 %s.queue && ./%s.queue' % \
+                               (self.executionpath,dirname,dirname,dirname,dirname,dirname,modelname,modelname)
+        issmssh(self.name,self.login,self.port,launchcommand)
+    # }}}
+    def Download(self,dirname,filelist):     # {{{
+
+        if m.ispc():
+            #do nothing
+            return
+
+        #copy files from cluster to current directory
+        directory='%s/%s/' % (self.executionpath,dirname)
+        issmscpin(self.name,self.login,self.port,directory,filelist)
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/levelset.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/levelset.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/classes/levelset.py	(revision 24115)
@@ -4,76 +4,80 @@
 from WriteData import WriteData
 
+
 class levelset(object):
-	"""
-	LEVELSET class definition
+    """
+    LEVELSET class definition
 
-	   Usage:
-	      levelset=levelset();
-	"""
+       Usage:
+          levelset=levelset();
+    """
 
-	def __init__(self): # {{{
+    def __init__(self):  # {{{
 
-		self.stabilization    = 0
-		self.spclevelset      = float('NaN')
-		self.reinit_frequency = 0
-                self.kill_icebergs    = 0
-		self.calving_max      = 0.
-		self.fe               = 'P1'
+        self.stabilization = 0
+        self.spclevelset = float('NaN')
+        self.reinit_frequency = 0
+        self.kill_icebergs = 0
+        self.calving_max = 0.
+        self.fe = 'P1'
 
-		#set defaults
-		self.setdefaultparameters()
+        #set defaults
+        self.setdefaultparameters()
 
-		#}}}
-	def __repr__(self): # {{{
-		string='   Level-set parameters:'
-		string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: streamline upwinding'))
-		string="%s\n%s"%(string,fielddisplay(self,'spclevelset','levelset constraints (NaN means no constraint)'))
-		string="%s\n%s"%(string,fielddisplay(self,'reinit_frequency','Amount of time steps after which the levelset function in re-initialized'))
-                string="%s\n%s"%(string,fielddisplay(self,'kill_icebergs','remove floating icebergs to prevent rigid body motions (1: true, 0: false)'))
-		string="%s\n%s"%(string,fielddisplay(self,'calving_max','maximum allowed calving rate (m/a)'))
-		string="%s\n%s"%(string,fielddisplay(self,'fe','Finite Element type: ''P1'' (default), or ''P2'''))
+    #}}}
+    def __repr__(self):  # {{{
+        string = '   Level-set parameters:'
+        string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', '0: no, 1: artificial_diffusivity, 2: streamline upwinding'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'spclevelset', 'levelset constraints (NaN means no constraint)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'reinit_frequency', 'Amount of time steps after which the levelset function in re-initialized'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'kill_icebergs', 'remove floating icebergs to prevent rigid body motions (1: true, 0: false)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'calving_max', 'maximum allowed calving rate (m/a)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'fe', 'Finite Element type: ''P1'' (default), or ''P2'''))
 
-		return string
-		#}}}
-	def extrude(self,md): # {{{
-		self.spclevelset=project3d(md,'vector',self.spclevelset,'type','node')
-		return self
-	#}}}
-	def setdefaultparameters(self): # {{{
+        return string
+    #}}}
 
-		#stabilization = 1 by default
-		self.stabilization		= 1
-		self.reinit_frequency = 5
-                self.kill_icebergs    = 1
-		self.calving_max      = 3000.
+    def extrude(self, md):  # {{{
+        self.spclevelset = project3d(md, 'vector', self.spclevelset, 'type', 'node')
+        return self
+    #}}}
 
-		#Linear elements by default
-		self.fe='P1'
+    def setdefaultparameters(self):  # {{{
 
-		return self
-	#}}}
-	def checkconsistency(self,md,solution,analyses):    # {{{
+        #stabilization = 1 by default
+        self.stabilization = 1
+        self.reinit_frequency = 5
+        self.kill_icebergs = 1
+        self.calving_max = 3000.
 
-		#Early return
-		if (solution!='TransientSolution') or (not md.transient.ismovingfront):
-			return md
+        #Linear elements by default
+        self.fe = 'P1'
 
-		md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1)
-		md = checkfield(md,'fieldname','levelset.stabilization','numel',[1],'values',[0,1,2]);
-                md = checkfield(md,'fieldname','levelset.kill_icebergs','numel',[1],'values',[0,1]);
-		md = checkfield(md,'fieldname','levelset.calving_max','numel',[1],'NaN',1,'Inf',1,'>',0);
-		md = checkfield(md,'fieldname','levelset.fe','values',['P1','P2']);
+        return self
+    #}}}
+    def checkconsistency(self, md, solution, analyses):    # {{{
 
-		return md
-	# }}}
-	def marshall(self,prefix,md,fid):    # {{{
+        #Early return
+        if (solution != 'TransientSolution') or (not md.transient.ismovingfront):
+            return md
 
-		yts=md.constants.yts;
+        md = checkfield(md, 'fieldname', 'levelset.spclevelset', 'Inf', 1, 'timeseries', 1)
+        md = checkfield(md, 'fieldname', 'levelset.stabilization', 'numel', [1], 'values', [0, 1, 2])
+        md = checkfield(md, 'fieldname', 'levelset.kill_icebergs', 'numel', [1], 'values', [0, 1])
+        md = checkfield(md, 'fieldname', 'levelset.calving_max', 'numel', [1], 'NaN', 1, 'Inf', 1, '>', 0)
+        md = checkfield(md, 'fieldname', 'levelset.fe', 'values', ['P1', 'P2'])
 
-		WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer');
-		WriteData(fid,prefix,'object',self,'fieldname','spclevelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
-		WriteData(fid,prefix,'object',self,'fieldname','reinit_frequency','format','Integer');
-                WriteData(fid,prefix,'object',self,'fieldname','kill_icebergs','format','Boolean');
-		WriteData(fid,prefix,'object',self,'fieldname','calving_max','format','Double','scale',1./yts);
-		WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
-	# }}}
+        return md
+    # }}}
+
+    def marshall(self, prefix, md, fid):    # {{{
+
+        yts = md.constants.yts
+
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'spclevelset', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'reinit_frequency', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'kill_icebergs', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'calving_max', 'format', 'Double', 'scale', 1. / yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'fe', 'format', 'String')
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/toolkits.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/toolkits.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/classes/toolkits.py	(revision 24115)
@@ -3,107 +3,110 @@
 from iluasmoptions import iluasmoptions
 from fielddisplay import fielddisplay
-from checkfield import checkfield
 from issmgslsolver import issmgslsolver
 from issmmumpssolver import issmmumpssolver
 
+
 class toolkits(object):
-	"""
-	TOOLKITS class definition
+    """
+    TOOLKITS class definition
 
-	   Usage:
-	      self=toolkits();
-	"""
+       Usage:
+          self=toolkits();
+    """
 
-	def __init__(self):    # {{{
-		#default toolkits
-		if IssmConfig('_HAVE_PETSC_')[0]:
-			#MUMPS is the default toolkits
-			if IssmConfig('_HAVE_MUMPS_')[0]:
-				self.DefaultAnalysis           = mumpsoptions()
-			else:
-				self.DefaultAnalysis           = iluasmoptions()
-		else:
-			if IssmConfig('_HAVE_MUMPS_')[0]:
-				self.DefaultAnalysis           = issmmumpssolver()
-			elif IssmConfig('_HAVE_GSL_')[0]:
-				self.DefaultAnalysis           = issmgslsolver()
-			else:
-				raise IOError("ToolkitsFile error: need at least Mumps or Gsl to define issm solver type")
+    def __init__(self):    # {{{
+        #default toolkits
+        if IssmConfig('_HAVE_PETSC_')[0]:
+            #MUMPS is the default toolkits
+            if IssmConfig('_HAVE_MUMPS_')[0]:
+                self.DefaultAnalysis = mumpsoptions()
+            else:
+                self.DefaultAnalysis = iluasmoptions()
+        else:
+            if IssmConfig('_HAVE_MUMPS_')[0]:
+                self.DefaultAnalysis = issmmumpssolver()
+            elif IssmConfig('_HAVE_GSL_')[0]:
+                self.DefaultAnalysis = issmgslsolver()
+            else:
+                raise IOError("ToolkitsFile error: need at least Mumps or Gsl to define issm solver type")
 
-                #Use same solver for Recovery mode
-                self.RecoveryAnalysis = self.DefaultAnalysis
+        #Use same solver for Recovery mode
+        self.RecoveryAnalysis = self.DefaultAnalysis
 
-                #The other properties are dynamic
-	# }}}
-	def __repr__(self):    # {{{
-		s ="List of toolkits options per analysis:\n\n"
-		for analysis in list(vars(self).keys()):
-			s+="%s\n" % fielddisplay(self,analysis,'')
+        #The other properties are dynamic
+    # }}}
+    def __repr__(self):    # {{{
+        s = "List of toolkits options per analysis:\n\n"
+        for analysis in list(vars(self).keys()):
+            s += "%s\n" % fielddisplay(self, analysis, '')
 
-		return s
-	# }}}
-	def addoptions(self,analysis,*args):    # {{{
-		# Usage example:
-		#    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
-		#    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
+            return s
+    # }}}
 
-		#Create dynamic property if property does not exist yet
-		if not hasattr(self,analysis):
-			setattr(self,analysis,None)
+    def addoptions(self, analysis, *args):    # {{{
+        # Usage example:
+        #    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
+        #    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
 
-		#Add toolkits options to analysis
-		if len(args)==1:
-			setattr(self,analysis,args[0])
+        #Create dynamic property if property does not exist yet
+        if not hasattr(self, analysis):
+            setattr(self, analysis, None)
 
-		return self
-	# }}}
-	def checkconsistency(self,md,solution,analyses):    # {{{
-		for analysis in list(vars(self).keys()):
-			if not getattr(self,analysis):
-				md.checkmessage("md.toolkits.%s is empty" % analysis)
+        #Add toolkits options to analysis
+        if len(args) == 1:
+            setattr(self, analysis, args[0])
 
-		return md
-	# }}}
-	def ToolkitsFile(self,filename):    # {{{
-		"""
-		TOOLKITSFILE- build toolkits file
+        return self
+    # }}}
 
-		   Build a Petsc compatible options file, from the toolkits model field  + return options string
-		   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
+    def checkconsistency(self, md, solution, analyses):    # {{{
+        for analysis in list(vars(self).keys()):
+            if not getattr(self, analysis):
+                md.checkmessage("md.toolkits.%s is empty" % analysis)
+
+        return md
+    # }}}
+
+    def ToolkitsFile(self, filename):    # {{{
+        """
+        TOOLKITSFILE- build toolkits file
+
+           Build a Petsc compatible options file, from the toolkits model field  + return options string
+           This file will also be used when the toolkit used is 'issm' instead of 'petsc'
 
 
-		   Usage:     ToolkitsFile(toolkits,filename);
-		"""
+           Usage:     ToolkitsFile(toolkits,filename);
+        """
 
-		#open file for writing
-		try:
-			fid=open(filename,'w')
-		except IOError as e:
-			raise IOError("ToolkitsFile error: could not open '%s' for writing." % filename)
+        #open file for writing
+        try:
+            fid = open(filename, 'w')
+        except IOError as e:
+            raise IOError("ToolkitsFile error: could not open '%s' for writing." % filename)
 
-		#write header
-		fid.write("%s%s%s\n" % ('%Toolkits options file: ',filename,' written from Python toolkits array'))
+        #write header
+        fid.write("%s%s%s\n" % ('%Toolkits options file: ', filename, ' written from Python toolkits array'))
 
-		#start writing options
-		for analysis in list(vars(self).keys()):
-			options=getattr(self,analysis)
+        #start writing options
+        for analysis in list(vars(self).keys()):
+            options = getattr(self, analysis)
 
-			#first write analysis:
-			fid.write("\n+%s\n" % analysis)    #append a + to recognize it's an analysis enum
-			#now, write options
-			for optionname,optionvalue in list(options.items()):
+            #first write analysis:
+            fid.write("\n+%s\n" % analysis)    #append a + to recognize it's an analysis enum
+            #now, write options
+            for optionname, optionvalue in list(options.items()):
 
-				if not optionvalue:
-					#this option has only one argument
-					fid.write("-%s\n" % optionname)
-				else:
-					#option with value. value can be string or scalar
-					if   isinstance(optionvalue,(bool,int,float)):
-						fid.write("-%s %g\n" % (optionname,optionvalue))
-					elif isinstance(optionvalue,str):
-						fid.write("-%s %s\n" % (optionname,optionvalue))
-					else:
-						raise TypeError("ToolkitsFile error: option '%s' is not well formatted." % optionname)
+                if not optionvalue:
+                    #this option has only one argument
+                    fid.write("-%s\n" % optionname)
+                else:
+                    #option with value. value can be string or scalar
+                    if isinstance(optionvalue, (bool, int, float)):
+                        fid.write("-%s %g\n" % (optionname, optionvalue))
+                    elif isinstance(optionvalue, str):
+                        fid.write("-%s %s\n" % (optionname, optionvalue))
+                    else:
+                        raise TypeError("ToolkitsFile error: option '%s' is not well formatted." % optionname)
 
-		fid.close()
-	# }}}
+        fid.close()
+    # }}}
Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py
===================================================================
--- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py	(revision 24115)
@@ -21,5 +21,4 @@
 
     TODO: - make time easily accessible
-          - make evolving geometry
 
     Basile de Fleurian:
@@ -191,9 +190,9 @@
                             Vxstruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'x'])
                             Vystruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'y'])
+                            treated_res += [field[:-1] + 'x', field[:-1] + 'y']
                             if dim == 3:
                                 Vzstruct = np.squeeze(spe_res_struct.__dict__[field[:-1] + 'z'])
-                                treated_res += [field[:-1] + 'x', field[:-1] + 'y', field[:-1] + 'z']
-                            elif dim == 2:
-                                treated_res += [field[:-1] + 'x', field[:-1] + 'y']
+                                treated_res += field[:-1] + 'z'
+
                         except KeyError:
                             fieldnames += field
@@ -211,5 +210,35 @@
 
                     elif field in tensors:
-                        print("nothing")
+                        try:
+                            Txxstruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'xx'])
+                            Txystruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'xy'])
+                            Tyystruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'yy'])
+                            treated_res += [field[:-2] + 'xx', field[:-2] + 'xy', field[:-2] + 'yy']
+                            if dim == 3:
+                                Tzzstruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'zz'])
+                                Txzstruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'xz'])
+                                Tyzstruct = np.squeeze(spe_res_struct.__dict__[field[:-2] + 'yz'])
+                                treated_res += [field[:-2] + 'zz', field[:-2] + 'xz', field[:-2] + 'yz']
+
+                        except KeyError:
+                            fieldnames += field
+                            tensors.remove(field)
+
+                        fid.write('TENSORS {} float \n'.format(field[:-2]))
+                        for node in range(0, num_of_points):
+                            Txx = cleanOutliers(Txxstruct[enveloppe_index[node]])
+                            Tyy = cleanOutliers(Tyystruct[enveloppe_index[node]])
+                            Txy = cleanOutliers(Txystruct[enveloppe_index[node]])
+                            if dim == 3:
+                                Tzz = cleanOutliers(Tzzstruct[enveloppe_index[node]])
+                                Txz = cleanOutliers(Txzstruct[enveloppe_index[node]])
+                                Tyz = cleanOutliers(Tyzstruct[enveloppe_index[node]])
+                                fid.write('{:f} {:f} {:f}\n'.format(Txx, Txy, Txz))
+                                fid.write('{:f} {:f} {:f}\n'.format(Txy, Tyy, Tyz))
+                                fid.write('{:f} {:f} {:f}\n'.format(Txz, Tyz, Tzz))
+                            elif dim == 2:
+                                fid.write('{:f} {:f} {:f}\n'.format(Txx, Txy, 0))
+                                fid.write('{:f} {:f} {:f}\n'.format(Txy, Tyy, 0))
+                                fid.write('{:f} {:f} {:f}\n'.format(0, 0, 0))
                     else:
                         if ((np.size(spe_res_struct.__dict__[field])) == every_nodes):
Index: /issm/trunk-jpl/src/m/geometry/slope.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/slope.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/geometry/slope.py	(revision 24115)
@@ -1,46 +1,46 @@
-import numpy as  np
-from GetNodalFunctionsCoeff import  GetNodalFunctionsCoeff
+import numpy as np
+from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
+from project3d import project3d
 
-def slope(md,*args):
-	"""
-	SLOPE - compute the surface slope
 
-	Usage:
-		sx,sy,s=slope(md)
-		sx,sy,s=slope(md,md.results.TransientSolution(1).Surface)
-	"""
+def slope(md, *args):
+    """
+    SLOPE - compute the surface slope
 
-	#load some variables (it is much faster if the variables are loaded from md once for all) 
-	if md.mesh.dimension()==2:
-		numberofelements=md.mesh.numberofelements
-		numberofnodes=md.mesh.numberofvertices
-		index=md.mesh.elements
-		x=md.mesh.x ; y=md.mesh.y
-	else:
-		numberofelements=md.mesh.numberofelements2d
-		numberofnodes=md.mesh.numberofvertices2d
-		index=md.mesh.elements2d
-		x=md.mesh.x2d; y=md.mesh.y2d
+    Usage:
+            sx,sy,s=slope(md)
+            sx,sy,s=slope(md,md.results.TransientSolution(1).Surface)
+    """
 
-	if len(args)==0:
-		surf=md.geometry.surface
-	elif len(args)==1:
-		surf=args[0]
-	else:
-		raise RuntimeError("slope.py usage error")
+    #load some variables (it is much faster if the variables are loaded from md once for all)
+    if md.mesh.dimension() == 2:
+        index = md.mesh.elements
+        x = md.mesh.x
+        y = md.mesh.y
+    else:
+        index = md.mesh.elements2d
+        x = md.mesh.x2d
+        y = md.mesh.y2d
 
-	#%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-	alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2]
+    if len(args) == 0:
+        surf = md.geometry.surface
+    elif len(args) == 1:
+        surf = args[0]
+    else:
+        raise RuntimeError("slope.py usage error")
 
-	summation=np.array([[1],[1],[1]])
-	sx=np.dot(surf[index-1,0]*alpha,summation).reshape(-1,)
-	sy=np.dot(surf[index-1,0]*beta,summation).reshape(-1,)
+    #%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
+    alpha, beta = GetNodalFunctionsCoeff(index, x, y)[0:2]
 
-	s=np.sqrt(sx**2+sy**2)
+    summation = np.array([[1], [1], [1]])
+    sx = np.dot(surf[index - 1, 0] * alpha, summation).reshape(-1,)
+    sy = np.dot(surf[index - 1, 0] * beta, summation).reshape(-1,)
 
-	if md.mesh.dimension()==3:
-		sx=project3d(md,'vector',sx,'type','element')
-		sy=project3d(md,'vector',sy,'type','element')
-		s=np.sqrt(sx**2+sy**2)
+    s = np.sqrt(sx**2 + sy**2)
 
-	return (sx,sy,s)
+    if md.mesh.dimension() == 3:
+        sx = project3d(md, 'vector', sx, 'type', 'element')
+        sy = project3d(md, 'vector', sy, 'type', 'element')
+        s = np.sqrt(sx**2 + sy**2)
+
+    return (sx, sy, s)
Index: /issm/trunk-jpl/src/m/interp/averaging.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/averaging.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/interp/averaging.py	(revision 24115)
@@ -1,96 +1,95 @@
-import numpy as  np
+import numpy as np
 from GetAreas import GetAreas
-import MatlabFuncs as m
 try:
-	from scipy.sparse import csc_matrix
+    from scipy.sparse import csc_matrix
 except ImportError:
-	print("could not import scipy, no averaging capabilities enabled")
+    print("could not import scipy, no averaging capabilities enabled")
 
-def averaging(md,data,iterations,layer=0):
-	'''
-	AVERAGING - smooths the input over the mesh
-	
-	   This routine takes a list over the elements or the nodes in input
-	   and return a list over the nodes.
-	   For each iterations it computes the average over each element (average 
-	   of the vertices values) and then computes the average over each node
-	   by taking the average of the element around a node weighted by the
-	   elements volume
-	   For 3d mesh, a last argument can be added to specify the layer to be averaged on.
-	
-	   Usage:
-	      smoothdata=averaging(md,data,iterations)
-	      smoothdata=averaging(md,data,iterations,layer)
-	
-	   Examples:
-	      velsmoothed=averaging(md,md.initialization.vel,4)
-	      pressure=averaging(md,md.initialization.pressure,0)
-	      temperature=averaging(md,md.initialization.temperature,1,1)
-	'''
 
-	if len(data)!=md.mesh.numberofelements and len(data)!=md.mesh.numberofvertices:
-		raise Exception('averaging error message: data not supported yet')
-	if md.mesh.dimension()==3 and layer!=0:
-		if layer<=0 or layer>md.mesh.numberoflayers:
-			raise ValueError('layer should be between 1 and md.mesh.numberoflayers')
-	else:
-		layer=0
-	
-	#initialization
-	if layer==0:
-		weights=np.zeros(md.mesh.numberofvertices,)
-		data=data.flatten(1)
-	else:
-		weights=np.zeros(md.mesh.numberofvertices2d,)
-		data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:]
-	
-	#load some variables (it is much faster if the variabes are loaded from md once for all)
-	if layer==0:
-		index=md.mesh.elements
-		numberofnodes=md.mesh.numberofvertices
-		numberofelements=md.mesh.numberofelements
-	else:
-		index=md.mesh.elements2d
-		numberofnodes=md.mesh.numberofvertices2d
-		numberofelements=md.mesh.numberofelements2d
+def averaging(md, data, iterations, layer=0):
+    '''
+    AVERAGING - smooths the input over the mesh
 
-	
-	#build some variables
-	if md.mesh.dimension()==3 and layer==0:
-		rep=6
-		areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z)
-	elif md.mesh.dimension()==2:
-		rep=3
-		areas=GetAreas(index,md.mesh.x,md.mesh.y)
-	else:
-		rep=3
-		areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d)
+       This routine takes a list over the elements or the nodes in input
+       and return a list over the nodes.
+       For each iterations it computes the average over each element (average
+       of the vertices values) and then computes the average over each node
+       by taking the average of the element around a node weighted by the
+       elements volume
+       For 3d mesh,  a last argument can be added to specify the layer to be averaged on.
 
-	index=index-1 # since python indexes starting from zero
-	line=index.flatten(1)
-	areas=np.vstack(areas).reshape(-1,)
-	summation=1./rep*np.ones(rep,)
-	linesize=rep*numberofelements
-	
-	#update weights that holds the volume of all the element holding the node i
-	weights=csc_matrix( (np.tile(areas,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
-	
-	#initialization
-	if len(data)==numberofelements:
-		average_node=csc_matrix( (np.tile(areas*data,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
-		average_node=average_node/weights
-		average_node = csc_matrix(average_node)
-	else:
-		average_node=csc_matrix(data.reshape(-1,1))
+       Usage:
+          smoothdata=averaging(md, data, iterations)
+          smoothdata=averaging(md, data, iterations, layer)
 
-	#loop over iteration
-	for i in np.arange(1,iterations+1):
-		average_el=np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements,rep),np.vstack(summation))).reshape(-1,)
-		average_node=csc_matrix( (np.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
-		average_node=average_node/weights
-		average_node=csc_matrix(average_node)
-	
-	#return output as a full matrix (C code do not like sparse matrices)
-	average=np.asarray(average_node.todense()).reshape(-1,)
+       Examples:
+          velsmoothed=averaging(md, md.initialization.vel, 4)
+          pressure=averaging(md, md.initialization.pressure, 0)
+          temperature=averaging(md, md.initialization.temperature, 1, 1)
+    '''
 
-	return average
+    if len(data) != md.mesh.numberofelements and len(data) != md.mesh.numberofvertices:
+        raise Exception('averaging error message: data not supported yet')
+    if md.mesh.dimension() == 3 and layer != 0:
+        if layer <= 0 or layer > md.mesh.numberoflayers:
+            raise ValueError('layer should be between 1 and md.mesh.numberoflayers')
+    else:
+        layer = 0
+
+    #initialization
+    if layer == 0:
+        weights = np.zeros(md.mesh.numberofvertices,)
+        data = data.flatten(1)
+    else:
+        weights = np.zeros(md.mesh.numberofvertices2d,)
+        data = data[(layer - 1) * md.mesh.numberofvertices2d + 1:layer * md.mesh.numberofvertices2d, :]
+
+    #load some variables (it is much faster if the variabes are loaded from md once for all)
+    if layer == 0:
+        index = md.mesh.elements
+        numberofnodes = md.mesh.numberofvertices
+        numberofelements = md.mesh.numberofelements
+    else:
+        index = md.mesh.elements2d
+        numberofnodes = md.mesh.numberofvertices2d
+        numberofelements = md.mesh.numberofelements2d
+
+    #build some variables
+    if md.mesh.dimension() == 3 and layer == 0:
+        rep = 6
+        areas = GetAreas(index, md.mesh.x, md.mesh.y, md.mesh.z)
+    elif md.mesh.dimension() == 2:
+        rep = 3
+        areas = GetAreas(index, md.mesh.x, md.mesh.y)
+    else:
+        rep = 3
+        areas = GetAreas(index, md.mesh.x2d, md.mesh.y2d)
+
+    index = index - 1  # since python indexes starting from zero
+    line = index.flatten(1)
+    areas = np.vstack(areas).reshape(-1,)
+    summation = 1. / rep * np.ones(rep,)
+    linesize = rep * numberofelements
+
+    #update weights that holds the volume of all the element holding the node i
+    weights = csc_matrix((np.tile(areas, (rep, 1)).reshape(-1,), (line, np.zeros(linesize,))), shape=(numberofnodes, 1))
+
+    #initialization
+    if len(data) == numberofelements:
+        average_node = csc_matrix((np.tile(areas * data, (rep, 1)).reshape(-1,), (line, np.zeros(linesize,))), shape=(numberofnodes, 1))
+        average_node = average_node / weights
+        average_node = csc_matrix(average_node)
+    else:
+        average_node = csc_matrix(data.reshape(-1, 1))
+
+    #loop over iteration
+    for i in np.arange(1, iterations + 1):
+        average_el = np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements, rep), np.vstack(summation))).reshape(-1,)
+        average_node = csc_matrix((np.tile(areas * average_el.reshape(-1), (rep, 1)).reshape(-1,), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
+        average_node = average_node / weights
+        average_node = csc_matrix(average_node)
+
+    #return output as a full matrix (C code do not like sparse matrices)
+    average = np.asarray(average_node.todense()).reshape(-1,)
+
+    return average
Index: /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py	(revision 24115)
@@ -1,58 +1,58 @@
 import numpy as np
 
-def GetNodalFunctionsCoeff(index,x,y):
-	"""
-	GETNODELFUNCTIONSCOEFF - compute nodal functions coefficients
 
-	   Compute the coefficients alpha beta and optionaly gamma of
-	   2d triangular elements. For each element, the nodal function
-	   is defined as:
-	   N(x,y)=sum(i=1:3) alpha_i * x + beta_i * y + gamma_i
+def GetNodalFunctionsCoeff(index, x, y):
+    """
+    GETNODELFUNCTIONSCOEFF - compute nodal functions coefficients
 
-	   Usage:
-	      [alpha beta]=GetNodalFunctionsCoeff(index,x,y);
-	      [alpha beta gamma]=GetNodalFunctionsCoeff(index,x,y);
+       Compute the coefficients alpha beta and optionaly gamma of
+       2d triangular elements. For each element,  the nodal function
+       is defined as:
+       N(x, y)=sum(i=1:3) alpha_i * x + beta_i * y + gamma_i
 
-	   Example:
-	      [alpha beta gamma]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y);
-	"""
+       Usage:
+          [alpha beta]=GetNodalFunctionsCoeff(index, x, y);
+          [alpha beta gamma]=GetNodalFunctionsCoeff(index, x, y);
 
-	#make columns out of x and y
-	x=x.reshape(-1)
-	y=y.reshape(-1)
+       Example:
+          [alpha beta gamma]=GetNodalFunctionsCoeff(md.mesh.elements, md.mesh.x, md.mesh.y);
+    """
 
-	#get nels and nods
-	nels=np.size(index,axis=0)
-	nods=np.size(x)
+    #make columns out of x and y
+    x = x.reshape(-1)
+    y = y.reshape(-1)
 
-	#some checks
-	if np.size(y)!=nods:
-		raise TypeError("GetNodalFunctionsCoeff error message: x and y do not have the same length.")
-	if np.max(index)>nods:
-		raise TypeError("GetNodalFunctionsCoeff error message: index should not have values above %d." % nods)
-	if np.size(index,axis=1)!=3:
-		raise TypeError("GetNodalFunctionsCoeff error message: only 2d meshes supported. index should have 3 columns.")
+    #get nels and nods
+    nels = np.size(index, axis=0)
+    nods = np.size(x)
 
-	#initialize output
-	alpha=np.zeros((nels,3))
-	beta=np.zeros((nels,3))
+    #some checks
+    if np.size(y) != nods:
+        raise TypeError("GetNodalFunctionsCoeff error message: x and y do not have the same length.")
+    if np.max(index) > nods:
+        raise TypeError("GetNodalFunctionsCoeff error message: index should not have values above {}.".format(nods))
+    if np.size(index, axis=1) != 3:
+        raise TypeError("GetNodalFunctionsCoeff error message: only 2d meshes supported. index should have 3 columns.")
 
-	#compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
-	x1=x[index[:,0]-1]
-	x2=x[index[:,1]-1]
-	x3=x[index[:,2]-1]
-	y1=y[index[:,0]-1]
-	y2=y[index[:,1]-1]
-	y3=y[index[:,2]-1]
-	invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
+    #initialize output
+    alpha = np.zeros((nels, 3))
+    beta = np.zeros((nels, 3))
 
-	#get alpha and beta
-	alpha=np.vstack(((invdet*(y2-y3)).reshape(-1,),(invdet*(y3-y1)).reshape(-1,),(invdet*(y1-y2)).reshape(-1,))).T
-	beta =np.vstack(((invdet*(x3-x2)).reshape(-1,),(invdet*(x1-x3)).reshape(-1,),(invdet*(x2-x1)).reshape(-1,))).T
+    #compute nodal functions coefficients N(x, y) = alpha x + beta y +gamma
+    x1 = x[index[:, 0] - 1]
+    x2 = x[index[:, 1] - 1]
+    x3 = x[index[:, 2] - 1]
+    y1 = y[index[:, 0] - 1]
+    y2 = y[index[:, 1] - 1]
+    y3 = y[index[:, 2] - 1]
+    invdet = 1. / (x1 * (y2 - y3) - x2 * (y1 - y3) + x3 * (y1 - y2))
 
-	#get gamma if requested
-	gamma=np.zeros((nels,3))
-	gamma=np.vstack(((invdet*(x2*y3-x3*y2)).reshape(-1,),(invdet*(y1*x3-y3*x1)).reshape(-1,),(invdet*(x1*y2-x2*y1)).reshape(-1,))).T
+    #get alpha and beta
+    alpha = np.vstack(((invdet * (y2 - y3)).reshape(-1,), (invdet * (y3 - y1)).reshape(-1,), (invdet * (y1 - y2)).reshape(-1,))).T
+    beta = np.vstack(((invdet * (x3 - x2)).reshape(-1,), (invdet * (x1 - x3)).reshape(-1,), (invdet * (x2 - x1)).reshape(-1,))).T
 
-	return alpha,beta,gamma
+    #get gamma if requested
+    gamma = np.zeros((nels, 3))
+    gamma = np.vstack(((invdet * (x2 * y3 - x3 * y2)).reshape(-1,), (invdet * (y1 * x3 - y3 * x1)).reshape(-1,), (invdet * (x1 * y2 - x2 * y1)).reshape(-1,))).T
 
+    return alpha, beta, gamma
Index: /issm/trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 24115)
@@ -1,4 +1,4 @@
 import os.path
-import numpy as  np
+import numpy as np
 from mesh2d import *
 from mesh2dvertical import *
@@ -15,610 +15,605 @@
 from ContourToNodes import ContourToNodes
 
-def bamg(md,*args):
-	"""
-	BAMG - mesh generation
-
-	   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
-
-	   - domain :            followed by an ARGUS file that prescribes the domain outline
-	   - holes :             followed by an ARGUS file that prescribes the holes
-     - subdomains :        followed by an ARGUS file that prescribes the list of
-	                         subdomains (that need to be inside domain)
-
-	   - hmin :              minimum edge length (default is 10^-100)
-	   - hmax :              maximum edge length (default is 10^100)
-	   - hVertices :         imposed edge length for each vertex (geometry or mesh)
-	   - hminVertices :      minimum edge length for each vertex (mesh)
-	   - hmaxVertices :      maximum edge length for each vertex (mesh)
-
-	   - anisomax :          maximum ratio between the smallest and largest edges (default is 10^30)
-	   - coeff :             coefficient applied to the metric (2-> twice as many elements, default is 1)
-	   - cutoff :            scalar used to compute the metric when metric type 2 or 3 are applied
-	   - err :               error used to generate the metric from a field
-	   - errg :              geometric error (default is 0.1)
-	   - field :             field of the model that will be used to compute the metric
-	                         to apply several fields, use one column per field
-	   - gradation :         maximum ratio between two adjacent edges
-	   - Hessiantype :       0 -> use double P2 projection (default)
-	                         1 -> use Green formula
-	   - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
-	   - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
-	   - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
-	   - metric :            matrix (numberofnodes x 3) used as a metric
-	   - Metrictype :        1 -> absolute error          c/(err coeff^2) * Abs(H)        (default)
-	                         2 -> relative error          c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))
-	                         3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)
-	   - nbjacoby :          correction used by Hessiantype=1 (default is 1)
-	   - nbsmooth :          number of metric smoothing procedure (default is 3)
-	   - omega :             relaxation parameter of the smoothing procedure (default is 1.8)
-	   - power :             power applied to the metric (default is 1)
-	   - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
-	   - verbose :           level of verbosity (default is 1)
-
-	   - rifts :             followed by an ARGUS file that prescribes the rifts
-	   - toltip :            tolerance to move tip on an existing point of the domain outline
-	   - tracks :            followed by an ARGUS file that prescribes the tracks that the mesh will stick to
-	   - RequiredVertices :  mesh vertices that are required. [x,y,ref]; ref is optional
-	   - tol :               if the distance between 2 points of the domain outline is less than tol, they
-	                         will be merged
-
-	   Examples:
-	      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
-	      md=bamg(md,'field',[md.inversion.vel_obs md.geometry.thickness],'hmax',20000,'hmin',1000);
-	      md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1);
-	"""
-
-	#process options
-	options=pairoptions(*args)
-#	options=deleteduplicates(options,1);
-
-	#initialize the structures required as input of Bamg
-	bamg_options=OrderedDict()
-	bamg_geometry=bamggeom()
-	bamg_mesh=bamgmesh()
-
-	subdomain_ref = 1
-	hole_ref = 1
-
-	# Bamg Geometry parameters {{{
-	if options.exist('domain'):
-		#Check that file exists
-		domainfile=options.getfieldvalue('domain')
-		if type(domainfile) == str:
-			if not os.path.exists(domainfile):
-				raise IOError("bamg error message: file '%s' not found" % domainfile)
-			domain=expread(domainfile)
-		else:
-			domain=domainfile
-
-		#Build geometry 
-		count=0
-		for i,domaini in enumerate(domain):
-			#Check that the domain is closed
-			if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
-				raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
-
-			#Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
-			if i:
-				flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0]
-				if np.any(np.logical_not(flags)):
-					raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
-
-			#Check orientation
-			nods=domaini['nods']-1#the domain are closed 1=end;
-
-			test = np.sum((domaini['x'][1:nods+1] - domaini['x'][0:nods])*(domaini['y'][1:nods+1] + domaini['y'][0:nods]))
-			if (i==0 and test>0) or (i>0 and test<0):
-				print('At least one contour was not correctly oriented and has been re-oriented')
-				domaini['x'] = np.flipud(domaini['x'])
-				domaini['y'] = np.flipud(domaini['y'])
-
-
-			#Add all points to bamg_geometry
-			nods=domaini['nods']-1    #the domain are closed 0=end
-			bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T))
-			bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
-			if i:
-				bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
-				subdomain_ref = subdomain_ref+1;
-			else:
-				bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,0]))
-			
-			#update counter
-			count+=nods
-
-		#Deal with domain holes
-		if options.exist('holes'):
-			holesfile=options.getfieldvalue('holes')
-			if type(holesfile) == str:
-				if not os.path.exists(holesfile):
-					raise IOError("bamg error message: file '%s' not found" % holesfile)
-				holes=expread(holesfile)
-			else:
-				holes=holesfile
-
-			#Build geometry 
-			for i,holei in enumerate(holes):
-				#Check that the hole is closed
-				if (holei['x'][0]!=holei['x'][-1] or holei['y'][0]!=holei['y'][-1]):
-					raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
-
-				#Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
-				flags=ContourToNodes(holei['x'],holei['y'],domainfile,0)[0]
-				if np.any(np.logical_not(flags)):
-					raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
-
-				#Check orientation
-				nods=holei['nods']-1#the hole are closed 1=end;
-				test = np.sum((holei['x'][1:nods+1] - holei['x'][0:nods])*(holei['y'][1:nods+1] + holei['y'][0:nods]))
-				if test<0:
-					print('At least one hole was not correctly oriented and has been re-oriented')
-					holei['x'] = np.flipud(holei['x'])
-					holei['y'] = np.flipud(holei['y'])
-
-				#Add all points to bamg_geometry
-				nods=holei['nods']-1    #the hole are closed 0=end
-				bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((holei['x'][0:nods],holei['y'][0:nods],np.ones((nods)))).T))
-				bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
-				#update counter
-				count+=nods
-
-		#And subdomains
-		if options.exist('subdomains'):
-			subdomainfile=options.getfieldvalue('subdomains')
-			if type(subdomainfile) == str:
-				if not os.path.exists(subdomainfile):
-					raise IOError("bamg error message: file '%s' not found" % subdomainfile)
-				subdomains=expread(subdomainfile)
-			else:
-				subdomains=subdomainfile
-
-			#Build geometry 
-			for i,subdomaini in enumerate(subdomains):
-				#Check that the subdomain is closed
-				if (subdomaini['x'][0]!=subdomaini['x'][-1] or subdomaini['y'][0]!=subdomaini['y'][-1]):
-					raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed")
-
-				#Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0
-				if i:
-					flags=ContourToNodes(subdomaini['x'],subdomaini['y'],domainfile,0)[0]
-					if np.any(np.logical_not(flags)):
-						raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain")
-
-				#Check orientation
-				nods=subdomaini['nods']-1#the subdomain are closed 1=end;
-
-				test = np.sum((subdomaini['x'][1:nods+1] - subdomaini['x'][0:nods])*(subdomaini['y'][1:nods+1] + subdomaini['y'][0:nods]))
-				if test>0:
-					print('At least one subcontour was not correctly oriented and has been re-oriented')
-					subdomaini['x'] = np.flipud(subdomaini['x'])
-					subdomaini['y'] = np.flipud(subdomaini['y'])
-
-				#Add all points to bamg_geometry
-				nods=subdomaini['nods']-1    #the subdomain are closed 0=end
-				bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((subdomaini['x'][0:nods],subdomaini['y'][0:nods],np.ones((nods)))).T))
-				bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
-				#update counter
-				count+=nods
-			
-		if options.getfieldvalue('vertical',0):
-			if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0):
-				raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0)))
-			if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0):
-				bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers')
-
-		#take care of rifts
-		if options.exist('rifts'):
-
-			#Check that file exists
-			riftfile=options.getfieldvalue('rifts')
-			if not os.path.exists(riftfile):
-				raise IOError("bamg error message: file '%s' not found" % riftfile)
-			rift=expread(riftfile)
-
-			for i,rifti in enumerate(rift):
-
-				#detect whether all points of the rift are inside the domain
-				flags=ContourToNodes(rifti['x'],rifti['y'],domain[0],0)[0]
-				if np.all(np.logical_not(flags)):
-					raise RuntimeError("one rift has all its points outside of the domain outline")
-
-				elif np.any(np.logical_not(flags)):
-					#We LOTS of work to do
-					print("Rift tip outside of or on the domain has been detected and is being processed...")
-
-					#check that only one point is outside (for now)
-					if np.sum(np.logical_not(flags).astype(int))!=1:
-						raise RuntimeError("bamg error message: only one point outside of the domain is supported yet")
-
-					#Move tip outside to the first position
-					if   not flags[0]:
-						#OK, first point is outside (do nothing),
-						pass
-					elif not flags[-1]:
-						rifti['x']=np.flipud(rifti['x'])
-						rifti['y']=np.flipud(rifti['y'])
-					else:
-						raise RuntimeError("bamg error message: only a rift tip can be outside of the domain")
-
-					#Get cordinate of intersection point
-					x1=rifti['x'][0]
-					y1=rifti['y'][0]
-					x2=rifti['x'][1]
-					y2=rifti['y'][1]
-					for j in range(0,np.size(domain[0]['x'])-1):
-						if SegIntersect(np.array([[x1,y1],[x2,y2]]),np.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])):
-
-							#Get position of the two nodes of the edge in domain
-							i1=j
-							i2=j+1
-
-							#rift is crossing edge [i1 i2] of the domain
-							#Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
-							x3=domain[0]['x'][i1]
-							y3=domain[0]['y'][i1]
-							x4=domain[0]['x'][i2]
-							y4=domain[0]['y'][i2]
-#							x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-#							y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-							x=np.linalg.det(np.array([[np.linalg.det(np.array([[x1,y1],[x2,y2]])),x1-x2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),x3-x4]]))/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))
-							y=np.linalg.det(np.array([[np.linalg.det(np.array([[x1,y1],[x2,y2]])),y1-y2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),y3-y4]]))/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))
-
-							segdis= sqrt((x4-x3)**2+(y4-y3)**2)
-							tipdis=np.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)])
-
-							if np.min(tipdis)/segdis < options.getfieldvalue('toltip',0):
-								print("moving tip-domain intersection point")
-
-								#Get position of the closer point
-								if tipdis[0]>tipdis[1]:
-									pos=i2
-								else:
-									pos=i1
-
-								#This point is only in Vertices (number pos).
-								#OK, now we can add our own rift
-								nods=rifti['nods']-1
-								bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1))))))
-								bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\
-									np.array([[pos,count+1,(1+i)]]),\
-									np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),(1+i)*np.ones((nods-1,1))))))
-								count+=nods
-
-								break
-
-							else:
-								#Add intersection point to Vertices
-								bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.array([[x,y,1]])))
-								count+=1
-
-								#Decompose the crossing edge into 2 subedges
-								pos=np.nonzero(np.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0]
-								if not pos:
-									raise RuntimeError("bamg error message: a problem occurred...")
-								bamg_geometry.Edges=np.vstack((bamg_geometry.Edges[0:pos-1,:],\
-									np.array([[bamg_geometry.Edges[pos,0],count                     ,bamg_geometry.Edges[pos,2]]]),\
-									np.array([[count                     ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\
-									bamg_geometry.Edges[pos+1:,:]))
-
-								#OK, now we can add our own rift
-								nods=rifti['nods']-1
-								bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.hstack((rifti['x'][1:].reshape(-1,),rifti['y'][1:].reshape(-1,),np.ones((nods,1))))))
-								bamg_geometry.Edges=np.vstack((bamg_geometry.Edges,\
-									np.array([[count,count+1,2]]),\
-									np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),(1+i)*np.ones((nods-1,1))))))
-								count+=nods
-
-								break
-
-				else:
-					nods=rifti['nods']-1
-					bamg_geometry.Vertices=np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:],rifti['y'][:],np.ones((nods+1,1))))
-					bamg_geometry.Edges   =np.vstack(bamg_geometry.Edges, np.hstack(np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),i*np.ones((nods,1))))
-					count=+nods+1
-
-		#Deal with tracks
-		if options.exist('tracks'):
-
-			#read tracks
-			track=options.getfieldvalue('tracks')
-			if all(isinstance(track,str)):
-				A=expread(track)
-				track=np.hstack((A.x.reshape(-1,),A.y.reshape(-1,)))
-			else:
-				track=float(track)    #for some reason, it is of class "single"
-			if np.size(track,axis=1)==2:
-				track=np.hstack((track,3.*np.ones((size(track,axis=0),1))))
-
-			#only keep those inside
-			flags=ContourToNodes(track[:,0],track[:,1],domainfile,0)[0]
-			track=track[np.nonzero(flags),:]
-
-			#Add all points to bamg_geometry
-			nods=np.size(track,axis=0)
-			bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,track))
-			bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.hstack((np.arange(count+1,count+nods).reshape(-1,),np.arange(count+2,count+nods+1).reshape(-1,),3.*np.ones((nods-1,1))))))
-
-			#update counter
-			count+=nods
-
-		#Deal with vertices that need to be kept by mesher
-		if options.exist('RequiredVertices'):
-
-			#recover RequiredVertices
-			requiredvertices=options.getfieldvalue('RequiredVertices')    #for some reason, it is of class "single"
-			if np.size(requiredvertices,axis=1)==2:
-				requiredvertices=np.hstack((requiredvertices,4.*np.ones((np.size(requiredvertices,axis=0),1))))
-
-			#only keep those inside
-			flags=ContourToNodes(requiredvertices[:,0],requiredvertices[:,1],domainfile,0)[0]
-			requiredvertices=requiredvertices[np.nonzero(flags)[0],:]
-			#Add all points to bamg_geometry
-			nods=np.size(requiredvertices,axis=0)
-			bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,requiredvertices))
-
-			#update counter
-			count+=nods
-
-		#process geom
-		#bamg_geometry=processgeometry(bamg_geometry,options.getfieldvalue('tol',float(nan)),domain[0])
-
-	elif isinstance(md.private.bamg,dict) and 'geometry' in md.private.bamg:
-		bamg_geometry=bamggeom(md.private.bamg['geometry'].__dict__) 
-	else:
-		#do nothing...
-		pass
-	#}}}
-	# Bamg Mesh parameters {{{
-	if not options.exist('domain') and md.mesh.numberofvertices and m.strcmp(md.mesh.elementtype(),'Tria'):
-
-		if isinstance(md.private.bamg,dict) and 'mesh' in md.private.bamg:
-			bamg_mesh=bamgmesh(md.private.bamg['mesh'].__dict__)
-		else:
-			bamg_mesh.Vertices=np.vstack((md.mesh.x,md.mesh.y,np.ones((md.mesh.numberofvertices)))).T
-			#bamg_mesh.Vertices=np.hstack((md.mesh.x.reshape(-1,),md.mesh.y.reshape(-1,),np.ones((md.mesh.numberofvertices,1))))
-			bamg_mesh.Triangles=np.hstack((md.mesh.elements,np.ones((md.mesh.numberofelements,1))))
-
-		if isinstance(md.rifts.riftstruct,dict):
-			raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg")
-	#}}}
-	# Bamg Options {{{
-	bamg_options['Crack']=options.getfieldvalue('Crack',0)
-	bamg_options['anisomax']=options.getfieldvalue('anisomax',10.**30)
-	bamg_options['coeff']=options.getfieldvalue('coeff',1.)
-	bamg_options['cutoff']=options.getfieldvalue('cutoff',10.**-5)
-	bamg_options['err']=options.getfieldvalue('err',np.array([[0.01]]))
-	bamg_options['errg']=options.getfieldvalue('errg',0.1)
-	bamg_options['field']=options.getfieldvalue('field',np.empty((0,1)))
-	bamg_options['gradation']=options.getfieldvalue('gradation',1.5)
-	bamg_options['Hessiantype']=options.getfieldvalue('Hessiantype',0)
-	bamg_options['hmin']=options.getfieldvalue('hmin',10.**-100)
-	bamg_options['hmax']=options.getfieldvalue('hmax',10.**100)
-	bamg_options['hminVertices']=options.getfieldvalue('hminVertices',np.empty((0,1)))
-	bamg_options['hmaxVertices']=options.getfieldvalue('hmaxVertices',np.empty((0,1)))
-	bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1)))
-	bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1)
-	bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6)
-	bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.)
-	bamg_options['metric']=options.getfieldvalue('metric',np.empty((0,1)))
-	bamg_options['Metrictype']=options.getfieldvalue('Metrictype',0)
-	bamg_options['nbjacobi']=options.getfieldvalue('nbjacobi',1)
-	bamg_options['nbsmooth']=options.getfieldvalue('nbsmooth',3)
-	bamg_options['omega']=options.getfieldvalue('omega',1.8)
-	bamg_options['power']=options.getfieldvalue('power',1.)
-	bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1)
-	bamg_options['verbose']=options.getfieldvalue('verbose',1)
-	#}}}
-
-	#call Bamg
-	bamgmesh_out,bamggeom_out=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options)
-
-	# plug results onto model
-	if options.getfieldvalue('vertical',0):
-		md.mesh=mesh2dvertical()
-		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-		#Fill in rest of fields:
-		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-		md.mesh.numberofvertices=np.size(md.mesh.x)
-		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
-
-		#Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
-		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,)
-		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1
-
-	elif options.getfieldvalue('3dsurface',0):
-		md.mesh=mesh3dsurface()
-		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-		md.mesh.z=md.mesh.x
-		md.mesh.z[:]=0
-		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-		#Fill in rest of fields:
-		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-		md.mesh.numberofvertices=np.size(md.mesh.x)
-		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
-
-	else:
-		md.mesh=mesh2d()
-		md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-		md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-		md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-		md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-		md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-		md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-		#Fill in rest of fields:
-		md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-		md.mesh.numberofvertices=np.size(md.mesh.x)
-		md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-		md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-		md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
-
-	#Bamg private fields
-	md.private.bamg=OrderedDict()
-	md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
-	md.private.bamg['geometry']=bamggeom(bamggeom_out)
-	md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
-	md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0
-	md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
-
-	#Check for orphan
-	if np.any(np.logical_not(np.in1d(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements.flat))):
-		raise RuntimeError("Output mesh has orphans. Check your Domain and/or RequiredVertices")
-
-	return md
-
-def processgeometry(geom,tol,outline):    # {{{
-
-	raise RuntimeError("bamg.py/processgeometry is not complete.")
-	#Deal with edges
-	print("Checking Edge crossing...")
-	i=0
-	while (i<np.size(geom.Edges,axis=0)):
-
-		#edge counter
-		i+=1
-
-		#Get coordinates
-		x1=geom.Vertices[geom.Edges[i,0],0]
-		y1=geom.Vertices[geom.Edges[i,0],1]
-		x2=geom.Vertices[geom.Edges[i,1],0]
-		y2=geom.Vertices[geom.Edges[i,1],1]
-		color1=geom.Edges[i,2]
-
-		j=i    #test edges located AFTER i only
-		while (j<np.size(geom.Edges,axis=0)):
-
-			#edge counter
-			j+=1
-
-			#Skip if the two edges already have a vertex in common
-			if any(m.ismember(geom.Edges[i,0:2],geom.Edges[j,0:2])):
-				continue
-
-			#Get coordinates
-			x3=geom.Vertices[geom.Edges[j,0],0]
-			y3=geom.Vertices[geom.Edges[j,0],1]
-			x4=geom.Vertices[geom.Edges[j,1],0]
-			y4=geom.Vertices[geom.Edges[j,1],1]
-			color2=geom.Edges[j,2]
-
-			#Check if the two edges are crossing one another
-			if SegIntersect(np.array([[x1,y1],[x2,y2]]),np.array([[x3,y3],[x4,y4]])):
-
-				#Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
-				x=np.linalg.det(np.array([np.linalg.det(np.array([[x1,y1],[x2,y2]])),x1-x2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),x3-x4])/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
-				y=np.linalg.det(np.array([np.linalg.det(np.array([[x1,y1],[x2,y2]])),y1-y2],[np.linalg.det(np.array([[x3,y3],[x4,y4]])),y3-y4])/np.linalg.det(np.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])))
-
-				#Add vertex to the list of vertices
-				geom.Vertices=np.vstack((geom.Vertices,[x,y,min(color1,color2)]))
-				id=np.size(geom.Vertices,axis=0)
-
-				#Update edges i and j
-				edgei=geom.Edges[i,:].copy()
-				edgej=geom.Edges[j,:].copy()
-				geom.Edges[i,:]    =[edgei(0),id      ,edgei(2)]
-				geom.Edges=np.vstack((geom.Edges,[id      ,edgei(1),edgei(2)]))
-				geom.Edges[j,:]    =[edgej(0),id      ,edgej(2)]
-				geom.Edges=np.vstack((geom.Edges,[id      ,edgej(1),edgej(2)]))
-
-				#update current edge second tip
-				x2=x
-				y2=y
-
-	#Check point outside
-	print("Checking for points outside the domain...")
-	i=0
-	num=0
-	while (i<np.size(geom.Vertices,axis=0)):
-
-		#vertex counter
-		i+=1
-
-		#Get coordinates
-		x=geom.Vertices[i,0]
-		y=geom.Vertices[i,1]
-		color=geom.Vertices[i,2]
-
-		#Check that the point is inside the domain
-		if color!=1 and not ContourToNodes(x,y,outline[0],1):
-
-			#Remove points from list of Vertices
-			num+=1
-			geom.Vertices[i,:]=[]
-
-			#update edges
-			posedges=np.nonzero(geom.Edges==i)
-			geom.Edges[posedges[0],:]=[]
-			posedges=np.nonzero(geom.Edges>i)
-			geom.Edges[posedges]=geom.Edges[posedges]-1
-
-			#update counter
-			i-=1
-
-	if num:
-		print(("WARNING: %d points outside the domain outline have been removed" % num))
-
-	"""
-	%Check point spacing
-	if ~isnan(tol),
-		disp('Checking point spacing...');
-		i=0;
-		while (i<size(geom.Vertices,1)),
-
-			%vertex counter
-			i=i+1;
-
-			%Get coordinates
-			x1=geom.Vertices(i,1);
-			y1=geom.Vertices(i,2);
-
-			j=i; %test edges located AFTER i only
-			while (j<size(geom.Vertices,1)),
-
-				%vertex counter
-				j=j+1;
-
-				%Get coordinates
-				x2=geom.Vertices(j,1);
-				y2=geom.Vertices(j,2);
-
-				%Check whether the two vertices are too close
-				if ((x2-x1)**2+(y2-y1)**2<tol**2)
-
-					%Remove points from list of Vertices
-					geom.Vertices(j,:)=[];
-
-					%update edges
-					posedges=find(m.ismember(geom.Edges,j));
-					geom.Edges(posedges)=i;
-					posedges=find(geom.Edges>j);
-					geom.Edges(posedges)=geom.Edges(posedges)-1;
-
-					%update counter
-					j=j-1;
-
-				end
-			end
-		end
-	end
-	%remove empty edges
-	geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
-	"""
-	return geom
+
+def bamg(md, *args):
+    """
+    BAMG - mesh generation
+
+       Available options (for more details see ISSM website http: / / issm.jpl.nasa.gov / ):
+
+       - domain :            followed by an ARGUS file that prescribes the domain outline
+       - holes :             followed by an ARGUS file that prescribes the holes
+       - subdomains :        followed by an ARGUS file that prescribes the list of
+                             subdomains (that need to be inside domain)
+
+       - hmin :              minimum edge length (default is 10^ - 100)
+       - hmax :              maximum edge length (default is 10^100)
+       - hVertices :         imposed edge length for each vertex (geometry or mesh)
+       - hminVertices :      minimum edge length for each vertex (mesh)
+       - hmaxVertices :      maximum edge length for each vertex (mesh)
+
+       - anisomax :          maximum ratio between the smallest and largest edges (default is 10^30)
+       - coeff :             coefficient applied to the metric (2 - > twice as many elements, default is 1)
+       - cutoff :            scalar used to compute the metric when metric type 2 or 3 are applied
+       - err :               error used to generate the metric from a field
+       - errg :              geometric error (default is 0.1)
+       - field :             field of the model that will be used to compute the metric
+                                   to apply several fields, use one column per field
+       - gradation :         maximum ratio between two adjacent edges
+       - Hessiantype :       0 - > use double P2 projection (default)
+                                   1 - > use Green formula
+       - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
+       - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
+       - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
+       - metric :            matrix (numberofnodes x 3) used as a metric
+       - Metrictype :        1 - > absolute error          c / (err coeff^2) * Abs(H)        (default)
+                                   2 - > relative error          c / (err coeff^2) * Abs(H) / max(s, cutoff * max(s))
+                                   3 - > rescaled absolute error c / (err coeff^2) * Abs(H) / (smax - smin)
+       - nbjacoby :          correction used by Hessiantype = 1 (default is 1)
+       - nbsmooth :          number of metric smoothing procedure (default is 3)
+       - omega :             relaxation parameter of the smoothing procedure (default is 1.8)
+       - power :             power applied to the metric (default is 1)
+       - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
+       - verbose :           level of verbosity (default is 1)
+
+       - rifts :             followed by an ARGUS file that prescribes the rifts
+       - toltip :            tolerance to move tip on an existing point of the domain outline
+       - tracks :            followed by an ARGUS file that prescribes the tracks that the mesh will stick to
+       - RequiredVertices :  mesh vertices that are required. [x, y, ref]; ref is optional
+       - tol :               if the distance between 2 points of the domain outline is less than tol, they
+                             will be merged
+
+       Examples:
+          md = bamg(md, 'domain', 'DomainOutline.exp', 'hmax', 3000)
+          md = bamg(md, 'field', [md.inversion.vel_obs md.geometry.thickness], 'hmax', 20000, 'hmin', 1000)
+          md = bamg(md, 'metric', A, 'hmin', 1000, 'hmax', 20000, 'gradation', 3, 'anisomax', 1)
+    """
+
+    #process options
+    options = pairoptions(*args)
+    #options = deleteduplicates(options, 1)
+
+    #initialize the structures required as input of Bamg
+    bamg_options = OrderedDict()
+    bamg_geometry = bamggeom()
+    bamg_mesh = bamgmesh()
+
+    subdomain_ref = 1
+    hole_ref = 1
+
+    # Bamg Geometry parameters {{{
+    if options.exist('domain'):
+        #Check that file exists
+        domainfile = options.getfieldvalue('domain')
+        if type(domainfile) == str:
+            if not os.path.exists(domainfile):
+                raise IOError("bamg error message: file '%s' not found" % domainfile)
+            domain = expread(domainfile)
+        else:
+            domain = domainfile
+
+        #Build geometry
+        count = 0
+        for i, domaini in enumerate(domain):
+            #Check that the domain is closed
+            if (domaini['x'][0] != domaini['x'][-1] or domaini['y'][0] != domaini['y'][-1]):
+                raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
+
+            #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
+            if i:
+                flags = ContourToNodes(domaini['x'], domaini['y'], domainfile, 0)[0]
+                if np.any(np.logical_not(flags)):
+                    raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
+
+            #Check orientation
+            nods = domaini['nods'] - 1  #the domain are closed 1 = end
+
+            test = np.sum((domaini['x'][1:nods + 1] - domaini['x'][0:nods]) * (domaini['y'][1:nods + 1] + domaini['y'][0:nods]))
+            if (i == 0 and test > 0) or (i > 0 and test < 0):
+                print('At least one contour was not correctly oriented and has been re - oriented')
+                domaini['x'] = np.flipud(domaini['x'])
+                domaini['y'] = np.flipud(domaini['y'])
+
+            #Add all points to bamg_geometry
+            nods = domaini['nods'] - 1    #the domain are closed 0 = end
+            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((domaini['x'][0:nods], domaini['y'][0:nods], np.ones((nods)))).T))
+            bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
+            if i:
+                bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, - subdomain_ref]))
+                subdomain_ref = subdomain_ref + 1
+            else:
+                bamg_geometry.SubDomains = np.vstack((bamg_geometry.SubDomains, [2, count + 1, 1, 0]))
+
+            #update counter
+            count += nods
+
+        #Deal with domain holes
+        if options.exist('holes'):
+            holesfile = options.getfieldvalue('holes')
+            if type(holesfile) == str:
+                if not os.path.exists(holesfile):
+                    raise IOError("bamg error message: file '%s' not found" % holesfile)
+                holes = expread(holesfile)
+            else:
+                holes = holesfile
+
+            #Build geometry
+            for i, holei in enumerate(holes):
+                #Check that the hole is closed
+                if (holei['x'][0] != holei['x'][-1] or holei['y'][0] != holei['y'][-1]):
+                    raise RuntimeError("bamg error message: all contours provided in ''hole'' should be closed")
+
+                #Checks that all holes are INSIDE the principle domain outline princial domain should be index 0
+                flags = ContourToNodes(holei['x'], holei['y'], domainfile, 0)[0]
+                if np.any(np.logical_not(flags)):
+                    raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
+
+                #Check orientation
+                nods = holei['nods'] - 1  #the hole are closed 1 = end
+                test = np.sum((holei['x'][1:nods + 1] - holei['x'][0:nods]) * (holei['y'][1:nods + 1] + holei['y'][0:nods]))
+                if test < 0:
+                    print('At least one hole was not correctly oriented and has been re - oriented')
+                    holei['x'] = np.flipud(holei['x'])
+                    holei['y'] = np.flipud(holei['y'])
+
+                #Add all points to bamg_geometry
+                nods = holei['nods'] - 1    #the hole are closed 0 = end
+                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((holei['x'][0:nods], holei['y'][0:nods], np.ones((nods)))).T))
+                bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
+                #update counter
+                count += nods
+
+        #And subdomains
+        if options.exist('subdomains'):
+            subdomainfile = options.getfieldvalue('subdomains')
+            if type(subdomainfile) == str:
+                if not os.path.exists(subdomainfile):
+                    raise IOError("bamg error message: file '{}' not found".format(subdomainfile))
+                subdomains = expread(subdomainfile)
+            else:
+                subdomains = subdomainfile
+
+            #Build geometry
+            for i, subdomaini in enumerate(subdomains):
+                #Check that the subdomain is closed
+                if (subdomaini['x'][0] != subdomaini['x'][-1] or subdomaini['y'][0] != subdomaini['y'][-1]):
+                    raise RuntimeError("bamg error message: all contours provided in ''subdomain'' should be closed")
+
+                #Checks that all subdomains are INSIDE the principle subdomain outline princial domain should be index 0
+                if i:
+                    flags = ContourToNodes(subdomaini['x'], subdomaini['y'], domainfile, 0)[0]
+                    if np.any(np.logical_not(flags)):
+                        raise RuntimeError("bamg error message: All subdomains should be strictly inside the principal subdomain")
+
+                #Check orientation
+                nods = subdomaini['nods'] - 1  #the subdomain are closed 1 = end
+
+                test = np.sum((subdomaini['x'][1:nods + 1] - subdomaini['x'][0:nods]) * (subdomaini['y'][1:nods + 1] + subdomaini['y'][0:nods]))
+                if test > 0:
+                    print('At least one subcontour was not correctly oriented and has been re - oriented')
+                    subdomaini['x'] = np.flipud(subdomaini['x'])
+                    subdomaini['y'] = np.flipud(subdomaini['y'])
+
+                #Add all points to bamg_geometry
+                nods = subdomaini['nods'] - 1    #the subdomain are closed 0 = end
+                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.vstack((subdomaini['x'][0:nods], subdomaini['y'][0:nods], np.ones((nods)))).T))
+                bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.vstack((np.arange(count + 1, count + nods + 1), np.hstack((np.arange(count + 2, count + nods + 1), count + 1)), 1. * np.ones((nods)))).T))
+                #update counter
+                count += nods
+
+        if options.getfieldvalue('vertical', 0):
+            if np.size(options.getfieldvalue('Markers', [])) != np.size(bamg_geometry.Edges, 0):
+                raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges, 0)))
+            if np.size(options.getfieldvalue('Markers', [])) == np.size(bamg_geometry.Edges, 0):
+                bamg_geometry.Edges[:, 2] = options.getfieldvalue('Markers')
+
+        #take care of rifts
+        if options.exist('rifts'):
+
+            #Check that file exists
+            riftfile = options.getfieldvalue('rifts')
+            if not os.path.exists(riftfile):
+                raise IOError("bamg error message: file '%s' not found" % riftfile)
+            rift = expread(riftfile)
+
+            for i, rifti in enumerate(rift):
+
+                #detect whether all points of the rift are inside the domain
+                flags = ContourToNodes(rifti['x'], rifti['y'], domain[0], 0)[0]
+                if np.all(np.logical_not(flags)):
+                    raise RuntimeError("one rift has all its points outside of the domain outline")
+
+                elif np.any(np.logical_not(flags)):
+                    #We LOTS of work to do
+                    print("Rift tip outside of or on the domain has been detected and is being processed...")
+
+                    #check that only one point is outside (for now)
+                    if np.sum(np.logical_not(flags).astype(int)) != 1:
+                        raise RuntimeError("bamg error message: only one point outside of the domain is supported yet")
+
+                    #Move tip outside to the first position
+                    if not flags[0]:
+                        #OK, first point is outside (do nothing),
+                        pass
+                    elif not flags[-1]:
+                        rifti['x'] = np.flipud(rifti['x'])
+                        rifti['y'] = np.flipud(rifti['y'])
+                    else:
+                        raise RuntimeError("bamg error message: only a rift tip can be outside of the domain")
+
+                    #Get cordinate of intersection point
+                    x1 = rifti['x'][0]
+                    y1 = rifti['y'][0]
+                    x2 = rifti['x'][1]
+                    y2 = rifti['y'][1]
+                    for j in range(0, np.size(domain[0]['x']) - 1):
+                        if SegIntersect(np.array([[x1, y1], [x2, y2]]), np.array([[domain[0]['x'][j], domain[0]['y'][j]], [domain[0]['x'][j + 1], domain[0]['y'][j + 1]]])):
+
+                            #Get position of the two nodes of the edge in domain
+                            i1 = j
+                            i2 = j + 1
+
+                            #rift is crossing edge [i1 i2] of the domain
+                            #Get coordinate of intersection point (http: / / mathworld.wolfram.com / Line - LineIntersection.html)
+                            x3 = domain[0]['x'][i1]
+                            y3 = domain[0]['y'][i1]
+                            x4 = domain[0]['x'][i2]
+                            y4 = domain[0]['y'][i2]
+                            #x = det([det([x1 y1; x2 y2])  x1 - x2;det([x3 y3; x4 y4])  x3 - x4]) / det([x1 - x2 y1 - y2;x3 - x4 y3 - y4])
+                            #y = det([det([x1 y1; x2 y2])  y1 - y2;det([x3 y3; x4 y4])  y3 - y4]) / det([x1 - x2 y1 - y2;x3 - x4 y3 - y4])
+                            x = np.linalg.det(np.array([[np.linalg.det(np.array([[x1, y1], [x2, y2]])), x1 - x2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), x3 - x4]])) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]]))
+                            y = np.linalg.det(np.array([[np.linalg.det(np.array([[x1, y1], [x2, y2]])), y1 - y2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), y3 - y4]])) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]]))
+
+                            segdis = sqrt((x4 - x3)**2 + (y4 - y3)**2)
+                            tipdis = np.array([sqrt((x - x3)**2 + (y - y3)**2), sqrt((x - x4)**2 + (y - y4)**2)])
+
+                            if np.min(tipdis) / segdis < options.getfieldvalue('toltip', 0):
+                                print("moving tip - domain intersection point")
+
+                                #Get position of the closer point
+                                if tipdis[0] > tipdis[1]:
+                                    pos = i2
+                                else:
+                                    pos = i1
+
+                                #This point is only in Vertices (number pos).
+                                #OK, now we can add our own rift
+                                nods = rifti['nods'] - 1
+                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(- 1, ), rifti['y'][1:].reshape(- 1, ), np.ones((nods, 1))))))
+                                bamg_geometry.Edges = np.vstack((bamg_geometry.Edges,
+                                                                 np.array([[pos, count + 1, (1 + i)]]),
+                                                                 np.hstack((np.arange(count + 1, count + nods).reshape(- 1, ), np.arange(count + 2, count + nods + 1).reshape(- 1, ), (1 + i) * np.ones((nods - 1, 1))))))
+                                count += nods
+
+                                break
+
+                            else:
+                                #Add intersection point to Vertices
+                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.array([[x, y, 1]])))
+                                count += 1
+
+                                #Decompose the crossing edge into 2 subedges
+                                pos = np.nonzero(np.logical_and(bamg_geometry.Edges[:, 0] == i1, bamg_geometry.Edges[:, 1] == i2))[0]
+                                if not pos:
+                                    raise RuntimeError("bamg error message: a problem occurred...")
+                                bamg_geometry.Edges = np.vstack((bamg_geometry.Edges[0:pos - 1, :],
+                                                                 np.array([[bamg_geometry.Edges[pos, 0], count, bamg_geometry.Edges[pos, 2]]]),
+                                                                 np.array([[count, bamg_geometry.Edges[pos, 1], bamg_geometry.Edges[pos, 2]]]),
+                                                                 bamg_geometry.Edges[pos + 1:, :]))
+
+                                #OK, now we can add our own rift
+                                nods = rifti['nods'] - 1
+                                bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, np.hstack((rifti['x'][1:].reshape(- 1, ), rifti['y'][1:].reshape(- 1, ), np.ones((nods, 1))))))
+                                bamg_geometry.Edges = np.vstack((bamg_geometry.Edges,
+                                                                 np.array([[count, count + 1, 2]]),
+                                                                 np.hstack((np.arange(count + 1, count + nods).reshape(- 1, ), np.arange(count + 2, count + nods + 1).reshape(- 1, ), (1 + i) * np.ones((nods - 1, 1))))))
+                                count += nods
+
+                                break
+
+                else:
+                    nods = rifti['nods'] - 1
+                    bamg_geometry.Vertices = np.vstack(bamg_geometry.Vertices, np.hstack(rifti['x'][:], rifti['y'][:], np.ones((nods + 1, 1))))
+                    bamg_geometry.Edges = np.vstack(bamg_geometry.Edges, np.hstack(np.arange(count + 1, count + nods).reshape(- 1, ), np.arange(count + 2, count + nods + 1).reshape(- 1, ), i * np.ones((nods, 1))))
+                    count += nods + 1
+
+        #Deal with tracks
+        if options.exist('tracks'):
+
+            #read tracks
+            track = options.getfieldvalue('tracks')
+            if all(isinstance(track, str)):
+                A = expread(track)
+                track = np.hstack((A.x.reshape(- 1, ), A.y.reshape(- 1, )))
+            else:
+                track = float(track)    #for some reason, it is of class "single"
+            if np.size(track, axis = 1) == 2:
+                track = np.hstack((track, 3. * np.ones((size(track, axis = 0), 1))))
+
+            #only keep those inside
+            flags = ContourToNodes(track[:, 0], track[:, 1], domainfile, 0)[0]
+            track = track[np.nonzero(flags), :]
+
+            #Add all points to bamg_geometry
+            nods = np.size(track, axis = 0)
+            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, track))
+            bamg_geometry.Edges = np.vstack((bamg_geometry.Edges, np.hstack((np.arange(count + 1, count + nods).reshape(-1, ), np.arange(count + 2, count + nods + 1).reshape(-1, ), 3. * np.ones((nods - 1, 1))))))
+
+            #update counter
+            count += nods
+
+        #Deal with vertices that need to be kept by mesher
+        if options.exist('RequiredVertices'):
+
+            #recover RequiredVertices
+            requiredvertices = options.getfieldvalue('RequiredVertices')    #for some reason, it is of class "single"
+            if np.size(requiredvertices, axis = 1) == 2:
+                requiredvertices = np.hstack((requiredvertices, 4. * np.ones((np.size(requiredvertices, axis = 0), 1))))
+
+            #only keep those inside
+            flags = ContourToNodes(requiredvertices[:, 0], requiredvertices[:, 1], domainfile, 0)[0]
+            requiredvertices = requiredvertices[np.nonzero(flags)[0], :]
+            #Add all points to bamg_geometry
+            nods = np.size(requiredvertices, axis = 0)
+            bamg_geometry.Vertices = np.vstack((bamg_geometry.Vertices, requiredvertices))
+
+            #update counter
+            count += nods
+
+        #process geom
+        #bamg_geometry = processgeometry(bamg_geometry, options.getfieldvalue('tol', float(nan)), domain[0])
+
+    elif isinstance(md.private.bamg, dict) and 'geometry' in md.private.bamg:
+        bamg_geometry = bamggeom(md.private.bamg['geometry'].__dict__)
+    else:
+        #do nothing...
+        pass
+    #}}}
+    # Bamg Mesh parameters {{{
+    if not options.exist('domain') and md.mesh.numberofvertices and m.strcmp(md.mesh.elementtype(), 'Tria'):
+        if isinstance(md.private.bamg, dict) and 'mesh' in md.private.bamg:
+            bamg_mesh = bamgmesh(md.private.bamg['mesh'].__dict__)
+        else:
+            bamg_mesh.Vertices = np.vstack((md.mesh.x, md.mesh.y, np.ones((md.mesh.numberofvertices)))).T
+            #bamg_mesh.Vertices = np.hstack((md.mesh.x.reshape(- 1, ), md.mesh.y.reshape(- 1, ), np.ones((md.mesh.numberofvertices, 1))))
+            bamg_mesh.Triangles = np.hstack((md.mesh.elements, np.ones((md.mesh.numberofelements, 1))))
+
+        if isinstance(md.rifts.riftstruct, dict):
+            raise TypeError("bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg")
+    #}}}
+    # Bamg Options {{{
+    bamg_options['Crack'] = options.getfieldvalue('Crack', 0)
+    bamg_options['anisomax'] = options.getfieldvalue('anisomax', 10.**18)
+    bamg_options['coeff'] = options.getfieldvalue('coeff', 1.)
+    bamg_options['cutoff'] = options.getfieldvalue('cutoff', 10.**-5)
+    bamg_options['err'] = options.getfieldvalue('err', np.array([[0.01]]))
+    bamg_options['errg'] = options.getfieldvalue('errg', 0.1)
+    bamg_options['field'] = options.getfieldvalue('field', np.empty((0, 1)))
+    bamg_options['gradation'] = options.getfieldvalue('gradation', 1.5)
+    bamg_options['Hessiantype'] = options.getfieldvalue('Hessiantype', 0)
+    bamg_options['hmin'] = options.getfieldvalue('hmin', 10.**-100)
+    bamg_options['hmax'] = options.getfieldvalue('hmax', 10.**100)
+    bamg_options['hminVertices'] = options.getfieldvalue('hminVertices', np.empty((0, 1)))
+    bamg_options['hmaxVertices'] = options.getfieldvalue('hmaxVertices', np.empty((0, 1)))
+    bamg_options['hVertices'] = options.getfieldvalue('hVertices', np.empty((0, 1)))
+    bamg_options['KeepVertices'] = options.getfieldvalue('KeepVertices', 1)
+    bamg_options['maxnbv'] = options.getfieldvalue('maxnbv', 10**6)
+    bamg_options['maxsubdiv'] = options.getfieldvalue('maxsubdiv', 10.)
+    bamg_options['metric'] = options.getfieldvalue('metric', np.empty((0, 1)))
+    bamg_options['Metrictype'] = options.getfieldvalue('Metrictype', 0)
+    bamg_options['nbjacobi'] = options.getfieldvalue('nbjacobi', 1)
+    bamg_options['nbsmooth'] = options.getfieldvalue('nbsmooth', 3)
+    bamg_options['omega'] = options.getfieldvalue('omega', 1.8)
+    bamg_options['power'] = options.getfieldvalue('power', 1.)
+    bamg_options['splitcorners'] = options.getfieldvalue('splitcorners', 1)
+    bamg_options['verbose'] = options.getfieldvalue('verbose', 1)
+    #}}}
+
+    #call Bamg
+    bamgmesh_out, bamggeom_out = BamgMesher(bamg_mesh.__dict__, bamg_geometry.__dict__, bamg_options)
+
+    # plug results onto model
+    if options.getfieldvalue('vertical', 0):
+        md.mesh = mesh2dvertical()
+        md.mesh.x = bamgmesh_out['Vertices'][:, 0].copy()
+        md.mesh.y = bamgmesh_out['Vertices'][:, 1].copy()
+        md.mesh.elements = bamgmesh_out['Triangles'][:, 0:3].astype(int)
+        md.mesh.edges = bamgmesh_out['IssmEdges'].astype(int)
+        md.mesh.segments = bamgmesh_out['IssmSegments'][:, 0:3].astype(int)
+        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
+
+        #Fill in rest of fields:
+        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
+        md.mesh.numberofvertices = np.size(md.mesh.x)
+        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
+        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
+        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
+
+        #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
+        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, )
+        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
+
+    elif options.getfieldvalue('3dsurface', 0):
+        md.mesh = mesh3dsurface()
+        md.mesh.x = bamgmesh_out['Vertices'][:, 0].copy()
+        md.mesh.y = bamgmesh_out['Vertices'][:, 1].copy()
+        md.mesh.z = md.mesh.x
+        md.mesh.z[:] = 0
+        md.mesh.elements = bamgmesh_out['Triangles'][:, 0:3].astype(int)
+        md.mesh.edges = bamgmesh_out['IssmEdges'].astype(int)
+        md.mesh.segments = bamgmesh_out['IssmSegments'][:, 0:3].astype(int)
+        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
+
+        #Fill in rest of fields:
+        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
+        md.mesh.numberofvertices = np.size(md.mesh.x)
+        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
+        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
+        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
+
+    else:
+        md.mesh = mesh2d()
+        md.mesh.x = bamgmesh_out['Vertices'][:, 0].copy()
+        md.mesh.y = bamgmesh_out['Vertices'][:, 1].copy()
+        md.mesh.elements = bamgmesh_out['Triangles'][:, 0:3].astype(int)
+        md.mesh.edges = bamgmesh_out['IssmEdges'].astype(int)
+        md.mesh.segments = bamgmesh_out['IssmSegments'][:, 0:3].astype(int)
+        md.mesh.segmentmarkers = bamgmesh_out['IssmSegments'][:, 3].astype(int)
+
+        #Fill in rest of fields:
+        md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
+        md.mesh.numberofvertices = np.size(md.mesh.x)
+        md.mesh.numberofedges = np.size(md.mesh.edges, axis=0)
+        md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
+        md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
+
+    #Bamg private fields
+    md.private.bamg = OrderedDict()
+    md.private.bamg['mesh'] = bamgmesh(bamgmesh_out)
+    md.private.bamg['geometry'] = bamggeom(bamggeom_out)
+    md.mesh.elementconnectivity = md.private.bamg['mesh'].ElementConnectivity
+    md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))] = 0
+    md.mesh.elementconnectivity = md.mesh.elementconnectivity.astype(int)
+
+    #Check for orphan
+    if np.any(np.logical_not(np.in1d(np.arange(1, md.mesh.numberofvertices + 1), md.mesh.elements.flat))):
+        raise RuntimeError("Output mesh has orphans. Check your Domain and / or RequiredVertices")
+
+    return md
+
+
+def processgeometry(geom, tol, outline):    # {{{
+
+    raise RuntimeError("bamg.py / processgeometry is not complete.")
+    #Deal with edges
+    print("Checking Edge crossing...")
+    i = 0
+    while (i < np.size(geom.Edges, axis=0)):
+        #edge counter
+        i += 1
+
+        #Get coordinates
+        x1 = geom.Vertices[geom.Edges[i, 0], 0]
+        y1 = geom.Vertices[geom.Edges[i, 0], 1]
+        x2 = geom.Vertices[geom.Edges[i, 1], 0]
+        y2 = geom.Vertices[geom.Edges[i, 1], 1]
+        color1 = geom.Edges[i, 2]
+
+        j = i    #test edges located AFTER i only
+        while (j < np.size(geom.Edges, axis=0)):
+            #edge counter
+            j += 1
+
+            #Skip if the two edges already have a vertex in common
+            if any(m.ismember(geom.Edges[i, 0:2], geom.Edges[j, 0:2])):
+                continue
+
+            #Get coordinates
+            x3 = geom.Vertices[geom.Edges[j, 0], 0]
+            y3 = geom.Vertices[geom.Edges[j, 0], 1]
+            x4 = geom.Vertices[geom.Edges[j, 1], 0]
+            y4 = geom.Vertices[geom.Edges[j, 1], 1]
+            color2 = geom.Edges[j, 2]
+
+            #Check if the two edges are crossing one another
+            if SegIntersect(np.array([[x1, y1], [x2, y2]]), np.array([[x3, y3], [x4, y4]])):
+
+                #Get coordinate of intersection point (http: / / mathworld.wolfram.com / Line - LineIntersection.html)
+                x = np.linalg.det(np.array([np.linalg.det(np.array([[x1, y1], [x2, y2]])), x1 - x2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), x3 - x4]) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]])))
+                y = np.linalg.det(np.array([np.linalg.det(np.array([[x1, y1], [x2, y2]])), y1 - y2], [np.linalg.det(np.array([[x3, y3], [x4, y4]])), y3 - y4]) / np.linalg.det(np.array([[x1 - x2, y1 - y2], [x3 - x4, y3 - y4]])))
+
+                #Add vertex to the list of vertices
+                geom.Vertices = np.vstack((geom.Vertices, [x, y, min(color1, color2)]))
+                id = np.size(geom.Vertices, axis=0)
+
+                #Update edges i and j
+                edgei = geom.Edges[i, :].copy()
+                edgej = geom.Edges[j, :].copy()
+                geom.Edges[i, :] = [edgei(0), id, edgei(2)]
+                geom.Edges = np.vstack((geom.Edges, [id, edgei(1), edgei(2)]))
+                geom.Edges[j, :] = [edgej(0), id, edgej(2)]
+                geom.Edges = np.vstack((geom.Edges, [id, edgej(1), edgej(2)]))
+
+                #update current edge second tip
+                x2 = x
+                y2 = y
+
+    #Check point outside
+    print("Checking for points outside the domain...")
+    i = 0
+    num = 0
+    while (i < np.size(geom.Vertices, axis=0)):
+        #vertex counter
+        i += 1
+
+        #Get coordinates
+        x = geom.Vertices[i, 0]
+        y = geom.Vertices[i, 1]
+        color = geom.Vertices[i, 2]
+
+        #Check that the point is inside the domain
+        if color != 1 and not ContourToNodes(x, y, outline[0], 1):
+            #Remove points from list of Vertices
+            num += 1
+            geom.Vertices[i, :]=[]
+
+            #update edges
+            posedges = np.nonzero(geom.Edges == i)
+            geom.Edges[posedges[0], :]=[]
+            posedges = np.nonzero(geom.Edges > i)
+            geom.Edges[posedges] = geom.Edges[posedges] - 1
+
+            #update counter
+            i -= 1
+
+    if num:
+        print(("WARNING: %d points outside the domain outline have been removed" % num))
+
+    """
+    %Check point spacing
+    if ~isnan(tol),
+            disp('Checking point spacing...')
+            i = 0
+            while (i < size(geom.Vertices, 1)),
+
+                    %vertex counter
+                    i = i + 1
+
+                    %Get coordinates
+                    x1 = geom.Vertices(i, 1)
+                    y1 = geom.Vertices(i, 2)
+
+                    j = i; %test edges located AFTER i only
+                    while (j < size(geom.Vertices, 1)),
+
+                            %vertex counter
+                            j = j + 1
+
+                            %Get coordinates
+                            x2 = geom.Vertices(j, 1)
+                            y2 = geom.Vertices(j, 2)
+
+                            %Check whether the two vertices are too close
+                            if ((x2 - x1)**2 + (y2 - y1)**2 < tol**2)
+
+                                    %Remove points from list of Vertices
+                                    geom.Vertices(j, :)=[]
+
+                                    %update edges
+                                    posedges = find(m.ismember(geom.Edges, j))
+                                    geom.Edges(posedges)=i
+                                    posedges = find(geom.Edges > j)
+                                    geom.Edges(posedges)=geom.Edges(posedges) - 1
+
+                                    %update counter
+                                    j = j - 1
+
+                            end
+                    end
+            end
+    end
+    %remove empty edges
+    geom.Edges(find(geom.Edges(:, 1) == geom.Edges(:, 2)), :)=[]
+    """
+    return geom
 # }}}
-
Index: /issm/trunk-jpl/src/m/mesh/roundmesh.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/roundmesh.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/mesh/roundmesh.py	(revision 24115)
@@ -5,55 +5,54 @@
 from triangle import triangle
 
-def roundmesh(md,radius,resolution):
-	"""
-	ROUNDMESH - create an unstructured round mesh 
+def roundmesh(md, radius, resolution):
+    """
+    ROUNDMESH - create an unstructured round mesh
 
-	   This script will generate a structured round mesh
-	   - radius     : specifies the radius of the circle in meters
-	   - resolution : specifies the resolution in meters
+       This script will generate a structured round mesh
+       - radius     : specifies the radius of the circle in meters
+       - resolution : specifies the resolution in meters
 
-	   Usage:
-	      md=roundmesh(md,radius,resolution)
-	"""
+       Usage:
+          md=roundmesh(md,radius,resolution)
+    """
 
-	#First we have to create the domain outline 
+    #First we have to create the domain outline
 
-	#Get number of points on the circle
-	pointsonedge=np.floor((2.*np.pi*radius) / resolution)+1 #+1 to close the outline
+    #Get number of points on the circle
+    pointsonedge = np.floor((2. * np.pi * radius) / resolution) + 1  #+1 to close the outline
 
-	#Calculate the cartesians coordinates of the points
-	theta=np.linspace(0.,2.*np.pi,pointsonedge)
-	x_list=roundsigfig(radius*np.cos(theta),12)
-	y_list=roundsigfig(radius*np.sin(theta),12)
-	A=OrderedDict()
-	A['x']=[x_list]
-	A['y']=[y_list]
-	A['density']=1.
-	expwrite(A,'RoundDomainOutline.exp')
+    #Calculate the cartesians coordinates of the points
+    theta = np.linspace(0., 2. * np.pi, pointsonedge)
+    x_list = roundsigfig(radius * np.cos(theta), 12)
+    y_list = roundsigfig(radius * np.sin(theta), 12)
+    A = OrderedDict()
+    A['x'] = [x_list]
+    A['y'] = [y_list]
+    A['density'] = 1.
+    print('now writing mesh')
+    expwrite(A, 'RoundDomainOutline.exp')
 
-	#Call Bamg
-	md=triangle(md,'RoundDomainOutline.exp',resolution)
-	#md=bamg(md,'domain','RoundDomainOutline.exp','hmin',resolution)
+    #Call Bamg
+    print('now meshing')
+    md = triangle(md, 'RoundDomainOutline.exp', resolution)
+    #md = bamg(md,'domain','RoundDomainOutline.exp','hmin',resolution)
 
-	#move the closest node to the center
-	pos=np.argmin(md.mesh.x**2+md.mesh.y**2)
-	md.mesh.x[pos]=0.
-	md.mesh.y[pos]=0.
+    #move the closest node to the center
+    pos = np.argmin(md.mesh.x**2 + md.mesh.y**2)
+    md.mesh.x[pos] = 0.
+    md.mesh.y[pos] = 0.
 
-	#delete domain
-	os.remove('RoundDomainOutline.exp')
+    #delete domain
+    os.remove('RoundDomainOutline.exp')
 
-	return md
+    return md
 
-def roundsigfig(x,n):
 
-	digits=np.ceil(np.log10(np.abs(x)))
-	x=x/10.**digits
-	x=np.round(x,decimals=n)
-	x=x*10.**digits
+def roundsigfig(x, n):
 
-	pos=np.nonzero(np.isnan(x))
-	x[pos]=0.
-
-	return x
-
+    nonzeros = np.where(x != 0)
+    digits = np.ceil(np.log10(np.abs(x[nonzeros])))
+    x[nonzeros] = x[nonzeros] / 10.**digits
+    x[nonzeros] = np.round(x[nonzeros], decimals=n)
+    x[nonzeros] = x[nonzeros] * 10.**digits
+    return x
Index: /issm/trunk-jpl/src/m/mesh/triangle.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 24115)
@@ -5,63 +5,66 @@
 from ElementConnectivity import ElementConnectivity
 from Triangle_python import Triangle_python
-import MatlabFuncs as m
 
-def triangle(md,domainname,*args):
-	"""
-	TRIANGLE - create model mesh using the triangle package
 
-	   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
-	   where md is a @model object, domainname is the name of an Argus domain outline file, 
-	   and resolution is a characteristic length for the mesh (same unit as the domain outline
-	   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
+def triangle(md, domainname, *args):
+    """
+    TRIANGLE - create model mesh using the triangle package
 
-	   Usage:
-	      md=triangle(md,domainname,resolution)
-	   or md=triangle(md,domainname, resolution, riftname)
+       This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
+       where md is a @model object, domainname is the name of an Argus domain outline file,
+       and resolution is a characteristic length for the mesh (same unit as the domain outline
+       unit). Riftname is an optional argument (Argus domain outline) describing rifts.
 
-	   Examples:
-	      md=triangle(md,'DomainOutline.exp',1000);
-	      md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp');
-	"""
+       Usage:
+          md=triangle(md,domainname,resolution)
+       or md=triangle(md,domainname, resolution, riftname)
 
-	#Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would 
-	#be made of 1000*1000 area squares). 
+       Examples:
+          md=triangle(md,'DomainOutline.exp',1000);
+          md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp');
+    """
 
-	if len(args)==1:
-		resolution=args[0]
-		riftname=''
-	if len(args)==2:
-		riftname=args[0]
-		resolution=args[1]
+    #Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would
+    #be made of 1000*1000 area squares).
 
-	#Check that mesh was not already run, and warn user: 
-	if md.mesh.numberofelements:
-		choice = eval(input('This model already has a mesh. Are you sure you want to go ahead? (y/n)'))
-		if not m.strcmp(choice,'y'):
-			print('no meshing done ... exiting')
-			return None
+    if len(args) == 1:
+        resolution = args[0]
+        riftname = ''
+    if len(args) == 2:
+        riftname = args[0]
+        resolution = args[1]
 
-	area = resolution**2
+    #Check that mesh was not already run, and warn user:
+    if md.mesh.numberofelements:
+        choice = eval(input('This model already has a mesh. Are you sure you want to go ahead? (y/n)'))
+        if choice not in ['y', 'n']:
+            print("bad answer try you should use 'y' or 'n' ... exiting")
+            return None
+        if choice == 'n':
+            print('no meshing done ... exiting')
+            return None
 
-	#Check that file exist (this is a very very common mistake)
-	if not os.path.exists(domainname):
-		raise IOError("file '%s' not found" % domainname)
+    area = resolution**2
 
-	#Mesh using Triangle
-	md.mesh=mesh2d()
-	md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=Triangle_python(domainname,riftname,area)
-	md.mesh.elements=md.mesh.elements.astype(int)
-	md.mesh.segments=md.mesh.segments.astype(int)
-	md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
+    #Check that file exist (this is a very very common mistake)
+    if not os.path.exists(domainname):
+        raise IOError("file '%s' not found" % domainname)
 
-	#Fill in rest of fields:
-	md.mesh.numberofelements = np.size(md.mesh.elements,axis=0)
-	md.mesh.numberofvertices = np.size(md.mesh.x)
-	md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
+    #Mesh using Triangle
+    md.mesh = mesh2d()
+    md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.segments, md.mesh.segmentmarkers = Triangle_python(domainname, riftname, area)
+    md.mesh.elements = md.mesh.elements.astype(int)
+    md.mesh.segments = md.mesh.segments.astype(int)
+    md.mesh.segmentmarkers = md.mesh.segmentmarkers.astype(int)
 
-	#Now, build the connectivity tables for this mesh.
-	md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
-	md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
+    #Fill in rest of fields:
+    md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
+    md.mesh.numberofvertices = np.size(md.mesh.x)
+    md.mesh.vertexonboundary = np.zeros(md.mesh.numberofvertices, bool)
+    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = True
 
-	return md
+    #Now, build the connectivity tables for this mesh.
+    md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)[0]
+    md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)[0]
+
+    return md
Index: /issm/trunk-jpl/src/m/modules/BamgMesher.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/BamgMesher.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/modules/BamgMesher.py	(revision 24115)
@@ -1,19 +1,20 @@
 from BamgMesher_python import BamgMesher_python
 
-def BamgMesher(bamgmesh,bamggeom,bamgoptions):
-	"""
-	BAMGMESHER
 
-	Usage:
-		bamgmesh,bamggeom = BamgMesher(bamgmesh,bamggeom,bamgoptions);
+def BamgMesher(bamgmesh, bamggeom, bamgoptions):
+    """
+    BAMGMESHER
 
-	bamgmesh: input bamg mesh
-	bamggeom: input bamg geometry for the mesh
-	bamgoptions: options for the bamg mesh
-	"""
-	
-	#Call mex module
-	bamgmesh, bamggeom = BamgMesher_python(bamgmesh, bamggeom, bamgoptions)
+    Usage:
+            bamgmesh,bamggeom = BamgMesher(bamgmesh,bamggeom,bamgoptions);
 
-	#return
-	return bamgmesh, bamggeom
+    bamgmesh: input bamg mesh
+    bamggeom: input bamg geometry for the mesh
+    bamgoptions: options for the bamg mesh
+    """
+
+    #Call mex module
+    bamgmesh, bamggeom = BamgMesher_python(bamgmesh, bamggeom, bamgoptions)
+
+    #return
+    return bamgmesh, bamggeom
Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 24115)
@@ -1,6 +1,5 @@
-from cmaptools import getcolormap,truncate_colormap
+from cmaptools import getcolormap, truncate_colormap
 from plot_quiver import plot_quiver
-from scipy.interpolate import griddata
-import numpy as  np
+import numpy as np
 try:
     import matplotlib as mpl
@@ -10,72 +9,73 @@
     from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 except ImportError:
-    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
-
-def plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options,fig,axgrid,gridindex):
+    print("could not import pylab,  matplotlib has not been installed,  no plotting capabilities enabled")
+
+
+def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, axgrid, gridindex):
     """
-    PLOT_UNIT - unit plot, display data
+    PLOT_UNIT - unit plot,  display data
 
     Usage:
-    plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-
-    See also: PLOTMODEL, PLOT_MANAGER
+    plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options)
+
+    See also: PLOTMODEL,  PLOT_MANAGER
     """
     #if we are plotting 3d replace the current axis
     if not is2d:
         axgrid[gridindex].axis('off')
-        ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
-        ax.set_axis_bgcolor((0.7,0.7,0.7))
-    else:
-        ax=axgrid[gridindex]
+        ax = inset_locator.inset_axes(axgrid[gridindex], width='100%', height='100%', loc=3, borderpad=0, axes_class=Axes3D)
+        ax.set_axis_bgcolor((0.7, 0.7, 0.7))
+    else:
+        ax = axgrid[gridindex]
 
     #edgecolor
-    edgecolor=options.getfieldvalue('edgecolor','None')
+    edgecolor = options.getfieldvalue('edgecolor', 'None')
 
     # colormap
     # {{{ give number of colorlevels and transparency
-    colorlevels=options.getfieldvalue('colorlevels',128)
-    alpha=options.getfieldvalue('alpha',1)
-    if alpha<1:
-        antialiased=True
-    else:
-        antialiased=False
+    colorlevels = options.getfieldvalue('colorlevels', 128)
+    alpha = options.getfieldvalue('alpha', 1)
+    if alpha < 1:
+        antialiased = True
+    else:
+        antialiased = False
     # }}}
     # {{{ define wich colormap to use
     try:
-        defaultmap=plt.cm.get_cmap('viridis',colorlevels)
+        defaultmap = plt.cm.get_cmap('viridis', colorlevels)
     except AttributeError:
         print("Viridis can't be found (probably too old Matplotlib) reverting to gnuplot colormap")
-        defaultmap=truncate_colormap('gnuplot2',0.1,0.9,colorlevels)
+        defaultmap = truncate_colormap('gnuplot2', 0.1, 0.9, colorlevels)
     if not options.exist('colormap'):
-        cmap=defaultmap
-    else:
-        cmap=getcolormap(options)
+        cmap = defaultmap
+    else:
+        cmap = getcolormap(options)
     if options.exist('cmap_set_over'):
-        over=options.getfieldvalue('cmap_set_over','0.5')
+        over = options.getfieldvalue('cmap_set_over', '0.5')
         cmap.set_over(over)
     if options.exist('cmap_set_under'):
-        under=options.getfieldvalue('cmap_set_under','0.5')
+        under = options.getfieldvalue('cmap_set_under', '0.5')
         cmap.set_under(under)
-    options.addfield('colormap',cmap)
-    # }}}
-    # {{{ if plotting only one of several layers reduce dataset, same for surface
-    if options.getfieldvalue('layer',0)>=1:
-        plotlayer=options.getfieldvalue('layer',0)
-        if datatype==1:
-            slicesize=np.shape(elements)[0]
-        elif datatype in [2,3]:
-            slicesize=len(x)
-        data=data[(plotlayer-1)*slicesize:plotlayer*slicesize]
+    options.addfield('colormap', cmap)
+    # }}}
+    # {{{ if plotting only one of several layers reduce dataset,  same for surface
+    if options.getfieldvalue('layer', 0) >= 1:
+        plotlayer = options.getfieldvalue('layer', 0)
+        if datatype == 1:
+            slicesize = np.shape(elements)[0]
+        elif datatype in [2, 3]:
+            slicesize = len(x)
+        data = data[(plotlayer - 1) * slicesize:plotlayer * slicesize]
     # }}}
     # {{{ Get the colormap limits
     if options.exist('clim'):
-        lims=options.getfieldvalue('clim',[np.amin(data),np.amax(data)])
+        lims = options.getfieldvalue('clim', [np.amin(data), np.amax(data)])
     elif options.exist('caxis'):
-        lims=options.getfieldvalue('caxis',[np.amin(data),np.amax(data)])
-    else:
-        if np.amin(data)==np.amax(data):
-            lims=[np.amin(data)-0.5,np.amax(data)+0.5]
-        else:
-            lims=[np.amin(data),np.amax(data)]
+        lims = options.getfieldvalue('caxis', [np.amin(data), np.amax(data)])
+    else:
+        if np.amin(data) == np.amax(data):
+            lims = [np.amin(data) * 0.9, np.amax(data) * 1.1]
+        else:
+            lims = [np.amin(data), np.amax(data)]
     # }}}
     # {{{ Set the spread of the colormap (default is normal
@@ -88,5 +88,5 @@
     else:
         norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
-    options.addfield('colornorm',norm)
+    options.addfield('colornorm', norm)
     # }}}
 
@@ -94,32 +94,32 @@
     # {{{ data are on elements
 
-    if datatype==1:
+    if datatype == 1:
         if is2d:
             if options.exist('mask'):
-                triangles=mpl.tri.Triangulation(x,y,elements,data.mask)
+                triangles = mpl.tri.Triangulation(x, y, elements, data.mask)
             else:
-                triangles=mpl.tri.Triangulation(x,y,elements)
-            tri=ax.tripcolor(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,edgecolors=edgecolor)
+                triangles = mpl.tri.Triangulation(x, y, elements)
+            tri = ax.tripcolor(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=alpha, edgecolors=edgecolor)
         else:
             #first deal with colormap
             loccmap = plt.cm.ScalarMappable(cmap=cmap)
-            loccmap.set_array([min(data),max(data)])
-            loccmap.set_clim(vmin=min(data),vmax=max(data))
+            loccmap.set_array([min(data), max(data)])
+            loccmap.set_clim(vmin=min(data), vmax=max(data))
 
             #dealing with prism sides
-            recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T
-            eltind=np.arange(0,np.shape(elements)[0])
-            recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T))
-            eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0])))
-            recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T))
-            eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0])))
+            recface = np.vstack((elements[:, 0], elements[:, 1], elements[:, 4], elements[:, 3])).T
+            eltind = np.arange(0, np.shape(elements)[0])
+            recface = np.vstack((recface, np.vstack((elements[:, 1], elements[:, 2], elements[:, 5], elements[:, 4])).T))
+            eltind = np.hstack((eltind, np.arange(0, np.shape(elements)[0])))
+            recface = np.vstack((recface, np.vstack((elements[:, 2], elements[:, 0], elements[:, 3], elements[:, 5])).T))
+            eltind = np.hstack((eltind, np.arange(0, np.shape(elements)[0])))
             tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
             _, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
-            recel= recface[idx[np.where(recur==1)]]
-            recindex=eltind[idx[np.where(recur==1)]]
-            for i,rectangle in enumerate(recel):
-                rec=list(zip(x[rectangle],y[rectangle],z[rectangle]))
-                pl3=Poly3DCollection([rec])
-                color=loccmap.to_rgba(data[recindex[i]])
+            recel = recface[idx[np.where(recur == 1)]]
+            recindex = eltind[idx[np.where(recur == 1)]]
+            for i, rectangle in enumerate(recel):
+                rec = list(zip(x[rectangle], y[rectangle], z[rectangle]))
+                pl3 = Poly3DCollection([rec])
+                color = loccmap.to_rgba(data[recindex[i]])
                 pl3.set_edgecolor(color)
                 pl3.set_color(color)
@@ -127,65 +127,65 @@
 
             #dealing with prism bases
-            triface=np.vstack((elements[:,0:3],elements[:,3:6]))
-            eltind=np.arange(0,np.shape(elements)[0])
-            eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0])))
+            triface = np.vstack((elements[:, 0:3], elements[:, 3:6]))
+            eltind = np.arange(0, np.shape(elements)[0])
+            eltind = np.hstack((eltind, np.arange(0, np.shape(elements)[0])))
             tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1])))
-            _, idx,recur = np.unique(tmp, return_index=True,return_counts=True)
+            _, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
             #we keep only top and bottom elements
-            triel= triface[idx[np.where(recur==1)]]
-            triindex=eltind[idx[np.where(recur==1)]]
-            for i,triangle in enumerate(triel):
-                tri=list(zip(x[triangle],y[triangle],z[triangle]))
-                pl3=Poly3DCollection([tri])
-                color=loccmap.to_rgba(data[triindex[i]])
-                pl3.set_edgecolor(color)
-                pl3.set_color(color)
-                ax.add_collection3d(pl3)
-    
-            ax.set_xlim([min(x),max(x)])
-            ax.set_ylim([min(y),max(y)])
-            ax.set_zlim([min(z),max(z)])
+            triel = triface[idx[np.where(recur == 1)]]
+            triindex = eltind[idx[np.where(recur == 1)]]
+            for i, triangle in enumerate(triel):
+                tri = list(zip(x[triangle], y[triangle], z[triangle]))
+                pl3 = Poly3DCollection([tri])
+                color = loccmap.to_rgba(data[triindex[i]])
+                pl3.set_edgecolor(color)
+                pl3.set_color(color)
+                ax.add_collection3d(pl3)
+
+            ax.set_xlim([min(x), max(x)])
+            ax.set_ylim([min(y), max(y)])
+            ax.set_zlim([min(z), max(z)])
 
             #raise ValueError('plot_unit error: 3D element plot not supported yet')
-        return 
+        return
 
     # }}}
     # {{{ data are on nodes
 
-    elif datatype==2:
+    elif datatype == 2:
         if is2d:
             if np.ma.is_masked(data):
-                EltMask=np.asarray([np.any(np.in1d(index,np.where(data.mask))) for index in elements])
-                triangles=mpl.tri.Triangulation(x,y,elements,EltMask)
+                EltMask = np.asarray([np.any(np.in1d(index, np.where(data.mask))) for index in elements])
+                triangles = mpl.tri.Triangulation(x, y, elements, EltMask)
             else:
-                triangles=mpl.tri.Triangulation(x,y,elements)
-                #tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha)
+                triangles = mpl.tri.Triangulation(x, y, elements)
+                #tri = ax.tricontourf(triangles, data, colorlevels, cmap = cmap, norm = norm, alpha = alpha)
             if options.exist('log'):
-                if alpha<1:#help with antialiasing
-                    tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased)
-                tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,antialiased=antialiased)
+                if alpha < 1:  #help with antialiasing
+                    tri = ax.tricontour(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=0.1, antialiased=antialiased)
+                tri = ax.tricontourf(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=alpha, antialiased=antialiased)
             else:
-                if alpha<1:#help with antialiasing
-                    tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased)
-                tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both',antialiased=antialiased)
+                if alpha < 1:  #help with antialiasing
+                    tri = ax.tricontour(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=0.1, antialiased=antialiased)
+                tri = ax.tricontourf(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=alpha, extend='both', antialiased=antialiased)
             if edgecolor != 'None':
-                ax.triplot(x,y,elements,color=edgecolor)
+                ax.triplot(x, y, elements, color=edgecolor)
         else:
             #first deal with the colormap
             loccmap = plt.cm.ScalarMappable(cmap=cmap)
-            loccmap.set_array([min(data),max(data)])
-            loccmap.set_clim(vmin=min(data),vmax=max(data))
+            loccmap.set_array([min(data), max(data)])
+            loccmap.set_clim(vmin=min(data), vmax=max(data))
 
             #deal with prism sides
-            recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T
-            recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T))
-            recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T))
+            recface = np.vstack((elements[:, 0], elements[:, 1], elements[:, 4], elements[:, 3])).T
+            recface = np.vstack((recface, np.vstack((elements[:, 1], elements[:, 2], elements[:, 5], elements[:, 4])).T))
+            recface = np.vstack((recface, np.vstack((elements[:, 2], elements[:, 0], elements[:, 3], elements[:, 5])).T))
             tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
             _, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
-            recel= recface[idx[np.where(recur==1)]]
+            recel = recface[idx[np.where(recur == 1)]]
             for rectangle in recel:
-                rec=list(zip(x[rectangle],y[rectangle],z[rectangle]))
-                pl3=Poly3DCollection([rec])
-                color=loccmap.to_rgba(np.mean(data[rectangle]))
+                rec = list(zip(x[rectangle], y[rectangle], z[rectangle]))
+                pl3 = Poly3DCollection([rec])
+                color = loccmap.to_rgba(np.mean(data[rectangle]))
                 pl3.set_edgecolor(color)
                 pl3.set_color(color)
@@ -193,20 +193,20 @@
 
             #deal with prism faces
-            triface=np.vstack((elements[:,0:3],elements[:,3:6]))
+            triface = np.vstack((elements[:, 0:3], elements[:, 3:6]))
             tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1])))
-            _, idx,recur = np.unique(tmp, return_index=True,return_counts=True)
+            _, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
             #we keep only top and bottom elements
-            triel= triface[idx[np.where(recur==1)]]
+            triel = triface[idx[np.where(recur == 1)]]
             for triangle in triel:
-                tri=list(zip(x[triangle],y[triangle],z[triangle]))
-                pl3=Poly3DCollection([tri])
-                color=loccmap.to_rgba(np.mean(data[triangle]))
-                pl3.set_edgecolor(color)
-                pl3.set_color(color)
-                ax.add_collection3d(pl3)
-
-            ax.set_xlim([min(x),max(x)])
-            ax.set_ylim([min(y),max(y)])
-            ax.set_zlim([min(z),max(z)])
+                tri = list(zip(x[triangle], y[triangle], z[triangle]))
+                pl3 = Poly3DCollection([tri])
+                color = loccmap.to_rgba(np.mean(data[triangle]))
+                pl3.set_edgecolor(color)
+                pl3.set_color(color)
+                ax.add_collection3d(pl3)
+
+            ax.set_xlim([min(x), max(x)])
+            ax.set_ylim([min(y), max(y)])
+            ax.set_zlim([min(z), max(z)])
             #raise ValueError('plot_unit error: 3D element plot not supported yet')
         return
@@ -214,15 +214,15 @@
     # }}}
     # {{{ plotting quiver
-    elif datatype==3:
+    elif datatype == 3:
         if is2d:
-            Q=plot_quiver(x,y,data,options,ax)
+            Q = plot_quiver(x, y, data, options, ax)
         else:
             raise ValueError('plot_unit error: 3D node plot not supported yet')
         return
-    
+
     # }}}
     # {{{ plotting P1 Patch (TODO)
 
-    elif datatype==4:
+    elif datatype == 4:
         print('plot_unit message: P1 patch plot not implemented yet')
         return
@@ -231,5 +231,5 @@
     # {{{ plotting P0 Patch (TODO)
 
-    elif datatype==5:
+    elif datatype == 5:
         print('plot_unit message: P0 patch plot not implemented yet')
         return
@@ -237,3 +237,3 @@
     # }}}
     else:
-        raise ValueError('datatype=%d not supported' % datatype)
+        raise ValueError('datatype = %d not supported' % datatype)
Index: /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py	(revision 24115)
@@ -7,84 +7,84 @@
 
 def loadresultsfromcluster(md,runtimename=False):
-	"""
-	LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster
- 
-	   Usage:
-	      md=loadresultsfromcluster(md,runtimename);
-	"""
+        """
+        LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster
 
-	#retrieve cluster, to be able to call its methods
-	cluster=md.cluster
+           Usage:
+              md=loadresultsfromcluster(md,runtimename);
+        """
 
-	if runtimename:
-		md.private.runtimename=runtimename
+        #retrieve cluster, to be able to call its methods
+        cluster=md.cluster
 
-	#Download outputs from the cluster
-	filelist=[md.miscellaneous.name+'.outlog',md.miscellaneous.name+'.errlog']
-	if md.qmu.isdakota:
-		filelist.append(md.miscellaneous.name+'.qmu.err')
-		filelist.append(md.miscellaneous.name+'.qmu.out')
-		if 'tabular_graphics_data' in fieldnames(md.qmu.params):
-			if md.qmu.params.tabular_graphics_data:
-				filelist.append('dakota_tabular.dat')
-	else:
-		filelist.append(md.miscellaneous.name+'.outbin')
-	cluster.Download(md.private.runtimename,filelist)
+        if runtimename:
+                md.private.runtimename=runtimename
 
-	#If we are here, no errors in the solution sequence, call loadresultsfromdisk.
-	if os.path.exists(md.miscellaneous.name+'.outbin'):
-		if os.path.getsize(md.miscellaneous.name+'.outbin')>0:
-			md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin')
-		else:
-			print(('WARNING, outbin file is empty for run '+md.miscellaneous.name))
-	elif not md.qmu.isdakota:
-		print(('WARNING, outbin file does not exist '+md.miscellaneous.name))
-		
-	#erase the log and output files
-	if md.qmu.isdakota:
-		#filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name)
-		filename = md.miscellaneous.name
+        #Download outputs from the cluster
+        filelist=[md.miscellaneous.name+'.outlog',md.miscellaneous.name+'.errlog']
+        if md.qmu.isdakota:
+                filelist.append(md.miscellaneous.name+'.qmu.err')
+                filelist.append(md.miscellaneous.name+'.qmu.out')
+                if 'tabular_graphics_data' in fieldnames(md.qmu.params):
+                        if md.qmu.params.tabular_graphics_data:
+                                filelist.append('dakota_tabular.dat')
+        else:
+                filelist.append(md.miscellaneous.name+'.outbin')
+        cluster.Download(md.private.runtimename,filelist)
 
-		# this will not work like normal as dakota doesn't generate outbin files,
-		#   instead calls postqmu to store results directly in the model
-		#   at md.results.dakota via dakota_out_parse
-		md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin')
-	else:
-		filename=md.miscellaneous.name
-		TryRem('.outbin',filename)
-		if not platform.system()=='Windows':
-			TryRem('.tar.gz',md.private.runtimename)
+        #If we are here, no errors in the solution sequence, call loadresultsfromdisk.
+        if os.path.exists(md.miscellaneous.name+'.outbin'):
+                if os.path.getsize(md.miscellaneous.name+'.outbin')>0:
+                        md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin')
+                else:
+                        print(('WARNING, outbin file is empty for run '+md.miscellaneous.name))
+        elif not md.qmu.isdakota:
+                print(('WARNING, outbin file does not exist '+md.miscellaneous.name))
 
-	TryRem('.errlog',filename)
-	TryRem('.outlog',filename)
-	
-	#erase input file if run was carried out on same platform.
-	hostname=socket.gethostname()
-	if hostname==cluster.name:
-		if md.qmu.isdakota:
-			#filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name)
-			filename = md.miscellaneous.name
-			TryRem('.queue',filename)
-		else:
-			filename=md.miscellaneous.name
-			TryRem('.toolkits',filename)
-			if not platform.system()=='Windows':
-				TryRem('.queue',filename)
-			else:
-				TryRem('.bat',filename)
+        #erase the log and output files
+        if md.qmu.isdakota:
+                #filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name)
+                filename = md.miscellaneous.name
 
-		# remove this for bin file debugging
-		TryRem('.bin',filename)
+                # this will not work like normal as dakota doesn't generate outbin files,
+                #   instead calls postqmu to store results directly in the model
+                #   at md.results.dakota via dakota_out_parse
+                md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin')
+        else:
+                filename=md.miscellaneous.name
+                TryRem('.outbin',filename)
+                if not platform.system()=='Windows':
+                        TryRem('.tar.gz',md.private.runtimename)
 
-	#cwd = os.getcwd().split('/')[-1]
-	if md.qmu.isdakota:
-		os.chdir('..')
-		#TryRem('',cwd)
+        TryRem('.errlog',filename)
+        TryRem('.outlog',filename)
 
-	return md
+        #erase input file if run was carried out on same platform.
+        hostname=socket.gethostname()
+        if hostname==cluster.name:
+                if md.qmu.isdakota:
+                        #filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name)
+                        filename = md.miscellaneous.name
+                        TryRem('.queue',filename)
+                else:
+                        filename=md.miscellaneous.name
+                        TryRem('.toolkits',filename)
+                        if not platform.system()=='Windows':
+                                TryRem('.queue',filename)
+                        else:
+                                TryRem('.bat',filename)
+
+                # remove this for bin file debugging
+                TryRem('.bin',filename)
+
+        #cwd = os.getcwd().split('/')[-1]
+        if md.qmu.isdakota:
+                os.chdir('..')
+                #TryRem('',cwd)
+
+        return md
 
 def TryRem(extension,filename):
-	try:
-		os.remove(filename+extension)
-	except OSError:
-		print(('WARNING, no '+extension+'  is present for run '+filename))
+        try:
+                os.remove(filename+extension)
+        except OSError:
+                print(('WARNING, no '+extension+'  is present for run '+filename))
Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 24115)
@@ -4,59 +4,60 @@
 from postqmu import postqmu
 
+
 def loadresultsfromdisk(md,filename):
-	"""
-	LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename"
+    """
+    LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename"
 
-	   Usage:
-	      md=loadresultsfromdisk(md=False,filename=False);
-	"""
+       Usage:
+          md=loadresultsfromdisk(md=False,filename=False);
+    """
 
-	#check number of inputs/outputs
-	if not md or not filename:
-		raise ValueError("loadresultsfromdisk: error message.")
+    #check number of inputs/outputs
+    if not md or not filename:
+        raise ValueError("loadresultsfromdisk: error message.")
 
-	if not md.qmu.isdakota:
+    if not md.qmu.isdakota:
 
-		#Check that file exists
-		if not os.path.exists(filename):
-			raise OSError("binary file '{}' not found.".format(filename))
+        #Check that file exists
+        if not os.path.exists(filename):
+            raise OSError("binary file '{}' not found.".format(filename))
 
-		#initialize md.results if not a structure yet
-		if not isinstance(md.results,results):
-			md.results=results()
+        #initialize md.results if not a structure yet
+        if not isinstance(md.results,results):
+            md.results=results()
 
-		#load results onto model
-		structure=parseresultsfromdisk(md,filename,not md.settings.io_gather)
-		if not len(structure):
-			raise RuntimeError("No result found in binary file '{}'. Check for solution crash.".format(filename))
+        #load results onto model
+        structure=parseresultsfromdisk(md,filename,not md.settings.io_gather)
+        if not len(structure):
+            raise RuntimeError("No result found in binary file '{}'. Check for solution crash.".format(filename))
 
-		setattr(md.results,structure[0].SolutionType,structure)
+        setattr(md.results,structure[0].SolutionType,structure)
 
-		#recover solution_type from results
-		md.private.solution=structure[0].SolutionType
+        #recover solution_type from results
+        md.private.solution=structure[0].SolutionType
 
-		#read log files onto fields
-		if os.path.exists(md.miscellaneous.name+'.errlog'):
-			with open(md.miscellaneous.name+'.errlog','r') as f:
-				setattr(getattr(md.results,structure[0].SolutionType)[0],'errlog',[line[:-1] for line in f])
-		else:
-			setattr(getattr(md.results,structure[0].SolutionType)[0],'errlog',[])
+        #read log files onto fields
+        if os.path.exists(md.miscellaneous.name+'.errlog'):
+                with open(md.miscellaneous.name+'.errlog','r') as f:
+                        setattr(getattr(md.results,structure[0].SolutionType)[0],'errlog',[line[:-1] for line in f])
+        else:
+                setattr(getattr(md.results,structure[0].SolutionType)[0],'errlog',[])
 
-		if os.path.exists(md.miscellaneous.name+'.outlog'):
-			with open(md.miscellaneous.name+'.outlog','r') as f:
-				setattr(getattr(md.results,structure[0].SolutionType)[0],'outlog',[line[:-1] for line in f])
-		else:
-			setattr(getattr(md.results,structure[0].SolutionType)[0],'outlog',[])
+        if os.path.exists(md.miscellaneous.name+'.outlog'):
+                with open(md.miscellaneous.name+'.outlog','r') as f:
+                        setattr(getattr(md.results,structure[0].SolutionType)[0],'outlog',[line[:-1] for line in f])
+        else:
+                setattr(getattr(md.results,structure[0].SolutionType)[0],'outlog',[])
 
-		if len(getattr(md.results,structure[0].SolutionType)[0].errlog):
-			print ("loadresultsfromcluster info message: error during solution. Check your errlog and outlog model fields.")
+        if len(getattr(md.results,structure[0].SolutionType)[0].errlog):
+                print ("loadresultsfromcluster info message: error during solution. Check your errlog and outlog model fields.")
 
-		#if only one solution, extract it from list for user friendliness
-		if len(structure) == 1 and not structure[0].SolutionType=='TransientSolution':
-			setattr(md.results,structure[0].SolutionType,structure[0])
+        #if only one solution, extract it from list for user friendliness
+        if len(structure) == 1 and not structure[0].SolutionType=='TransientSolution':
+                setattr(md.results,structure[0].SolutionType,structure[0])
 
-	#post processes qmu results if necessary
-	else:
-		md=postqmu(md,filename)
+    #post processes qmu results if necessary
+    else:
+        md=postqmu(md,filename)
 
-	return md
+    return md
Index: /issm/trunk-jpl/src/m/solve/marshall.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/marshall.py	(revision 24114)
+++ /issm/trunk-jpl/src/m/solve/marshall.py	(revision 24115)
@@ -1,44 +1,44 @@
 from WriteData import WriteData
 
+
 def marshall(md):
-	"""
-	MARSHALL - outputs a compatible binary file from @model md, for certain solution type.
+    """
+    MARSHALL - outputs a compatible binary file from @model md, for certain solution type.
 
-	   The routine creates a compatible binary file from @model md
-	   This binary file will be used for parallel runs in JPL-package
+       The routine creates a compatible binary file from @model md
+       This binary file will be used for parallel runs in JPL-package
 
-	   Usage:
-	      marshall(md)
-	"""
-	if md.verbose.solution:
-		print(("marshalling file '%s.bin'." % md.miscellaneous.name))
+       Usage:
+          marshall(md)
+    """
+    if md.verbose.solution:
+        print(("marshalling file '%s.bin'." % md.miscellaneous.name))
 
-	#open file for binary writing
-	try:
-		fid=open(md.miscellaneous.name+'.bin','wb')
-	except IOError as e:
-		raise IOError("marshall error message: could not open '%s.bin' file for binary writing." % md.miscellaneous.name)
+    #open file for binary writing
+    try:
+        fid = open(md.miscellaneous.name + '.bin', 'wb')
+    except IOError as e:
+        raise IOError("marshall error message: could not open '%s.bin' file for binary writing." % md.miscellaneous.name)
 
-	for field in md.properties():
+    for field in md.properties():
 
-		#Some properties do not need to be marshalled
-		if field in ['results','radaroverlay','toolkits','cluster','private']:
-			continue
+        #Some properties do not need to be marshalled
+        if field in ['results', 'radaroverlay', 'toolkits', 'cluster', 'private']:
+            continue
 
-		#Check that current field is an object
-		if not hasattr(getattr(md,field),'marshall'):
-			raise TypeError("field '%s' is not an object." % field)
+        #Check that current field is an object
+        if not hasattr(getattr(md, field), 'marshall'):
+            raise TypeError("field '{}' is not an object.".format(field))
 
-		#Marshall current object
-		#print "marshalling %s ..." % field
-		exec("md.{}.marshall('md.{}',md,fid)".format(field,field))
+        #Marshall current object
+        #print "marshalling %s ..." % field
+        exec("md.{}.marshall('md.{}',md,fid)".format(field, field))
 
-	#Last, write "md.EOF" to make sure that the binary file is not corrupt
-	WriteData(fid,'XXX','name','md.EOF','data',True,'format','Boolean');
+    #Last, write "md.EOF" to make sure that the binary file is not corrupt
+    WriteData(fid, 'XXX', 'name', 'md.EOF', 'data', True, 'format', 'Boolean')
 
-	#close file
-	try:
-		fid.close()
-	except IOError as e:
-		raise IOError("marshall error message: could not close file '%s.bin'." % md.miscellaneous.name)
-
+    #close file
+    try:
+        fid.close()
+    except IOError as e:
+        raise IOError("marshall error message: could not close file '%s.bin'." % md.miscellaneous.name)