Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 17723) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 17724) @@ -160,5 +160,5 @@ end - %Start with changing alle the fields from the 3d mesh + %Start with changing all the fields from the 3d mesh %drag is limited to nodes that are on the bedrock. @@ -244,11 +244,12 @@ %Initialize with the 2d mesh - md.mesh=mesh2d(); - md.mesh.x=md.mesh.x2d; - md.mesh.y=md.mesh.y2d; - md.mesh.z=zeros(size(md.mesh.x2d)); - md.mesh.numberofvertices=md.mesh.numberofvertices2d; - md.mesh.numberofelements=md.mesh.numberofelements2d; - md.mesh.elements=md.mesh.elements2d; + mesh=mesh2d(); + mesh.x=md.mesh.x2d; + mesh.y=md.mesh.y2d; + mesh.z=zeros(size(md.mesh.x2d)); + mesh.numberofvertices=md.mesh.numberofvertices2d; + mesh.numberofelements=md.mesh.numberofelements2d; + mesh.elements=md.mesh.elements2d; + md.mesh=mesh; %Keep a trace of lower and upper nodes Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 17723) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 17724) @@ -42,8 +42,8 @@ from miscellaneous import miscellaneous from private import private -from EnumDefinitions import * from mumpsoptions import mumpsoptions from iluasmoptions import iluasmoptions from project3d import project3d +from project2d import project2d from FlagElements import FlagElements from NodeConnectivity import NodeConnectivity @@ -51,4 +51,5 @@ from contourenvelope import contourenvelope import MatlabFuncs as m +from DepthAverage import DepthAverage #}}} @@ -175,5 +176,5 @@ # }}} def checkmessage(self,string): # {{{ - print ("model not consistent: %s" % string) + print "model not consistent: ", string self.private.isconsistent=False return self @@ -328,5 +329,5 @@ #Edges if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'): - if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs... + if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some npy.nans... #renumber first two columns pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0] @@ -394,6 +395,6 @@ md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2] else: - md2.stressbalance.spcvx[nodestoflag2]=float('NaN') - md2.stressbalance.spcvy[nodestoflag2]=float('NaN') + md2.stressbalance.spcvx[nodestoflag2]=npy.nan + md2.stressbalance.spcvy[nodestoflag2]=npy.nan print "\n!! extract warning: spc values should be checked !!\n\n" #put 0 for vz @@ -640,7 +641,7 @@ 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',float('NaN')) + 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=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.thermal.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1)) if hasattr(md.mesh,'vertexonsurface'): pos=numpy.nonzero(md.mesh.vertexonsurface)[0] @@ -703,2 +704,124 @@ return md # }}} + def collapse(md): #{{{ + ''' + collapses a 3d mesh into a 2d mesh + + This routine collapses a 3d model into a 2d model and collapses all + the fileds of the 3d model by taking their depth-averaged values + + Usage: + md=collapse(md) + ''' + + #Check that the model is really a 3d model + if md.mesh.domaintype().lower() != '3d': + raise StandardError("only a 3D model can be collapsed") + + #drag is limited to nodes that are on the bedrock. + md.friction.coefficient=project2d(md,md.friction.coefficient,1) + + #p and q (same deal, except for element that are on the bedrock: ) + md.friction.p=project2d(md,md.friction.p,1) + md.friction.q=project2d(md,md.friction.q,1) + + #observations + if not numpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) + if not numpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) + if not numpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) + if not numpy.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) + if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) + if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers) + if not numpy.isnan(md.surfaceforcings.mass_balance).all(): + md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers) + + if not numpy.isnan(md.balancethickness.thickening_rate).all(): md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers) + + #results + if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) + if not numpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) + if not numpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) + if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) + if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) + + #gia + if not numpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) + if not numpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) + + #elementstype + if not numpy.isnan(md.flowequation.element_equation).all(): + md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1) + md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1) + md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1) + md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1) + md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1) + + #boundary conditions + md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers) + md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers) + md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers) + md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers) + md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers) + md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers) + md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) + md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1) + + #materials + md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B) + md.materials.rheology_n=project2d(md,md.materials.rheology_n,1) + + #damage: + md.damage.D=DepthAverage(md,md.damage.D) + + #special for thermal modeling: + md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1) + md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1) #bedrock only gets geothermal flux + + #update of connectivity matrix + md.mesh.average_vertex_connectivity=25 + + #Collapse the mesh + nodes2d=md.mesh.numberofvertices2d + elements2d=md.mesh.numberofelements2d + + #parameters + md.geometry.surface=project2d(md,md.geometry.surface,1) + md.geometry.thickness=project2d(md,md.geometry.thickness,1) + md.geometry.base=project2d(md,md.geometry.base,1) + md.geometry.bed=project2d(md,md.geometry.bed,1) + md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) + md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) + md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1) + md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1) + + #lat long + if md.mesh.lat.size==md.mesh.numberofvertices: md.mesh.lat=project2d(md,md.mesh.lat,1) + if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1) + + #Initialize with the 2d mesh + mesh=mesh2d() + mesh.x=md.mesh.x2d + mesh.y=md.mesh.y2d + mesh.z=numpy.zeros_like(md.mesh.x2d) + mesh.numberofvertices=md.mesh.numberofvertices2d + mesh.numberofelements=md.mesh.numberofelements2d + mesh.elements=md.mesh.elements2d + md.mesh=mesh + + #Keep a trace of lower and upper nodes + md.mesh.lowervertex=numpy.nan + md.mesh.uppervertex=numpy.nan + md.mesh.lowerelements=numpy.nan + md.mesh.upperelements=numpy.nan + + #Remove old mesh + md.mesh.x2d=numpy.nan + md.mesh.y2d=numpy.nan + md.mesh.elements2d=numpy.nan + md.mesh.numberofelements2d=md.mesh.numberofelements + md.mesh.numberofvertices2d=md.mesh.numberofvertices + md.mesh.numberoflayers=0 + + return md + +#}}}