Changeset 17724


Ignore:
Timestamp:
04/14/14 15:18:49 (11 years ago)
Author:
cborstad
Message:

CHG: added collapse function to python model, fixed bug when initializing 2d mesh wiped out old mesh before assignment

Location:
issm/trunk-jpl/src/m/classes
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r17687 r17724  
    160160                        end
    161161
    162                         %Start with changing alle the fields from the 3d mesh
     162                        %Start with changing all the fields from the 3d mesh
    163163
    164164                        %drag is limited to nodes that are on the bedrock.
    … …  
    244244
    245245                        %Initialize with the 2d mesh
    246                         md.mesh=mesh2d();
    247                         md.mesh.x=md.mesh.x2d;
    248                         md.mesh.y=md.mesh.y2d;
    249                         md.mesh.z=zeros(size(md.mesh.x2d));
    250                         md.mesh.numberofvertices=md.mesh.numberofvertices2d;
    251                         md.mesh.numberofelements=md.mesh.numberofelements2d;
    252                         md.mesh.elements=md.mesh.elements2d;
     246                        mesh=mesh2d();
     247                        mesh.x=md.mesh.x2d;
     248                        mesh.y=md.mesh.y2d;
     249                        mesh.z=zeros(size(md.mesh.x2d));
     250                        mesh.numberofvertices=md.mesh.numberofvertices2d;
     251                        mesh.numberofelements=md.mesh.numberofelements2d;
     252                        mesh.elements=md.mesh.elements2d;
     253                        md.mesh=mesh;
    253254
    254255                        %Keep a trace of lower and upper nodes

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

    r17714 r17724  
    4242from miscellaneous import miscellaneous
    4343from private import private
    44 from EnumDefinitions import *
    4544from mumpsoptions import mumpsoptions
    4645from iluasmoptions import iluasmoptions
    4746from project3d import project3d
     47from project2d import project2d
    4848from FlagElements import FlagElements
    4949from NodeConnectivity import NodeConnectivity
    … …  
    5151from contourenvelope import contourenvelope
    5252import MatlabFuncs as m
     53from DepthAverage import DepthAverage
    5354#}}}
    5455
    … …  
    175176        # }}}
    176177        def checkmessage(self,string):    # {{{
    177                 print ("model not consistent: %s" % string)
     178                print "model not consistent: ", string
    178179                self.private.isconsistent=False
    179180                return self
    … …  
    328329                #Edges
    329330                if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
    330                         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...
     331                        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...
    331332                                #renumber first two columns
    332333                                pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
    … …  
    394395                                md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
    395396                        else:
    396                                 md2.stressbalance.spcvx[nodestoflag2]=float('NaN')
    397                                 md2.stressbalance.spcvy[nodestoflag2]=float('NaN')
     397                                md2.stressbalance.spcvx[nodestoflag2]=npy.nan
     398                                md2.stressbalance.spcvy[nodestoflag2]=npy.nan
    398399                                print "\n!! extract warning: spc values should be checked !!\n\n"
    399400                        #put 0 for vz
    … …  
    640641                md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node')
    641642                md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node')
    642                 md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',float('NaN'))
     643                md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
    643644                if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
    644                         md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
     645                        md.thermal.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
    645646                        if hasattr(md.mesh,'vertexonsurface'):
    646647                                pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
    … …  
    703704                return md
    704705                # }}}
     706        def collapse(md): #{{{
     707                '''
     708                collapses a 3d mesh into a 2d mesh
     709                       
     710                This routine collapses a 3d model into a 2d model and collapses all
     711                the fileds of the 3d model by taking their depth-averaged values
     712                       
     713                Usage:
     714                        md=collapse(md)
     715                '''     
     716
     717                #Check that the model is really a 3d model
     718                if md.mesh.domaintype().lower() != '3d':
     719                        raise StandardError("only a 3D model can be collapsed")
     720               
     721                #drag is limited to nodes that are on the bedrock.
     722                md.friction.coefficient=project2d(md,md.friction.coefficient,1)
     723
     724                #p and q (same deal, except for element that are on the bedrock: )
     725                md.friction.p=project2d(md,md.friction.p,1)
     726                md.friction.q=project2d(md,md.friction.q,1)
     727
     728                #observations
     729                if not numpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)
     730                if not numpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)
     731                if not numpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)
     732                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)
     733                if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers)
     734                if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers)
     735                if not numpy.isnan(md.surfaceforcings.mass_balance).all():
     736                        md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers)
     737               
     738                if not numpy.isnan(md.balancethickness.thickening_rate).all(): md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers)
     739
     740                #results
     741                if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)
     742                if not numpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy)
     743                if not numpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz)
     744                if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)
     745                if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
     746
     747                #gia
     748                if not numpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)
     749                if not numpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)
     750
     751                #elementstype
     752                if not numpy.isnan(md.flowequation.element_equation).all():
     753                        md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1)
     754                        md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1)
     755                        md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1)
     756                        md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1)
     757                        md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1)
     758
     759                #boundary conditions
     760                md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers)
     761                md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers)
     762                md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers)
     763                md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers)
     764                md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers)
     765                md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers)
     766                md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)
     767                md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1)
     768
     769                #materials
     770                md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B)
     771                md.materials.rheology_n=project2d(md,md.materials.rheology_n,1)
     772               
     773                #damage:
     774                md.damage.D=DepthAverage(md,md.damage.D)
     775
     776                #special for thermal modeling:
     777                md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1)
     778                md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1) #bedrock only gets geothermal flux
     779
     780                #update of connectivity matrix
     781                md.mesh.average_vertex_connectivity=25
     782
     783                #Collapse the mesh
     784                nodes2d=md.mesh.numberofvertices2d
     785                elements2d=md.mesh.numberofelements2d
     786
     787                #parameters
     788                md.geometry.surface=project2d(md,md.geometry.surface,1)
     789                md.geometry.thickness=project2d(md,md.geometry.thickness,1)
     790                md.geometry.base=project2d(md,md.geometry.base,1)
     791                md.geometry.bed=project2d(md,md.geometry.bed,1)
     792                md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)
     793                md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)
     794                md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1)
     795                md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1)
     796
     797                #lat long
     798                if md.mesh.lat.size==md.mesh.numberofvertices:  md.mesh.lat=project2d(md,md.mesh.lat,1)
     799                if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1)
     800
     801                #Initialize with the 2d mesh
     802                mesh=mesh2d()
     803                mesh.x=md.mesh.x2d
     804                mesh.y=md.mesh.y2d
     805                mesh.z=numpy.zeros_like(md.mesh.x2d)
     806                mesh.numberofvertices=md.mesh.numberofvertices2d
     807                mesh.numberofelements=md.mesh.numberofelements2d
     808                mesh.elements=md.mesh.elements2d
     809                md.mesh=mesh
     810
     811                #Keep a trace of lower and upper nodes
     812                md.mesh.lowervertex=numpy.nan
     813                md.mesh.uppervertex=numpy.nan
     814                md.mesh.lowerelements=numpy.nan
     815                md.mesh.upperelements=numpy.nan
     816
     817                #Remove old mesh
     818                md.mesh.x2d=numpy.nan
     819                md.mesh.y2d=numpy.nan
     820                md.mesh.elements2d=numpy.nan
     821                md.mesh.numberofelements2d=md.mesh.numberofelements
     822                md.mesh.numberofvertices2d=md.mesh.numberofvertices
     823                md.mesh.numberoflayers=0
     824
     825                return md
     826
     827#}}}
Note: See TracChangeset for help on using the changeset viewer.