Changeset 19048

Timestamp:

01/29/15 08:44:39 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: extrude is now a method of each class so that model.m does not need to know what class is being used

Location:

issm/trunk-jpl/src/m

Files:

• classes/basalforcings.m (modified) (1 diff)
• classes/basalforcings.py (modified) (2 diffs)
• classes/calving.m (modified) (1 diff)
• classes/calving.py (modified) (2 diffs)
• classes/calvinglevermann.m (modified) (1 diff)
• classes/calvinglevermann.py (modified) (1 diff)
• classes/calvingpi.m (modified) (1 diff)
• classes/damage.m (modified) (2 diffs)
• classes/damage.py (modified) (2 diffs)
• classes/flowequation.m (modified) (1 diff)
• classes/flowequation.py (modified) (2 diffs)
• classes/friction.m (modified) (1 diff)
• classes/friction.py (modified) (2 diffs)
• classes/frictionhydro.m (modified) (1 diff)
• classes/frictionwaterlayer.m (modified) (1 diff)
• classes/frictionweertman.m (modified) (2 diffs)
• classes/frictionweertman.py (modified) (1 diff)
• classes/geometry.m (modified) (1 diff)
• classes/geometry.py (modified) (2 diffs)
• classes/gia.m (modified) (1 diff)
• classes/gia.py (modified) (2 diffs)
• classes/hydrologydc.m (modified) (11 diffs)
• classes/hydrologydc.py (modified) (4 diffs)
• classes/hydrologyshreve.m (modified) (1 diff)
• classes/hydrologyshreve.py (modified) (1 diff)
• classes/initialization.m (modified) (1 diff)
• classes/initialization.py (modified) (2 diffs)
• classes/inversion.m (modified) (1 diff)
• classes/inversion.py (modified) (2 diffs)
• classes/inversionvalidation.m (modified) (1 diff)
• classes/m1qn3inversion.m (modified) (1 diff)
• classes/m1qn3inversion.py (modified) (2 diffs)
• classes/mask.m (modified) (1 diff)
• classes/mask.py (modified) (2 diffs)
• classes/masstransport.m (modified) (1 diff)
• classes/masstransport.py (modified) (2 diffs)
• classes/matdamageice.m (modified) (1 diff)
• classes/matdamageice.py (modified) (2 diffs)
• classes/matice.m (modified) (1 diff)
• classes/matice.py (modified) (2 diffs)
• classes/model.m (modified) (2 diffs)
• classes/model.py (modified) (2 diffs)
• classes/qmu.m (modified) (1 diff)
• classes/qmu.py (modified) (3 diffs)
• classes/stressbalance.m (modified) (1 diff)
• classes/stressbalance.py (modified) (2 diffs)
• classes/taoinversion.m (modified) (1 diff)
• classes/thermal.m (modified) (1 diff)
• classes/thermal.py (modified) (2 diffs)
• dev/issmversion.m (modified) (1 diff)

issm/trunk-jpl/src/m/classes/basalforcings.m

@@ -11 +11 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- basalforcings -->');
-                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="melting_rate" type="',            class(self.melting_rate),'" default="',              num2str(self.melting_rate),'">',              '     <section name="basalforcings" />','     <help> basal melting rate (positive if melting) [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="geothermalflux" type="',          class(self.geothermalflux),'" default="',            num2str(self.geothermalflux),'">',            '     <section name="basalforcings" />','     <help> geothermal heat flux [W/m^2] </help>','</parameter>');
-             
-        end % }}}
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- basalforcings -->');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="melting_rate" type="',            class(self.melting_rate),'" default="',              num2str(self.melting_rate),'">',              '     <section name="basalforcings" />','     <help> basal melting rate (positive if melting) [m/yr] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="geothermalflux" type="',          class(self.geothermalflux),'" default="',            num2str(self.geothermalflux),'">',            '     <section name="basalforcings" />','     <help> geothermal heat flux [W/m^2] </help>','</parameter>');
+
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); 
+                        self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1); 
+                        self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
+                end % }}}
                 function self = basalforcings(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/basalforcings.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -30 +31 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1)
+                self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1)
+                self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1)    #bedrock only gets geothermal flux
+                return self
+        #}}}
         def initialize(self,md): # {{{
 

issm/trunk-jpl/src/m/classes/calving.m

@@ -29 +29 @@
                         end
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.calvingrate=project3d(md,'vector',self.calvingrate,'type','node');
+                        self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
+                end % }}}
                 function self = setdefaultparameters(self) % {{{
 

issm/trunk-jpl/src/m/classes/calving.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -31 +32 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.calvingrate=project3d(md,'vector',self.calvingrate,'type','node')
+                self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
 

issm/trunk-jpl/src/m/classes/calvinglevermann.m

@@ -28 +28 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.coeff=project3d(md,'vector',self.coeff,'type','node');
+                        self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
                 end % }}}
                 function self = setdefaultparameters(self) % {{{

issm/trunk-jpl/src/m/classes/calvinglevermann.py

@@ -30 +30 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.coeff=project3d(md,'vector',self.coeff,'type','node')
+                self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
 

issm/trunk-jpl/src/m/classes/calvingpi.m

@@ -28 +28 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.coeff=project3d(md,'vector',self.coeff,'type','node');
+                        self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
                 end % }}}
                 function self = setdefaultparameters(self) % {{{

issm/trunk-jpl/src/m/classes/damage.m

@@ -85 +85 @@
 
                 end % }}}
-        
                 function self = damage(varargin) % {{{
                         switch nargin
@@ -103 +102 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.D=project3d(md,'vector',self.D,'type','node');
+                        self.spcdamage=project3d(md,'vector',self.spcdamage,'type','node');
                 end % }}}
                 function self = setdefaultparameters(self) % {{{

issm/trunk-jpl/src/m/classes/damage.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -74 +75 @@
                 return s
         # }}}
+        def extrude(self,md): # {{{
+                self.D=project3d(md,'vector',self.D,'type','node')
+                self.spcdamage=project3d(md,'vector',self.spcdamage,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self):    # {{{
 

issm/trunk-jpl/src/m/classes/flowequation.m

@@ -107 +107 @@
                         fprintf(fid,'%s\n%s\n','</frame>');
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.element_equation=project3d(md,'vector',self.element_equation,'type','element');
+                        self.vertex_equation=project3d(md,'vector',self.vertex_equation,'type','node');
+                        self.borderSSA=project3d(md,'vector',self.borderSSA,'type','node');
+                        self.borderHO=project3d(md,'vector',self.borderHO,'type','node');
+                        self.borderFS=project3d(md,'vector',self.borderFS,'type','node');
+                end % }}}
                 function self = flowequation(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/flowequation.py

@@ -1 +1 @@
 import numpy
 import copy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -58 +59 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.element_equation=project3d(md,'vector',self.element_equation,'type','element')
+                self.vertex_equation=project3d(md,'vector',self.vertex_equation,'type','node')
+                self.borderSSA=project3d(md,'vector',self.borderSSA,'type','node')
+                self.borderHO=project3d(md,'vector',self.borderHO,'type','node')
+                self.borderFS=project3d(md,'vector',self.borderFS,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
 

issm/trunk-jpl/src/m/classes/friction.m

@@ -11 +11 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="friction" />');   
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',          class(self.coefficient),'" default="',          convert2str(self.coefficient),'">',              '     <section name="friction" />','     <help> friction coefficient [SI] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="friction" />','     <help> p exponent </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="friction" />','     <help> q exponent </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="friction" />');   
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',      class(self.coefficient),'" default="',          convert2str(self.coefficient),'">',              '     <section name="friction" />','     <help> friction coefficient [SI] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="friction" />','     <help> p exponent </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="friction" />','     <help> q exponent </help>','</parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
+                        self.p=project3d(md,'vector',self.p,'type','element');
+                        self.q=project3d(md,'vector',self.q,'type','element');
+                end % }}}
                 function self = friction(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/friction.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -29 +30 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1)
