Changeset 1755


Ignore:
Timestamp:
08/18/09 15:16:26 (15 years ago)
Author:
Mathieu Morlighem
Message:

Boundary conditions new names

Location:
issm/trunk/src/m/classes
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/classes/@model/model.m

    r1666 r1755  
    133133        %Boundary conditions
    134134        md.gridonboundary=NaN;
    135         %Diagnostic
    136         md.segmentonneumann_diag=NaN;
    137         md.neumannvalues_diag=NaN;
    138         md.gridondirichlet_diag=NaN;
    139         md.dirichletvalues_diag=NaN;
    140        
    141         %Diagnostic stokes
    142         md.segmentonneumann_diag_stokes=NaN;
    143 
    144         %Thermal
     135        md.pressureload=NaN;
     136        md.pressureload_stokes=NaN;
     137        md.spcvelocity=NaN;
     138        md.spctemperature=NaN;
     139        md.spcthickness=NaN;
    145140        md.min_thermal_constraints=0;
    146         md.gridondirichlet_thermal=NaN;
    147         md.dirichletvalues_thermal=NaN;
    148 
    149         %Transient
    150         md.segmentonneumann_prog=NaN;
    151         md.neumannvalues_prog=NaN;
    152         md.segmentonneumann_prog2=NaN;
    153         md.neumannvalues_prog2=NaN;
    154         md.gridondirichlet_prog=NaN;
    155         md.dirichletvalues_prog=NaN;
    156141
    157142        %Observations
  • issm/trunk/src/m/classes/public/BasinConstrain.m

    r1338 r1755  
    4848
    4949%all elements outside the constraint domain are equivalent to water. all grids outside are spc'd.
    50 md.gridondirichlet_diag(gridnotondomain)=1;
    51 md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);
    52 md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);
     50md.spcvelocity(gridnotondomain,1:2)=1;
     51md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain);
     52md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain);
    5353md.elementonwater(elementnotondomain)=1;
    5454
     
    5757numpos=unique(md.elements(pos,:));
    5858grids=setdiff(1:1:md.numberofgrids,numpos);
    59 md.gridondirichlet_diag(grids)=1;
    60 md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
    61 md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
     59md.spcvelocity(grids,1:2)=1;
     60md.spcvelocity(grids,4)=md.vx_obs(grids);
     61md.spcvelocity(grids,5)=md.vy_obs(grids);
    6262
    6363%make sure icefronts that are completely spc'd are taken out:
    64 free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
    65 md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);
     64free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2);
     65md.pressureload=md.pressureload(free_segments,:);
  • issm/trunk/src/m/classes/public/BasinConstrain2.m

    r1317 r1755  
    4848
    4949%all elements outside the constraint domain are equivalent to water. all grids outside are spc'd.
    50 md.gridondirichlet_diag(gridnotondomain)=1;
    51 md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);
    52 md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);
     50md.spcvelocity(gridnotondomain,1:2)=1;
     51md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain);
     52md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain);
    5353md.elementonwater(elementnotondomain)=1;
    5454
     
    5757numpos=unique(md.elements(pos,:));
    5858grids=setdiff(1:1:md.numberofgrids,numpos);
    59 md.gridondirichlet_diag(grids)=1;
    60 md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
    61 md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
     59md.spcvelocity(grids,1:2)=1;
     60md.spcvelocity(grids,4)=md.vx_obs(grids);
     61md.spcvelocity(grids,5)=md.vy_obs(grids);
    6262
    6363
    6464%make sure icefronts that are completely spc'd are taken out:
    65 free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
    66 md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);
     65free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2);
     66md.pressureload=md.pressureload(free_segments,:);
  • issm/trunk/src/m/classes/public/BasinConstrainShelf.m

    r1317 r1755  
    4848
    4949%all elements outside the constraint domain are equivalent to water. all grids outside are spc'd.
    50 md.gridondirichlet_diag(gridnotondomain)=1;
    51 md.dirichletvalues_diag(gridnotondomain,1)=md.vx_obs(gridnotondomain);
    52 md.dirichletvalues_diag(gridnotondomain,2)=md.vy_obs(gridnotondomain);
     50md.spcvelocity(gridnotondomain,1:2)=1;
     51md.spcvelocity(gridnotondomain,4)=md.vx_obs(gridnotondomain);
     52md.spcvelocity(gridnotondomain,5)=md.vy_obs(gridnotondomain);
    5353md.elementonwater(elementnotondomain)=1;
    5454
     
    5757numpos=unique(md.elements(pos,:));
    5858grids=setdiff(1:1:md.numberofgrids,numpos);
    59 md.gridondirichlet_diag(grids)=1;
    60 md.dirichletvalues_diag(grids,1)=md.vx_obs(grids);
    61 md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
     59md.spcvelocity(grids,1:2)=1;
     60md.spcvelocity(grids,4)=md.vx_obs(grids);
     61md.spcvelocity(grids,5)=md.vy_obs(grids);
    6262
    6363
    6464%make sure any grid with NaN velocity is spc'd:
    6565pos=find(isnan(md.vel_obs_raw));
    66 md.gridondirichlet_diag(pos)=1;
     66md.spcvelocity(pos,1:2)=1;
    6767%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.
    68 md.dirichletvalues_diag(pos,1)=md.vx_obs(pos);
    69 md.dirichletvalues_diag(pos,2)=md.vy_obs(pos);
     68md.spcvelocity(pos,4)=md.vx_obs(pos);
     69md.spcvelocity(pos,5)=md.vy_obs(pos);
    7070
    7171%iceshelves: any grid on icesheet is spc'd
    7272pos=find(md.gridonicesheet);
    73 md.gridondirichlet_diag(pos)=1;
    74 md.dirichletvalues_diag(pos,1)=md.vx_obs(pos);
    75 md.dirichletvalues_diag(pos,2)=md.vy_obs(pos);
     73md.spcvelocity(pos,1:2)=1;
     74md.spcvelocity(pos,4)=md.vx_obs(pos);
     75md.spcvelocity(pos,5)=md.vy_obs(pos);
    7676
    7777%make sure icefronts that are completely spc'd are taken out:
    78 free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);
    79 md.segmentonneumann_diag=md.segmentonneumann_diag(free_segments,:);
     78free_segments=find(sum(md.spcvelocity(md.segmentonneumann_diag(:,1:2),1:2),2)~=2);
     79md.pressureload=md.pressureload(free_segments,:);
  • issm/trunk/src/m/classes/public/collapse.m

    r1 r1755  
    6161
    6262%boundary conditions
    63 md.gridondirichlet_diag=project2d(md,md.gridondirichlet_diag,md.numlayers);
    64 dirichletvalues_diag2d(:,1)=project2d(md,md.dirichletvalues_diag(:,1),md.numlayers);
    65 dirichletvalues_diag2d(:,2)=project2d(md,md.dirichletvalues_diag(:,2),md.numlayers);
    66 md.dirichletvalues_diag=dirichletvalues_diag2d;
     63md.spcvelocity=project2d(md,md.spcvelocity,md.numlayers);
     64md.spcthickness=project2d(md,md.spcthickness,md.numlayers);
     65md.spctemperature=project2d(md,md.spctemperature,md.numlayers);
    6766
    6867%Extrusion of Neumann BC
    6968%in 2d, segmentonnumann is: [grid1 grid2 element]
    70 numberofneumann2d=size(md.segmentonneumann_diag,1)/md.numlayers;
    71 md.segmentonneumann_diag=[md.segmentonneumann_diag(1:numberofneumann2d,1:2) md.segmentonneumann_diag(1:numberofneumann2d,5)]; %Add two columns on the first layer
    72 
    73 %Prognostic
    74 md.gridondirichlet_prog=project2d(md,md.gridondirichlet_prog,md.numlayers);
    75 md.dirichletvalues_prog=project2d(md,md.dirichletvalues_prog,md.numlayers);
    76 %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,5))];
    77 %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,5))];
     69numberofneumann2d=size(md.pressureload,1)/md.numlayers;
     70md.pressureload=[md.pressureload(1:numberofneumann2d,1:2) md.pressureload(1:numberofneumann2d,5)]; %Add two columns on the first layer
    7871
    7972%materials
     
    8578md.observed_temperature=DepthAverage(md,md.observed_temperature);
    8679md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux
    87 md.gridondirichlet_thermal=project2d(md,md.gridondirichlet_thermal,md.numlayers); %surface temperature
    88 md.dirichletvalues_thermal=project2d(md,md.dirichletvalues_thermal,md.numlayers); %surface temperature
    89 
    90 %NaN the values that are not on an spc'd temperature.
    91 pos=find(~md.gridondirichlet_thermal);
    92 md.dirichletvalues_thermal(pos)=NaN;
    9380
    9481%update of connectivity matrix
  • issm/trunk/src/m/classes/public/display/displaybc.m

    r1277 r1755  
    1919
    2020disp(sprintf('\n      diagnostic:'));
    21 fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list');
    22 fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]');
    23 fielddisplay(md,'segmentonneumann_diag','segments on ice front list');
    24 fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]');
     21fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)');
     22fielddisplay(md,'pressureload','segments on ice front list');
    2523
    2624disp(sprintf('\n      prognostic:'));
    27 fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list');
    28 fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]');
    29 fielddisplay(md,'segmentonneumann_prog','segments on ice front list');
     25fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)');
    3026fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]');
    3127
    3228disp(sprintf('\n      thermal:'));
    33 fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list');
    34 fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]');
     29fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)');
    3530fielddisplay(md,'melting','melting rate [m/a]');
  • issm/trunk/src/m/classes/public/display/displaydiagnostic.m

    r1666 r1755  
    1818
    1919disp(sprintf('\n      boundary conditions:'));
    20 fielddisplay(md,'gridondirichlet_diag','grid on dirichlet flags list');
    21 fielddisplay(md,'dirichletvalues_diag','values of the dirichlet [m/a]');
    22 fielddisplay(md,'segmentonneumann_diag','segments on ice front list');
    23 fielddisplay(md,'neumannvalues_diag','values of the Neumann [N]');
     20fielddisplay(md,'spcvelocity','constraints flag list (first 3 columns) and values [m/yr] (last 3 columns)');
     21fielddisplay(md,'pressureload','segments on ice front list');
    2422
    2523disp(sprintf('\n      %s','Penalties:'));
  • issm/trunk/src/m/classes/public/display/displayprognostic.m

    r1277 r1755  
    1717
    1818disp(sprintf('\n      boundary conditions:'));
    19 fielddisplay(md,'gridondirichlet_prog','grid on dirichlet flags list');
    20 fielddisplay(md,'dirichletvalues_prog','values of the dirichlet [m]');
     19fielddisplay(md,'spcthickness','constraints flag list (first column) and values (second column)');
    2120fielddisplay(md,'segmentonneumann_prog','segments on ice front list');
    2221fielddisplay(md,'neumannvalues_prog','values of the Neumann [m/a]');
  • issm/trunk/src/m/classes/public/display/displaythermal.m

    r1277 r1755  
    1212
    1313disp(sprintf('\n      boundary conditions:'));
    14 fielddisplay(md,'gridondirichlet_thermal','grid on dirichlet flags list');
    15 fielddisplay(md,'dirichletvalues_thermal','values of the dirichlet [m]');
     14fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)');
    1615fielddisplay(md,'melting','melting rate [m/a]');
    1716
  • issm/trunk/src/m/classes/public/extrude.m

    r1651 r1755  
    145145
    146146%boundary conditions
    147 md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node');
    148 md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node');
     147md.spcvelocity=project3d(md,md.spcvelovity,'node');
     148md.spctemperature=project3d(md,md.spcvelovity,'node');
     149md.spcthickness=project3d(md,md.spcvelovity,'node');
    149150
    150151%Extrusion of Neumann BC
    151152%in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element]
    152 oldsegmentonneumann_diag=md.segmentonneumann_diag;
    153 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
    154 segmentonneumann_diag=[];
     153oldpressureload=md.pressureload;
     154pressureload_layer1=[oldpressureload(:,1:2)  oldpressureload(:,2)+md.numberofgrids2d  oldpressureload(:,1)+md.numberofgrids2d  oldpressureload(:,3)]; %Add two columns on the first layer
     155pressureload=[];
    155156for i=1:numlayers-1,
    156         segmentonneumann_diag=[segmentonneumann_diag ;segmentonneumann_diag_layer1(:,1:4)+(i-1)*md.numberofgrids2d segmentonneumann_diag_layer1(:,5)+(i-1)*md.numberofelements2d ];
     157        pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.numberofgrids2d pressureload_layer1(:,5)+(i-1)*md.numberofelements2d ];
    157158end
    158159
    159160%plug into md
    160 md.segmentonneumann_diag=segmentonneumann_diag;
    161 md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node');
    162 md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node');
    163 %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))];
    164 %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))];
     161md.pressureload=pressureload;
    165162
    166163%materials
     
    184181md.observed_temperature=project3d(md,md.observed_temperature,'node');
    185182md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux
    186 md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers); %surface temperature
    187 md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers); %surface temperature
    188 
    189 %NaN the values that are not on an spc'd temperature.
    190 pos=find(~md.gridondirichlet_thermal);
    191 md.dirichletvalues_thermal(pos)=NaN;
    192183
    193184%increase connectivity if less than 25:
    194185if md.connectivity<25,
    195         md.connectivity=48;
     186        md.connectivity=70;
    196187end
    197188
  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r1737 r1755  
    6666
    6767%NO NAN
    68 fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
     68fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
    6969        'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
    7070        'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
     
    7979
    8080%FIELDS >= 0
    81 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
     81fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
    8282        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
    8383        'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
     
    9696
    9797%FIELDS ~=0
    98 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
     98fields={'numberofelements','numberofgrids','elements','drag_type',...
    9999        'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
    100100        'sparsity','deltaH','DeltaH','timeacc','timedec'};
     
    118118
    119119%SIZE NUMBEROFGRIDS
    120 fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
     120fields={'x','y','z','B','drag','scpvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
    121121for i=1:length(fields),
    122122        if length(md.(fields{i}))~=md.numberofgrids,
     
    194194
    195195        %SINGULAR
    196         if ~any(md.gridondirichlet_diag),
     196        if ~any(sum(md.spcvelocity(:,1:3),2)==3),
    197197                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
    198                 bool=0;return;
    199         end
    200 
    201         %DIRICHLET VALUES
    202         if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
    203                 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
    204198                bool=0;return;
    205199        end
     
    208202        if any(md.thickness<=0),
    209203                pos=find(md.thickness<=0);
    210                 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
     204                if any(find(md.spcthickness(pos,1)==0)),
    211205                        disp(['model ' md.name ' has some grids with 0 thickness']);
    212206                        bool=0; return;
     
    217211%PROGNOSTIC
    218212if strcmp(md.analysis_type,'prognostic'),
    219 
    220         %DIRICHLET VALUES
    221         if isempty(md.gridondirichlet_prog) | isempty(md.dirichletvalues_prog),
    222                 disp(['model ' md.name ' is not well posed. Missing dirichlet values for prognostic run']);
    223                 bool=0;return;
    224         end
    225213
    226214        %VELOCITIES
     
    339327
    340328        %SINGULAR
    341         if ~any(md.gridondirichlet_diag),
     329        if ~any(sum(md.spcvelocity(:,1:3),2)==3),
    342330                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
    343                 bool=0;return;
    344         end
    345 
    346         %DIRICHLET VALUES
    347         if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
    348                 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
    349331                bool=0;return;
    350332        end
     
    353335        if any(md.thickness<=0),
    354336                pos=find(md.thickness<=0);
    355                 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
     337                if any(find(md.spcthickness(pos,1)==0)),
    356338                        disp(['model ' md.name ' has some grids with 0 thickness']);
    357339                        bool=0; return;
  • issm/trunk/src/m/classes/public/marshall.m

    r1666 r1755  
    7777WriteData(fid,md.gridoniceshelf,'Mat','gridoniceshelf');
    7878
    79 WriteData(fid,md.segmentonneumann_diag,'Mat','segmentonneumann_diag');
    80 WriteData(fid,md.segmentonneumann_diag_stokes,'Mat','segmentonneumann_diag_stokes');
    81 WriteData(fid,md.neumannvalues_diag,'Mat','neumannvalues_diag');
    82 WriteData(fid,md.gridondirichlet_diag,'Mat','gridondirichlet_diag');
    83 WriteData(fid,md.dirichletvalues_diag,'Mat','dirichletvalues_diag');
     79WriteData(fid,md.spcvelocity,'Mat','spcvelocity');
     80WriteData(fid,md.spctemperature,'Mat','spctemperature');
     81WriteData(fid,md.spcthickness,'Mat','spcthickness');
    8482
    85 WriteData(fid,md.segmentonneumann_prog,'Mat','segmentonneumann_prog');
    86 WriteData(fid,md.neumannvalues_prog,'Mat','neumannvalues_prog');
    87 WriteData(fid,md.gridondirichlet_prog,'Mat','gridondirichlet_prog');
    88 WriteData(fid,md.dirichletvalues_prog,'Mat','dirichletvalues_prog');
    89 WriteData(fid,md.segmentonneumann_prog2,'Mat','segmentonneumann_prog2');
    90 WriteData(fid,md.neumannvalues_prog2,'Mat','neumannvalues_prog2');
     83WriteData(fid,md.pressureload,'Mat','pressureload');
     84WriteData(fid,md.pressureload_stokes,'Mat','pressureload_stokes');
    9185
    92 WriteData(fid,md.gridondirichlet_thermal,'Mat','gridondirichlet_thermal');
    93 WriteData(fid,md.dirichletvalues_thermal,'Mat','dirichletvalues_thermal');
    9486WriteData(fid,md.geothermalflux,'Mat','geothermalflux');
    9587WriteData(fid,md.melting,'Mat','melting');
  • issm/trunk/src/m/classes/public/modelextract.m

    r1730 r1755  
    188188        gridstoflag1=intersect(orphans_grid,pos_grid);
    189189        gridstoflag2=Pgrid(gridstoflag1);
    190         if ~isnan(md1.gridondirichlet_diag),
    191                 md2.gridondirichlet_diag(gridstoflag2)=1;
     190        if ~isnan(md1.spcvelocity),
     191                md2.spcvelocity(gridstoflag2,1:2)=1;
    192192                if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
    193                         md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2);
    194                         md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
     193                        md2.spcvelocity(gridstoflag2,4)=md2.vx_obs(gridstoflag2);
     194                        md2.spcvelocity(gridstoflag2,5)=md2.vy_obs(gridstoflag2);
    195195                else
    196                         md2.dirichletvalues_diag=zeros(numberofgrids2,2);
    197                         disp(' ')
    198                         disp('!! modelextract warning: dirichlet values should be checked !!')
    199                         disp(' ')
    200                 end
    201         end
    202         if ~isnan(md1.gridondirichlet_thermal),
    203                 md2.gridondirichlet_thermal(gridstoflag2)=1;
     196                        md2.spcvelocity(gridstoflag2,4:5)=zeros(length(gridstoflag2)),2);
     197                        disp(' ')
     198                        disp('!! modelextract warning: spc values should be checked !!')
     199                        disp(' ')
     200                end
     201        end
     202        if ~isnan(md1.spctemperature),
     203                md2.spctemperature(gridstoflag2,1)=1;
    204204                if ~isnan(md1.observed_temperature)
    205                         md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2);
     205                        md2.spctemperature(gridstoflag2,2)=md2.observed_temperature(gridstoflag2);
    206206                else
    207                         md2.dirichletvalues_thermal=zeros(numberofgrids2,2);
    208                         disp(' ')
    209                         disp('!! modelextract warning: dirichlet values should be checked !!')
     207                        md2.spctemperature(gridstoflag2,2)=zeros(length(gridstoflag2),2);
     208                        disp(' ')
     209                        disp('!! modelextract warning: spc values should be checked !!')
    210210                        disp(' ')
    211211                end
     
    214214        %Recreate border stokes
    215215        if md1.isstokes,
    216                 pos=find(md2.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
     216                pos=find(sum(md2.spcvelocity(:,1:3),2));                 %find all the grids on the boundary of the domain without icefront
    217217                md2.gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront
    218218                stokes_elements=find(md2.elements_type(:,2)==StokesFormulationEnum()); %find the elements on the stokes domain
     
    223223
    224224        %Diagnostic
    225         if ~isnan(md2.segmentonneumann_diag)
    226                 md2.segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1));
    227                 md2.segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2));
    228                 md2.segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end));
     225        if ~isnan(md2.pressureload)
     226                md2.pressureload(:,1)=Pgrid(md1.pressureload(:,1));
     227                md2.pressureload(:,2)=Pgrid(md1.pressureload(:,2));
     228                md2.pressureload(:,end)=Pelem(md1.pressureload(:,end));
    229229                if strcmpi(md1.type,'3d')
    230                         md2.segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3));
    231                         md2.segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4));
    232                 end
    233                 md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);
    234         end
    235         if ~isnan(md2.segmentonneumann_diag_stokes)
    236                 md2.segmentonneumann_diag_stokes(:,1)=Pgrid(md1.segmentonneumann_diag_stokes(:,1));
    237                 md2.segmentonneumann_diag_stokes(:,2)=Pgrid(md1.segmentonneumann_diag_stokes(:,2));
    238                 md2.segmentonneumann_diag_stokes(:,end)=Pelem(md1.segmentonneumann_diag_stokes(:,end));
     230                        md2.pressureload(:,3)=Pgrid(md1.pressureload(:,3));
     231                        md2.pressureload(:,4)=Pgrid(md1.pressureload(:,4));
     232                end
     233                md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2)),:);
     234        end
     235        if ~isnan(md2.pressureload_stokes)
     236                md2.pressureload_stokes(:,1)=Pgrid(md1.pressureload_stokes(:,1));
     237                md2.pressureload_stokes(:,2)=Pgrid(md1.pressureload_stokes(:,2));
     238                md2.pressureload_stokes(:,end)=Pelem(md1.pressureload_stokes(:,end));
    239239                if strcmpi(md1.type,'3d')
    240                         md2.segmentonneumann_diag_stokes(:,3)=Pgrid(md1.segmentonneumann_diag_stokes(:,3));
    241                         md2.segmentonneumann_diag_stokes(:,4)=Pgrid(md1.segmentonneumann_diag_stokes(:,4));
    242                 end
    243                 md2.segmentonneumann_diag_stokes=md2.segmentonneumann_diag_stokes(find(md2.segmentonneumann_diag_stokes(:,1) & md2.segmentonneumann_diag_stokes(:,2)),:);
    244         end
    245         md2.neumannvalues_diag=NaN;
    246 
    247         %Transient
    248         if ~isnan(md2.segmentonneumann_prog)
    249                 md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1));
    250                 md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2));
    251                 md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end));
    252                 md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:);
    253         end
    254         md2.neumannvalues_prog=NaN;
    255         if ~isnan(md2.segmentonneumann_prog2)
    256                 md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1));
    257                 md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2));
    258                 md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end));
    259                 md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:);
     240                        md2.pressureload_stokes(:,3)=Pgrid(md1.pressureload_stokes(:,3));
     241                        md2.pressureload_stokes(:,4)=Pgrid(md1.pressureload_stokes(:,4));
     242                end
     243                md2.pressureload_stokes=md2.pressureload_stokes(find(md2.pressureload_stokes(:,1) & md2.pressureload_stokes(:,2)),:);
    260244        end
    261245
  • issm/trunk/src/m/classes/public/setelementstype.m

    r1682 r1755  
    6161gridonstokes=ones(md.numberofgrids,1);
    6262gridonstokes(md.elements(nonstokesflag,:))=0;      %non stokes grids
    63 pos=find(md.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
     63pos=find(sum(md.spcvelocity(:,1:3),2)==3);         %find all the grids on the boundary of the domain without icefront
    6464gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront
    6565
Note: See TracChangeset for help on using the changeset viewer.