Changeset 13150

Timestamp:

08/23/12 15:49:08 (13 years ago)

Author:

jschierm

Message:

NEW: Python versions of model.extrude and project3d (and corresponding Matlab changes).

Location:

issm/trunk-jpl/src/m

Files:

• classes/model/model.m (modified) (1 diff)
• classes/model/model.py (modified) (3 diffs)
• extrusion/project3d.m (modified) (3 diffs)
• miscellaneous/fielddisplay.py (modified) (1 diff)

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

@@ -543 +543 @@
                          %Extrude the mesh
                          if nargin==2, %list of coefficients
-                                 list=varargin{1};
-                                 if any(list<0) | any(list>1),
+                                 clist=varargin{1};
+                                 if any(clist<0) | any(clist>1),
                                          error('extrusioncoefficients must be between 0 and 1');
                                  end
-                                 extrusionlist=sort(unique([list(:);0;1]));
+                                 extrusionlist=sort(unique([clist(:);0;1]));
                                  numlayers=length(extrusionlist);
                          elseif nargin==3, %one polynomial law

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

@@ -1 +1 @@
 #module imports {{{
+import numpy
 from mesh import mesh
 from mask import mask
@@ -37 +38 @@
 from mumpsoptions import *
 from iluasmoptions import *
+from project3d import *
 #}}}
 
@@ -167 +169 @@
         # }}}
 