+                self.p=project3d(md,'vector',self.p,'type','element')
+                self.q=project3d(md,'vector',self.q,'type','element')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 return self

issm/trunk-jpl/src/m/classes/frictionhydro.m

@@ -31 +31 @@
                         md = checkfield(md,'fieldname','friction.As','NaN',1,'size',[md.mesh.numberofelements 1]);
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.q=project3d(md,'vector',self.q,'type','element');
+                        self.C=project3d(md,'vector',self.C,'type','element');
+                        self.As=project3d(md,'vector',self.As,'type','element');
+                        self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
+                end % }}}
                 function disp(self) % {{{
                         disp(sprintf('Effective Pressure based friction law described in Gagliardini 2007'));

issm/trunk-jpl/src/m/classes/frictionwaterlayer.m

@@ -36 +36 @@
 
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
+                        self.p=project3d(md,'vector',self.p,'type','element');
+                        self.q=project3d(md,'vector',self.q,'type','element');
+                        self.water_layer=project3d(md,'vector',self.water_layer,'type','node','layer',1);
+                end % }}}
                 function disp(self) % {{{
                         disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(bed+water_layer), r=q/p and s=1/p)'));

issm/trunk-jpl/src/m/classes/frictionweertman.m

@@ -10 +10 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="frictionweertman" />');   
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',          class(self.coefficient),'" default="',          convert2str(self.coefficient),'">',              '     <section name="frictionweertman" />','     <help> frictionweertman coefficient [SI] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="frictionweertman" />','     <help> p exponent </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="frictionweertman" />','     <help> q exponent </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="frictionweertman" />');   
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',      class(self.coefficient),'" default="',          convert2str(self.coefficient),'">',              '     <section name="frictionweertman" />','     <help> frictionweertman coefficient [SI] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="frictionweertman" />','     <help> p exponent </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="frictionweertman" />','     <help> q exponent </help>','</parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
                 function self = frictionweertman(varargin) % {{{
                         switch nargin
@@ -26 +26 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        md.friction.C=project3d(md,'vector',md.friction.C,'type','node','layer',1);
+                        md.friction.m=project3d(md,'vector',md.friction.m,'type','element');
                 end % }}}
                 function self = setdefaultparameters(self) % {{{

issm/trunk-jpl/src/m/classes/frictionweertman.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield

issm/trunk-jpl/src/m/classes/geometry.m

@@ -33 +33 @@
         end
         methods
-        %function createxml(self,fid)
-        function createxml(self, fid)% {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- geometry -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Geometry parameters">','<section name="geometry" />'); 
+                function createxml(self, fid)% {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- geometry -->');
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Geometry parameters">','<section name="geometry" />'); 
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="surface" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> surface elevation [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thickness" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> ice thickness [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bed" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> bed elevation [m] </help>','</parameter>');
-            %fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bathymetry" type="',class(self.bathymetry),'" default="',convert2str(self.bathymetry),'">','     <section name="geometry" />','     <help> bathymetry elevation [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="hydrostatic_ratio" type="',class(self.hydrostatic_ratio),'" default="',convert2str(self.hydrostatic_ratio),'">','     <section name="geometry" />','     <help> coefficient for ice shelves; thickness correction: hydrostatic_ratio H_obs+ (1-hydrostatic_ratio) H_hydro </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thickness" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> ice thickness [m] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bed" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> bed elevation [m] </help>','</parameter>');
+                        %fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bathymetry" type="',class(self.bathymetry),'" default="',convert2str(self.bathymetry),'">','     <section name="geometry" />','     <help> bathymetry elevation [m] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="hydrostatic_ratio" type="',class(self.hydrostatic_ratio),'" default="',convert2str(self.hydrostatic_ratio),'">','     <section name="geometry" />','     <help> coefficient for ice shelves; thickness correction: hydrostatic_ratio H_obs+ (1-hydrostatic_ratio) H_hydro </help>','</parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.surface=project3d(md,'vector',self.surface,'type','node');
+                        self.thickness=project3d(md,'vector',self.thickness,'type','node');
+                        self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node');
+                        self.base=project3d(md,'vector',self.base,'type','node');
+                        self.bed=project3d(md,'vector',self.bed,'type','node');
+                end % }}}
                 function self = geometry(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/geometry.py

@@ +1 @@
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -32 +33 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.surface=project3d(md,'vector',self.surface,'type','node')
+                self.thickness=project3d(md,'vector',self.thickness,'type','node')
+                self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node')
+                self.base=project3d(md,'vector',self.base,'type','node')
+                self.bed=project3d(md,'vector',self.bed,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 return self

issm/trunk-jpl/src/m/classes/gia.m

@@ -11 +11 @@
         end
         methods
-         function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- gia -->\n');            
-                    
-            % gia solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="gia parameters">','<section name="gia" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mantle_viscosity" type="',class(self.mantle_viscosity),'" default="',convert2str(self.mantle_viscosity),'">','     <section name="gia" />','     <help> mantle viscosity[Pa s] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lithosphere_thickness" type="',class(self.lithosphere_thickness),'" default="',convert2str(self.lithosphere_thickness),'">','     <section name="gia" />','     <help> lithosphere thickness (km) </help>','  </parameter>');
-            %cross_section_shape drop-down (1 o r 2)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="cross_section_shape" type="alternative" optional="false">','     <section name="gia" />','     <help> 1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore </help>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
-            
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '<!-- gia -->\n');            
+
+                        % gia solution parameters
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="gia parameters">','<section name="gia" />');                    
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mantle_viscosity" type="',class(self.mantle_viscosity),'" default="',convert2str(self.mantle_viscosity),'">','     <section name="gia" />','     <help> mantle viscosity[Pa s] </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lithosphere_thickness" type="',class(self.lithosphere_thickness),'" default="',convert2str(self.lithosphere_thickness),'">','     <section name="gia" />','     <help> lithosphere thickness (km) </help>','  </parameter>');
+                        %cross_section_shape drop-down (1 o r 2)
+                        fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="cross_section_shape" type="alternative" optional="false">','     <section name="gia" />','     <help> 1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore </help>');
+                        fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
+                        fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+
+                        fprintf(fid,'%s\n%s\n','</frame>');    
+
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node');
+                        self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node');
+                end % }}}
                 function self = gia(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/gia.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -30 +31 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node')
+                self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
 

issm/trunk-jpl/src/m/classes/hydrologydc.m

@@ -37 +37 @@
   end
         methods
+                function self = extrude(self,md) % {{{
+                        self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1);
+                        self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1);
+                        self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1);
+                        self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1);
+                        if(self.isefficientlayer==1);
+                                self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1);
+                        end
+                end % }}}
                 % {{{ function self = hydrologydc(varargin) 
                 function self = hydrologydc(varargin) 
@@ -44 +53 @@
                                 otherwise
                                         error('constructor not supported');
-            end
+                        end
                 end 
                 % }}}
