Changeset 13857

Timestamp:

10/30/12 15:03:19 (12 years ago)

Author:

jschierm

Message:

NEW: Added python model.extract and contourenvelope (plus matlab cosmetic changes).

Location:

issm/trunk-jpl/src/m

Files:

• boundaryconditions/SetMarineIceSheetBC.py (modified) (4 diffs)
• classes/model/model.m (modified) (7 diffs)
• classes/model/model.py (modified) (3 diffs)
• classes/verbose.py (modified) (1 diff)
• mesh/bamg.py (modified) (5 diffs)
• parameterization/contourenvelope.m (modified) (5 diffs)
• parameterization/contourenvelope.py (added)

issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py

@@ -38 +38 @@
         pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]
         if not numpy.size(pos):
-                print("SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually.")
+                print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually."
 
         md.diagnostic.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
@@ -50 +50 @@
         #Dirichlet Values
         if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
-                print("      boundary conditions for diagnostic model: spc set as observed velocities")
+                print "      boundary conditions for diagnostic model: spc set as observed velocities"
                 md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
                 md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos]
         else:
-                print("      boundary conditions for diagnostic model: spc set as zero")
+                print "      boundary conditions for diagnostic model: spc set as zero"
 
         md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2))
@@ -83 +83 @@
         if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
                 md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
-                print("      no surfaceforcings.precipitation specified: values set as zero")
+                print "      no surfaceforcings.precipitation specified: values set as zero"
         if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
                 md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
-                print("      no surfaceforcings.mass_balance specified: values set as zero")
+                print "      no surfaceforcings.mass_balance specified: values set as zero"
         if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
                 md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-                print("      no basalforcings.melting_rate specified: values set as zero")
+                print "      no basalforcings.melting_rate specified: values set as zero"
         if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
                 md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
-                print("      no balancethickness.thickening_rate specified: values set as zero")
+                print "      no balancethickness.thickening_rate specified: values set as zero"
 
         md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
@@ -104 +104 @@
                 if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
                         md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
-                        md.basalforcings.geothermalflux[numpy.nonzero(md.mask.vertexongroundedice)]=50.*10.**-3; #50mW/m2
+                        md.basalforcings.geothermalflux[numpy.nonzero(md.mask.vertexongroundedice)]=50.*10.**-3    #50mW/m2
         else:
-                print("      no thermal boundary conditions created: no observed temperature found");
+                print "      no thermal boundary conditions created: no observed temperature found"
 
         return md

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

@@ -252 +252 @@
                         %   an empty string '' will be considered as an empty domain
                         %   a string 'all' will be considered as the entire domain
-                        %   add an argument 0 if you do not want the elements to be checked (faster)
                         %
                         %   Usage:
@@ -272 +271 @@
                         end
 
-                        %get check option
-                        if (nargin==3 & varargin{1}==0),
-                                checkoutline=0;
-                        else
-                                checkoutline=1;
-                        end
-
                         %get elements that are inside area
                         flag_elem=FlagElements(md1,area);
@@ -311 +303 @@
                         Pnode(pos_node)=[1:numberofvertices2]';
 
-                        %renumber the elements (some node won't exist anymore)
+                        %renumber the elements (some nodes won't exist anymore)
                         elements_1=md1.mesh.elements;
                         elements_2=elements_1(pos_elem,:);
@@ -323 +315 @@
                         end
 
-                        %OK, now create the new model !
-
-                        %take every fields from model
+                        %OK, now create the new model!
+
+                        %take every field from model
                         md2=md1;
 
@@ -347 +339 @@
                                                 elseif (fieldsize(1)==numberofvertices1+1)
                                                         md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
-                                                        %size = number of elements * n
+                                                %size = number of elements * n
                                                 elseif fieldsize(1)==numberofelements1
                                                         md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
@@ -358 +350 @@
                                         elseif (fieldsize(1)==numberofvertices1+1)
                                                 md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
-                                                %size = number of elements * n
+                                        %size = number of elements * n
                                         elseif fieldsize(1)==numberofelements1
                                                 md2.(model_fields{i})=field(pos_elem,:);
@@ -413 +405 @@
                                 %renumber first two columns
                                 pos=find(md2.mesh.edges(:,4)~=-1);
-                                md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1)); 
-                                md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2)); 
+                                md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1));
+                                md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2));
                                 md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
                                 md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
                                 %remove edges when the 2 vertices are not in the domain.
                                 md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
