Index: issm/trunk-jpl/src/m/classes/giamme.py
===================================================================
--- issm/trunk-jpl/src/m/classes/giamme.py	(revision 25164)
+++ issm/trunk-jpl/src/m/classes/giamme.py	(revision 25165)
@@ -42,5 +42,5 @@
 
     def checkconsistency(self, md, solution, analyses): # {{{
-        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.istidallove == 0):
+        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
             return md
 
Index: issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- issm/trunk-jpl/src/m/classes/model.py	(revision 25164)
+++ issm/trunk-jpl/src/m/classes/model.py	(revision 25165)
@@ -688,6 +688,5 @@
 
         md.materials.extrude(md)
-        if md.damage.isdamage == 1:
-            md.damage.extrude(md)
+        md.damage.extrude(md)
         md.mask.extrude(md)
         md.qmu.extrude(md)
Index: issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py
===================================================================
--- issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py	(revision 25164)
+++ issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py	(revision 25165)
@@ -1,53 +1,15 @@
-def AnalysisConfiguration(solutiontype):  #{{{
-    """
-    ANALYSISCONFIGURATION - return type of analyses, number of analyses
-
-            Usage:
-                    [analyses] = AnalysisConfiguration(solutiontype)
-    """
-
-    if solutiontype == 'StressbalanceSolution':
-        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis']
-    elif solutiontype == 'SteadystateSolution':
-        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis']
-    elif solutiontype == 'ThermalSolution':
-        analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
-    elif solutiontype == 'MasstransportSolution':
-        analyses = ['MasstransportAnalysis']
-    elif solutiontype == 'BalancethicknessSolution':
-        analyses = ['BalancethicknessAnalysis']
-    elif solutiontype == 'SurfaceSlopeSolution':
-        analyses = ['L2ProjectionBaseAnalysis']
-    elif solutiontype == 'BalancevelocitySolution':
-        analyses = ['BalancevelocityAnalysis']
-    elif solutiontype == 'BedSlopeSolution':
-        analyses = ['L2ProjectionBaseAnalysis']
-    elif solutiontype == 'GiaSolution':
-        analyses = ['GiaIvinsAnalysis']
-    elif solutiontype == 'LoveSolution':
-        analyses = ['LoveAnalysis']
-    elif solutiontype == 'TransientSolution':
-        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
-    elif solutiontype == 'HydrologySolution':
-        analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
-    elif 'DamageEvolutionSolution':
-        analyses = ['DamageEvolutionAnalysis']
-
-    else:
-        raise TypeError("solution type: '%s' not supported yet!" % solutiontype)
-    return analyses
-    #}}}
-
-
-def ismodelselfconsistent(md):
-    """
+def ismodelselfconsistent(md): #{{{
+    '''
     ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
 
        Usage:
           ismodelselfconsistent(md),
-    """
+    '''
 
     #initialize consistency as true
     md.private.isconsistent = True
+
+    # print(md.mesh.z)
+    # exit()
 
     #Get solution and associated analyses
@@ -64,5 +26,5 @@
         #Check that current field is an object
         if not hasattr(getattr(md, field), 'checkconsistency'):
-            md.checkmessage("field '%s' is not an object." % field)
+            md.checkmessage('field {} is not an object.'.format(field))
 
         #Check consistency of the object
@@ -72,2 +34,51 @@
     if not md.private.isconsistent:
         raise RuntimeError('Model not consistent, see messages above.')
+#}}}
+
+def AnalysisConfiguration(solutiontype): #{{{
+    '''
+    ANALYSISCONFIGURATION - return type of analyses, number of analyses
+
+        Usage:
+                [analyses] = AnalysisConfiguration(solutiontype)
+    '''
+
+    if solutiontype == 'StressbalanceSolution':
+        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis']
+    elif solutiontype == 'SteadystateSolution':
+        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis']
+    elif solutiontype == 'ThermalSolution':
+        analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
+    elif solutiontype == 'MasstransportSolution':
+        analyses = ['MasstransportAnalysis']
+    elif solutiontype == 'BalancethicknessSolution':
+        analyses = ['BalancethicknessAnalysis']
+    elif solutiontype == 'Balancethickness2Solution':
+        analyses = ['Balancethickness2Analysis']
+    elif solutiontype == 'BalancethicknessSoftSolution':
+        analyses = ['BalancethicknessAnalysis']
+    elif solutiontype == 'BalancevelocitySolution':
+        analyses = ['BalancevelocityAnalysis']
+    elif solutiontype == 'SurfaceSlopeSolution':
+        analyses = ['L2ProjectionBaseAnalysis']
+    elif solutiontype == 'BedSlopeSolution':
+        analyses = ['L2ProjectionBaseAnalysis']
+    elif solutiontype == 'GiaSolution':
+        analyses = ['GiaIvinsAnalysis']
+    elif solutiontype == 'LoveSolution':
+        analyses = ['LoveAnalysis']
+    elif solutiontype == 'EsaSolution':
+        analyses = ['EsaAnalysis']
+    elif solutiontype == 'TransientSolution':
+        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
+    elif solutiontype == 'SealevelriseSolution':
+        analyses = ['SealevelriseAnalysis']
+    elif solutiontype == 'HydrologySolution':
+        analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
+    elif 'DamageEvolutionSolution':
+        analyses = ['DamageEvolutionAnalysis']
+    else:
+        raise TypeError('solution type: {} not supported yet!'.format(solutiontype))
+
+    return analyses
+#}}}
Index: issm/trunk-jpl/src/m/extrusion/project3d.m
===================================================================
--- issm/trunk-jpl/src/m/extrusion/project3d.m	(revision 25164)
+++ issm/trunk-jpl/src/m/extrusion/project3d.m	(revision 25165)
@@ -38,8 +38,7 @@
 
 elseif strcmpi(type,'node'),
-
 	%Initialize 3d vector
 	if size(vector2d,1)==md.mesh.numberofvertices2d
-		projected_vector=paddingvalue*ones(md.mesh.numberofvertices,  size(vector2d,2));
+		projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2));
 	elseif size(vector2d,1)==md.mesh.numberofvertices2d+1
 		projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2));
Index: issm/trunk-jpl/src/m/extrusion/project3d.py
===================================================================
--- issm/trunk-jpl/src/m/extrusion/project3d.py	(revision 25164)
+++ issm/trunk-jpl/src/m/extrusion/project3d.py	(revision 25165)
@@ -1,26 +1,32 @@
 import numpy as np
+
 from pairoptions import pairoptions
 
 
 def project3d(md, *args):
-    """
+    '''
     PROJECT3D - vertically project a vector from 2d mesh
 