@@ -57 +66 @@
                 function self = setdefaultparameters(self) 
 
-                %Parameters from de Fleurian 2014
+                        %Parameters from de Fleurian 2014
                         self.water_compressibility    = 5.04e-10;
                         self.isefficientlayer         = 1;
@@ -86 +95 @@
                 % {{{ function md = checkconsistency(self,md,solution,analyses) 
                 function md = checkconsistency(self,md,solution,analyses) 
-                %Early return
+                        %Early return
                         if ~ismember(HydrologyDCInefficientAnalysisEnum(),analyses) & ~ismember(HydrologyDCEfficientAnalysisEnum(),analyses),
                                 return;
@@ -101 +110 @@
                         if self.sedimentlimit_flag==1,
                                 md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1);
-            end
+                        end
                         if self.transfer_flag==1,
                                 md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1);
-            end
+                        end
                         md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'timeseries',1);
 
@@ -127 +136 @@
                                         md = checkmessage(md,'Colapsing thickness for EPL larger than initial thickness');
                                 end 
-            end
+                        end
                 end 
                 % }}}
@@ -147 +156 @@
                         if self.sedimentlimit_flag==1,
                                 fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]');
-            end
+                        end
                         fielddisplay(self,'transfer_flag',['what kind of transfer method is applied between the layers']);
                         disp(sprintf('%55s  0: no transfer',' '));
@@ -153 +162 @@
                         if self.transfer_flag==1,
                                 fielddisplay(self,'leakage_factor','user defined leakage factor [m]');
-            end
+                        end
                         fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]');
                         disp(sprintf('   - for the sediment layer'));
@@ -174 +183 @@
                                 fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]');
                                 fielddisplay(self,'eplflip_lock','lock the epl activation to avoid fli-floping (default is 0, no stabilization)');
-            end
+                        end
 
                 end 
@@ -191 +200 @@
                         if self.sedimentlimit_flag==1,
                                 WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double');
-            end
+                        end
                         if self.transfer_flag==1,
                                 WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double');
-            end
+                        end
                         WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1)
 
@@ -214 +223 @@
                                 WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double');
                                 WriteData(fid,'object',self,'fieldname','eplflip_lock','format','Integer');
-            end
-                end 
-% }}}
-  end
+                        end
+                end 
+                % }}}
+        end
 end

issm/trunk-jpl/src/m/classes/hydrologydc.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -47 +48 @@
                 self.setdefaultparameters()
         #}}}
-
         def __repr__(self): # {{{
                 string='   hydrology Dual Porous Continuum Equivalent parameters:'
@@ -95 +95 @@
                 return string
 #}}}
+        def extrude(self,md): # {{{
+                self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1)
+                self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1)
+                self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1)
+                self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1)
+                self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1)
+                if self.isefficientlayer==1 :
+                        self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1)
+                return self
+        #}}}
         def setdefaultparameters(self): #{{{ 
 
@@ -125 +135 @@
                 return self
         # }}}
-
         def initialize(self,md): # {{{
                 if numpy.all(numpy.isnan(self.basal_moulin_input)):

issm/trunk-jpl/src/m/classes/hydrologyshreve.m

@@ -10 +10 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Hydrology -->');
-            
-            % Convergence criteria            
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Hydrologyshreve solution parameters">','<section name="hydrologyshreve" />');                    
-            
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- Hydrology -->');
+
+                        % Convergence criteria            
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Hydrologyshreve solution parameters">','<section name="hydrologyshreve" />');                    
+
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcwatercolumn" type="',      class(self.spcwatercolumn),'" default="',        convert2str(self.spcwatercolumn),'">',      '     <section name="hydrologyshreve" />','     <help> water thickness constraints (NaN means no constraint) [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="stabilization" type="',         class(self.stabilization),'" default="',           convert2str(self.stabilization),'">',       '     <section name="hydrologyshreve" />','     <help> artificial diffusivity (default is 1). can be more than 1 to increase diffusivity. </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="stabilization" type="',         class(self.stabilization),'" default="',           convert2str(self.stabilization),'">',   '     <section name="hydrologyshreve" />','     <help> artificial diffusivity (default is 1). can be more than 1 to increase diffusivity. </help>','</parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
+                function self = extrude(self,md) % {{{
+                end % }}}
                 function self = hydrologyshreve(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/hydrologyshreve.py

@@ -27 +27 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 

issm/trunk-jpl/src/m/classes/initialization.m

@@ -35 +35 @@
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="watercolumn" type="',class(self.watercolumn),'" default="',self.watercolumn,'">','     <section name="initialization" />','     <help> thickness of subglacial water [m] </help>','</parameter>');
                         fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.vx=project3d(md,'vector',self.vx,'type','node');
+                        self.vy=project3d(md,'vector',self.vy,'type','node');
+                        self.vz=project3d(md,'vector',self.vz,'type','node');
+                        self.vel=project3d(md,'vector',self.vel,'type','node');
+                        self.temperature=project3d(md,'vector',self.temperature,'type','node');
+                        self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
+                        self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
+                        self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
+                        self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
+                        self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
+
+                        %Lithostatic pressure by default
+                        self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
                 end % }}}
                 function self = initialization(varargin) % {{{

issm/trunk-jpl/src/m/classes/initialization.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -48 +49 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.vx=project3d(md,'vector',self.vx,'type','node')
+                self.vy=project3d(md,'vector',self.vy,'type','node')
+                self.vz=project3d(md,'vector',self.vz,'type','node')
+                self.vel=project3d(md,'vector',self.vel,'type','node')
+                self.temperature=project3d(md,'vector',self.temperature,'type','node')
+                self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node')
+                self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node')
+                self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1)
+                self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1)
+                self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1)
+
+                #Lithostatic pressure by default
+                self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 return self

