Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 1754) +++ /issm/trunk/src/m/classes/@model/model.m (revision 1755) @@ -133,25 +133,10 @@ %Boundary conditions md.gridonboundary=NaN; - %Diagnostic - md.segmentonneumann_diag=NaN; - md.neumannvalues_diag=NaN; - md.gridondirichlet_diag=NaN; - md.dirichletvalues_diag=NaN; - - %Diagnostic stokes - md.segmentonneumann_diag_stokes=NaN; - - %Thermal + md.pressureload=NaN; + md.pressureload_stokes=NaN; + md.spcvelocity=NaN; + md.spctemperature=NaN; + md.spcthickness=NaN; md.min_thermal_constraints=0; - md.gridondirichlet_thermal=NaN; - md.dirichletvalues_thermal=NaN; - - %Transient - md.segmentonneumann_prog=NaN; - md.neumannvalues_prog=NaN; - md.segmentonneumann_prog2=NaN; - md.neumannvalues_prog2=NaN; - md.gridondirichlet_prog=NaN; - md.dirichletvalues_prog=NaN; %Observations Index: /issm/trunk/src/m/classes/public/BasinConstrain.m =================================================================== --- /issm/trunk/src/m/classes/public/BasinConstrain.m (revision 1754) +++ /issm/trunk/src/m/classes/public/BasinConstrain.m (revision 1755) @@ -48,7 +48,7 @@ %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. -md.gridondirichlet_diag(gridnotondomain)=1; -md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain); -md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,1:2)=1; +md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); md.elementonwater(elementnotondomain)=1; @@ -57,9 +57,9 @@ numpos=unique(md.elements(pos,:)); grids=setdiff(1:1:md.numberofgrids,numpos); -md.gridondirichlet_diag(grids)=1; -md.dirichletvalues_diag(grids,1)=md.vx_obs(grids); -md.dirichletvalues_diag(grids,2)=md.vy_obs(grids); +md.spcvelocity(grids,1:2)=1; +md.spcvelocity(grids,4)=md.vx_obs(grids); +md.spcvelocity(grids,5)=md.vy_obs(grids); %make sure icefronts that are completely spc'd are taken out: -free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2); -md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:); +free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); +md.pressureload=md.pressureload(free_segments,:); Index: /issm/trunk/src/m/classes/public/BasinConstrain2.m =================================================================== --- /issm/trunk/src/m/classes/public/BasinConstrain2.m (revision 1754) +++ /issm/trunk/src/m/classes/public/BasinConstrain2.m (revision 1755) @@ -48,7 +48,7 @@ %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. -md.gridondirichlet_diag(gridnotondomain)=1; -md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain); -md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,1:2)=1; +md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); md.elementonwater(elementnotondomain)=1; @@ -57,10 +57,10 @@ numpos=unique(md.elements(pos,:)); grids=setdiff(1:1:md.numberofgrids,numpos); -md.gridondirichlet_diag(grids)=1; -md.dirichletvalues_diag(grids,1)=md.vx_obs(grids); -md.dirichletvalues_diag(grids,2)=md.vy_obs(grids); +md.spcvelocity(grids,1:2)=1; +md.spcvelocity(grids,4)=md.vx_obs(grids); +md.spcvelocity(grids,5)=md.vy_obs(grids); %make sure icefronts that are completely spc'd are taken out: -free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2); -md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:); +free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); +md.pressureload=md.pressureload(free_segments,:); Index: /issm/trunk/src/m/classes/public/BasinConstrainShelf.m =================================================================== --- /issm/trunk/src/m/classes/public/BasinConstrainShelf.m (revision 1754) +++ /issm/trunk/src/m/classes/public/BasinConstrainShelf.m (revision 1755) @@ -48,7 +48,7 @@ %all elements outside the constraint domain are equivalent to water. all grids outside are spc'd. -md.gridondirichlet_diag(gridnotondomain)=1; -md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain); -md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,1:2)=1; +md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain); +md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain); md.elementonwater(elementnotondomain)=1; @@ -57,23 +57,23 @@ numpos=unique(md.elements(pos,:)); grids=setdiff(1:1:md.numberofgrids,numpos); -md.gridondirichlet_diag(grids)=1; -md.dirichletvalues_diag(grids,1)=md.vx_obs(grids); -md.dirichletvalues_diag(grids,2)=md.vy_obs(grids); +md.spcvelocity(grids,1:2)=1; +md.spcvelocity(grids,4)=md.vx_obs(grids); +md.spcvelocity(grids,5)=md.vy_obs(grids); %make sure any grid with NaN velocity is spc'd: pos=find(isnan(md.vel_obs_raw)); -md.gridondirichlet_diag(pos)=1; +md.spcvelocity(pos,1:2)=1; %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. -md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); -md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); +md.spcvelocity(pos,4)=md.vx_obs(pos); +md.spcvelocity(pos,5)=md.vy_obs(pos); %iceshelves: any grid on icesheet is spc'd pos=find(md.gridonicesheet); -md.gridondirichlet_diag(pos)=1; -md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); -md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); +md.spcvelocity(pos,1:2)=1; +md.spcvelocity(pos,4)=md.vx_obs(pos); +md.spcvelocity(pos,5)=md.vy_obs(pos); %make sure icefronts that are completely spc'd are taken out: -free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2); -md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:); +free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2); +md.pressureload=md.pressureload(free_segments,:); Index: /issm/trunk/src/m/classes/public/collapse.m =================================================================== --- /issm/trunk/src/m/classes/public/collapse.m (revision 1754) +++ /issm/trunk/src/m/classes/public/collapse.m (revision 1755) @@ -61,19 +61,12 @@ %boundary conditions -md.gridondirichlet_diag=project2d(md,md.gridondirichlet_diag,md.numlayers); -dirichletvalues_diag2d(:,1)=project2d(md,md.dirichletvalues_diag(:,1),md.numlayers); -dirichletvalues_diag2d(:,2)=project2d(md,md.dirichletvalues_diag(:,2),md.numlayers); -md.dirichletvalues_diag=dirichletvalues_diag2d; +md.spcvelocity=project2d(md,md.spcvelocity,md.numlayers); +md.spcthickness=project2d(md,md.spcthickness,md.numlayers); +md.spctemperature=project2d(md,md.spctemperature,md.numlayers); %Extrusion of Neumann BC %in 2d, segmentonnumann is: [grid1 grid2 element] -numberofneumann2d=size(md.segmentonneumann_diag,1)/md.numlayers; -md.segmentonneumann_diag=[md.segmentonneumann_diag(1:numberofneumann2d,1:2) md.segmentonneumann_diag(1:numberofneumann2d,5)]; %Add two columns on the first layer - -%Prognostic -md.gridondirichlet_prog=project2d(md,md.gridondirichlet_prog,md.numlayers); -md.dirichletvalues_prog=project2d(md,md.dirichletvalues_prog,md.numlayers); -%md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,5))]; -%md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,5))]; +numberofneumann2d=size(md.pressureload,1)/md.numlayers; +md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5)]; %Add two columns on the first layer %materials @@ -85,10 +78,4 @@ md.observed_temperature=DepthAverage(md,md.observed_temperature); md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux -md.gridondirichlet_thermal=project2d(md,md.gridondirichlet_thermal,md.numlayers); %surface temperature -md.dirichletvalues_thermal=project2d(md,md.dirichletvalues_thermal,md.numlayers); %surface temperature - -%NaN the values that are not on an spc'd temperature. -pos=find(~md.gridondirichlet_thermal); -md.dirichletvalues_thermal(pos)=NaN; %update of connectivity matrix Index: /issm/trunk/src/m/classes/public/display/displaybc.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaybc.m (revision 1754) +++ /issm/trunk/src/m/classes/public/display/displaybc.m (revision 1755) @@ -19,17 +19,12 @@ disp(sprintf('\n diagnostic:')); -fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]'); -fielddisplay(md,'segmentonneumann_diag','segments on ice front list'); -fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]'); +fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)'); +fielddisplay(md,'pressureload','segments on ice front list'); disp(sprintf('\n prognostic:')); -fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]'); -fielddisplay(md,'segmentonneumann_prog','segments on ice front list'); +fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)'); fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]'); disp(sprintf('\n thermal:')); -fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]'); +fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)'); fielddisplay(md,'melting','melting rate [m/a]'); Index: /issm/trunk/src/m/classes/public/display/displaydiagnostic.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 1754) +++ /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 1755) @@ -18,8 +18,6 @@ disp(sprintf('\n boundary conditions:')); -fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]'); -fielddisplay(md,'segmentonneumann_diag','segments on ice front list'); -fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]'); +fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)'); +fielddisplay(md,'pressureload','segments on ice front list'); disp(sprintf('\n %s','Penalties:')); Index: /issm/trunk/src/m/classes/public/display/displayprognostic.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displayprognostic.m (revision 1754) +++ /issm/trunk/src/m/classes/public/display/displayprognostic.m (revision 1755) @@ -17,6 +17,5 @@ disp(sprintf('\n boundary conditions:')); -fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]'); +fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)'); fielddisplay(md,'segmentonneumann_prog','segments on ice front list'); fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]'); Index: /issm/trunk/src/m/classes/public/display/displaythermal.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaythermal.m (revision 1754) +++ /issm/trunk/src/m/classes/public/display/displaythermal.m (revision 1755) @@ -12,6 +12,5 @@ disp(sprintf('\n boundary conditions:')); -fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list'); -fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]'); +fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)'); fielddisplay(md,'melting','melting rate [m/a]'); Index: /issm/trunk/src/m/classes/public/extrude.m =================================================================== --- /issm/trunk/src/m/classes/public/extrude.m (revision 1754) +++ /issm/trunk/src/m/classes/public/extrude.m (revision 1755) @@ -145,22 +145,19 @@ %boundary conditions -md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node'); -md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node'); +md.spcvelocity=project3d(md,md.spcvelovity,'node'); +md.spctemperature=project3d(md,md.spcvelovity,'node'); +md.spcthickness=project3d(md,md.spcvelovity,'node'); %Extrusion of Neumann BC %in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element] -oldsegmentonneumann_diag=md.segmentonneumann_diag; -segmentonneumann_diag_layer1=[oldsegmentonneumann_diag(:,1:2) oldsegmentonneumann_diag(:,2)+md.numberofgrids2d oldsegmentonneumann_diag(:,1)+md.numberofgrids2d oldsegmentonneumann_diag(:,3)]; %Add two columns on the first layer -segmentonneumann_diag=[]; +oldpressureload=md.pressureload; +pressureload_layer1=[oldpressureload(:,1:2) oldpressureload(:,2)+md.numberofgrids2d oldpressureload(:,1)+md.numberofgrids2d oldpressureload(:,3)]; %Add two columns on the first layer +pressureload=[]; for i=1:numlayers-1, - segmentonneumann_diag=[segmentonneumann_diag ;segmentonneumann_diag_layer1(:,1:4)+(i-1)*md.numberofgrids2d segmentonneumann_diag_layer1(:,5)+(i-1)*md.numberofelements2d ]; + pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.numberofgrids2d pressureload_layer1(:,5)+(i-1)*md.numberofelements2d ]; end %plug into md -md.segmentonneumann_diag=segmentonneumann_diag; -md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node'); -md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node'); -%md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))]; -%md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))]; +md.pressureload=pressureload; %materials @@ -184,14 +181,8 @@ md.observed_temperature=project3d(md,md.observed_temperature,'node'); md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux -md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers); %surface temperature -md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers); %surface temperature - -%NaN the values that are not on an spc'd temperature. -pos=find(~md.gridondirichlet_thermal); -md.dirichletvalues_thermal(pos)=NaN; %increase connectivity if less than 25: if md.connectivity<25, - md.connectivity=48; + md.connectivity=70; end Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 1754) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 1755) @@ -66,5 +66,5 @@ %NO NAN -fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... +fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',... 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',... 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; @@ -79,5 +79,5 @@ %FIELDS >= 0 -fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... +fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',... 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',... 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; @@ -96,5 +96,5 @@ %FIELDS ~=0 -fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',... +fields={'numberofelements','numberofgrids','elements','drag_type',... 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',... 'sparsity','deltaH','DeltaH','timeacc','timedec'}; @@ -118,5 +118,5 @@ %SIZE NUMBEROFGRIDS -fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'}; +fields={'x','y','z','B','drag','scpvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'}; for i=1:length(fields), if length(md.(fields{i}))~=md.numberofgrids, @@ -194,12 +194,6 @@ %SINGULAR - if ~any(md.gridondirichlet_diag), + if ~any(sum(md.spcvelocity(:,1:3),2)==3), disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) - bool=0;return; - end - - %DIRICHLET VALUES - if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag), - disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']); bool=0;return; end @@ -208,5 +202,5 @@ if any(md.thickness<=0), pos=find(md.thickness<=0); - if ~isempty(find(md.gridondirichlet_diag(pos)==0)), + if any(find(md.spcthickness(pos,1)==0)), disp(['model ' md.name ' has some grids with 0 thickness']); bool=0; return; @@ -217,10 +211,4 @@ %PROGNOSTIC if strcmp(md.analysis_type,'prognostic'), - - %DIRICHLET VALUES - if isempty(md.gridondirichlet_prog) | isempty(md.dirichletvalues_prog), - disp(['model ' md.name ' is not well posed. Missing dirichlet values for prognostic run']); - bool=0;return; - end %VELOCITIES @@ -339,12 +327,6 @@ %SINGULAR - if ~any(md.gridondirichlet_diag), + if ~any(sum(md.spcvelocity(:,1:3),2)==3), disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!']) - bool=0;return; - end - - %DIRICHLET VALUES - if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag), - disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']); bool=0;return; end @@ -353,5 +335,5 @@ if any(md.thickness<=0), pos=find(md.thickness<=0); - if ~isempty(find(md.gridondirichlet_diag(pos)==0)), + if any(find(md.spcthickness(pos,1)==0)), disp(['model ' md.name ' has some grids with 0 thickness']); bool=0; return; Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 1754) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 1755) @@ -77,19 +77,11 @@ WriteData(fid,md.gridoniceshelf,'Mat','gridoniceshelf'); -WriteData(fid,md.segmentonneumann_diag,'Mat','segmentonneumann_diag'); -WriteData(fid,md.segmentonneumann_diag_stokes,'Mat','segmentonneumann_diag_stokes'); -WriteData(fid,md.neumannvalues_diag,'Mat','neumannvalues_diag'); -WriteData(fid,md.gridondirichlet_diag,'Mat','gridondirichlet_diag'); -WriteData(fid,md.dirichletvalues_diag,'Mat','dirichletvalues_diag'); +WriteData(fid,md.spcvelocity,'Mat','spcvelocity'); +WriteData(fid,md.spctemperature,'Mat','spctemperature'); +WriteData(fid,md.spcthickness,'Mat','spcthickness'); -WriteData(fid,md.segmentonneumann_prog,'Mat','segmentonneumann_prog'); -WriteData(fid,md.neumannvalues_prog,'Mat','neumannvalues_prog'); -WriteData(fid,md.gridondirichlet_prog,'Mat','gridondirichlet_prog'); -WriteData(fid,md.dirichletvalues_prog,'Mat','dirichletvalues_prog'); -WriteData(fid,md.segmentonneumann_prog2,'Mat','segmentonneumann_prog2'); -WriteData(fid,md.neumannvalues_prog2,'Mat','neumannvalues_prog2'); +WriteData(fid,md.pressureload,'Mat','pressureload'); +WriteData(fid,md.pressureload_stokes,'Mat','pressureload_stokes'); -WriteData(fid,md.gridondirichlet_thermal,'Mat','gridondirichlet_thermal'); -WriteData(fid,md.dirichletvalues_thermal,'Mat','dirichletvalues_thermal'); WriteData(fid,md.geothermalflux,'Mat','geothermalflux'); WriteData(fid,md.melting,'Mat','melting'); Index: /issm/trunk/src/m/classes/public/modelextract.m =================================================================== --- /issm/trunk/src/m/classes/public/modelextract.m (revision 1754) +++ /issm/trunk/src/m/classes/public/modelextract.m (revision 1755) @@ -188,24 +188,24 @@ gridstoflag1=intersect(orphans_grid,pos_grid); gridstoflag2=Pgrid(gridstoflag1); - if ~isnan(md1.gridondirichlet_diag), - md2.gridondirichlet_diag(gridstoflag2)=1; + if ~isnan(md1.spcvelocity), + md2.spcvelocity(gridstoflag2,1:2)=1; if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs) - md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); - md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2); + md2.spcvelocity(gridstoflag2,4)=md2.vx_obs(gridstoflag2); + md2.spcvelocity(gridstoflag2,5)=md2.vy_obs(gridstoflag2); else - md2.dirichletvalues_diag=zeros(numberofgrids2,2); - disp(' ') - disp('!! modelextract warning: dirichlet values should be checked !!') - disp(' ') - end - end - if ~isnan(md1.gridondirichlet_thermal), - md2.gridondirichlet_thermal(gridstoflag2)=1; + md2.spcvelocity(gridstoflag2,4:5)=zeros(length(gridstoflag2)),2); + disp(' ') + disp('!! modelextract warning: spc values should be checked !!') + disp(' ') + end + end + if ~isnan(md1.spctemperature), + md2.spctemperature(gridstoflag2,1)=1; if ~isnan(md1.observed_temperature) - md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2); + md2.spctemperature(gridstoflag2,2)=md2.observed_temperature(gridstoflag2); else - md2.dirichletvalues_thermal=zeros(numberofgrids2,2); - disp(' ') - disp('!! modelextract warning: dirichlet values should be checked !!') + md2.spctemperature(gridstoflag2,2)=zeros(length(gridstoflag2),2); + disp(' ') + disp('!! modelextract warning: spc values should be checked !!') disp(' ') end @@ -214,5 +214,5 @@ %Recreate border stokes if md1.isstokes, - pos=find(md2.gridondirichlet_diag); %find all the grids on the boundary of the domain without icefront + pos=find(sum(md2.spcvelocity(:,1:3),2)); %find all the grids on the boundary of the domain without icefront md2.gridonstokes(pos)=0; %we have to constrain all the boundary of the domain without icefront stokes_elements=find(md2.elements_type(:,2)==StokesFormulationEnum()); %find the elements on the stokes domain @@ -223,39 +223,23 @@ %Diagnostic - if ~isnan(md2.segmentonneumann_diag) - md2.segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1)); - md2.segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2)); - md2.segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end)); + if ~isnan(md2.pressureload) + md2.pressureload(:,1)=Pgrid(md1.pressureload(:,1)); + md2.pressureload(:,2)=Pgrid(md1.pressureload(:,2)); + md2.pressureload(:,end)=Pelem(md1.pressureload(:,end)); if strcmpi(md1.type,'3d') - md2.segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3)); - md2.segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4)); - end - md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:); - end - if ~isnan(md2.segmentonneumann_diag_stokes) - md2.segmentonneumann_diag_stokes(:,1)=Pgrid(md1.segmentonneumann_diag_stokes(:,1)); - md2.segmentonneumann_diag_stokes(:,2)=Pgrid(md1.segmentonneumann_diag_stokes(:,2)); - md2.segmentonneumann_diag_stokes(:,end)=Pelem(md1.segmentonneumann_diag_stokes(:,end)); + md2.pressureload(:,3)=Pgrid(md1.pressureload(:,3)); + md2.pressureload(:,4)=Pgrid(md1.pressureload(:,4)); + end + md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2)),:); + end + if ~isnan(md2.pressureload_stokes) + md2.pressureload_stokes(:,1)=Pgrid(md1.pressureload_stokes(:,1)); + md2.pressureload_stokes(:,2)=Pgrid(md1.pressureload_stokes(:,2)); + md2.pressureload_stokes(:,end)=Pelem(md1.pressureload_stokes(:,end)); if strcmpi(md1.type,'3d') - md2.segmentonneumann_diag_stokes(:,3)=Pgrid(md1.segmentonneumann_diag_stokes(:,3)); - md2.segmentonneumann_diag_stokes(:,4)=Pgrid(md1.segmentonneumann_diag_stokes(:,4)); - end - md2.segmentonneumann_diag_stokes=md2.segmentonneumann_diag_stokes(find(md2.segmentonneumann_diag_stokes(:,1) & md2.segmentonneumann_diag_stokes(:,2)),:); - end - md2.neumannvalues_diag=NaN; - - %Transient - if ~isnan(md2.segmentonneumann_prog) - md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1)); - md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2)); - md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end)); - md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:); - end - md2.neumannvalues_prog=NaN; - if ~isnan(md2.segmentonneumann_prog2) - md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1)); - md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2)); - md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end)); - md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:); + md2.pressureload_stokes(:,3)=Pgrid(md1.pressureload_stokes(:,3)); + md2.pressureload_stokes(:,4)=Pgrid(md1.pressureload_stokes(:,4)); + end + md2.pressureload_stokes=md2.pressureload_stokes(find(md2.pressureload_stokes(:,1) & md2.pressureload_stokes(:,2)),:); end Index: /issm/trunk/src/m/classes/public/setelementstype.m =================================================================== --- /issm/trunk/src/m/classes/public/setelementstype.m (revision 1754) +++ /issm/trunk/src/m/classes/public/setelementstype.m (revision 1755) @@ -61,5 +61,5 @@ gridonstokes=ones(md.numberofgrids,1); gridonstokes(md.elements(nonstokesflag,:))=0; %non stokes grids -pos=find(md.gridondirichlet_diag); %find all the grids on the boundary of the domain without icefront +pos=find(sum(md.spcvelocity(:,1:3),2)==3); %find all the grids on the boundary of the domain without icefront gridonstokes(pos)=0; %we have to constrain all the boundary of the domain without icefront