-                                %Replace all zeros by -1 in the last two columns;
+                                %Replace all zeros by -1 in the last two columns
                                 pos=find(md2.mesh.edges(:,3)==0);
                                 md2.mesh.edges(pos,3)=-1;

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

@@ -1 +1 @@
 #module imports {{{
 import numpy
+import copy
 from mesh import mesh
 from mask import mask
@@ -39 +40 @@
 from iluasmoptions import *
 from project3d import *
+from FlagElements import *
+from NodeConnectivity import *
+from ElementConnectivity import *
+from contourenvelope import *
 #}}}
 
@@ -170 +175 @@
         # }}}
 
-        """
-        function md2 = extract(md,area) % {{{
-                %extract - extract a model according to an Argus contour or flag list
-                %
-                %   This routine extracts a submodel from a bigger model with respect to a given contour
-                %   md must be followed by the corresponding exp file or flags list
-                %   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
-                %   If user wants every element outside the domain to be 
-                %   extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp');
-                %   an empty string '' will be considered as an empty domain
-                %   a string 'all' will be considered as the entire domain
-                %   add an argument 0 if you do not want the elements to be checked (faster)
-                %
-                %   Usage:
-                %      md2=extract(md,area);
-                %
-                %   Examples:
-                %      md2=extract(md,'Domain.exp');
-                %      md2=extract(md,md.mask.elementonfloatingice);
-                %
-                %   See also: EXTRUDE, COLLAPSE
-
-                %copy model
-                md1=md;
-
-                %some checks
-                if ((nargin~=2) | (nargout~=1)),
-                        help extract
-                        error('extract error message: bad usage');
-                end
-
-                %get check option
-                if (nargin==3 & varargin{1}==0),
-                        checkoutline=0;
-                else
-                        checkoutline=1;
-                end
-
-                %get elements that are inside area
-                flag_elem=FlagElements(md1,area);
-                if ~any(flag_elem),
-                        error('extracted model is empty');
-                end
-
-                %kick out all elements with 3 dirichlets
-                spc_elem=find(~flag_elem);
-                spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
-                flag=ones(md1.mesh.numberofvertices,1);
-                flag(spc_node)=0;
-                pos=find(sum(flag(md1.mesh.elements),2)==0);
-                flag_elem(pos)=0;
-
-                %extracted elements and nodes lists
-                pos_elem=find(flag_elem);
-                pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));
-
-                %keep track of some fields
-                numberofvertices1=md1.mesh.numberofvertices;
-                numberofelements1=md1.mesh.numberofelements;
-                numberofvertices2=length(pos_node);
-                numberofelements2=length(pos_elem);
-                flag_node=zeros(numberofvertices1,1);
-                flag_node(pos_node)=1;
-
-                %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
-                Pelem=zeros(numberofelements1,1);
-                Pelem(pos_elem)=[1:numberofelements2]';
-                Pnode=zeros(numberofvertices1,1);
-                Pnode(pos_node)=[1:numberofvertices2]';
-
-                %renumber the elements (some node won't exist anymore)
-                elements_1=md1.mesh.elements;
-                elements_2=elements_1(pos_elem,:);
-                elements_2(:,1)=Pnode(elements_2(:,1));
-                elements_2(:,2)=Pnode(elements_2(:,2));
-                elements_2(:,3)=Pnode(elements_2(:,3));
-                if md1.mesh.dimension==3,
-                        elements_2(:,4)=Pnode(elements_2(:,4));
-                        elements_2(:,5)=Pnode(elements_2(:,5));
-                        elements_2(:,6)=Pnode(elements_2(:,6));
-                end
-
-                %OK, now create the new model !
-
-                %take every fields from model
-                md2=md1;
-
-                %automatically modify fields
-
-                %loop over model fields
-                model_fields=fields(md1);
-                for i=1:length(model_fields),
-                        %get field
-                        field=md1.(model_fields{i});
-                        fieldsize=size(field);
-                        if isobject(field), %recursive call
-                                object_fields=fields(md1.(model_fields{i}));
-                                for j=1:length(object_fields),
-                                        %get field
-                                        field=md1.(model_fields{i}).(object_fields{j});
-                                        fieldsize=size(field);
-                                        %size = number of nodes * n
-                                        if fieldsize(1)==numberofvertices1
-                                                md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
-                                        elseif (fieldsize(1)==numberofvertices1+1)
-                                                md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
-                                                %size = number of elements * n
-                                        elseif fieldsize(1)==numberofelements1
-                                                md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
-                                        end
-                                end
-                        else
-                                %size = number of nodes * n
-                                if fieldsize(1)==numberofvertices1
-                                        md2.(model_fields{i})=field(pos_node,:);
-                                elseif (fieldsize(1)==numberofvertices1+1)
-                                        md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
-                                        %size = number of elements * n
-                                elseif fieldsize(1)==numberofelements1
-                                        md2.(model_fields{i})=field(pos_elem,:);
-                                end
-                        end
-                end
-
-                %modify some specific fields
-
-                %Mesh
-                md2.mesh.numberofelements=numberofelements2;
-                md2.mesh.numberofvertices=numberofvertices2;
-                md2.mesh.elements=elements_2;
-
-                %mesh.uppervertex mesh.lowervertex
-                if md1.mesh.dimension==3
-                        md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
-                        pos=find(~isnan(md2.mesh.uppervertex));
-                        md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));
-
-                        md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
-                        pos=find(~isnan(md2.mesh.lowervertex));
-                        md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));
-
-                        md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
-                        pos=find(~isnan(md2.mesh.upperelements));
-                        md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));
-
-                        md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
-                        pos=find(~isnan(md2.mesh.lowerelements));
-                        md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
-                end
-
-                %Initial 2d mesh 
-                if md1.mesh.dimension==3
-                        flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
-                        pos_elem_2d=find(flag_elem_2d);
-                        flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
-                        pos_node_2d=find(flag_node_2d);
-
-                        md2.mesh.numberofelements2d=length(pos_elem_2d);
-                        md2.mesh.numberofvertices2d=length(pos_node_2d);
-                        md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
-                        md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
-                        md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
-                        md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));
-
-                        md2.mesh.x2d=md1.mesh.x(pos_node_2d);
-                        md2.mesh.y2d=md1.mesh.y(pos_node_2d);
-                end
-
-                %Edges
-                if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
-                        %renumber first two columns
-                        pos=find(md2.mesh.edges(:,4)~=-1);
-                        md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1)); 
-                        md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2)); 
-                        md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
-                        md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
-                        %remove edges when the 2 vertices are not in the domain.
-                        md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
-                        %Replace all zeros by -1 in the last two columns;
-                        pos=find(md2.mesh.edges(:,3)==0);
-                        md2.mesh.edges(pos,3)=-1;
-                        pos=find(md2.mesh.edges(:,4)==0);
-                        md2.mesh.edges(pos,4)=-1;
-                        %Invert -1 on the third column with last column (Also invert first two columns!!)
-                        pos=find(md2.mesh.edges(:,3)==-1);
-                        md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
-                        md2.mesh.edges(pos,4)=-1;
-                        values=md2.mesh.edges(pos,2);
-                        md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
-                        md2.mesh.edges(pos,1)=values;
-                        %Finally remove edges that do not belong to any element
-                        pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
-                        md2.mesh.edges(pos,:)=[];
-                end
-
-                %Penalties
-                if ~isnan(md2.diagnostic.vertex_pairing),
-                        for i=1:size(md1.diagnostic.vertex_pairing,1);
-                                md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:));
-                        end
-                        md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:);
-                end
-                if ~isnan(md2.prognostic.vertex_pairing),
-                        for i=1:size(md1.prognostic.vertex_pairing,1);
-                                md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:));
-                        end
-                        md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:);
-                end
-
-                %recreate segments
-                if md1.mesh.dimension==2
-                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
-                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
-                        md2.mesh.segments=contourenvelope(md2);
-                        md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
-                else
-                        %First do the connectivity for the contourenvelope in 2d
-                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
-                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
-                        md2.mesh.segments=contourenvelope(md2);
-                        md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
-                        md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
-                        %Then do it for 3d as usual
-                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
-                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
-                end
-
-                %Boundary conditions: Dirichlets on new boundary
-                %Catch the elements that have not been extracted
-                orphans_elem=find(~flag_elem);
-                orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
-                %Figure out which node are on the boundary between md2 and md1
-                nodestoflag1=intersect(orphans_node,pos_node);
-                nodestoflag2=Pnode(nodestoflag1);
-                if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2,
-                        if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
-                                md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 
-                                md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
-                        else
-                                md2.diagnostic.spcvx(nodestoflag2)=NaN;
-                                md2.diagnostic.spcvy(nodestoflag2)=NaN;
-                                disp(' ')
-                                disp('!! extract warning: spc values should be checked !!')
-                                disp(' ')
-                        end
-                        %put 0 for vz
-                        md2.diagnostic.spcvz(nodestoflag2)=0;
-                end
-                if ~isnan(md1.thermal.spctemperature),
-                        md2.thermal.spctemperature(nodestoflag2,1)=1;
-                end
-
-                %Diagnostic
-                if ~isnan(md2.diagnostic.icefront)
-                        md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); 
-                        md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); 
-                        md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1));
-                        if md1.mesh.dimension==3
-                                md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); 
-                                md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); 
-                        end
-                        md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:);
-                end
-
-                %Results fields
-                if isstruct(md1.results),
-                        md2.results=struct();
-                        solutionfields=fields(md1.results);
-                        for i=1:length(solutionfields),
-                                %get subfields
-                                solutionsubfields=fields(md1.results.(solutionfields{i}));
-                                for j=1:length(solutionsubfields),
-                                        field=md1.results.(solutionfields{i}).(solutionsubfields{j});
-                                        if length(field)==numberofvertices1,
-                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node);
-                                        elseif length(field)==numberofelements1,
-                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem);
-                                        else
-                                                md2.results.(solutionfields{i}).(solutionsubfields{j})=field;
-                                        end
-                                end
-                        end
-                end
-
-                %Keep track of pos_node and pos_elem
-                md2.mesh.extractedvertices=pos_node;
-                md2.mesh.extractedelements=pos_elem;
-        end % }}}
-        """
+        def extract(md,area):    # {{{
+                """
+                extract - extract a model according to an Argus contour or flag list
+
+                   This routine extracts a submodel from a bigger model with respect to a given contour
+                   md must be followed by the corresponding exp file or flags list
+                   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
+                   If user wants every element outside the domain to be 
+                   extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp');
+                   an empty string '' will be considered as an empty domain
+                   a string 'all' will be considered as the entire domain
+
+                   Usage:
+                      md2=extract(md,area);
+
+                   Examples:
+                      md2=extract(md,'Domain.exp');
+                      md2=extract(md,md.mask.elementonfloatingice);
+
+                   See also: EXTRUDE, COLLAPSE
+                """
+
+                #copy model
+                md1=copy.deepcopy(md)
+
+                #get elements that are inside area
+                flag_elem=FlagElements(md1,area)
+                if not numpy.any(flag_elem):
+                        raise RuntimeError("extracted model is empty")
+
+                #kick out all elements with 3 dirichlets
+                spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
+                spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1
+                flag=numpy.ones(md1.mesh.numberofvertices)
+                flag[spc_node]=0
+                pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0]
+                flag_elem[pos]=0
+
+                #extracted elements and nodes lists
+                pos_elem=numpy.nonzero(flag_elem)[0]
+                pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1
+
+                #keep track of some fields
+                numberofvertices1=md1.mesh.numberofvertices
+                numberofelements1=md1.mesh.numberofelements
+                numberofvertices2=numpy.size(pos_node)
+                numberofelements2=numpy.size(pos_elem)
+                flag_node=numpy.zeros(numberofvertices1)
+                flag_node[pos_node]=1
+
+                #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
+                Pelem=numpy.zeros(numberofelements1)
+                Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
+                Pnode=numpy.zeros(numberofvertices1)
+                Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
+
+                #renumber the elements (some node won't exist anymore)
+                elements_1=copy.deepcopy(md1.mesh.elements)
+                elements_2=elements_1[pos_elem,:]
+                elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1]
+                elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1]
+                elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1]
+                if md1.mesh.dimension==3:
+                        elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1]
+                        elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1]
+                        elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1]
+
+                #OK, now create the new model!
+
+                #take every field from model
+                md2=copy.deepcopy(md1)
+
+                #automatically modify fields
+
+                #loop over model fields
+                model_fields=vars(md1)
+                for fieldi in model_fields:
+                        #get field
+                        field=getattr(md1,fieldi)
+                        fieldsize=numpy.shape(field)
+                        if hasattr(field,'__dict__') and not ismember(fieldi,['results'])[0]:    #recursive call
+                                object_fields=vars(field)
+                                for fieldj in object_fields:
+                                        #get field
+                                        field=getattr(getattr(md1,fieldi),fieldj)
+                                        fieldsize=numpy.shape(field)
+                                        if len(fieldsize):
+                                                #size = number of nodes * n
+                                                if   fieldsize[0]==numberofvertices1:
+                                                        setattr(getattr(md2,fieldi),fieldj,field[pos_node,:])
+                                                elif fieldsize[0]==numberofvertices1+1:
+                                                        setattr(getattr(md2,fieldi),fieldj,numpy.vstack((field[pos_node,:],field[-1,:])))
+                                                #size = number of elements * n
+                                                elif fieldsize[0]==numberofelements1:
+                                                        setattr(getattr(md2,fieldi),fieldj,field[pos_elem,:])
+                        else:
+                                if len(fieldsize):
+                                        #size = number of nodes * n
+                                        if   fieldsize[0]==numberofvertices1:
+                                                setattr(md2,fieldi,field[pos_node,:])
+                                        elif fieldsize[0]==numberofvertices1+1:
+                                                setattr(md2,fieldi,numpy.hstack((field[pos_node,:],field[-1,:])))
+                                        #size = number of elements * n
+                                        elif fieldsize[0]==numberofelements1:
+                                                setattr(md2,fieldi,field[pos_elem,:])
+
+                #modify some specific fields
+
+                #Mesh
+                md2.mesh.numberofelements=numberofelements2
+                md2.mesh.numberofvertices=numberofvertices2
+                md2.mesh.elements=elements_2
+
+                #mesh.uppervertex mesh.lowervertex
+                if md1.mesh.dimension==3:
+                        md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
+                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.uppervertex)))[0]
+                        md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1]
+
+                        md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node]
+                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowervertex)))[0]
+                        md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1]
+
+                        md2.mesh.upperelements=md1.mesh.upperelements[pos_elem]
+                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.upperelements)))[0]
+                        md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1]
+
+                        md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem]
+                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowerelements)))[0]
+                        md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1]
+
+                #Initial 2d mesh 
+                if md1.mesh.dimension==3:
+                        flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)]
+                        pos_elem_2d=numpy.nonzero(flag_elem_2d)[0]
+                        flag_node_2d=flag_node[numpy.arange(0,md1.mesh.numberofvertices2d)]
+                        pos_node_2d=numpy.nonzero(flag_node_2d)[0]
+
+                        md2.mesh.numberofelements2d=numpy.size(pos_elem_2d)
+                        md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
+                        md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
+                        md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1]
+                        md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1]
+                        md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1]
+
+                        md2.mesh.x2d=md1.mesh.x[pos_node_2d]
+                        md2.mesh.y2d=md1.mesh.y[pos_node_2d]
+
+                #Edges
+                if len(numpy.shape(md2.mesh.edges))>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
+                        #renumber first two columns
+                        pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
+                        md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0].astype(int)-1]
+                        md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1]
+                        md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1]
+                        md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1]
+                        #remove edges when the 2 vertices are not in the domain.
+                        md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
+                        #Replace all zeros by -1 in the last two columns
+                        pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0]
+                        md2.mesh.edges[pos,2]=-1
+                        pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0]
+                        md2.mesh.edges[pos,3]=-1
+                        #Invert -1 on the third column with last column (Also invert first two columns!!)
+                        pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0]
+                        md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3]
+                        md2.mesh.edges[pos,3]=-1
+                        values=md2.mesh.edges[pos,1]
+                        md2.mesh.edges[pos,1]=md2.mesh.edges[pos,0]
+                        md2.mesh.edges[pos,0]=values
+                        #Finally remove edges that do not belong to any element
+                        pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
+                        md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0)
+
+                #Penalties
+                if numpy.any(numpy.logical_not(numpy.isnan(md2.diagnostic.vertex_pairing))):
+                        for i in xrange(numpy.size(md1.diagnostic.vertex_pairing,axis=0)):
+                                md2.diagnostic.vertex_pairing[i,:]=Pnode[md1.diagnostic.vertex_pairing[i,:]]
+                        md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing[numpy.nonzero(md2.diagnostic.vertex_pairing[:,0])[0],:]
+                if numpy.any(numpy.logical_not(numpy.isnan(md2.prognostic.vertex_pairing))):
+                        for i in xrange(numpy.size(md1.prognostic.vertex_pairing,axis=0)):
+                                md2.prognostic.vertex_pairing[i,:]=Pnode[md1.prognostic.vertex_pairing[i,:]]
+                        md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing[numpy.nonzero(md2.prognostic.vertex_pairing[:,0])[0],:]
+
+                #recreate segments
+                if md1.mesh.dimension==2:
+                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
+                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
+                        md2.mesh.segments=contourenvelope(md2)
+                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2)
+                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+                else:
+                        #First do the connectivity for the contourenvelope in 2d
+                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d)
+                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
+                        md2.mesh.segments=contourenvelope(md2)
+                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers)
+                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+                        md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
+                        #Then do it for 3d as usual
+                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
+                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
+
+                #Boundary conditions: Dirichlets on new boundary
+                #Catch the elements that have not been extracted
+                orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
+                orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1
+                #Figure out which node are on the boundary between md2 and md1
+                nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
+                nodestoflag2=Pnode[nodestoflag1].astype(int)-1
+                if numpy.size(md1.diagnostic.spcvx)>1 and numpy.size(md1.diagnostic.spcvy)>2 and numpy.size(md1.diagnostic.spcvz)>2:
+                        if numpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1:
+                                md2.diagnostic.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] 
+                                md2.diagnostic.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
+                        else:
+                                md2.diagnostic.spcvx[nodestoflag2]=float('NaN')
+                                md2.diagnostic.spcvy[nodestoflag2]=float('NaN')
+                                print "\n!! extract warning: spc values should be checked !!\n\n"
+                        #put 0 for vz
+                        md2.diagnostic.spcvz[nodestoflag2]=0
+                if numpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))):
+                        md2.thermal.spctemperature[nodestoflag2,0]=1
+
+                #Diagnostic
+                if numpy.any(numpy.logical_not(numpy.isnan(md2.diagnostic.icefront))):
+                        md2.diagnostic.icefront[:,0]=Pnode[md1.diagnostic.icefront[:,0].astype(int)-1]
+                        md2.diagnostic.icefront[:,1]=Pnode[md1.diagnostic.icefront[:,1].astype(int)-1]
+                        md2.diagnostic.icefront[:,-2]=Pelem[md1.diagnostic.icefront[:,-2].astype(int)-1]
+                        if md1.mesh.dimension==3:
+                                md2.diagnostic.icefront[:,2]=Pnode[md1.diagnostic.icefront[:,2].astype(int)-1]
+                                md2.diagnostic.icefront[:,3]=Pnode[md1.diagnostic.icefront[:,3].astype(int)-1]
+                        md2.diagnostic.icefront=md2.diagnostic.icefront[numpy.nonzero(numpy.logical_and(numpy.logical_and(md2.diagnostic.icefront[:,0],md2.diagnostic.icefront[:,1]),md2.diagnostic.icefront[:,-1]))[0],:]
+
+                #Results fields
+                if md1.results:
+                        md2.results=OrderedDict()
+                        for solutionfield in md1.results.iterkeys():
+                                #get time step
+                                for i in m1.results[solutionfield].iterkeys():
+                                        #get subfields
+                                        for solutionsubfield,field in m1.results[solutionfield][i].iteritems():
+                                                if   numpy.size(field)==numberofvertices1:
+                                                        md2.results[solutionfield][i][solutionsubfield]=field[pos_node]
+                                                elif numpy.size(field)==numberofelements1:
+                                                        md2.results[solutionfield][i][solutionsubfield]=field[pos_elem]
+                                                else:
+                                                        md2.results[solutionfield][i][solutionsubfield]=field
+
+                #Keep track of pos_node and pos_elem
+                md2.mesh.extractedvertices=pos_node.astype(float)+1
+                md2.mesh.extractedelements=pos_elem.astype(float)+1
+
+                return md2
+        # }}}
 
         def extrude(md,*args):    # {{{

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