issm/trunk-jpl/src/m/classes/inversion.m

@@ -84 +84 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+                        self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+                        self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+                        self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+                        if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+                        if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+                        if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
                 end % }}}
                 function self = setdefaultparameters(self) % {{{

issm/trunk-jpl/src/m/classes/inversion.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -72 +73 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
+                self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
+                self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
+                self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
+                if not numpy.any(numpy.isnan(self.cost_functions_coefficients)):
+                        self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
+                if not numpy.any(numpy.isnan(self.min_parameters)):
+                        self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
+                if not numpy.any(numpy.isnan(self.max_parameters)):
+                        self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 

issm/trunk-jpl/src/m/classes/inversionvalidation.m

@@ -22 +22 @@
         end
         methods
+                function self = extrude(self,md) % {{{
+                        self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+                        self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+                        self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+                        self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+                        if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+                        if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+                        if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+                end % }}}
                 function self = inversionvalidation(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/m1qn3inversion.m

@@ -27 +27 @@
         end
         methods
+                function self = extrude(self,md) % {{{
+                        self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+                        self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+                        self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+                        self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+                        if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+                        if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+                        if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+                end % }}}
                 function self = m1qn3inversion(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/m1qn3inversion.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -94 +95 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
+                self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
+                self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
+                self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
+                if not numpy.any(numpy.isnan(self.cost_functions_coefficients)):
+                        self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
+                if not numpy.any(numpy.isnan(self.min_parameters)):
+                        self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
+                if not numpy.any(numpy.isnan(self.max_parameters)):
+                        self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 

issm/trunk-jpl/src/m/classes/mask.m

@@ -22 +22 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- mask -->');
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- mask -->');
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="groundedice_levelset" type="',class(self.groundedice_levelset),'" default="',self.groundedice_levelset,'">','     <section name="mask" />','     <help> is ice grounded ? grounded ice if > 0, grounding line position if = 0, floating ice if &amp;lt; 0 </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="ice_levelset" type="',class(self.ice_levelset),'" default="',self.ice_levelset,'">','     <section name="mask" />','     <help> presence of ice if > 0, icefront position if = 0, no ice if &amp;lt; 0 </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="ice_levelset" type="',class(self.ice_levelset),'" default="',self.ice_levelset,'">','     <section name="mask" />','     <help> presence of ice if > 0, icefront position if = 0, no ice if &amp;lt; 0 </help>','</parameter>');
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node');
+                        self.ice_levelset=project3d(md,'vector',self.ice_levelset,'type','node');
                 end % }}}
                 function self = mask(varargin) % {{{

issm/trunk-jpl/src/m/classes/mask.py

@@ -1 +1 @@
 import numpy
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -28 +29 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.ice_levelset=project3d(md,'vector',self.ice_levelset,'type','node')
+                self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 return self

issm/trunk-jpl/src/m/classes/masstransport.m

@@ -58 +58 @@
                                         error('constructor not supported');
                         end
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.spcthickness=project3d(md,'vector',self.spcthickness,'type','node');
                 end % }}}
                 function list = defaultoutputs(self,md) % {{{

issm/trunk-jpl/src/m/classes/masstransport.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -38 +39 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.spcthickness=project3d(md,'vector',self.spcthickness,'type','node')
+                return self
+        #}}}
         def defaultoutputs(self,md): # {{{
 

issm/trunk-jpl/src/m/classes/matdamageice.m

@@ -63 +63 @@
 
 
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
+                        self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
                 end % }}}
                 function self = matdamageice(varargin) % {{{

issm/trunk-jpl/src/m/classes/matdamageice.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import MaterialsEnum, MatdamageiceEnum, MaterialsRheologyLawEnum, MaterialsRhoSeawaterEnum
 from StringToEnum import StringToEnum
@@ -63 +64 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+                self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 #ice density (kg/m^3)

issm/trunk-jpl/src/m/classes/matice.m

@@ -70 +70 @@
 
 
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
+                        self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
                 end % }}}
                 function self = matice(varargin) % {{{

issm/trunk-jpl/src/m/classes/matice.py

@@ -1 +1 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -63 +64 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+                self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 #ice density (kg/m^3)

issm/trunk-jpl/src/m/classes/model.m

@@ -750 +750 @@
                         %Ok, now deal with the other fields from the 2d mesh:
 
+                        %bedinfo and surface info
+                        md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
+                        md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
+                        md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
+
                         %lat long
                         md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
                         md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
 
-                        %drag coefficient is limited to nodes that are on the bedrock.
-                        if isa(md.friction,'friction'),
-                                md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
-                                md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
-                                md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-                        elseif isa(md.friction,'frictionhydro'),
-                                md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-                                md.friction.C=project3d(md,'vector',md.friction.C,'type','element');
-                                md.friction.As=project3d(md,'vector',md.friction.As,'type','element');
-                                md.friction.effective_pressure=project3d(md,'vector',md.friction.effective_pressure,'type','node','layer',1);
-                        elseif isa(md.friction,'frictionwaterlayer'),
-                                md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
-                                md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
-                                md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-                                md.friction.water_layer=project3d(md,'vector',md.friction.water_layer,'type','node','layer',1);
-                        elseif isa(md.friction,'frictionweertman'),
-                                md.friction.C=project3d(md,'vector',md.friction.C,'type','node','layer',1);
-                                md.friction.m=project3d(md,'vector',md.friction.m,'type','element');
-            end
-         
-                        %observations
-                        md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node');
-                        md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node');
-                        md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node');
-                        md.inversion.thickness_obs=project3d(md,'vector',md.inversion.thickness_obs,'type','node');
+                        md.geometry=extrude(md.geometry,md);
+                        md.friction  = extrude(md.friction,md);
+                        md.inversion = extrude(md.inversion,md);
                         md.surfaceforcings = extrude(md.surfaceforcings,md);
-                        md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
-
-                        %results
-                        if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end;
-                        if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end;
-                        if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end;
-                        if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end;
-                        if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end;
-                        if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end;
-      if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);end;
-      if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1);end;
-      if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1);end;
-      if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1);end;
-
-                        %bedinfo and surface info
-                        md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
-                        md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
-
-                        %elementstype
-                        if ~isnan(md.flowequation.element_equation)
-                                oldelements_type=md.flowequation.element_equation;
-                                md.flowequation.element_equation=zeros(number_el3d,1);
-                                md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element');
-                        end
-
-                        %verticestype
-                        if ~isnan(md.flowequation.vertex_equation)
-                                oldvertices_type=md.flowequation.vertex_equation;
-                                md.flowequation.vertex_equation=zeros(number_nodes3d,1);
-                                md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node');
-                        end
-                        md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node');
-                        md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node');
-                        md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node');
-
-                        %boundary conditions
-                        md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node');
-                        md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node');
-                        md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node');
-                        md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
-                        if (length(md.initialization.temperature)==md.mesh.numberofvertices),
-                                md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
-                                if isprop(md.mesh,'vertexonsurface'),
-                                        pos=find(md.mesh.vertexonsurface);
-                                        md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-                                end
-                        end
-                        md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node');
-                        md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
-                        md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node');
-                        md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
-                        md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
-
-                        % Calving variables
-                        if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end;
-                        if ~isnan(md.calving.meltingrate),md.calving.meltingrate=project3d(md,'vector',md.calving.meltingrate,'type','node');end;
-                        if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
-                        if isa(md.calving,'calvingpi'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
-
-                        % Hydrologydc variables
-                        if isa(md.hydrology,'hydrologydc');
-                                md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1);
-                                md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1);
-                                md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1);
-                                md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1);
-                                if(md.hydrology.isefficientlayer==1);
-                                        md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1);
-                    end
-            end
+                        md.initialization = extrude(md.initialization,md);
+
+                        md.flowequation.extrude(md);
+                        md.stressbalance=extrude(md.stressbalance,md);
+                        md.thermal.extrude(md);
+                        md.masstransport.extrude(md);
+                        md.calving=extrude(md.calving,md);
+                        md.hydrology = extrude(md.hydrology,md);
 
                         %connectivity