-       vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
-       This vector can be a node vector of size (md.mesh.numberofvertices2d, N / A) or an
-       element vector of size (md.mesh.numberofelements2d, N / A).
-       arguments:
-          'vector': 2d vector
-          'type': 'element' or 'node'.
-       options:
-          'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
-                 value of the 2d vector.
-          'padding': default to 0 (value adopted by other 3d layers not being projected
+        vertically project a vector from 2d mesh (split in noncoll and coll 
+        areas) into a 3d mesh.
+        This vector can be a node vector of size (md.mesh.numberofvertices2d, 
+        N/A) or an element vector of size (md.mesh.numberofelements2d, N/A).
 
-       Examples:
-          extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
-          extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0)
-          extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node')
-    """
+        arguments:
+            'vector': 2d vector
+            'type': 'element' or 'node'
+
+        options:
+            'layer'     a layer number where vector should keep its values. If 
+                        not specified, all layers adopt the value of the 2d 
+                        vector.
+            'padding':  default to 0 (value adopted by other 3d layers not 
+                        being projected.
+
+        Examples:
+            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
+            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0)
+            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node')
+    '''
 
     #some regular checks
@@ -37,10 +43,9 @@
     paddingvalue = options.getfieldvalue('padding', 0)  #0 by default
 
-    vector1d = False
-    if isinstance(vector2d, np.ndarray) and np.ndim(vector2d) == 1:
-        vector1d = True
-        vector2d = vector2d.reshape(-1, )
+    #Handle special case where vector2d is single element (differs from representation in MATLAB)
+    if isinstance(vector2d, (bool, int, float)):
+        projected_vector = vector2d
 
-    if isinstance(vector2d, (bool, int, float)) or np.size(vector2d) == 1:
+    if np.size(vector2d) == 1:
         projected_vector = vector2d
 
@@ -56,5 +61,5 @@
             else:
                 raise TypeError("vector length not supported")
-    #Fill in
+            #Fill in
             if layer == 0:
                 for i in range(md.mesh.numberoflayers):
@@ -71,5 +76,5 @@
             else:
                 raise TypeError("vector length not supported")
-    #Fill in
+            #Fill in
             if layer == 0:
                 for i in range(md.mesh.numberoflayers):
@@ -89,5 +94,5 @@
             else:
                 raise TypeError("vector length not supported")
-    #Fill in
+            #Fill in
             if layer == 0:
                 for i in range(md.mesh.numberoflayers - 1):
@@ -104,5 +109,5 @@
             else:
                 raise TypeError("vector length not supported")
-    #Fill in
+            #Fill in
             if layer == 0:
                 for i in range(md.mesh.numberoflayers - 1):
@@ -114,6 +119,3 @@
         raise TypeError("project3d error message: unknown projection type")
 
-    if vector1d:
-        projected_vector = projected_vector.reshape(-1, )
-
     return projected_vector
Index: issm/trunk-jpl/src/m/partition/partitioner.m
===================================================================
--- issm/trunk-jpl/src/m/partition/partitioner.m	(revision 25164)
+++ issm/trunk-jpl/src/m/partition/partitioner.m	(revision 25165)
@@ -14,4 +14,5 @@
 %   Usage:
 %      partitionvector=partitioner(md,'package','chaco','npart',100,'weighting','on');
+%      [partitionvector,md]=partitioner(md,'package','chaco','npart',100,'weighting','on');
 %
 
@@ -36,5 +37,5 @@
 	%partitioning essentially happens in 2D. So partition in 2D, then 
 	%extrude the partition vector vertically. 
-	md3d=md; %save  for later
+	md3d=md; %save for later
 	md.mesh.elements=md.mesh.elements2d;
 	md.mesh.x=md.mesh.x2d;
@@ -78,5 +79,5 @@
 		%  partition into nparts
 		if isa(md.mesh,'mesh2d'),
-			part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis.
+			part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x,md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis.
 		else
 			part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[])'+1; %index partitions from 1 up. like metis.
Index: issm/trunk-jpl/src/m/partition/partitioner.py
===================================================================
--- issm/trunk-jpl/src/m/partition/partitioner.py	(revision 25164)
+++ issm/trunk-jpl/src/m/partition/partitioner.py	(revision 25165)
@@ -1,7 +1,8 @@
+import copy
+
 import numpy as np
 
 from adjacency import *
 from Chaco import *
-import copy
 from mesh2d import *
 from MeshPartition import *
@@ -10,41 +11,40 @@
 
 
-def partitioner(md, *varargin):
+def partitioner(md, *args):
     '''
     PARTITIONER - partition mesh
 
-    List of options to partitioner:
+        List of options to partitioner:
+            package: 'chaco', 'metis', or 'scotch'
+            npart: number of partitions
+            weighting: 'on' or 'off': default off
+            section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
+            recomputeadjacency: 'on' by default (set to 'off' to compute existing one)
+            type: 'node' or 'element' partition vector (default to 'node')
+            Output: partitionvector: the partition vector
 
-    package: 'chaco', 'metis', or 'scotch'
-    npart: number of partitions
-    weighting: 'on' or 'off': default off
-    section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
-    recomputeadjacency: 'on' by default (set to 'off' to compute existing one)
-    type: 'node' or 'element' partition vector (default to 'node')
-    Output: partitionvector: the partition vector
-
-    Usage:
-        partitionvector = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on')
+        Usage:
+            partitionvector, md = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on')
     '''
 
     #get options:
-    options = pairoptions(*varargin)
+    options = pairoptions(*args)
 
-    #get options:
+    #set defaults
+    package = options.getfieldvalue('package', 'chaco')
+    npart = options.getfieldvalue('npart', 10)
+    weighting = options.getfieldvalue('weighting', 'on')
     section = options.getfieldvalue('section', 1)
-    weighting = options.getfieldvalue('weighting', 'on')
-    package = options.getfieldvalue('package', 'chaco')  #default to chaco
-    npart = options.getfieldvalue('npart', 10)  # default to 10
-    recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on')  # default to on
-    vectortype = options.getfieldvalue('type', 'node')  #default to node
+    recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on')
+    vectortype = options.getfieldvalue('type', 'node')
 
     # Python only: short-circuit
     if vectortype == 'element' and not package == 'linear':
-        raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package)
+        raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package))
 
-    if(md.mesh.dimension() == 3):
+    if md.mesh.dimension() == 3:
         #partitioning essentially happens in 2D. So partition in 2D, then
         #extrude the partition vector vertically.
-        md3d = copy.deepcopy(md)
+        md3d = copy.deepcopy(md) # save for later
         md.mesh.elements = md.mesh.elements2d
         md.mesh.x = md.mesh.x2d
@@ -64,8 +64,8 @@
     if package == 'chaco':
         if vectortype == 'element':
-            raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package)
+            raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package))
         else:
-            #  default method (from chaco.m)
-            method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321])
+            # default method (from chaco.m)
+            method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321]).conj().transpose()
             method[0] = 3  #  global method (3 = inertial (geometric))
             method[2] = 0  #  vertex weights (0 = off, 1 = on)
@@ -81,12 +81,11 @@
                 weights = []
 
-            method = method.reshape(-1, 1)  # transpose to 1x10 instead of 10
+            #method = method.reshape(-1, 1)  # transpose to 1x10 instead of 10
 
             #  partition into nparts
             if isinstance(md.mesh, mesh2d):
-                part = np.array(Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, np.zeros((md.mesh.numberofvertices, )), method, npart, np.array([]))).T + 1  #index partitions from 1 up. like metis.
+                part = Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, np.zeros((md.mesh.numberofvertices, 1)), method, npart, np.array([]))[0].conj().transpose() + 1 #index partitions from 1 up. like metis.
             else:
-                part = np.array(Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, md.mesh.z, method, npart, np.array([]))).T + 1  #index partitions from 1 up. like metis.
-
+                part = Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, md.mesh.z, method, npart, np.array([]))[0].conj().transpose() + 1 #index partitions from 1 up. like metis.
     elif package == 'scotch':
         if vectortype == 'element':
@@ -116,5 +115,5 @@
 
     else:
-        raise RuntimeError('partitioner error message: could not find %s partitioner' % package)
+        raise RuntimeError('partitioner error message: could not find {} partitioner'.format(package))
 
     #extrude if we are in 3D:
@@ -122,11 +121,12 @@
         md3d.qmu.vertex_weight = md.qmu.vertex_weight
         md3d.qmu.adjacency = md.qmu.adjacency
-        md = md3d
+        md = copy.deepcopy(md3d)
         if vectortype == 'element':
-            part = project3d(md, 'vector', np.squeeze(part), 'type', 'element')
+            part = project3d(md, 'vector', part.conj().transpose(), 'type', 'element')
         else:
-            part = project3d(md, 'vector', np.squeeze(part), 'type', 'node')
+            part = project3d(md, 'vector', part.conj().transpose(), 'type', 'node')
 
-        part = part.reshape(-1, 1)
+    if part.shape[0] == 1:
+        part = part.conj().transpose()
 
     # Output