@@ -61 +61 @@
                         #Cast to logicals
                         listproperties=vars(self)
-                        for [fieldname,fieldvalue] in listproperties.iteritems():
+                        for fieldname,fieldvalue in listproperties.iteritems():
                                 if isinstance(fieldvalue,bool) or isinstance(fieldvalue,(int,long,float)):
                                         setattr(self,fieldname,bool(fieldvalue))

issm/trunk-jpl/src/m/mesh/bamg.py

@@ -123 +123 @@
                                 elif numpy.any(numpy.logical_not(flags)):
                                         #We LOTS of work to do
-                                        print("Rift tip outside of or on the domain has been detected and is being processed...")
+                                        print "Rift tip outside of or on the domain has been detected and is being processed..."
 
                                         #check that only one point is outside (for now)
@@ -166 +166 @@
 
                                                         if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0):
-                                                                print("moving tip-domain intersection point")
+                                                                print "moving tip-domain intersection point"
 
                                                                 #Get position of the closer point
@@ -350 +350 @@
         raise RuntimeError("bamg.py/processgeometry is not complete.")
         #Deal with edges
-        print("Checking Edge crossing...")
+        print "Checking Edge crossing..."
         i=0
         while (i<numpy.size(geom.Edges,axis=0)):
@@ -405 +405 @@
 
         #Check point outside
