Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 13149) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 13150) @@ -543,9 +543,9 @@ %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 Index: /issm/trunk-jpl/src/m/classes/model/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.py (revision 13149) +++ /issm/trunk-jpl/src/m/classes/model/model.py (revision 13150) @@ -1,3 +1,4 @@ #module imports {{{ +import numpy from mesh import mesh from mask import mask @@ -37,4 +38,5 @@ from mumpsoptions import * from iluasmoptions import * +from project3d import * #}}} @@ -167,2 +169,249 @@ # }}} + 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 + # }}} + Index: /issm/trunk-jpl/src/m/extrusion/project3d.m =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.m (revision 13149) +++ /issm/trunk-jpl/src/m/extrusion/project3d.m (revision 13150) @@ -36,4 +36,5 @@ if length(vector2d)==1, projected_vector=vector2d; + elseif strcmpi(type,'node'), @@ -57,4 +58,5 @@ projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d; end + elseif strcmpi(type,'element'), @@ -70,12 +72,13 @@ 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'); Index: /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py (revision 13149) +++ /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py (revision 13150) @@ -78,10 +78,10 @@ #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