Ignore:
Timestamp:
01/29/15 08:44:39 (10 years ago)
Author:
Mathieu Morlighem
Message:

NEW: extrude is now a method of each class so that model.m does not need to know what class is being used

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19040 r19048  
    750750                        %Ok, now deal with the other fields from the 2d mesh:
    751751
     752                        %bedinfo and surface info
     753                        md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
     754                        md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
     755                        md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
     756
    752757                        %lat long
    753758                        md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
    754759                        md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
    755760
    756                         %drag coefficient is limited to nodes that are on the bedrock.
    757                         if isa(md.friction,'friction'),
    758                                 md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
    759                                 md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
    760                                 md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
    761                         elseif isa(md.friction,'frictionhydro'),
    762                                 md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
    763                                 md.friction.C=project3d(md,'vector',md.friction.C,'type','element');
    764                                 md.friction.As=project3d(md,'vector',md.friction.As,'type','element');
    765                                 md.friction.effective_pressure=project3d(md,'vector',md.friction.effective_pressure,'type','node','layer',1);
    766                         elseif isa(md.friction,'frictionwaterlayer'),
    767                                 md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
    768                                 md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
    769                                 md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
    770                                 md.friction.water_layer=project3d(md,'vector',md.friction.water_layer,'type','node','layer',1);
    771                         elseif isa(md.friction,'frictionweertman'),
    772                                 md.friction.C=project3d(md,'vector',md.friction.C,'type','node','layer',1);
    773                                 md.friction.m=project3d(md,'vector',md.friction.m,'type','element');
    774             end
    775          
    776                         %observations
    777                         md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node');
    778                         md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node');
    779                         md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node');
    780                         md.inversion.thickness_obs=project3d(md,'vector',md.inversion.thickness_obs,'type','node');
     761                        md.geometry=extrude(md.geometry,md);
     762                        md.friction  = extrude(md.friction,md);
     763                        md.inversion = extrude(md.inversion,md);
    781764                        md.surfaceforcings = extrude(md.surfaceforcings,md);
    782                         md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
    783 
    784                         %results
    785                         if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end;
    786                         if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end;
    787                         if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end;
    788                         if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end;
    789                         if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end;
    790                         if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end;
    791       if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);end;
    792       if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1);end;
    793       if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1);end;
    794       if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1);end;
    795 
    796                         %bedinfo and surface info
    797                         md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
    798                         md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
    799 
    800                         %elementstype
    801                         if ~isnan(md.flowequation.element_equation)
    802                                 oldelements_type=md.flowequation.element_equation;
    803                                 md.flowequation.element_equation=zeros(number_el3d,1);
    804                                 md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element');
    805                         end
    806 
    807                         %verticestype
    808                         if ~isnan(md.flowequation.vertex_equation)
    809                                 oldvertices_type=md.flowequation.vertex_equation;
    810                                 md.flowequation.vertex_equation=zeros(number_nodes3d,1);
    811                                 md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node');
    812                         end
    813                         md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node');
    814                         md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node');
    815                         md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node');
    816 
    817                         %boundary conditions
    818                         md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node');
    819                         md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node');
    820                         md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node');
    821                         md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
    822                         if (length(md.initialization.temperature)==md.mesh.numberofvertices),
    823                                 md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
    824                                 if isprop(md.mesh,'vertexonsurface'),
    825                                         pos=find(md.mesh.vertexonsurface);
    826                                         md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
    827                                 end
    828                         end
    829                         md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node');
    830                         md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
    831                         md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node');
    832                         md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
    833                         md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
    834 
    835                         % Calving variables
    836                         if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end;
    837                         if ~isnan(md.calving.meltingrate),md.calving.meltingrate=project3d(md,'vector',md.calving.meltingrate,'type','node');end;
    838                         if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
    839                         if isa(md.calving,'calvingpi'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
    840 
    841                         % Hydrologydc variables
    842                         if isa(md.hydrology,'hydrologydc');
    843                                 md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1);
    844                                 md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1);
    845                                 md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1);
    846                                 md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1);
    847                                 if(md.hydrology.isefficientlayer==1);
    848                                         md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1);
    849                     end
    850             end
     765                        md.initialization = extrude(md.initialization,md);
     766
     767                        md.flowequation.extrude(md);
     768                        md.stressbalance=extrude(md.stressbalance,md);
     769                        md.thermal.extrude(md);
     770                        md.masstransport.extrude(md);
     771                        md.calving=extrude(md.calving,md);
     772                        md.hydrology = extrude(md.hydrology,md);
    851773
    852774                        %connectivity
    … …  
    861783                        end
    862784
    863                         %materials
    864                         md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
    865                         md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
    866                                
    867                         %damage
    868                         md.damage.D=project3d(md,'vector',md.damage.D,'type','node');
    869 
    870                         %parameters
    871                         md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node');
    872                         md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node');
    873                         md.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node');
    874                         md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_thickness,'type','node');
    875                         md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node');
    876                         md.geometry.base=project3d(md,'vector',md.geometry.base,'type','node');
    877                         md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node');
    878                         md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
    879                         md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node');
    880                         md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node');
    881                         if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end;
    882                         if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end;
    883                         if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end;
    884                         if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end
    885 
    886                         %Put lithostatic pressure if there is an existing pressure
    887                         if ~isnan(md.initialization.pressure),
    888                                 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
    889             end
    890 
    891                         %special for thermal modeling:
    892                         md.basalforcings.groundedice_melting_rate=project3d(md,'vector',md.basalforcings.groundedice_melting_rate,'type','node','layer',1);
    893                         md.basalforcings.floatingice_melting_rate=project3d(md,'vector',md.basalforcings.floatingice_melting_rate,'type','node','layer',1);
    894                         if ~isnan(md.basalforcings.geothermalflux)
    895                                 md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
    896                         end
     785                        md.materials=extrude(md.materials,md);
     786                        md.damage=extrude(md.damage,md);
     787                        md.mask=extrude(md.mask,md);
     788                        md.qmu=extrude(md.qmu,md);
     789                        md.basalforcings=extrude(md.basalforcings,md);
    897790
    898791                        %increase connectivity if less than 25:
Note: See TracChangeset for help on using the changeset viewer.