+        def extrude(md,*args):    # {{{
+                """
+                EXTRUDE - vertically extrude a 2d mesh
+
+                   vertically extrude a 2d mesh and create corresponding 3d mesh.
+                   The vertical distribution can:
+                    - follow a polynomial law
+                    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
+                    - be discribed by a list of coefficients (between 0 and 1)
+ 
+
+                   Usage:
+                      md=extrude(md,numlayers,extrusionexponent);
+                      md=extrude(md,numlayers,lowerexponent,upperexponent);
+                      md=extrude(md,listofcoefficients);
+
+                   Example:
+                      md=extrude(md,8,3);
+                      md=extrude(md,8,3,2);
+                      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
+
+                   See also: MODELEXTRACT, COLLAPSE
+                """
+
+                #some checks on list of arguments
+                if len(args)>3 or len(args)<1:
+                        raise RuntimeError("extrude error message")
+
+                #Extrude the mesh
+                if   len(args)==1:    #list of coefficients
+                        clist=args[0]
+                        if any(clist<0) or any(clist>1):
+                                raise TypeError("extrusioncoefficients must be between 0 and 1")
+                        clist.extend([0.,1.])
+                        clist.sort()
+                        extrusionlist=list(set(clist))
+                        numlayers=len(extrusionlist)
+
+                elif len(args)==2:    #one polynomial law
+                        if args[1]<=0:
+                                raise TypeError("extrusionexponent must be >=0")
+                        numlayers=args[0]
+                        extrusionlist=(numpy.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1]
+
+                elif len(args)==3:    #two polynomial laws
+                        numlayers=args[0]
+                        lowerexp=args[1]
+                        upperexp=args[2]
+
+                        if args[1]<=0 or args[2]<=0:
+                                raise TypeError("lower and upper extrusionexponents must be >=0")
+
+                        lowerextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2.
+                        upperextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2.
+                        extrusionlist=numpy.unique(numpy.concatenate((lowerextrusionlist,1.-upperextrusionlist)))
+
+                if numlayers<2:
+                        raise TypeError("number of layers should be at least 2")
+                if md.mesh.dimension==3:
+                        raise TypeError("Cannot extrude a 3d mesh (extrude cannot be called more than once)")
+
+                #Initialize with the 2d mesh
+                x3d=numpy.empty((0))
+                y3d=numpy.empty((0))
+                z3d=numpy.empty((0))    #the lower node is on the bed
+                thickness3d=md.geometry.thickness    #thickness and bed for these nodes
+                bed3d=md.geometry.bed
+
+                #Create the new layers
+                for i in xrange(numlayers):
+                        x3d=numpy.concatenate((x3d,md.mesh.x))
+                        y3d=numpy.concatenate((y3d,md.mesh.y))
+                        #nodes are distributed between bed and surface accordingly to the given exponent
+                        z3d=numpy.concatenate((z3d,bed3d+thickness3d*extrusionlist[i]))
+                number_nodes3d=numpy.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
+
+                #Extrude elements 
+                elements3d=numpy.empty((0,6))
+                for i in xrange(numlayers-1):
+                        elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices))))    #Create the elements of the 3d mesh for the non extruded part
+                number_el3d=numpy.size(elements3d,axis=0)    #number of 3d nodes for the non extruded part of the mesh
+
+                #Keep a trace of lower and upper nodes
+                mesh.lowervertex=float('NaN')*numpy.ones(number_nodes3d)
+                mesh.uppervertex=float('NaN')*numpy.ones(number_nodes3d)
+                mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
+                mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
+                md.mesh.lowervertex=mesh.lowervertex
+                md.mesh.uppervertex=mesh.uppervertex
+
+                #same for lower and upper elements
+                mesh.lowerelements=float('NaN')*numpy.ones(number_el3d)
+                mesh.upperelements=float('NaN')*numpy.ones(number_el3d)
+                mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
+                mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
+                md.mesh.lowerelements=mesh.lowerelements
+                md.mesh.upperelements=mesh.upperelements
+
+                #Save old mesh 
+                md.mesh.x2d=md.mesh.x
+                md.mesh.y2d=md.mesh.y
+                md.mesh.elements2d=md.mesh.elements
+                md.mesh.numberofelements2d=md.mesh.numberofelements
+                md.mesh.numberofvertices2d=md.mesh.numberofvertices
+
+                #Update mesh type
+                md.mesh.dimension=3
+
+                #Build global 3d mesh 
+                md.mesh.elements=elements3d
+                md.mesh.x=x3d
+                md.mesh.y=y3d
+                md.mesh.z=z3d
+                md.mesh.numberofelements=number_el3d
+                md.mesh.numberofvertices=number_nodes3d
+                md.mesh.numberoflayers=numlayers
+
+                #Ok, now deal with the other fields from the 2d mesh:
+
+                #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.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node')
+                md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node')
+                md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node')
+                md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'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')
+
+                #bedinfo and surface info
+                md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',1)
+                md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',md.mesh.numberoflayers-1)
+                md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',1)
+                md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'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)
+                        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)
+                        md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
+
+                md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node')
+                md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node')
+                md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node')
+
+                #boundary conditions
+                md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node')
+                md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node')
+                md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node')
+                md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',float('NaN'))
+                md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node')
+                md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node')
+                md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node')
+
+                #in 3d, pressureload: [node1 node2 node3 node4 element]
+                pressureload_layer1=numpy.hstack((md.diagnostic.icefront[:,0:2],md.diagnostic.icefront[:,1:2]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,0:1]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,2:4]))    #Add two columns on the first layer 
+                pressureload=numpy.empty((0,6))
+                for i in xrange(numlayers-1):
+                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+i*md.mesh.numberofvertices2d,pressureload_layer1[:,4:5]+i*md.mesh.numberofelements2d,pressureload_layer1[:,5:6]))))
+                md.diagnostic.icefront=pressureload
+
+                #connectivity
+                md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
+                md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=float('NaN')
+                for i in xrange(1,numlayers):
+                        md.mesh.elementconnectivity[i*md.mesh.numberofelements2d+1:(i+1)*md.mesh.numberofelements2d,:] \
+                                =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d+1:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
+                md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=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')
+
+                #parameters
+                md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node')
+                md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node')
+                md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node')
+                md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node')
+                md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node')
+                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
+                md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element')
+                md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node')
+                md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element')
+                md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node')
+                md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element')
+                md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'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');end;
+                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')
+                if(md.surfaceforcings.isdelta18o):
+                        md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node')
+                if(md.surfaceforcings.isdelta18o):
+                        md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node')
+                if(md.surfaceforcings.isdelta18o):
+                        md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'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)
+
+                #special for thermal modeling:
+                md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.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
+
+                #increase connectivity if less than 25:
+                if md.mesh.average_vertex_connectivity<=25:
+                        md.mesh.average_vertex_connectivity=100
+
+                return md
+                # }}}
+

issm/trunk-jpl/src/m/extrusion/project3d.m

@@ -36 +36 @@
 if length(vector2d)==1,
         projected_vector=vector2d;
+
 elseif strcmpi(type,'node'),
 
@@ -57 +58 @@
                 projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d;
         end
+
 elseif strcmpi(type,'element'),
 
@@ -70 +72 @@
         end
 
+        %Fill in
         if layer==0,
                 for i=1:(md.mesh.numberoflayers-1),
                         projected_vector( ((i-1)*md.mesh.numberofelements2d+1):(i*md.mesh.numberofelements2d),:)=vector2d;
                 end
-
         else
                 projected_vector( ((layer-1)*md.mesh.numberofelements2d+1):(layer*md.mesh.numberofelements2d),:)=vector2d;
         end
+
 else
         error('project3d error message: unknown projection type');

issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py

@@ -78 +78 @@
 
         #initialization
-        if   isinstance(field,tuple):
+        if   isinstance(field,list):
+                sbeg='['
+                send=']'
+        elif isinstance(field,tuple):
                 sbeg='('
                 send=')'
-        elif isinstance(field,list):
-                sbeg='['
-                send=']'
         string=sbeg