@@ -861 +783 @@
                         end
 
-                        %materials
-                        md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
-                        md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
-                                
-                        %damage
-                        md.damage.D=project3d(md,'vector',md.damage.D,'type','node');
-
-                        %parameters
-                        md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node');
-                        md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node');
-                        md.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node');
-                        md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_thickness,'type','node');
-                        md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node');
-                        md.geometry.base=project3d(md,'vector',md.geometry.base,'type','node');
-                        md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node');
-                        md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
-                        md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node');
-                        md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node');
-                        if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end;
-                        if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end;
-                        if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end;
-                        if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end
-
-                        %Put lithostatic pressure if there is an existing pressure
-                        if ~isnan(md.initialization.pressure),
-                                md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
-            end
-
-                        %special for thermal modeling:
-                        md.basalforcings.groundedice_melting_rate=project3d(md,'vector',md.basalforcings.groundedice_melting_rate,'type','node','layer',1); 
-                        md.basalforcings.floatingice_melting_rate=project3d(md,'vector',md.basalforcings.floatingice_melting_rate,'type','node','layer',1); 
-                        if ~isnan(md.basalforcings.geothermalflux)
-                                md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
-                        end
+                        md.materials=extrude(md.materials,md);
+                        md.damage=extrude(md.damage,md);
+                        md.mask=extrude(md.mask,md);
+                        md.qmu=extrude(md.qmu,md);
+                        md.basalforcings=extrude(md.basalforcings,md);
 
                         %increase connectivity if less than 25:

issm/trunk-jpl/src/m/classes/model.py

@@ -595 +595 @@
                 #Ok, now deal with the other fields from the 2d mesh:
 
+                #bedinfo and surface info
+                md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
+                md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
+                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
+
                 #lat long
                 md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node')
                 md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node')
 
-                #drag coefficient is limited to nodes that are on the bedrock.
-                md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1)
-
-                #p and q (same deal, except for element that are on the bedrock: )
-                md.friction.p=project3d(md,'vector',md.friction.p,'type','element')
-                md.friction.q=project3d(md,'vector',md.friction.q,'type','element')
-
-                #observations
-                md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node')
-                md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node')
-                md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node')
-                md.inversion.thickness_obs=project3d(md,'vector',md.inversion.thickness_obs,'type','node')
+                md.geometry.extrude(md)
+                md.friction.extrude(md)
+                md.inversion.extrude(md)
                 md.surfaceforcings.extrude(md)
-                md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node')
-
-                #results
-                if not numpy.any(numpy.isnan(md.initialization.vx)):
-                        md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.vy)):
-                        md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.vz)):
-                        md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.vel)):
-                        md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.temperature)):
-                        md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.waterfraction)):
-                        md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.watercolumn)):
-                        md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.sediment_head)):
-                        md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1)
-                if not numpy.any(numpy.isnan(md.initialization.epl_head)):
-                        md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1)
-                if not numpy.any(numpy.isnan(md.initialization.epl_thickness)):
-                        md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1)
-
-                #bedinfo and surface info
-                md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
-                md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
-
-                #elementstype
-                if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
-                        oldelements_type=md.flowequation.element_equation
-                        md.flowequation.element_equation=numpy.zeros(number_el3d,int)
-                        md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
-
-                #verticestype
-                if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
-                        oldvertices_type=md.flowequation.vertex_equation
-                        md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int)
-                        md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
-
-                md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node')
-                md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node')
-                md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node')
-
-                #boundary conditions
-                md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node')
-                md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node')
-                md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node')
-                md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
-                if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-                        md.thermal.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
-                        if hasattr(md.mesh,'vertexonsurface'):
-                                pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
-                                md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-                md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node')
-                md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node')
-                md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node')
-                md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node')
-                md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
+                md.initialization.extrude(md)
+                md.flowequation.extrude(md)
+
+                md.stressbalance.extrude(md)
+                md.thermal.extrude(md)
+                md.masstransport.extrude(md)
 
                 # Calving variables
-                if isinstance(md.calving,calving):
-                        md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node')
-                if isinstance(md.calving,calvinglevermann):
-                        md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node')
-                #if isinstance(md.calving,calvingpi):
-                #       md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node')
-                if not numpy.any(numpy.isnan(md.calving.meltingrate)):
-                        md.calving.meltingrate=project3d(md,'vector',md.calving.meltingrate,'type','node')
-
-                # Hydrologydc variables
-                if hasattr(md.hydrology,'hydrologydc'):
-                        md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1)
-                        md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
-                        md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1)
-                        md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1)
-                        md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1)
-                        if md.hydrology.isefficientlayer==1 :
-                                md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
+                md.hydrology.extrude(md)
+                md.calving.extrude(md)
 
                 #connectivity
@@ -701 +628 @@
                                 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
 