-        print("Checking for points outside the domain...")
+        print "Checking for points outside the domain..."
         i=0
         num=0
@@ -435 +435 @@
 
         if num:
-                print("WARNING: %d points outside the domain outline have been removed" % num)
+                print "WARNING: %d points outside the domain outline have been removed" % num
 
         """

issm/trunk-jpl/src/m/parameterization/contourenvelope.m

@@ -20 +20 @@
         if ischar(flags),
                 file=flags;
-                file=varargin{1};
                 if ~exist(file),
                         error(['contourenvelope error message: file ' file ' not found']);
@@ -29 +28 @@
                 isfile=0;
         else
-                error('contourenvelope error message:  second argument should a file or an elements flag');
+                error('contourenvelope error message:  second argument should be a file or an elements flag');
         end
 end
@@ -70 +69 @@
         else
                 %get flag list of elements and nodes inside the contour
-                nodein=zeros(mesh.numberofvertices,1); 
-                elemin=zeros(mesh.numberofelements,1); 
+                nodein=zeros(mesh.numberofvertices,1);
+                elemin=zeros(mesh.numberofelements,1);
 
-                pos=find(flags); 
+                pos=find(flags);
                 elemin(pos)=1;
                 nodein(mesh.elements(pos,:))=1;
@@ -95 +94 @@
 pos=find(elementonboundary);
 num_segments=length(pos);
-segments=zeros(num_segments,3);
+segments=zeros(num_segments*3,3);
 count=1;
 
@@ -138 +137 @@
         end
 end
+segments=segments(1:count-1,:);
+