Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 13004) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 13005) @@ -91,4 +91,667 @@ end %}}} + function md = collapse(md)% {{{ + %COLLAPSE - 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) + % + % See also: EXTRUDE, MODELEXTRACT + + %Check that the model is really a 3d model + if ~md.mesh.dimension==3, + error('collapse error message: only 3d mesh can be collapsed') + end + + %Start with changing alle the fields from the 3d mesh + + %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 ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; + if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; + if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; + if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; + if ~isnan(md.surfaceforcings.mass_balance), + md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); + end; + if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; + + %results + if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; + if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; + if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; + if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; + if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; + + %bedinfo and surface info + md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); + md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); + md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); + md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); + + %elementstype + if ~isnan(md.flowequation.element_equation) + md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); + md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); + md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); + md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); + md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); + end + + %boundary conditions + md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); + md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); + md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); + md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); + md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); + md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); + + %Extrusion of Neumann BC + if ~isnan(md.diagnostic.icefront), + numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); + md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer + end + + %materials + md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); + md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); + + %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.bed=project2d(md,md.geometry.bed,1); + md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); + md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); + md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); + md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); + md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); + md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); + md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); + md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); + md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); + + %lat long + if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end + if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end + + %Initialize with the 2d mesh + 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; + + %Keep a trace of lower and upper nodes + md.mesh.lowervertex=NaN; + md.mesh.uppervertex=NaN; + md.mesh.lowerelements=NaN; + md.mesh.upperelements=NaN; + + %Remove old mesh + md.mesh.x2d=NaN; + md.mesh.y2d=NaN; + md.mesh.elements2d=NaN; + md.mesh.numberofelements2d=md.mesh.numberofelements; + md.mesh.numberofvertices2d=md.mesh.numberofvertices; + md.mesh.numberoflayers=0; + + %Update mesh type + md.mesh.dimension=2; + end % }}} + function md2 = extract(md,area) % {{{ + %extract - extract a model according to an Argus contour or flag list + % + % 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 + % extract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); + % an empty string '' will be considered as an empty domain + % a string 'all' will be considered as the entire domain + % add an argument 0 if you do not want the elements to be checked (faster) + % + % Usage: + % md2=extract(md,area); + % + % Examples: + % md2=extract(md,'Domain.exp'); + % md2=extract(md,md.mask.elementonfloatingice); + % + % See also: EXTRUDE, COLLAPSE + + %copy model + md1=md; + + %some checks + if ((nargin~=2) | (nargout~=1)), + help extract + error('extract error message: bad usage'); + end + + %get check option + if (nargin==3 & varargin{1}==0), + checkoutline=0; + else + checkoutline=1; + end + + %get elements that are inside area + flag_elem=FlagElements(md1,area); + if ~any(flag_elem), + error('extracted model is empty'); + end + + %kick out all elements with 3 dirichlets + spc_elem=find(~flag_elem); + spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); + flag=ones(md1.mesh.numberofvertices,1); + flag(spc_node)=0; + pos=find(sum(flag(md1.mesh.elements),2)==0); + flag_elem(pos)=0; + + %extracted elements and nodes lists + pos_elem=find(flag_elem); + pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); + + %keep track of some fields + numberofvertices1=md1.mesh.numberofvertices; + numberofelements1=md1.mesh.numberofelements; + numberofvertices2=length(pos_node); + numberofelements2=length(pos_elem); + flag_node=zeros(numberofvertices1,1); + flag_node(pos_node)=1; + + %Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) + Pelem=zeros(numberofelements1,1); + Pelem(pos_elem)=[1:numberofelements2]'; + Pnode=zeros(numberofvertices1,1); + Pnode(pos_node)=[1:numberofvertices2]'; + + %renumber the elements (some node won't exist anymore) + elements_1=md1.mesh.elements; + elements_2=elements_1(pos_elem,:); + elements_2(:,1)=Pnode(elements_2(:,1)); + elements_2(:,2)=Pnode(elements_2(:,2)); + elements_2(:,3)=Pnode(elements_2(:,3)); + if md1.mesh.dimension==3, + elements_2(:,4)=Pnode(elements_2(:,4)); + elements_2(:,5)=Pnode(elements_2(:,5)); + elements_2(:,6)=Pnode(elements_2(:,6)); + end + + %OK, now create the new model ! + + %take every fields from model + md2=md1; + + %automatically modify fields + + %loop over model fields + model_fields=fields(md1); + for i=1:length(model_fields), + %get field + field=md1.(model_fields{i}); + fieldsize=size(field); + if isobject(field), %recursive call + object_fields=fields(md1.(model_fields{i})); + for j=1:length(object_fields), + %get field + field=md1.(model_fields{i}).(object_fields{j}); + fieldsize=size(field); + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); + end + end + else + %size = number of nodes * n + if fieldsize(1)==numberofvertices1 + md2.(model_fields{i})=field(pos_node,:); + elseif (fieldsize(1)==numberofvertices1+1) + md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; + %size = number of elements * n + elseif fieldsize(1)==numberofelements1 + md2.(model_fields{i})=field(pos_elem,:); + end + end + end + + %modify some specific fields + + %Mesh + md2.mesh.numberofelements=numberofelements2; + md2.mesh.numberofvertices=numberofvertices2; + md2.mesh.elements=elements_2; + + %mesh.uppervertex mesh.lowervertex + if md1.mesh.dimension==3 + md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); + pos=find(~isnan(md2.mesh.uppervertex)); + md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); + + md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); + pos=find(~isnan(md2.mesh.lowervertex)); + md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); + + md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); + pos=find(~isnan(md2.mesh.upperelements)); + md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); + + md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); + pos=find(~isnan(md2.mesh.lowerelements)); + md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); + end + + %Initial 2d mesh + if md1.mesh.dimension==3 + flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); + pos_elem_2d=find(flag_elem_2d); + flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); + pos_node_2d=find(flag_node_2d); + + md2.mesh.numberofelements2d=length(pos_elem_2d); + md2.mesh.numberofvertices2d=length(pos_node_2d); + md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); + md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); + md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); + md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); + + md2.mesh.x2d=md1.mesh.x(pos_node_2d); + md2.mesh.y2d=md1.mesh.y(pos_node_2d); + end + + %Edges + if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... + %renumber first two columns + pos=find(md2.mesh.edges(:,4)~=-1); + md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); + md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); + md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); + md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); + %remove edges when the 2 vertices are not in the domain. + md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); + %Replace all zeros by -1 in the last two columns; + pos=find(md2.mesh.edges(:,3)==0); + md2.mesh.edges(pos,3)=-1; + pos=find(md2.mesh.edges(:,4)==0); + md2.mesh.edges(pos,4)=-1; + %Invert -1 on the third column with last column (Also invert first two columns!!) + pos=find(md2.mesh.edges(:,3)==-1); + md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); + md2.mesh.edges(pos,4)=-1; + values=md2.mesh.edges(pos,2); + md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); + md2.mesh.edges(pos,1)=values; + %Finally remove edges that do not belong to any element + pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); + md2.mesh.edges(pos,:)=[]; + end + + %Penalties + if ~isnan(md2.diagnostic.vertex_pairing), + for i=1:size(md1.diagnostic.vertex_pairing,1); + md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); + end + md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); + end + if ~isnan(md2.prognostic.vertex_pairing), + for i=1:size(md1.prognostic.vertex_pairing,1); + md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); + end + md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); + end + + %recreate segments + if md1.mesh.dimension==2 + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + else + %First do the connectivity for the contourenvelope in 2d + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); + md2.mesh.segments=contourenvelope(md2); + md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; + md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); + %Then do it for 3d as usual + md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); + md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); + end + + %Boundary conditions: Dirichlets on new boundary + %Catch the elements that have not been extracted + orphans_elem=find(~flag_elem); + orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; + %Figure out which node are on the boundary between md2 and md1 + nodestoflag1=intersect(orphans_node,pos_node); + nodestoflag2=Pnode(nodestoflag1); + if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, + if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 + md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); + md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); + else + md2.diagnostic.spcvx(nodestoflag2)=NaN; + md2.diagnostic.spcvy(nodestoflag2)=NaN; + disp(' ') + disp('!! extract warning: spc values should be checked !!') + disp(' ') + end + %put 0 for vz + md2.diagnostic.spcvz(nodestoflag2)=0; + end + if ~isnan(md1.thermal.spctemperature), + md2.thermal.spctemperature(nodestoflag2,1)=1; + end + + %Diagnostic + if ~isnan(md2.diagnostic.icefront) + md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); + md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); + md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); + if md1.mesh.dimension==3 + md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); + md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); + end + md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); + end + + %Results fields + if isstruct(md1.results), + md2.results=struct(); + solutionfields=fields(md1.results); + for i=1:length(solutionfields), + %get subfields + solutionsubfields=fields(md1.results.(solutionfields{i})); + for j=1:length(solutionsubfields), + field=md1.results.(solutionfields{i}).(solutionsubfields{j}); + if length(field)==numberofvertices1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); + elseif length(field)==numberofelements1, + md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); + else + md2.results.(solutionfields{i}).(solutionsubfields{j})=field; + end + end + end + end + + %Keep track of pos_node and pos_elem + md2.mesh.extractedvertices=pos_node; + md2.mesh.extractedelements=pos_elem; + end % }}} + function md = extrude(md,varargin) % {{{ + %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 ((nargin>4) | (nargin<2) | (nargout~=1)), + help extrude; + error('extrude error message'); + end + + %Extrude the mesh + if nargin==2, %list of coefficients + list=varargin{1}; + if any(list<0) | any(list>1), + error('extrusioncoefficients must be between 0 and 1'); + end + extrusionlist=sort(unique([list(:);0;1])); + numlayers=length(extrusionlist); + elseif nargin==3, %one polynomial law + if varargin{2}<=0, + help extrude; + error('extrusionexponent must be >=0'); + end + numlayers=varargin{1}; + extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; + elseif nargin==4, %two polynomial laws + numlayers=varargin{1}; + lowerexp=varargin{2}; + upperexp=varargin{3}; + + if varargin{2}<=0 | varargin{3}<=0, + help extrude; + error('lower and upper extrusionexponents must be >=0'); + end + + lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; + upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; + extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); + + end + + if numlayers<2, + error('number of layers should be at least 2'); + end + if md.mesh.dimension==3, + error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); + end + + %Initialize with the 2d mesh + x3d=[]; + y3d=[]; + z3d=[]; %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=1:numlayers, + x3d=[x3d; md.mesh.x]; + y3d=[y3d; md.mesh.y]; + %nodes are distributed between bed and surface accordingly to the given exponent + z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; + end + number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh + + %Extrude elements + elements3d=[]; + for i=1:numlayers-1, + elements3d=[elements3d;[md.mesh.elements+(i-1)*md.mesh.numberofvertices md.mesh.elements+i*md.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part + end + number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh + + %Keep a trace of lower and upper nodes + mesh.lowervertex=NaN*ones(number_nodes3d,1); + mesh.uppervertex=NaN*ones(number_nodes3d,1); + mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; + mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; + md.mesh.lowervertex=mesh.lowervertex; + md.mesh.uppervertex=mesh.uppervertex; + + %same for lower and upper elements + mesh.lowerelements=NaN*ones(number_el3d,1); + mesh.upperelements=NaN*ones(number_el3d,1); + mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; + mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; + 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 ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; + if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; + if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; + if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; + if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; + if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; + + %bedinfo and surface info + md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); + md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); + md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); + md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); + + %elementstype + if ~isnan(md.flowequation.element_equation) + oldelements_type=md.flowequation.element_equation; + md.flowequation.element_equation=zeros(number_el3d,1); + md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); + end + + %verticestype + if ~isnan(md.flowequation.vertex_equation) + oldvertices_type=md.flowequation.vertex_equation; + md.flowequation.vertex_equation=zeros(number_nodes3d,1); + md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); + end + 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',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=[md.diagnostic.icefront(:,1:2) md.diagnostic.icefront(:,2)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,1)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,3:4)]; %Add two columns on the first layer + pressureload=[]; + for i=1:numlayers-1, + pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; + end + md.diagnostic.icefront=pressureload; + + %connectivity + md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); + md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; + for i=2:numlayers-1, + md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... + =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; + end + md.mesh.elementconnectivity(find(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 ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; + if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; + if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; + if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end + if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end + + %Put lithostatic pressure if there is an existing pressure + if ~isnan(md.initialization.pressure), + md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); + end + + %special for thermal modeling: + md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); + if ~isnan(md.basalforcings.geothermalflux) + md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux + end + + %increase connectivity if less than 25: + if md.mesh.average_vertex_connectivity<=25, + md.mesh.average_vertex_connectivity=100; + end + end % }}} function md = structtomodel(md,structmd) % {{{ Index: /issm/trunk-jpl/src/m/interp/averaging.m =================================================================== --- /issm/trunk-jpl/src/m/interp/averaging.m (revision 13005) +++ /issm/trunk-jpl/src/m/interp/averaging.m (revision 13005) @@ -0,0 +1,90 @@ +function average=averaging(md,data,iterations,varargin) +%AVERAGING - smooths the input over the mesh +% +% This routine takes a list over the elements or the nodes in input +% and return a list over the nodes. +% For each iterations it computes the average over each element (average +% of the vertices values) and then computes the average over each node +% by taking the average of the element around a node weighted by the +% elements volume +% For 3d mesh, a last argument can be added to specify the layer to be averaged on. +% +% Usage: +% smoothdata=averaging(md,data,iterations) +% smoothdata=averaging(md,data,iterations,layer) +% +% Examples: +% velsmoothed=averaging(md,md.initialization.vel,4); +% pressure=averaging(md,md.initialization.pressure,0); +% temperature=averaging(md,md.initialization.temperature,1,1); + +if ((nargin~=4) & (nargin~=3)), + error('averaging error message'); +end +if (length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices), + error('averaging error message: data not supported yet'); +end +if md.mesh.dimension==3 & nargin==4, + if varargin{1}<=0 | varargin{1}>md.mesh.numberoflayers, + error('layer should be between 1 and md.mesh.numberoflayers'); + end + layer=varargin{1}; +else + layer=0; +end + +%initialization +if layer==0, + weights=zeros(md.mesh.numberofvertices,1); + data=data(:); +else + weights=zeros(md.mesh.numberofvertices2d,1); + data=data((layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:); +end + +%load some variables (it is much faster if the variabes are loaded from md once for all) +if layer==0, + index=md.mesh.elements; + numberofnodes=md.mesh.numberofvertices; + numberofelements=md.mesh.numberofelements; +else + index=md.mesh.elements2d; + numberofnodes=md.mesh.numberofvertices2d; + numberofelements=md.mesh.numberofelements2d; +end + +%build some variables +line=index(:); +if md.mesh.dimension==3 & layer==0, + rep=6; + areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z); +elseif md.mesh.dimension==2, + rep=3; + areas=GetAreas(index,md.mesh.x,md.mesh.y); +else + rep=3; + areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d); +end +summation=1/rep*ones(rep,1); +linesize=rep*numberofelements; + +%update weights that holds the volume of all the element holding the node i +weights=sparse(line,ones(linesize,1),repmat(areas,rep,1),numberofnodes,1); + +%initialization +if length(data)==numberofelements + average_node=sparse(line,ones(linesize,1),repmat(areas.*data,rep,1),numberofnodes,1); + average_node=average_node./weights; +else + average_node=data; +end + +%loop over iteration +for i=1:iterations + average_el=average_node(index)*summation; + average_node=sparse(line,ones(linesize,1),repmat(areas.*average_el,rep,1),numberofnodes,1); + average_node=average_node./weights; +end + +%return output as a full matrix (C code do not like sparse matrices) +average=full(average_node); Index: /issm/trunk-jpl/src/m/meca/basalstress.m =================================================================== --- /issm/trunk-jpl/src/m/meca/basalstress.m (revision 13005) +++ /issm/trunk-jpl/src/m/meca/basalstress.m (revision 13005) @@ -0,0 +1,22 @@ +function [bx by b]=basalstress(md) +%BASALSTRESS - compute basal stress from basal drag and geometric information. +% +% Usage: +% [bx by b]=basalstress(md); +% +% See also: plot_basaldrag + + +%compute exponents +s=averaging(md,1./md.friction.p,0); +r=averaging(md,md.friction.q./md.friction.p,0); + +%compute horizontal velocity +ub=sqrt(md.initialization.vx.^2+md.initialization.vy.^2)/md.constants.yts; +ubx=md.initialization.vx/md.constants.yts; +uby=md.initialization.vy/md.constants.yts; + +%compute basal drag +bx=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*ubx.^s; +by=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*uby.^s; +b=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*ub.^s; Index: /issm/trunk-jpl/src/m/meca/mechanicalproperties.m =================================================================== --- /issm/trunk-jpl/src/m/meca/mechanicalproperties.m (revision 13005) +++ /issm/trunk-jpl/src/m/meca/mechanicalproperties.m (revision 13005) @@ -0,0 +1,118 @@ +function md=mechanicalproperties(md,vx,vy) +%MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity +% +% this routine computes the components of the stress tensor +% strain rate tensor and their respective principal directions. +% the results are in the model md: md.results +% +% Usage: +% md=mechanicalproperties(md,vx,vy) +% +% Example: +% md=mechanicalproperties(md,md.initialization.vx,md.initialization.vy); +% md=mechanicalproperties(md,md.inversion.vx_obs,md.inversion.vy_obs); + +%some checks +if length(vx)~=md.mesh.numberofvertices | length(vy)~=md.mesh.numberofvertices, + error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!']) +end +if ~(md.mesh.dimension==2) + error('only 2d model supported yet'); +end +if any(md.flowequation.element_equation~=2), + disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements'); +end + +%initialization +numberofelements=md.mesh.numberofelements; +index=md.mesh.elements; +summation=[1;1;1]; +directionsstress=zeros(numberofelements,4); +directionsstrain=zeros(numberofelements,4); +valuesstress=zeros(numberofelements,2); +valuesstrain=zeros(numberofelements,2); + +%compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma +[alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y); + +%compute shear +vxlist=vx(index)/md.constants.yts; +vylist=vy(index)/md.constants.yts; +ux=(vxlist.*alpha)*summation; +uy=(vxlist.*beta)*summation; +vx=(vylist.*alpha)*summation; +vy=(vylist.*beta)*summation; +uyvx=(vx+uy)./2; +clear vxlist vylist + +%compute viscosity +nu=zeros(numberofelements,1); +B_bar=md.materials.rheology_B(index)*summation/3; +power=(md.materials.rheology_n-1)./(2*md.materials.rheology_n); +second_inv=(ux.^2+vy.^2+((uy+vx).^2)/4+ux.*vy); +%some corrections +location=find(second_inv~=0); +nu(location)=B_bar(location)./(second_inv(location).^power(location)); +location=find(second_inv==0 & power~=0); +nu(location)=10^18; %arbitrary maximum viscosity to apply where there is no effective shear +location=find(second_inv==0 & power==0); +nu(location)=B_bar(location); +clear B_bar location second_inv power + +%compute stress +tau_xx=nu.*ux; +tau_yy=nu.*vy; +tau_xy=nu.*uyvx; + +%compute principal properties of stress +for i=1:numberofelements, + + %compute stress and strainrate matrices + stress=[tau_xx(i) tau_xy(i) + tau_xy(i) tau_yy(i)]; + strain=[ux(i) uyvx(i) + uyvx(i) vy(i)]; + + %eigen values and vectors + [directions,value]=eig(stress); + valuesstress(i,:)=[value(1,1) value(2,2)]; + directionsstress(i,:)=directions(:)'; + [directions,value]=eig(strain); + valuesstrain(i,:)=[value(1,1) value(2,2)]; + directionsstrain(i,:)=directions(:)'; +end + +%plug onto the model +%NB: Matlab sorts the eigen value in increasing order, we want the reverse +stress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); +stress.xx=tau_xx; +stress.yy=tau_yy; +stress.xy=tau_xy; +stress.principalvalue2=valuesstress(:,1); +stress.principalaxis2=directionsstress(:,1:2); +stress.principalvalue1=valuesstress(:,2); +stress.principalaxis1=directionsstress(:,3:4); +stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); +md.results.stress=stress; + +strainrate=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); +strainrate.xx=ux; +strainrate.yy=vy; +strainrate.xy=uyvx; +strainrate.principalvalue2=valuesstrain(:,1)*(365.25*24*3600); %strain rate in 1/a instead of 1/s +strainrate.principalaxis2=directionsstrain(:,1:2); +strainrate.principalvalue1=valuesstrain(:,2)*(365.25*24*3600); %strain rate in 1/a instead of 1/s +strainrate.principalaxis1=directionsstrain(:,3:4); +strainrate.effectivevalue=1/sqrt(2)*sqrt(strainrate.xx.^2+strainrate.yy.^2+2*strainrate.xy.^2); +md.results.strainrate=strainrate; + +deviatoricstress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); +deviatoricstress.xx=tau_xx; +deviatoricstress.yy=tau_yy; +deviatoricstress.xy=tau_xy; +deviatoricstress.principalvalue2=valuesstress(:,1); +deviatoricstress.principalaxis2=directionsstress(:,1:2); +deviatoricstress.principalvalue1=valuesstress(:,2); +deviatoricstress.principalaxis1=directionsstress(:,3:4); +deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); +md.results.deviatoricstress=deviatoricstress; Index: sm/trunk-jpl/src/m/model/AnalysisConfiguration.m =================================================================== --- /issm/trunk-jpl/src/m/model/AnalysisConfiguration.m (revision 13004) +++ (revision ) @@ -1,58 +1,0 @@ -function [analyses,numanalyses]=AnalysisConfiguration(solutiontype), -%ANALYSISCONFIGURATION - return type of analyses, number of analyses -% -% Usage: -% [analyses, numanalyses]=AnalysisConfiguration(solutiontype); - - - -switch solutiontype, - - case DiagnosticSolutionEnum, - numanalyses=5; - analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum]; - - case SteadystateSolutionEnum, - numanalyses=7; - analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum]; - - case ThermalSolutionEnum, - numanalyses=2; - analyses=[ThermalAnalysisEnum;MeltingAnalysisEnum]; - - case EnthalpySolutionEnum, - numanalyses=1; - analyses=[EnthalpyAnalysisEnum]; - - case PrognosticSolutionEnum, - numanalyses=1; - analyses=[PrognosticAnalysisEnum]; - - case BalancethicknessSolutionEnum, - numanalyses=1; - analyses=[BalancethicknessAnalysisEnum]; - - case SurfaceSlopeSolutionEnum, - numanalyses=1; - analyses=[SurfaceSlopeAnalysisEnum]; - - case BedSlopeSolutionEnum, - numanalyses=1; - analyses=[BedSlopeAnalysisEnum]; - - case TransientSolutionEnum, - numanalyses=9; - analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum]; - - case FlaimSolutionEnum, - numanalyses=1; - analyses=[FlaimAnalysisEnum]; - - case HydrologySolutionEnum, - numanalyses=3; - analyses=[BedSlopeAnalysisEnum;SurfaceSlopeAnalysisEnum;HydrologyAnalysisEnum]; - - otherwise - error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!'); - -end Index: sm/trunk-jpl/src/m/model/AnalysisConfiguration.py =================================================================== --- /issm/trunk-jpl/src/m/model/AnalysisConfiguration.py (revision 13004) +++ (revision ) @@ -1,59 +1,0 @@ -from EnumDefinitions import * - -def AnalysisConfiguration(solutiontype): - """ - ANALYSISCONFIGURATION - return type of analyses, number of analyses - - Usage: - [analyses, numanalyses]=AnalysisConfiguration(solutiontype); - """ - - if solutiontype == DiagnosticSolutionEnum: - numanalyses=5 - analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum] - - elif solutiontype == SteadystateSolutionEnum: - numanalyses=7 - analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum] - - elif solutiontype == ThermalSolutionEnum: - numanalyses=2 - analyses=[ThermalAnalysisEnum,MeltingAnalysisEnum] - - elif solutiontype == EnthalpySolutionEnum: - numanalyses=1 - analyses=[EnthalpyAnalysisEnum] - - elif solutiontype == PrognosticSolutionEnum: - numanalyses=1 - analyses=[PrognosticAnalysisEnum] - - elif solutiontype == BalancethicknessSolutionEnum: - numanalyses=1 - analyses=[BalancethicknessAnalysisEnum] - - elif solutiontype == SurfaceSlopeSolutionEnum: - numanalyses=1 - analyses=[SurfaceSlopeAnalysisEnum] - - elif solutiontype == BedSlopeSolutionEnum: - numanalyses=1 - analyses=[BedSlopeAnalysisEnum] - - elif solutiontype == TransientSolutionEnum: - numanalyses=9 - analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum,EnthalpyAnalysisEnum,PrognosticAnalysisEnum] - - elif solutiontype == FlaimSolutionEnum: - numanalyses=1 - analyses=[FlaimAnalysisEnum] - - elif solutiontype == HydrologySolutionEnum: - numanalyses=3 - analyses=[BedSlopeAnalysisEnum,SurfaceSlopeAnalysisEnum,HydrologyAnalysisEnum] - - else: - raise TypeError("solution type: '%s' not supported yet!" % EnumToString(solutiontype)) - - return analyses,numanalyses - Index: sm/trunk-jpl/src/m/model/BasinConstrain.m =================================================================== --- /issm/trunk-jpl/src/m/model/BasinConstrain.m (revision 13004) +++ (revision ) @@ -1,63 +1,0 @@ -function md=BasinConstrain(md,domain); -%BASINCONSTRAIN - constrain basin -% -% Constrain basin using a constraint domain outline, -% to dirichlet boundary conditions. -% constraindomain is an Argus domain outline file enclosing -% the geographical area of interest. -% -% Usage: -% md=BasinConstrain(md,constraindomain) -% -% Example: -% md=BasinConstrain(md,'DomainOutline.exp'); -% md=BasinConstrain(md,'~Iceshelves.exp'); - -%now, flag nodes and elements outside the domain outline. -if ischar(domain), - if isempty(domain), - elementondomain=zeros(md.mesh.numberofelements,1); - vertexondomain=zeros(md.mesh.numberofvertices,1); - invert=0; - elseif strcmpi(domain,'all') - elementondomain=ones(md.mesh.numberofelements,1); - vertexondomain=ones(md.mesh.numberofvertices,1); - invert=0; - else - %make sure that we actually don't want the elements outside the domain outline! - if strcmpi(domain(1),'~'), - domain=domain(2:end); - invert=1; - else - invert=0; - end - %ok, flag elements and nodes - [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); - end - if invert, - vertexondomain=~vertexondomain; - elementondomain=~elementondomain; - end -else - error('BasinConstrain error message: domain type not supported yet'); -end - -%list of elements and nodes not on domain -vertexnotondomain=find(~vertexondomain); -elementnotondomain=find(~elementondomain); - -%all elements outside the constraint domain are equivalent to water. all nodes outside are spc'd. -md.diagnostic.spcvx(vertexnotondomain)=md.inversion.vx_obs(vertexnotondomain); -md.diagnostic.spcvy(vertexnotondomain)=md.inversion.vy_obs(vertexnotondomain); -md.mask.elementonwater(elementnotondomain)=1; - -%now, make sure all elements on water have nodes that are spc'd, otherwise, we'll get a singular problem. -pos=find(~md.mask.elementonwater); -numpos=unique(md.mesh.elements(pos,:)); -nodes=setdiff(1:1:md.mesh.numberofvertices,numpos); -md.diagnostic.spcvx(nodes)=md.inversion.vx_obs(nodes); -md.diagnostic.spcvy(nodes)=md.inversion.vy_obs(nodes); - -%make sure icefronts that are completely spc'd are taken out: -free_segments=find((~isnan(md.diagnostic.spcvx(md.diagnostic.icefront(:,1:2))) + ~isnan(md.diagnostic.spcvy(md.diagnostic.icefront(:,1:2))))~=2); -md.diagnostic.icefront=md.diagnostic.icefront(free_segments,:); Index: sm/trunk-jpl/src/m/model/BasinConstrainShelf.m =================================================================== --- /issm/trunk-jpl/src/m/model/BasinConstrainShelf.m (revision 13004) +++ (revision ) @@ -1,74 +1,0 @@ -function md=BasinConstrainShelf(md,domain); -%BASINCONSTRAIN - constrain basin -% -% Constrain basin using a constraint domain outline, -% to dirichlet boundary conditions. -% constraindomain is an Argus domain outline file enclosing -% the geographical area of interest. -% -% Usage: -% md=BasinConstrain(md,constraindomain) -% -% Example: -% md=BasinConstrain(md,'DomainOutline.exp'); -% md=BasinConstrain(md,'~Iceshelves.exp'); - -%now, flag nodes and elements outside the domain outline. -if ischar(domain), - if isempty(domain), - elementondomain=zeros(md.mesh.numberofelements,1); - vertexondomain=zeros(md.mesh.numberofvertices,1); - invert=0; - elseif strcmpi(domain,'all') - elementondomain=ones(md.mesh.numberofelements,1); - vertexondomain=ones(md.mesh.numberofvertices,1); - invert=0; - else - %make sure that we actually don't want the elements outside the domain outline! - if strcmpi(domain(1),'~'), - domain=domain(2:end); - invert=1; - else - invert=0; - end - %ok, flag elements and nodes - [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); - end - if invert, - vertexondomain=~vertexondomain; - elementondomain=~elementondomain; - end -else - error('BasinConstrain error message: domain type not supported yet'); -end - -%list of elements and nodes not on domain -vertexnotondomain=find(~vertexondomain); -elementnotondomain=find(~elementondomain); - -%all elements outside the constraint domain are equivalent to water. all nodes outside are spc'd. -md.diagnostic.spcvx(vertexnotondomain)=md.inversion.vx_obs(vertexnotondomain); -md.diagnostic.spcvy(vertexnotondomain)=md.inversion.vy_obs(vertexnotondomain); -md.mask.elementonwater(elementnotondomain)=1; - -%now, make sure all elements on water have nodes that are spc'd, otherwise, we'll get a singular problem. -pos=find(~md.mask.elementonwater); -numpos=unique(md.mesh.elements(pos,:)); -nodes=setdiff(1:1:md.mesh.numberofvertices,numpos); -md.diagnostic.spcvx(nodes)=md.inversion.vx_obs(nodes); -md.diagnostic.spcvy(nodes)=md.inversion.vy_obs(nodes); - -%make sure any node with NaN velocity is spc'd: -%we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand. -pos=find(isnan(md.inversion.vel_obs_raw)); -md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos); -md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos); - -%iceshelves: any vertex on floating ice is spc'd -pos=find(md.mask.vertexongroundedice); -md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos); -md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos); - -%make sure icefronts that are completely spc'd are taken out: -free_segments=find((~isnan(md.diagnostic.spcvx(md.diagnostic.icefront(:,1:2))) + ~isnan(md.diagnostic.spcvy(md.diagnostic.icefront(:,1:2))) )~=2); -md.diagnostic.icefront=md.diagnostic.icefront(free_segments,:); Index: sm/trunk-jpl/src/m/model/addnote.m =================================================================== --- /issm/trunk-jpl/src/m/model/addnote.m (revision 13004) +++ (revision ) @@ -1,31 +1,0 @@ -function md=addnote(md,string) -%ADDNOTE - add a note to the existing model notes field -% -% Usage: -% md=addnote(md,string); -% -% Example: -% md=addnote(md,'Pine Island, Geometry of 2007'); - -if (nargin~=2) & (nargout~=1), - help addnote -end - -if ~ischar(string), - error('addnote error message: second input argument should be a string'); -end -notes=md.miscellaneous.notes; - -if ischar(notes), - newnotes=cell(2,1); - newnotes(1)={notes}; - newnotes(2)={string}; -else - newnotes=cell(length(notes)+1,1); - for i=1:length(notes), - newnotes(i)=notes(i); - end - newnotes(length(notes)+1)={string}; -end - -md.miscellaneous.notes=newnotes; Index: sm/trunk-jpl/src/m/model/addnote.py =================================================================== --- /issm/trunk-jpl/src/m/model/addnote.py (revision 13004) +++ (revision ) @@ -1,25 +1,0 @@ -def addnote(md,string): - """ - ADDNOTE - add a note to the existing model notes field - - Usage: - md=addnote(md,string); - - Example: - md=addnote(md,'Pine Island, Geometry of 2007'); - """ - - if not isinstance(string,str): - raise TypeError('addnote error message: second input argument should be a string') - - notes=md.miscellaneous.notes - - if isinstance(notes,str): - newnotes=[notes,string] - else: - newnotes=notes.append(string) - - md.miscellaneous.notes=newnotes - - return md - Index: sm/trunk-jpl/src/m/model/averageconnectivity.m =================================================================== --- /issm/trunk-jpl/src/m/model/averageconnectivity.m (revision 13004) +++ (revision ) @@ -1,11 +1,0 @@ -function conn=averageconnectivity(md) -%AVERAGECONNECTIVITY - computes the average connectivity of a model -% -% Usage: -% connectivity=averageconnectivity(md); - -nnz=0; -for i=1:md.mesh.numberofvertices, - nnz=nnz+length(find(md.mesh.elements==i)); -end -conn=nnz/md.mesh.numberofvertices; Index: sm/trunk-jpl/src/m/model/averaging.m =================================================================== --- /issm/trunk-jpl/src/m/model/averaging.m (revision 13004) +++ (revision ) @@ -1,90 +1,0 @@ -function average=averaging(md,data,iterations,varargin) -%AVERAGING - smooths the input over the mesh -% -% This routine takes a list over the elements or the nodes in input -% and return a list over the nodes. -% For each iterations it computes the average over each element (average -% of the vertices values) and then computes the average over each node -% by taking the average of the element around a node weighted by the -% elements volume -% For 3d mesh, a last argument can be added to specify the layer to be averaged on. -% -% Usage: -% smoothdata=averaging(md,data,iterations) -% smoothdata=averaging(md,data,iterations,layer) -% -% Examples: -% velsmoothed=averaging(md,md.initialization.vel,4); -% pressure=averaging(md,md.initialization.pressure,0); -% temperature=averaging(md,md.initialization.temperature,1,1); - -if ((nargin~=4) & (nargin~=3)), - error('averaging error message'); -end -if (length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices), - error('averaging error message: data not supported yet'); -end -if md.mesh.dimension==3 & nargin==4, - if varargin{1}<=0 | varargin{1}>md.mesh.numberoflayers, - error('layer should be between 1 and md.mesh.numberoflayers'); - end - layer=varargin{1}; -else - layer=0; -end - -%initialization -if layer==0, - weights=zeros(md.mesh.numberofvertices,1); - data=data(:); -else - weights=zeros(md.mesh.numberofvertices2d,1); - data=data((layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:); -end - -%load some variables (it is much faster if the variabes are loaded from md once for all) -if layer==0, - index=md.mesh.elements; - numberofnodes=md.mesh.numberofvertices; - numberofelements=md.mesh.numberofelements; -else - index=md.mesh.elements2d; - numberofnodes=md.mesh.numberofvertices2d; - numberofelements=md.mesh.numberofelements2d; -end - -%build some variables -line=index(:); -if md.mesh.dimension==3 & layer==0, - rep=6; - areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z); -elseif md.mesh.dimension==2, - rep=3; - areas=GetAreas(index,md.mesh.x,md.mesh.y); -else - rep=3; - areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d); -end -summation=1/rep*ones(rep,1); -linesize=rep*numberofelements; - -%update weights that holds the volume of all the element holding the node i -weights=sparse(line,ones(linesize,1),repmat(areas,rep,1),numberofnodes,1); - -%initialization -if length(data)==numberofelements - average_node=sparse(line,ones(linesize,1),repmat(areas.*data,rep,1),numberofnodes,1); - average_node=average_node./weights; -else - average_node=data; -end - -%loop over iteration -for i=1:iterations - average_el=average_node(index)*summation; - average_node=sparse(line,ones(linesize,1),repmat(areas.*average_el,rep,1),numberofnodes,1); - average_node=average_node./weights; -end - -%return output as a full matrix (C code do not like sparse matrices) -average=full(average_node); Index: sm/trunk-jpl/src/m/model/basalstress.m =================================================================== --- /issm/trunk-jpl/src/m/model/basalstress.m (revision 13004) +++ (revision ) @@ -1,22 +1,0 @@ -function [bx by b]=basalstress(md) -%BASALSTRESS - compute basal stress from basal drag and geometric information. -% -% Usage: -% [bx by b]=basalstress(md); -% -% See also: plot_basaldrag - - -%compute exponents -s=averaging(md,1./md.friction.p,0); -r=averaging(md,md.friction.q./md.friction.p,0); - -%compute horizontal velocity -ub=sqrt(md.initialization.vx.^2+md.initialization.vy.^2)/md.constants.yts; -ubx=md.initialization.vx/md.constants.yts; -uby=md.initialization.vy/md.constants.yts; - -%compute basal drag -bx=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*ubx.^s; -by=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*uby.^s; -b=(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)).^r.*(md.friction.coefficient).^2.*ub.^s; Index: sm/trunk-jpl/src/m/model/basevert.m =================================================================== --- /issm/trunk-jpl/src/m/model/basevert.m (revision 13004) +++ (revision ) @@ -1,35 +1,0 @@ -function wb=basevert(md) -%BASEVERT - computes the basal vertical velcities -% -% This routine computes the basal vertical velocities of ice shelves -% for 2d models only using the following formula: -% wb=rho_ice/rho_water*div(thickness*vel_horiz)+vel_horiz.grad(base) -% -% Usage: -% wb=basevert(md); - -alpha=zeros(md.mesh.numberofelements,3); -beta=zeros(md.mesh.numberofelements,3); -gamma=zeros(md.mesh.numberofelements,3); - -for n=1:md.mesh.numberofelements - X=inv([md.mesh.x(md.mesh.elements(n,:)) md.mesh.y(md.mesh.elements(n,:)) ones(3,1)]); - alpha(n,:)=X(1,:); - beta(n,:)=X(2,:); - gamma(n,:)=X(3,:); -end - -hu=md.geometry.thickness.*md.initialization.vx; -hv=md.geometry.thickness.*md.initialization.vy; - -summation=[1;1;1]; -hux=(hu(md.mesh.elements).*alpha)*summation; -hvy=(hv(md.mesh.elements).*beta)*summation; - -uelem=md.initialization.vx(md.mesh.elements)*summation/3; -velem=md.initialization.vy(md.mesh.elements)*summation/3; - -dbdx=(md.geometry.bed(md.mesh.elements).*alpha)*summation; -dbdy=(md.geometry.bed(md.mesh.elements).*beta)*summation; - -wb=-md.materials.rho_ice/md.materials.rho_water*(hux+hvy)+uelem.*dbdx+velem.*dbdy; Index: sm/trunk-jpl/src/m/model/bedslope.m =================================================================== --- /issm/trunk-jpl/src/m/model/bedslope.m (revision 13004) +++ (revision ) @@ -1,32 +1,0 @@ -function [bx,by,b]=bedslope(md) -%BEDSLOPE - compute the bed slope -% -% Usage: -% [bx,by,s]=bedslope(md) - -%load some variables (it is much faster if the variab;es are loaded from md once for all) -if (md.mesh.dimension==2), - numberofelements=md.mesh.numberofelements; - numberofnodes=md.mesh.numberofvertices; - index=md.mesh.elements; - x=md.mesh.x; y=md.mesh.y; -else - numberofelements=md.mesh.numberofelements2d; - numberofnodes=md.mesh.numberofvertices2d; - index=md.mesh.elements2d; - x=md.mesh.x2d; y=md.mesh.y2d; -end - -%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma -[alpha beta]=GetNodalFunctionsCoeff(index,x,y); - -summation=[1;1;1]; -bx=(md.geometry.bed(index).*alpha)*summation; -by=(md.geometry.bed(index).*beta)*summation; -b=sqrt(bx.^2+by.^2); - -if md.mesh.dimension==3, - bx=project3d(md,bx,'element'); - by=project3d(md,by,'element'); - b=sqrt(bx.^2+by.^2); -end Index: sm/trunk-jpl/src/m/model/collapse.m =================================================================== --- /issm/trunk-jpl/src/m/model/collapse.m (revision 13004) +++ (revision ) @@ -1,131 +1,0 @@ -function md=collapse(md) -%COLLAPSE - 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) -% -% See also: EXTRUDE, MODELEXTRACT - -%Check that the model is really a 3d model -if ~md.mesh.dimension==3, - error('collapse error message: only 3d mesh can be collapsed') -end - -%Start with changing alle the fields from the 3d mesh - -%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 ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end; -if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end; -if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; -if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; -if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; -if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; -if ~isnan(md.surfaceforcings.mass_balance), - md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); -end; -if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end; - -%results -if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end; -if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end; -if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end; -if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end; -if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end; - -%bedinfo and surface info -md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1); -md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1); -md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1); -md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1); - -%elementstype -if ~isnan(md.flowequation.element_equation) - md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); - md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); - md.flowequation.bordermacayeal=project2d(md,md.flowequation.bordermacayeal,1); - md.flowequation.borderpattyn=project2d(md,md.flowequation.borderpattyn,1); - md.flowequation.borderstokes=project2d(md,md.flowequation.borderstokes,1); -end - -%boundary conditions -md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers); -md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers); -md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers); -md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers); -md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); -md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); - -%Extrusion of Neumann BC -if ~isnan(md.diagnostic.icefront), - numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); - md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer -end - -%materials -md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); -md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); - -%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.bed=project2d(md,md.geometry.bed,1); -md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1); -md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); -md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); -md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1); -md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1); -md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1); -md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1); -md.mask.elementonwater=project2d(md,md.mask.elementonwater,1); -md.mask.vertexonwater=project2d(md,md.mask.vertexonwater,1); - -%lat long -if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end -if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end - -%Initialize with the 2d mesh -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; - -%Keep a trace of lower and upper nodes -md.mesh.lowervertex=NaN; -md.mesh.uppervertex=NaN; -md.mesh.lowerelements=NaN; -md.mesh.upperelements=NaN; - -%Remove old mesh -md.mesh.x2d=NaN; -md.mesh.y2d=NaN; -md.mesh.elements2d=NaN; -md.mesh.numberofelements2d=md.mesh.numberofelements; -md.mesh.numberofvertices2d=md.mesh.numberofvertices; -md.mesh.numberoflayers=0; - -%Update mesh type -md.mesh.dimension=2; Index: sm/trunk-jpl/src/m/model/extrude.m =================================================================== --- /issm/trunk-jpl/src/m/model/extrude.m (revision 13004) +++ (revision ) @@ -1,243 +1,0 @@ -function md=extrude(md,varargin) -%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 ((nargin>4) | (nargin<2) | (nargout~=1)), - help extrude; - error('extrude error message'); -end - -%Extrude the mesh -if nargin==2, %list of coefficients - list=varargin{1}; - if any(list<0) | any(list>1), - error('extrusioncoefficients must be between 0 and 1'); - end - extrusionlist=sort(unique([list(:);0;1])); - numlayers=length(extrusionlist); -elseif nargin==3, %one polynomial law - if varargin{2}<=0, - help extrude; - error('extrusionexponent must be >=0'); - end - numlayers=varargin{1}; - extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2}; -elseif nargin==4, %two polynomial laws - numlayers=varargin{1}; - lowerexp=varargin{2}; - upperexp=varargin{3}; - - if varargin{2}<=0 | varargin{3}<=0, - help extrude; - error('lower and upper extrusionexponents must be >=0'); - end - - lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; - upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; - extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); - -end - -if numlayers<2, - error('number of layers should be at least 2'); -end -if md.mesh.dimension==3, - error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); -end - -%Initialize with the 2d mesh -x3d=[]; -y3d=[]; -z3d=[]; %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=1:numlayers, - x3d=[x3d; md.mesh.x]; - y3d=[y3d; md.mesh.y]; - %nodes are distributed between bed and surface accordingly to the given exponent - z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; -end -number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh - -%Extrude elements -elements3d=[]; -for i=1:numlayers-1, - elements3d=[elements3d;[md.mesh.elements+(i-1)*md.mesh.numberofvertices md.mesh.elements+i*md.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part -end -number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh - -%Keep a trace of lower and upper nodes -mesh.lowervertex=NaN*ones(number_nodes3d,1); -mesh.uppervertex=NaN*ones(number_nodes3d,1); -mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; -mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; -md.mesh.lowervertex=mesh.lowervertex; -md.mesh.uppervertex=mesh.uppervertex; - -%same for lower and upper elements -mesh.lowerelements=NaN*ones(number_el3d,1); -mesh.upperelements=NaN*ones(number_el3d,1); -mesh.lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; -mesh.upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; -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 ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end; -if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end; -if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end; -if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end; -if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; -if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; - -%bedinfo and surface info -md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1); -md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1); -md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); -md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); - -%elementstype -if ~isnan(md.flowequation.element_equation) - oldelements_type=md.flowequation.element_equation; - md.flowequation.element_equation=zeros(number_el3d,1); - md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element'); -end - -%verticestype -if ~isnan(md.flowequation.vertex_equation) - oldvertices_type=md.flowequation.vertex_equation; - md.flowequation.vertex_equation=zeros(number_nodes3d,1); - md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node'); -end -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',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=[md.diagnostic.icefront(:,1:2) md.diagnostic.icefront(:,2)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,1)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,3:4)]; %Add two columns on the first layer -pressureload=[]; -for i=1:numlayers-1, - pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; -end -md.diagnostic.icefront=pressureload; - -%connectivity -md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); -md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; -for i=2:numlayers-1, - md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... - =md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; -end -md.mesh.elementconnectivity(find(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 ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end; -if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end; -if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end; -if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end -if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end -if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end -if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end - -%Put lithostatic pressure if there is an existing pressure -if ~isnan(md.initialization.pressure), - md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); -end - -%special for thermal modeling: -md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); -if ~isnan(md.basalforcings.geothermalflux) - md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux -end - -%increase connectivity if less than 25: -if md.mesh.average_vertex_connectivity<=25, - md.mesh.average_vertex_connectivity=100; -end Index: sm/trunk-jpl/src/m/model/mechanicalproperties.m =================================================================== --- /issm/trunk-jpl/src/m/model/mechanicalproperties.m (revision 13004) +++ (revision ) @@ -1,118 +1,0 @@ -function md=mechanicalproperties(md,vx,vy) -%MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity -% -% this routine computes the components of the stress tensor -% strain rate tensor and their respective principal directions. -% the results are in the model md: md.results -% -% Usage: -% md=mechanicalproperties(md,vx,vy) -% -% Example: -% md=mechanicalproperties(md,md.initialization.vx,md.initialization.vy); -% md=mechanicalproperties(md,md.inversion.vx_obs,md.inversion.vy_obs); - -%some checks -if length(vx)~=md.mesh.numberofvertices | length(vy)~=md.mesh.numberofvertices, - error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!']) -end -if ~(md.mesh.dimension==2) - error('only 2d model supported yet'); -end -if any(md.flowequation.element_equation~=2), - disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements'); -end - -%initialization -numberofelements=md.mesh.numberofelements; -index=md.mesh.elements; -summation=[1;1;1]; -directionsstress=zeros(numberofelements,4); -directionsstrain=zeros(numberofelements,4); -valuesstress=zeros(numberofelements,2); -valuesstrain=zeros(numberofelements,2); - -%compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma -[alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y); - -%compute shear -vxlist=vx(index)/md.constants.yts; -vylist=vy(index)/md.constants.yts; -ux=(vxlist.*alpha)*summation; -uy=(vxlist.*beta)*summation; -vx=(vylist.*alpha)*summation; -vy=(vylist.*beta)*summation; -uyvx=(vx+uy)./2; -clear vxlist vylist - -%compute viscosity -nu=zeros(numberofelements,1); -B_bar=md.materials.rheology_B(index)*summation/3; -power=(md.materials.rheology_n-1)./(2*md.materials.rheology_n); -second_inv=(ux.^2+vy.^2+((uy+vx).^2)/4+ux.*vy); -%some corrections -location=find(second_inv~=0); -nu(location)=B_bar(location)./(second_inv(location).^power(location)); -location=find(second_inv==0 & power~=0); -nu(location)=10^18; %arbitrary maximum viscosity to apply where there is no effective shear -location=find(second_inv==0 & power==0); -nu(location)=B_bar(location); -clear B_bar location second_inv power - -%compute stress -tau_xx=nu.*ux; -tau_yy=nu.*vy; -tau_xy=nu.*uyvx; - -%compute principal properties of stress -for i=1:numberofelements, - - %compute stress and strainrate matrices - stress=[tau_xx(i) tau_xy(i) - tau_xy(i) tau_yy(i)]; - strain=[ux(i) uyvx(i) - uyvx(i) vy(i)]; - - %eigen values and vectors - [directions,value]=eig(stress); - valuesstress(i,:)=[value(1,1) value(2,2)]; - directionsstress(i,:)=directions(:)'; - [directions,value]=eig(strain); - valuesstrain(i,:)=[value(1,1) value(2,2)]; - directionsstrain(i,:)=directions(:)'; -end - -%plug onto the model -%NB: Matlab sorts the eigen value in increasing order, we want the reverse -stress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); -stress.xx=tau_xx; -stress.yy=tau_yy; -stress.xy=tau_xy; -stress.principalvalue2=valuesstress(:,1); -stress.principalaxis2=directionsstress(:,1:2); -stress.principalvalue1=valuesstress(:,2); -stress.principalaxis1=directionsstress(:,3:4); -stress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); -md.results.stress=stress; - -strainrate=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); -strainrate.xx=ux; -strainrate.yy=vy; -strainrate.xy=uyvx; -strainrate.principalvalue2=valuesstrain(:,1)*(365.25*24*3600); %strain rate in 1/a instead of 1/s -strainrate.principalaxis2=directionsstrain(:,1:2); -strainrate.principalvalue1=valuesstrain(:,2)*(365.25*24*3600); %strain rate in 1/a instead of 1/s -strainrate.principalaxis1=directionsstrain(:,3:4); -strainrate.effectivevalue=1/sqrt(2)*sqrt(strainrate.xx.^2+strainrate.yy.^2+2*strainrate.xy.^2); -md.results.strainrate=strainrate; - -deviatoricstress=struct('xx',[],'yy',[],'xy',[],'principalvalue1',[],'principalaxis1',[],'principalvalue2',[],'principalaxis2',[],'effectivevalue',[]); -deviatoricstress.xx=tau_xx; -deviatoricstress.yy=tau_yy; -deviatoricstress.xy=tau_xy; -deviatoricstress.principalvalue2=valuesstress(:,1); -deviatoricstress.principalaxis2=directionsstress(:,1:2); -deviatoricstress.principalvalue1=valuesstress(:,2); -deviatoricstress.principalaxis1=directionsstress(:,3:4); -deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2); -md.results.deviatoricstress=deviatoricstress; Index: sm/trunk-jpl/src/m/model/modelextract.m =================================================================== --- /issm/trunk-jpl/src/m/model/modelextract.m (revision 13004) +++ (revision ) @@ -1,283 +1,0 @@ -function md2=modelextract(md1,area) -%modelextract - extract a model according to an Argus contour or flag list -% -% 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 -% modelextract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp'); -% an empty string '' will be considered as an empty domain -% a string 'all' will be considered as the entire domain -% add an argument 0 if you do not want the elements to be checked (faster) -% -% Usage: -% md2=modelextract(md1,area); -% -% Examples: -% md2=modelextract(md,'Domain.exp'); -% md2=modelextract(md,md.mask.elementonfloatingice); -% -% See also: EXTRUDE, COLLAPSE - -%some checks -if ((nargin~=2) | (nargout~=1)), - help modelextract - error('modelextract error message: bad usage'); -end - -%get check option -if (nargin==3 & varargin{1}==0), - checkoutline=0; -else - checkoutline=1; -end - -%get elements that are inside area -flag_elem=FlagElements(md1,area); -if ~any(flag_elem), - error('extracted model is empty'); -end - -%kick out all elements with 3 dirichlets -spc_elem=find(~flag_elem); -spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); -flag=ones(md1.mesh.numberofvertices,1); -flag(spc_node)=0; -pos=find(sum(flag(md1.mesh.elements),2)==0); -flag_elem(pos)=0; - -%extracted elements and nodes lists -pos_elem=find(flag_elem); -pos_node=sort(unique(md1.mesh.elements(pos_elem,:))); - -%keep track of some fields -numberofvertices1=md1.mesh.numberofvertices; -numberofelements1=md1.mesh.numberofelements; -numberofvertices2=length(pos_node); -numberofelements2=length(pos_elem); -flag_node=zeros(numberofvertices1,1); -flag_node(pos_node)=1; - -%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements) -Pelem=zeros(numberofelements1,1); -Pelem(pos_elem)=[1:numberofelements2]'; -Pnode=zeros(numberofvertices1,1); -Pnode(pos_node)=[1:numberofvertices2]'; - -%renumber the elements (some node won't exist anymore) -elements_1=md1.mesh.elements; -elements_2=elements_1(pos_elem,:); -elements_2(:,1)=Pnode(elements_2(:,1)); -elements_2(:,2)=Pnode(elements_2(:,2)); -elements_2(:,3)=Pnode(elements_2(:,3)); -if md1.mesh.dimension==3, - elements_2(:,4)=Pnode(elements_2(:,4)); - elements_2(:,5)=Pnode(elements_2(:,5)); - elements_2(:,6)=Pnode(elements_2(:,6)); -end - -%OK, now create the new model ! - - %take every fields from model - md2=md1; - - %automatically modify fields - - %loop over model fields - model_fields=fields(md1); - for i=1:length(model_fields), - %get field - field=md1.(model_fields{i}); - fieldsize=size(field); - if isobject(field), %recursive call - object_fields=fields(md1.(model_fields{i})); - for j=1:length(object_fields), - %get field - field=md1.(model_fields{i}).(object_fields{j}); - fieldsize=size(field); - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); - end - end - else - %size = number of nodes * n - if fieldsize(1)==numberofvertices1 - md2.(model_fields{i})=field(pos_node,:); - elseif (fieldsize(1)==numberofvertices1+1) - md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n - elseif fieldsize(1)==numberofelements1 - md2.(model_fields{i})=field(pos_elem,:); - end - end - end - - %modify some specific fields - - %Mesh - md2.mesh.numberofelements=numberofelements2; - md2.mesh.numberofvertices=numberofvertices2; - md2.mesh.elements=elements_2; - - %mesh.uppervertex mesh.lowervertex - if md1.mesh.dimension==3 - md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); - pos=find(~isnan(md2.mesh.uppervertex)); - md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos)); - - md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node); - pos=find(~isnan(md2.mesh.lowervertex)); - md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos)); - - md2.mesh.upperelements=md1.mesh.upperelements(pos_elem); - pos=find(~isnan(md2.mesh.upperelements)); - md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos)); - - md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem); - pos=find(~isnan(md2.mesh.lowerelements)); - md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos)); - end - - %Initial 2d mesh - if md1.mesh.dimension==3 - flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d); - pos_elem_2d=find(flag_elem_2d); - flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d); - pos_node_2d=find(flag_node_2d); - - md2.mesh.numberofelements2d=length(pos_elem_2d); - md2.mesh.numberofvertices2d=length(pos_node_2d); - md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:); - md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1)); - md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2)); - md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3)); - - md2.mesh.x2d=md1.mesh.x(pos_node_2d); - md2.mesh.y2d=md1.mesh.y(pos_node_2d); - end - - %Edges - if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs... - %renumber first two columns - pos=find(md2.mesh.edges(:,4)~=-1); - md2.mesh.edges(: ,1)=Pnode(md2.mesh.edges(:,1)); - md2.mesh.edges(: ,2)=Pnode(md2.mesh.edges(:,2)); - md2.mesh.edges(: ,3)=Pelem(md2.mesh.edges(:,3)); - md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4)); - %remove edges when the 2 vertices are not in the domain. - md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:); - %Replace all zeros by -1 in the last two columns; - pos=find(md2.mesh.edges(:,3)==0); - md2.mesh.edges(pos,3)=-1; - pos=find(md2.mesh.edges(:,4)==0); - md2.mesh.edges(pos,4)=-1; - %Invert -1 on the third column with last column (Also invert first two columns!!) - pos=find(md2.mesh.edges(:,3)==-1); - md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4); - md2.mesh.edges(pos,4)=-1; - values=md2.mesh.edges(pos,2); - md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1); - md2.mesh.edges(pos,1)=values; - %Finally remove edges that do not belong to any element - pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1); - md2.mesh.edges(pos,:)=[]; - end - - %Penalties - if ~isnan(md2.diagnostic.vertex_pairing), - for i=1:size(md1.diagnostic.vertex_pairing,1); - md2.diagnostic.vertex_pairing(i,:)=Pnode(md1.diagnostic.vertex_pairing(i,:)); - end - md2.diagnostic.vertex_pairing=md2.diagnostic.vertex_pairing(find(md2.diagnostic.vertex_pairing(:,1)),:); - end - if ~isnan(md2.prognostic.vertex_pairing), - for i=1:size(md1.prognostic.vertex_pairing,1); - md2.prognostic.vertex_pairing(i,:)=Pnode(md1.prognostic.vertex_pairing(i,:)); - end - md2.prognostic.vertex_pairing=md2.prognostic.vertex_pairing(find(md2.prognostic.vertex_pairing(:,1)),:); - end - - %recreate segments - if md1.mesh.dimension==2 - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - else - %First do the connectivity for the contourenvelope in 2d - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity); - md2.mesh.segments=contourenvelope(md2); - md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1; - md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1); - %Then do it for 3d as usual - md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices); - md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity); - end - - %Boundary conditions: Dirichlets on new boundary - %Catch the elements that have not been extracted - orphans_elem=find(~flag_elem); - orphans_node=unique(md1.mesh.elements(orphans_elem,:))'; - %Figure out which node are on the boundary between md2 and md1 - nodestoflag1=intersect(orphans_node,pos_node); - nodestoflag2=Pnode(nodestoflag1); - if numel(md1.diagnostic.spcvx)>1 & numel(md1.diagnostic.spcvy)>2 & numel(md1.diagnostic.spcvz)>2, - if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1 - md2.diagnostic.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); - md2.diagnostic.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2); - else - md2.diagnostic.spcvx(nodestoflag2)=NaN; - md2.diagnostic.spcvy(nodestoflag2)=NaN; - disp(' ') - disp('!! modelextract warning: spc values should be checked !!') - disp(' ') - end - %put 0 for vz - md2.diagnostic.spcvz(nodestoflag2)=0; - end - if ~isnan(md1.thermal.spctemperature), - md2.thermal.spctemperature(nodestoflag2,1)=1; - end - - %Diagnostic - if ~isnan(md2.diagnostic.icefront) - md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); - md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); - md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); - if md1.mesh.dimension==3 - md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); - md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); - end - md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); - end - - %Results fields - if isstruct(md1.results), - md2.results=struct(); - solutionfields=fields(md1.results); - for i=1:length(solutionfields), - %get subfields - solutionsubfields=fields(md1.results.(solutionfields{i})); - for j=1:length(solutionsubfields), - field=md1.results.(solutionfields{i}).(solutionsubfields{j}); - if length(field)==numberofvertices1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node); - elseif length(field)==numberofelements1, - md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem); - else - md2.results.(solutionfields{i}).(solutionsubfields{j})=field; - end - end - end - end - -%Keep track of pos_node and pos_elem -md2.mesh.extractedvertices=pos_node; -md2.mesh.extractedelements=pos_elem; Index: /issm/trunk-jpl/src/m/model/regional/BasinConstrain.m =================================================================== --- /issm/trunk-jpl/src/m/model/regional/BasinConstrain.m (revision 13005) +++ /issm/trunk-jpl/src/m/model/regional/BasinConstrain.m (revision 13005) @@ -0,0 +1,63 @@ +function md=BasinConstrain(md,domain); +%BASINCONSTRAIN - constrain basin +% +% Constrain basin using a constraint domain outline, +% to dirichlet boundary conditions. +% constraindomain is an Argus domain outline file enclosing +% the geographical area of interest. +% +% Usage: +% md=BasinConstrain(md,constraindomain) +% +% Example: +% md=BasinConstrain(md,'DomainOutline.exp'); +% md=BasinConstrain(md,'~Iceshelves.exp'); + +%now, flag nodes and elements outside the domain outline. +if ischar(domain), + if isempty(domain), + elementondomain=zeros(md.mesh.numberofelements,1); + vertexondomain=zeros(md.mesh.numberofvertices,1); + invert=0; + elseif strcmpi(domain,'all') + elementondomain=ones(md.mesh.numberofelements,1); + vertexondomain=ones(md.mesh.numberofvertices,1); + invert=0; + else + %make sure that we actually don't want the elements outside the domain outline! + if strcmpi(domain(1),'~'), + domain=domain(2:end); + invert=1; + else + invert=0; + end + %ok, flag elements and nodes + [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); + end + if invert, + vertexondomain=~vertexondomain; + elementondomain=~elementondomain; + end +else + error('BasinConstrain error message: domain type not supported yet'); +end + +%list of elements and nodes not on domain +vertexnotondomain=find(~vertexondomain); +elementnotondomain=find(~elementondomain); + +%all elements outside the constraint domain are equivalent to water. all nodes outside are spc'd. +md.diagnostic.spcvx(vertexnotondomain)=md.inversion.vx_obs(vertexnotondomain); +md.diagnostic.spcvy(vertexnotondomain)=md.inversion.vy_obs(vertexnotondomain); +md.mask.elementonwater(elementnotondomain)=1; + +%now, make sure all elements on water have nodes that are spc'd, otherwise, we'll get a singular problem. +pos=find(~md.mask.elementonwater); +numpos=unique(md.mesh.elements(pos,:)); +nodes=setdiff(1:1:md.mesh.numberofvertices,numpos); +md.diagnostic.spcvx(nodes)=md.inversion.vx_obs(nodes); +md.diagnostic.spcvy(nodes)=md.inversion.vy_obs(nodes); + +%make sure icefronts that are completely spc'd are taken out: +free_segments=find((~isnan(md.diagnostic.spcvx(md.diagnostic.icefront(:,1:2))) + ~isnan(md.diagnostic.spcvy(md.diagnostic.icefront(:,1:2))))~=2); +md.diagnostic.icefront=md.diagnostic.icefront(free_segments,:); Index: /issm/trunk-jpl/src/m/model/regional/BasinConstrainShelf.m =================================================================== --- /issm/trunk-jpl/src/m/model/regional/BasinConstrainShelf.m (revision 13005) +++ /issm/trunk-jpl/src/m/model/regional/BasinConstrainShelf.m (revision 13005) @@ -0,0 +1,74 @@ +function md=BasinConstrainShelf(md,domain); +%BASINCONSTRAIN - constrain basin +% +% Constrain basin using a constraint domain outline, +% to dirichlet boundary conditions. +% constraindomain is an Argus domain outline file enclosing +% the geographical area of interest. +% +% Usage: +% md=BasinConstrain(md,constraindomain) +% +% Example: +% md=BasinConstrain(md,'DomainOutline.exp'); +% md=BasinConstrain(md,'~Iceshelves.exp'); + +%now, flag nodes and elements outside the domain outline. +if ischar(domain), + if isempty(domain), + elementondomain=zeros(md.mesh.numberofelements,1); + vertexondomain=zeros(md.mesh.numberofvertices,1); + invert=0; + elseif strcmpi(domain,'all') + elementondomain=ones(md.mesh.numberofelements,1); + vertexondomain=ones(md.mesh.numberofvertices,1); + invert=0; + else + %make sure that we actually don't want the elements outside the domain outline! + if strcmpi(domain(1),'~'), + domain=domain(2:end); + invert=1; + else + invert=0; + end + %ok, flag elements and nodes + [vertexondomain elementondomain]=ContourToMesh(md.mesh.elements(:,1:3),md.mesh.x,md.mesh.y,domain,'element and node',2); + end + if invert, + vertexondomain=~vertexondomain; + elementondomain=~elementondomain; + end +else + error('BasinConstrain error message: domain type not supported yet'); +end + +%list of elements and nodes not on domain +vertexnotondomain=find(~vertexondomain); +elementnotondomain=find(~elementondomain); + +%all elements outside the constraint domain are equivalent to water. all nodes outside are spc'd. +md.diagnostic.spcvx(vertexnotondomain)=md.inversion.vx_obs(vertexnotondomain); +md.diagnostic.spcvy(vertexnotondomain)=md.inversion.vy_obs(vertexnotondomain); +md.mask.elementonwater(elementnotondomain)=1; + +%now, make sure all elements on water have nodes that are spc'd, otherwise, we'll get a singular problem. +pos=find(~md.mask.elementonwater); +numpos=unique(md.mesh.elements(pos,:)); +nodes=setdiff(1:1:md.mesh.numberofvertices,numpos); +md.diagnostic.spcvx(nodes)=md.inversion.vx_obs(nodes); +md.diagnostic.spcvy(nodes)=md.inversion.vy_obs(nodes); + +%make sure any node with NaN velocity is spc'd: +%we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand. +pos=find(isnan(md.inversion.vel_obs_raw)); +md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos); +md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos); + +%iceshelves: any vertex on floating ice is spc'd +pos=find(md.mask.vertexongroundedice); +md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos); +md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos); + +%make sure icefronts that are completely spc'd are taken out: +free_segments=find((~isnan(md.diagnostic.spcvx(md.diagnostic.icefront(:,1:2))) + ~isnan(md.diagnostic.spcvy(md.diagnostic.icefront(:,1:2))) )~=2); +md.diagnostic.icefront=md.diagnostic.icefront(free_segments,:); Index: /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.m (revision 13004) +++ /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.m (revision 13005) @@ -39,2 +39,62 @@ error('Model not consistent, see messages above'); end +end + +function [analyses,numanalyses]=AnalysisConfiguration(solutiontype), % {{{ +%ANALYSISCONFIGURATION - return type of analyses, number of analyses +% +% Usage: +% [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + + + +switch solutiontype, + + case DiagnosticSolutionEnum, + numanalyses=5; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum]; + + case SteadystateSolutionEnum, + numanalyses=7; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case ThermalSolutionEnum, + numanalyses=2; + analyses=[ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case EnthalpySolutionEnum, + numanalyses=1; + analyses=[EnthalpyAnalysisEnum]; + + case PrognosticSolutionEnum, + numanalyses=1; + analyses=[PrognosticAnalysisEnum]; + + case BalancethicknessSolutionEnum, + numanalyses=1; + analyses=[BalancethicknessAnalysisEnum]; + + case SurfaceSlopeSolutionEnum, + numanalyses=1; + analyses=[SurfaceSlopeAnalysisEnum]; + + case BedSlopeSolutionEnum, + numanalyses=1; + analyses=[BedSlopeAnalysisEnum]; + + case TransientSolutionEnum, + numanalyses=9; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum]; + + case FlaimSolutionEnum, + numanalyses=1; + analyses=[FlaimAnalysisEnum]; + + case HydrologySolutionEnum, + numanalyses=3; + analyses=[BedSlopeAnalysisEnum;SurfaceSlopeAnalysisEnum;HydrologyAnalysisEnum]; + + otherwise + error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!'); + +end % }}} Index: /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.py (revision 13004) +++ /issm/trunk-jpl/src/m/model/solve/ismodelselfconsistent.py (revision 13005) @@ -1,3 +1,62 @@ from AnalysisConfiguration import * +from EnumDefinitions import * + +def AnalysisConfiguration(solutiontype): #{{{ + """ + ANALYSISCONFIGURATION - return type of analyses, number of analyses + + Usage: + [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + """ + + if solutiontype == DiagnosticSolutionEnum: + numanalyses=5 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum] + + elif solutiontype == SteadystateSolutionEnum: + numanalyses=7 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == ThermalSolutionEnum: + numanalyses=2 + analyses=[ThermalAnalysisEnum,MeltingAnalysisEnum] + + elif solutiontype == EnthalpySolutionEnum: + numanalyses=1 + analyses=[EnthalpyAnalysisEnum] + + elif solutiontype == PrognosticSolutionEnum: + numanalyses=1 + analyses=[PrognosticAnalysisEnum] + + elif solutiontype == BalancethicknessSolutionEnum: + numanalyses=1 + analyses=[BalancethicknessAnalysisEnum] + + elif solutiontype == SurfaceSlopeSolutionEnum: + numanalyses=1 + analyses=[SurfaceSlopeAnalysisEnum] + + elif solutiontype == BedSlopeSolutionEnum: + numanalyses=1 + analyses=[BedSlopeAnalysisEnum] + + elif solutiontype == TransientSolutionEnum: + numanalyses=9 + analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticHutterAnalysisEnum,SurfaceSlopeAnalysisEnum,BedSlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum,EnthalpyAnalysisEnum,PrognosticAnalysisEnum] + + elif solutiontype == FlaimSolutionEnum: + numanalyses=1 + analyses=[FlaimAnalysisEnum] + + elif solutiontype == HydrologySolutionEnum: + numanalyses=3 + analyses=[BedSlopeAnalysisEnum,SurfaceSlopeAnalysisEnum,HydrologyAnalysisEnum] + + else: + raise TypeError("solution type: '%s' not supported yet!" % EnumToString(solutiontype)) + + return analyses,numanalyses +#}}} def ismodelselfconsistent(md):