-                #materials
-                md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node')
-                md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element')
-
-                #damage
-                md.damage.D=project3d(md,'vector',md.damage.D,'type','node')
-
-                #parameters
-                md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node')
-                md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node')
-                md.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node')
-                md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_thickness,'type','node')
-                md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node')
-                md.geometry.base=project3d(md,'vector',md.geometry.base,'type','node')
-                md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node')
-                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
-                md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node')
-                md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node')
-                if not numpy.any(numpy.isnan(md.inversion.cost_functions_coefficients)):
-                        md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node')
-                if not numpy.any(numpy.isnan(md.inversion.min_parameters)):
-                        md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node')
-                if not numpy.any(numpy.isnan(md.inversion.max_parameters)):
-                        md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node')
-                if not numpy.any(numpy.isnan(md.qmu.partition)):
-                        md.qmu.partition=project3d(md,'vector',numpy.transpose(md.qmu.partition),'type','node')
-
-                #Put lithostatic pressure if there is an existing pressure
-                if not numpy.any(numpy.isnan(md.initialization.pressure)):
-                        md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
-
-                #special for thermal modeling:
-                md.basalforcings.groundedice_melting_rate=project3d(md,'vector',md.basalforcings.groundedice_melting_rate,'type','node','layer',1)
-                md.basalforcings.floatingice_melting_rate=project3d(md,'vector',md.basalforcings.floatingice_melting_rate,'type','node','layer',1)
-                if not numpy.any(numpy.isnan(md.basalforcings.geothermalflux)):
-                        md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1)    #bedrock only gets geothermal flux
+                md.materials.extrude(md)
+                md.damage.extrude(md)
+                md.gia.extrude(md)
+                md.mask.extrude(md)
+                md.qmu.extrude(md)
+                md.basalforcings.extrude(md)
 
                 #increase connectivity if less than 25:

issm/trunk-jpl/src/m/classes/qmu.m

