Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 23430) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 23431) @@ -16,5 +16,5 @@ from basalforcings import basalforcings from matice import matice -from levelset import levelset +from levelset import levelset from calving import calving from fourierlove import fourierlove @@ -75,5 +75,5 @@ # classtype=model.properties - + # for classe in dict.keys(classtype): # print classe @@ -225,6 +225,6 @@ 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 + 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: '~HO.exp') an empty string '' will be considered as an empty domain @@ -349,5 +349,5 @@ md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1] - #Initial 2d mesh + #Initial 2d mesh if md1.mesh.__class__.__name__=='mesh3dprisms': flag_elem_2d=flag_elem[np.arange(0,md1.mesh.numberofelements2d)] @@ -431,5 +431,5 @@ if np.size(md1.stressbalance.spcvx)>1 and np.size(md1.stressbalance.spcvy)>2 and np.size(md1.stressbalance.spcvz)>2: if np.size(md1.inversion.vx_obs)>1 and np.size(md1.inversion.vy_obs)>1: - md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] + md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2] else: @@ -508,5 +508,5 @@ - 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: @@ -596,5 +596,5 @@ number_nodes3d=np.size(x3d) #number of 3d nodes for the non extruded part of the mesh - #Extrude elements + #Extrude elements elements3d=np.empty((0,6),int) for i in xrange(numlayers-1): @@ -618,5 +618,5 @@ md.mesh.upperelements=upperelements - #Save old mesh + #Save old mesh md.mesh.x2d=md.mesh.x md.mesh.y2d=md.mesh.y @@ -625,5 +625,5 @@ md.mesh.numberofvertices2d=md.mesh.numberofvertices - #Build global 3d mesh + #Build global 3d mesh md.mesh.elements=elements3d md.mesh.x=x3d @@ -672,5 +672,6 @@ md.materials.extrude(md) - md.damage.extrude(md) + if md.damage.isdamage==1: + md.damage.extrude(md) md.gia.extrude(md) md.mask.extrude(md) @@ -688,83 +689,83 @@ ''' 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") - + #dealing with the friction law #drag is limited to nodes that are on the bedrock. - if hasattr(md.friction,'coefficient'): + if hasattr(md.friction,'coefficient'): md.friction.coefficient=project2d(md,md.friction.coefficient,1) #p and q (same deal, except for element that are on the bedrock: ) - if hasattr(md.friction,'p'): + if hasattr(md.friction,'p'): md.friction.p=project2d(md,md.friction.p,1) - if hasattr(md.friction,'q'): + if hasattr(md.friction,'q'): md.friction.q=project2d(md,md.friction.q,1) - - if hasattr(md.friction,'coefficientcoulomb'): + + if hasattr(md.friction,'coefficientcoulomb'): md.friction.coefficientcoulomb=project2d(md,md.friction.coefficientcoulomb,1) - if hasattr(md.friction,'C'): + if hasattr(md.friction,'C'): md.friction.C=project2d(md,md.friction.C,1) - if hasattr(md.friction,'As'): + if hasattr(md.friction,'As'): md.friction.As=project2d(md,md.friction.As,1) if hasattr(md.friction,'effective_pressure') and not np.isnan(md.friction.effective_pressure).all(): md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1) - if hasattr(md.friction,'water_layer'): + if hasattr(md.friction,'water_layer'): md.friction.water_layer=project2d(md,md.friction.water_layer,1) - if hasattr(md.friction,'m'): + if hasattr(md.friction,'m'): md.friction.m=project2d(md,md.friction.m,1) #observations - if not np.isnan(md.inversion.vx_obs).all(): - md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) - if not np.isnan(md.inversion.vy_obs).all(): - md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) - if not np.isnan(md.inversion.vel_obs).all(): - md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) - if not np.isnan(md.inversion.cost_functions_coefficients).all(): - md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) + if not np.isnan(md.inversion.vx_obs).all(): + md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.vy_obs).all(): + md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.vel_obs).all(): + md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers) + if not np.isnan(md.inversion.cost_functions_coefficients).all(): + md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers) if isinstance(md.inversion.min_parameters,np.ndarray): - if md.inversion.min_parameters.size>1: - md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) + if md.inversion.min_parameters.size>1: + md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers) if isinstance(md.inversion.max_parameters,np.ndarray): - if md.inversion.max_parameters.size>1: - md.inversion.max_parameters=project2d(md,md.inversion.max_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 np.isnan(md.smb.mass_balance).all(): - md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) - + md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers) + #results - if not np.isnan(md.initialization.vx).all(): + if not np.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx) - if not np.isnan(md.initialization.vy).all(): + if not np.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy) - if not np.isnan(md.initialization.vz).all(): + if not np.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz) - if not np.isnan(md.initialization.vel).all(): + if not np.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel) - if not np.isnan(md.initialization.temperature).all(): + if not np.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature) - if not np.isnan(md.initialization.pressure).all(): + if not np.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1) - if not np.isnan(md.initialization.sediment_head).all(): + if not np.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1) - if not np.isnan(md.initialization.epl_head).all(): + if not np.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1) - if not np.isnan(md.initialization.epl_thickness).all(): + if not np.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1) #giaivins - if not np.isnan(md.gia.mantle_viscosity).all(): - md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) - if not np.isnan(md.gia.lithosphere_thickness).all(): - md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) + if not np.isnan(md.gia.mantle_viscosity).all(): + md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1) + if not np.isnan(md.gia.lithosphere_thickness).all(): + md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1) #elementstype @@ -777,10 +778,11 @@ # Hydrologydc variables - if hasattr(md.hydrology,'hydrologydc'): + if type(md.hydrology is 'hydrologydc'): md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1) - md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1) md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1) md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1) + md.hydrology.mask_thawed_node=project2d(md,md.hydrology.mask_thawed_node,1) if md.hydrology.isefficientlayer == 1: + md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1) md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1) @@ -793,5 +795,5 @@ md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers) md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1) - if not np.isnan(md.damage.spcdamage).all(): + if not np.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1) @@ -799,12 +801,12 @@ md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B) md.materials.rheology_n=project2d(md,md.materials.rheology_n,1) - - #damage: + + #damage: if md.damage.isdamage: md.damage.D=DepthAverage(md,md.damage.D) #special for thermal modeling: - md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1) - md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1) + md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1) + md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1) md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1) #bedrock only gets geothermal flux @@ -848,18 +850,18 @@ mesh.numberofelements=md.mesh.numberofelements2d mesh.elements=md.mesh.elements2d - if not np.isnan(md.mesh.vertexonboundary).all(): + if not np.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1) - if not np.isnan(md.mesh.elementconnectivity).all(): + if not np.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1) if isinstance(md.mesh.lat,np.ndarray): - if md.mesh.lat.size==md.mesh.numberofvertices: - mesh.lat=project2d(md,md.mesh.lat,1) + if md.mesh.lat.size==md.mesh.numberofvertices: + mesh.lat=project2d(md,md.mesh.lat,1) if isinstance(md.mesh.long,np.ndarray): - if md.mesh.long.size==md.mesh.numberofvertices: - md.mesh.long=project2d(md,md.mesh.long,1) + if md.mesh.long.size==md.mesh.numberofvertices: + md.mesh.long=project2d(md,md.mesh.long,1) mesh.epsg=md.mesh.epsg if isinstance(md.mesh.scale_factor,np.ndarray): - if md.mesh.scale_factor.size==md.mesh.numberofvertices: - md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1) + if md.mesh.scale_factor.size==md.mesh.numberofvertices: + md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1) md.mesh=mesh md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices)[0]