@@ -24 +24 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- qmu -->\n');            
-                    
-            % qmu solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="qmu parameters">','<section name="qmu" />');                    
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdakota" type="',class(self.isdakota),'" default="',convert2str(self.isdakota),'">','     <section name="qmu" />','     <help> is qmu analysis activated? </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variables" type="',class(self.variables),'" default="',convert2str(self.variables),'">','     <section name="qmu" />','     <help> (arrays of each variable class) </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responses" type="',class(self.responses),'" default="',convert2str(self.responses),'">','     <section name="qmu" />','     <help> (arrays of each response class) </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofresponses" type="',class(self.numberofresponses),'" default="',convert2str(self.numberofresponses),'">','     <section name="qmu" />','     <help> number of responses </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="params" type="',class(self.params),'" default="',convert2str(self.params),'">','     <section name="qmu" />','     <help> (array of method-independent parameters)  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="results" type="',class(self.results),'" default="',convert2str(self.results),'">','     <section name="qmu" />','     <help> (information from dakota files) </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="partition" type="',class(self.partition),'" default="',convert2str(self.partition),'">','     <section name="qmu" />','     <help> user provided mesh partitioning, defaults to metis if not specified </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofpartitions" type="',class(self.numberofpartitions),'" default="',convert2str(self.numberofpartitions),'">','     <section name="qmu" />','     <help> number of partitions for semi-discrete qmu  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variabledescriptors" type="',class(self.variabledescriptors),'" default="',convert2str(self.variabledescriptors),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-            
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responsedescriptors" type="',class(self.responsedescriptors),'" default="',convert2str(self.responsedescriptors),'">','     <section name="qmu" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="method" type="',class(self.method),'" default="',convert2str(self.method),'">','     <section name="qmu" />','     <help> array of dakota_method class </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profile_directory" type="',class(self.mass_flux_profile_directory),'" default="',convert2str(self.mass_flux_profile_directory),'">','     <section name="qmu" />','     <help> directory for mass flux profiles </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profiles" type="',class(self.mass_flux_profiles),'" default="',convert2str(self.mass_flux_profiles),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_segments" type="',class(self.mass_flux_segments),'" default="',convert2str(self.mass_flux_segments),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="adjacency" type="',class(self.adjacency),'" default="',convert2str(self.adjacency),'">','     <section name="qmu" />','     <help> additional outputs requested </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertex_weight" type="',class(self.vertex_weight),'" default="',convert2str(self.vertex_weight),'">','     <section name="qmu" />','     <help> weight applied to each mesh vertex </help>','  </parameter>');
-            
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}        
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '<!-- qmu -->\n');            
+
+                        % qmu solution parameters
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="qmu parameters">','<section name="qmu" />');                    
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdakota" type="',class(self.isdakota),'" default="',convert2str(self.isdakota),'">','     <section name="qmu" />','     <help> is qmu analysis activated? </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variables" type="',class(self.variables),'" default="',convert2str(self.variables),'">','     <section name="qmu" />','     <help> (arrays of each variable class) </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responses" type="',class(self.responses),'" default="',convert2str(self.responses),'">','     <section name="qmu" />','     <help> (arrays of each response class) </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofresponses" type="',class(self.numberofresponses),'" default="',convert2str(self.numberofresponses),'">','     <section name="qmu" />','     <help> number of responses </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="params" type="',class(self.params),'" default="',convert2str(self.params),'">','     <section name="qmu" />','     <help> (array of method-independent parameters)  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="results" type="',class(self.results),'" default="',convert2str(self.results),'">','     <section name="qmu" />','     <help> (information from dakota files) </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="partition" type="',class(self.partition),'" default="',convert2str(self.partition),'">','     <section name="qmu" />','     <help> user provided mesh partitioning, defaults to metis if not specified </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofpartitions" type="',class(self.numberofpartitions),'" default="',convert2str(self.numberofpartitions),'">','     <section name="qmu" />','     <help> number of partitions for semi-discrete qmu  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variabledescriptors" type="',class(self.variabledescriptors),'" default="',convert2str(self.variabledescriptors),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responsedescriptors" type="',class(self.responsedescriptors),'" default="',convert2str(self.responsedescriptors),'">','     <section name="qmu" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="method" type="',class(self.method),'" default="',convert2str(self.method),'">','     <section name="qmu" />','     <help> array of dakota_method class </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profile_directory" type="',class(self.mass_flux_profile_directory),'" default="',convert2str(self.mass_flux_profile_directory),'">','     <section name="qmu" />','     <help> directory for mass flux profiles </help>','  </parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profiles" type="',class(self.mass_flux_profiles),'" default="',convert2str(self.mass_flux_profiles),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_segments" type="',class(self.mass_flux_segments),'" default="',convert2str(self.mass_flux_segments),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="adjacency" type="',class(self.adjacency),'" default="',convert2str(self.adjacency),'">','     <section name="qmu" />','     <help> additional outputs requested </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertex_weight" type="',class(self.vertex_weight),'" default="',convert2str(self.vertex_weight),'">','     <section name="qmu" />','     <help> weight applied to each mesh vertex </help>','  </parameter>');
+
+                        fprintf(fid,'%s\n%s\n','</frame>');    
+
+                end % }}}        
+                function self = extrude(self,md) % {{{
+                        self.partition=project3d(md,'vector',self.partition','type','node');
+                end % }}}
                 function self = qmu(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/qmu.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from collections import OrderedDict
 from fielddisplay import fielddisplay
@@ -37 +38 @@
 
                 #}}}
+        def __repr__(self):    # {{{
+                s ='   qmu parameters:\n'
+
+                s+="%s\n" % fielddisplay(self,'isdakota','is qmu analysis activated?')
+                for i,variable in enumerate(self.variables.iteritems()):
+                        s+="         variables%s:  (arrays of each variable class)\n" % \
+                                        string_dim(self.variables,i)
+                        fnames=vars(variable)
+                        maxlen=0
+                        for fname in fnames:
+                                maxlen=max(maxlen,len(fname))
+
+                        for fname in fnames:
+                                s+="'            %-*s:    [%ix%i]    '%s'\n" % \
+                                                (maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
+
+                for i,response in enumerate(self.responses.iteritems()):
+                        s+="         responses%s:  (arrays of each response class)\n" % \
+                                        string_dim(self.responses,i)
+                        fnames=vars(response)
+                        maxlen=0
+                        for fname in fnames:
+                                maxlen=max(maxlen,len(fname))
+
+                        for fname in fnames:
+                                s+="            %-*s:    [%ix%i]    '%s'\n" % \
+                                                (maxlen+1,fname,size(getattr(response,fname)),type(getattr(response,fname)))
+
+                s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses') 
+
+                for i,method in enumerate(self.method.iteritems()):
+                        if isinstance(method,'dakota_method'):
+                                s+="            method%s :    '%s'\n" % \
+                                                (string_dim(method,i),method.method)
+
+                for i,param in enumerate(self.params.iteritems()):
+                        s+="         params%s:  (array of method-independent parameters)\n" % \
+                                        string_dim(self.params,i)
+                        fnames=vars(param)
+                        maxlen=0
+                        for fname in fnames:
+                                maxlen=max(maxlen,len(fname))
+
+                        for fname in fnames:
+                                s+="            %-*s: %s\n" % \
+                                                (maxlen+1,fname,any2str(getattr(param,fname)))
+
+                for i,result in enumerate(self.results.iteritems()):
+                        s+="         results%s:  (information from dakota files)\n" % \
+                                        string_dim(self.results,i)
+                        fnames=vars(result)
+                        maxlen=0
+                        for fname in fnames:
+                                maxlen=max(maxlen,len(fname))
+
+                        for fname in fnames:
+                                s+="            %-*s:    [%ix%i]    '%s'\n" % \
+                                                (maxlen+1,fname,size(getattr(result,fname)),type(getattr(result,fname)))
+
+                s+="%s\n" % fielddisplay(self,'partition','user provided mesh partitioning, defaults to metis if not specified') 
+                s+="%s\n" % fielddisplay(self,'numberofpartitions','number of partitions for semi-discrete qmu') 
+                s+="%s\n" % fielddisplay(self,'variabledescriptors','')
+                s+="%s\n" % fielddisplay(self,'responsedescriptors','')
+                s+="%s\n" % fielddisplay(self,'method','array of dakota_method class')
+                s+="%s\n" % fielddisplay(self,'mass_flux_profile_directory','directory for mass flux profiles')
+                s+="%s\n" % fielddisplay(self,'mass_flux_profiles','list of mass_flux profiles')
+                s+="%s\n" % fielddisplay(self,'mass_flux_segments','')
+                s+="%s\n" % fielddisplay(self,'adjacency','')
+                s+="%s\n" % fielddisplay(self,'vertex_weight','weight applied to each mesh vertex')
+
+                return s
+        # }}}
+        def extrude(self,md): # {{{
+                self.partition=project3d(md,'vector',numpy.transpose(self.partition),'type','node')
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 return self
@@ -62 +139 @@
                 return md
         # }}}
-        def __repr__(self):    # {{{
-                s ='   qmu parameters:\n'
-
-                s+="%s\n" % fielddisplay(self,'isdakota','is qmu analysis activated?')
-                for i,variable in enumerate(self.variables.iteritems()):
-                        s+="         variables%s:  (arrays of each variable class)\n" % \
-                                string_dim(self.variables,i)
-                        fnames=vars(variable)
-                        maxlen=0
-                        for fname in fnames:
-                                maxlen=max(maxlen,len(fname))
-
-                        for fname in fnames:
-                                s+="'            %-*s:    [%ix%i]    '%s'\n" % \
-                                        (maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
-
-                for i,response in enumerate(self.responses.iteritems()):
-                        s+="         responses%s:  (arrays of each response class)\n" % \
-                                string_dim(self.responses,i)
-                        fnames=vars(response)
-                        maxlen=0
-                        for fname in fnames:
-                                maxlen=max(maxlen,len(fname))
-
-                        for fname in fnames:
-                                s+="            %-*s:    [%ix%i]    '%s'\n" % \
-                                        (maxlen+1,fname,size(getattr(response,fname)),type(getattr(response,fname)))
-
-                s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses') 
-
-                for i,method in enumerate(self.method.iteritems()):
-                        if isinstance(method,'dakota_method'):
-                                s+="            method%s :    '%s'\n" % \
-                                        (string_dim(method,i),method.method)
-
-                for i,param in enumerate(self.params.iteritems()):
-                        s+="         params%s:  (array of method-independent parameters)\n" % \
-                                string_dim(self.params,i)
-                        fnames=vars(param)
-                        maxlen=0
-                        for fname in fnames:
-                                maxlen=max(maxlen,len(fname))
-
-                        for fname in fnames:
-                                s+="            %-*s: %s\n" % \
-                                        (maxlen+1,fname,any2str(getattr(param,fname)))
-
-                for i,result in enumerate(self.results.iteritems()):
-                        s+="         results%s:  (information from dakota files)\n" % \
-                                string_dim(self.results,i)
-                        fnames=vars(result)
-                        maxlen=0
-                        for fname in fnames:
-                                maxlen=max(maxlen,len(fname))
-
-                        for fname in fnames:
-                                s+="            %-*s:    [%ix%i]    '%s'\n" % \
-                                        (maxlen+1,fname,size(getattr(result,fname)),type(getattr(result,fname)))
-
-                s+="%s\n" % fielddisplay(self,'partition','user provided mesh partitioning, defaults to metis if not specified') 
-                s+="%s\n" % fielddisplay(self,'numberofpartitions','number of partitions for semi-discrete qmu') 
-                s+="%s\n" % fielddisplay(self,'variabledescriptors','')
-                s+="%s\n" % fielddisplay(self,'responsedescriptors','')
-                s+="%s\n" % fielddisplay(self,'method','array of dakota_method class')
-                s+="%s\n" % fielddisplay(self,'mass_flux_profile_directory','directory for mass flux profiles')
-                s+="%s\n" % fielddisplay(self,'mass_flux_profiles','list of mass_flux profiles')
-                s+="%s\n" % fielddisplay(self,'mass_flux_segments','')
-                s+="%s\n" % fielddisplay(self,'adjacency','')
-                s+="%s\n" % fielddisplay(self,'vertex_weight','weight applied to each mesh vertex')
-
-                return s
-        # }}}
         def marshall(self,md,fid):    # {{{
                 WriteData(fid,'object',self,'fieldname','isdakota','format','Boolean')

issm/trunk-jpl/src/m/classes/stressbalance.m

@@ -70 +70 @@
                         fprintf(fid,'%s\n%s\n','</frame>'); 
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.spcvx=project3d(md,'vector',self.spcvx,'type','node');
+                        self.spcvy=project3d(md,'vector',self.spcvy,'type','node');
+                        self.spcvz=project3d(md,'vector',self.spcvz,'type','node');
+                        self.referential=project3d(md,'vector',self.referential,'type','node');
+                        self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node');
+
+                end % }}}
                 function self = stressbalance(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/stressbalance.py

@@ -2 +2 @@
 import sys
 import copy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -75 +76 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.spcvx=project3d(md,'vector',self.spcvx,'type','node')
+                self.spcvy=project3d(md,'vector',self.spcvy,'type','node')
+                self.spcvz=project3d(md,'vector',self.spcvz,'type','node')
+                self.referential=project3d(md,'vector',self.referential,'type','node')
+                self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node')
+
+                return self
+        #}}}
         def setdefaultparameters(self): # {{{
                 #maximum of non-linear iterations.

issm/trunk-jpl/src/m/classes/taoinversion.m

@@ -29 +29 @@
         end
         methods
+                function self = extrude(self,md) % {{{
+                        self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+                        self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+                        self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+                        self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+                        if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+                        if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+                        if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+                end % }}}
                 function self = taoinversion(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/thermal.m

@@ -18 +18 @@
         end
         methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- thermal -->\n');            
-                    
-            % thermal solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Thermal solution parameters">','<section name="thermal" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spctemperature" type="',class(self.spctemperature),'" default="',convert2str(self.spctemperature),'">','     <section name="thermal" />','     <help> temperature constraints (NaN means no constraint) [K] </help>','  </parameter>');
-                
-            % penalty_threshold drop-down (0, 1, or 2)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="penalty_threshold" type="alternative"  optional="false">','     <section name="thermal" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+                function createxml(self,fid) % {{{
+                        fprintf(fid, '<!-- thermal -->\n');            
 
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(self.stabilization),'" default="',convert2str(self.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
-                                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="reltol" type="',class(self.reltol),'" default="',convert2str(self.reltol),'">','     <section name="steadystate" />','     <help> relative tolerance criterion [K] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(self.maxiter),'" default="',convert2str(self.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(self.penalty_lock),'" default="',convert2str(self.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(self.penalty_factor),'" default="',convert2str(self.penalty_factor),'">','     <section name="thermal" />','     <help> scaling exponent (default is 3) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isenthalpy" type="',class(self.isenthalpy),'" default="',convert2str(self.isenthalpy),'">','     <section name="thermal" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdynamicbasalspc" type="',class(self.isdynamicbasalspc),'" default="',convert2str(self.isdynamicbasalspc),'">','     <section name="thermal" />','     <help> enable dynamic setting of basal forcing. recommended for enthalpy formulation (default is 0)  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(self.requested_outputs),'" default="',convert2str(self.requested_outputs),'">','     <section name="thermal" />','     <help> additional outputs requested </help>','  </parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}
+                        % thermal solution parameters
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Thermal solution parameters">','<section name="thermal" />');                    
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spctemperature" type="',class(self.spctemperature),'" default="',convert2str(self.spctemperature),'">','     <section name="thermal" />','     <help> temperature constraints (NaN means no constraint) [K] </help>','  </parameter>');
+
+                        % penalty_threshold drop-down (0, 1, or 2)
+                        fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="penalty_threshold" type="alternative"  optional="false">','     <section name="thermal" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>');
+                        fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
+                        fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
+                        fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(self.stabilization),'" default="',convert2str(self.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="reltol" type="',class(self.reltol),'" default="',convert2str(self.reltol),'">','     <section name="steadystate" />','     <help> relative tolerance criterion [K] </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(self.maxiter),'" default="',convert2str(self.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(self.penalty_lock),'" default="',convert2str(self.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(self.penalty_factor),'" default="',convert2str(self.penalty_factor),'">','     <section name="thermal" />','     <help> scaling exponent (default is 3) </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isenthalpy" type="',class(self.isenthalpy),'" default="',convert2str(self.isenthalpy),'">','     <section name="thermal" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdynamicbasalspc" type="',class(self.isdynamicbasalspc),'" default="',convert2str(self.isdynamicbasalspc),'">','     <section name="thermal" />','     <help> enable dynamic setting of basal forcing. recommended for enthalpy formulation (default is 0)  </help>','  </parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(self.requested_outputs),'" default="',convert2str(self.requested_outputs),'">','     <section name="thermal" />','     <help> additional outputs requested </help>','  </parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');    
+
+                end % }}}
+                function self = extrude(self,md) % {{{
+                        self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
+                        if (length(md.initialization.temperature)==md.mesh.numberofvertices),
+                                self.spctemperature=NaN(md.mesh.numberofvertices,1);
+                                pos=find(md.mesh.vertexonsurface);
+                                self.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+                        end
+                end % }}}
                 function self = thermal(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/thermal.py

@@ -1 +1 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -43 +44 @@
                 return string
                 #}}}
+        def extrude(self,md): # {{{
+                self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
+                if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+                        self.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
+                        pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+                        self.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
+                return self
+        #}}}
         def defaultoutputs(self,md): # {{{
 

issm/trunk-jpl/src/m/dev/issmversion.m

@@ -17 +17 @@
 disp(['Compiled on ' IssmConfig('HOST_VENDOR') ' ' IssmConfig('HOST_OS') ' ' IssmConfig('HOST_ARCH') ' by ' IssmConfig('USER_NAME')]);
 disp([' ']);
-disp(['Copyright (c) 2009-2014 California Institute of Technology']);
+disp(['Copyright (c) 2009-2015 California Institute of Technology']);
 disp([' ']);
 disp(['    to get started type: issmdoc']);