0001 function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorDiagnosticHoriz(md);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 global cluster
0014
0015 if cluster,
0016
0017 element_partitioning=MeshPartition(md,numlabs);
0018 else
0019
0020 element_partitioning=ones(md.numberofelements,1);
0021 labindex=1;
0022 end
0023
0024
0025 if strcmpi(md.type,'2d'),
0026 numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum()));
0027 else
0028 numberofelements_horiz=length(find(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()));
0029 end
0030 elements=struct('element',cell(numberofelements_horiz,1));
0031 materials=struct('material',cell(numberofelements_horiz+1,1));
0032 mygrids=zeros(md.numberofgrids,1);
0033
0034
0035 if strcmpi(md.type,'2d'),
0036 el3pos=find((md.elements_type(:,1)==macayealenum()) & (element_partitioning==labindex));
0037 if ~isempty(el3pos),
0038 pos=[1:length(el3pos)]';
0039 [elements(pos).element]=deal(triaelem);
0040
0041 elements(pos)=SetStructureField(elements(pos),'element','type','triaelem');
0042 elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0043 elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(el3pos,:));
0044 elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(el3pos,1:3)));
0045 elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(el3pos,1:3)));
0046 elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(el3pos,1:3)));
0047 elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
0048 elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(el3pos,1:3)));
0049 elements(pos)=SetStructureField(elements(pos),'element','p',md.p(el3pos));
0050 elements(pos)=SetStructureField(elements(pos),'element','q',md.q(el3pos));
0051 elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(el3pos));
0052 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
0053 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
0054 elements(pos)=SetStructureField(elements(pos),'element','acceleration',md.acceleration);
0055 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0056
0057
0058 [materials(pos).material]=deal(matice);
0059
0060 materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0061 materials(pos)=SetStructureField(materials(pos),'material','B',md.B(md.elements(el3pos,1:3))*[1;1;1]/3);
0062 materials(pos)=SetStructureField(materials(pos),'material','n',md.n(el3pos));
0063 end
0064 else
0065
0066
0067
0068 el6pos=find((md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum()) & (element_partitioning==labindex));
0069
0070 if ~isempty(el6pos),
0071 pos=[1:length(el6pos)]';
0072 [elements(pos).element]=deal(pentaelem);
0073
0074 elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem');
0075 elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0076 elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(el6pos,1:6));
0077 elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(el6pos,1:6)));
0078 elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(el6pos,1:6)));
0079 elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(el6pos,1:6)));
0080 elements(pos)=SetStructureField(elements(pos),'element','friction_type',md.drag_type);
0081 elements(pos)=SetStructureField(elements(pos),'element','k',md.drag(md.elements(el6pos,1:6)));
0082 elements(pos)=SetStructureField(elements(pos),'element','p',md.p(el6pos));
0083 elements(pos)=SetStructureField(elements(pos),'element','q',md.q(el6pos));
0084 elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(el6pos));
0085 elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(el6pos));
0086 elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(el6pos));
0087 elements(pos)=SetStructureField(elements(pos),'element','meanvel',md.meanvel);
0088 elements(pos)=SetStructureField(elements(pos),'element','epsvel',md.epsvel);
0089 elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);
0090 elements(pos)=SetStructureField(elements(pos),'element','collapse',0);
0091 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0092
0093 [materials(pos).material]=deal(matice);
0094
0095 materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0096 materials(pos)=SetStructureField(materials(pos),'material','B',md.B(md.elements(el6pos,1:6))*[1;1;1;1;1;1]/6);
0097 materials(pos)=SetStructureField(materials(pos),'material','n',md.n(el6pos));
0098
0099
0100
0101 pos=find(ismember(el6pos,find(md.elements_type(:,1)==macayealenum())));
0102 elements(pos)=SetStructureField(elements(pos),'element','collapse',ones(length(pos),1));
0103
0104 end
0105
0106 end
0107
0108
0109
0110 [materials(end).constants]=deal(matpar);
0111 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0112 materials(end)=SetStructureField(materials(end),'constants','viscosity_overshoot',md.viscosity_overshoot);
0113 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0114 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0115 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0116 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0117 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0118 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0119 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0120 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0121 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0122
0123 if cluster,
0124
0125 mygrids(md.elements(el3pos,:))=1;
0126 mygrids(md.elements(el6pos,:))=1;
0127 end
0128
0129
0130 if md.acceleration & strcmpi(md.type,'2d'),
0131
0132 elements(end+1).element=acceleratedtriaelem;
0133 elements(end).element.type='acceleratedtriaelem';
0134 elements(end).element.x=md.x;
0135 elements(end).element.y=md.y;
0136 elements(end).element.index=md.elements;
0137 elements(end).element.nods=md.numberofgrids;
0138 elements(end).element.nel=md.numberofelements;
0139 elements(end).element.thickness=md.thickness;
0140 elements(end).element.surface=md.surface;
0141 elements(end).element.bed=md.bed;
0142 elements(end).element.thickness_el=(md.thickness(md.elements))*[1;1;1]/3;
0143 elements(end).element.B_bar=(md.B(md.elements))*[1;1;1]/3;
0144 elements(end).element.glen_coeff=md.n;
0145 elements(end).element.friction_type=md.drag_type;
0146 elements(end).element.drag=md.drag;
0147 elements(end).element.p=md.p;
0148 elements(end).element.q=md.q;
0149 elements(end).element.meanvel=md.meanvel;
0150 elements(end).element.epsvel=md.epsvel;
0151 if isempty(md.segmentonneumann_diag),
0152 elements(end).element.index_icefront=[];
0153 else
0154 elements(end).element.index_icefront=md.segmentonneumann_diag(:,1:2);
0155 end
0156 elements(end).element.gridoniceshelf=md.gridoniceshelf;
0157 elements(end).element.elementonicesheet=md.elementonicesheet;
0158 elements(end).element.matid=1;
0159 end
0160
0161 if cluster,
0162
0163
0164 bordergrids=double(gplus(mygrids)>1);
0165 else
0166 bordergrids=zeros(md.numberofgrids,1);
0167 end
0168
0169
0170 grids=struct('grid',cell(md.numberofgrids,1));
0171
0172 pos=[1:md.numberofgrids]';
0173 [grids(pos).grid]=deal(node);
0174 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0175 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0176 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0177 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0178 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0179 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 if strcmpi(md.type,'3d')
0193 grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456');
0194 else
0195 grids(pos)=SetStructureField(grids(pos),'grid','gridset','3456');
0196 end
0197
0198
0199 if strcmpi(md.type,'3d'),
0200 for n=1:length(elements),
0201 if md.elements_type(el6pos(n),1)==macayealenum(),
0202
0203
0204 if elements(n).element.onbed,
0205 for j=1:3,
0206 grids(elements(n).element.g(j)).grid.gridset='3456';
0207 end
0208 end
0209 else
0210 for j=1:length(elements(n).element.g),
0211 grids(elements(n).element.g(j)).grid.gridset='3456';
0212 end
0213 end
0214 end
0215 end
0216
0217
0218
0219
0220 loads=struct('load',cell(length(md.segmentonneumann_diag),1));
0221
0222 for i=1:size(md.segmentonneumann_diag,1),
0223
0224 if (element_partitioning(md.segmentonneumann_diag(i,end))~=labindex),
0225 continue;
0226 end
0227
0228 if strcmpi(md.type,'3d'),
0229 if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
0230 loads(i).load=icefront;
0231 loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
0232 loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0233 loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0234 loads(i).load.rho_water=md.rho_water;
0235 loads(i).load.type='segment';
0236
0237 elseif md.elements_type(md.segmentonneumann_diag(i,end))==pattynenum(),
0238
0239 loads(i).load=icefront;
0240 loads(i).load.eid=find(el6pos==md.segmentonneumann_diag(i,end));
0241 loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0242 loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0243 loads(i).load.g(3)=md.segmentonneumann_diag(i,3);
0244 loads(i).load.g(4)=md.segmentonneumann_diag(i,4);
0245 loads(i).load.type='quad';
0246
0247 else
0248 if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
0249 error('diagnostic error message: unsupported element type');
0250 end
0251 end
0252 else
0253 if md.elements_type(md.segmentonneumann_diag(i,end))==macayealenum(),
0254 loads(i).load=icefront;
0255 loads(i).load.eid=find(el3pos==md.segmentonneumann_diag(i,end));
0256 loads(i).load.g(1)=md.segmentonneumann_diag(i,1);
0257 loads(i).load.g(2)=md.segmentonneumann_diag(i,2);
0258 loads(i).load.rho_water=md.rho_water;
0259 loads(i).load.type='segment';
0260
0261 else
0262 if ~(md.elements_type(md.segmentonneumann_diag(i,end))==hutterenum()),
0263 error('diagnostic error message: unsupported element type');
0264 end
0265 end
0266
0267 end
0268
0269 end
0270
0271
0272 constraints=struct('constraint',cell(2*length(find(md.gridondirichlet_diag)),1));
0273
0274 count=1;
0275 for i=1:md.numberofgrids,
0276
0277 if md.gridonhutter(i),
0278
0279
0280 constraints(count).constraint=spc;
0281 constraints(count).constraint.grid=i;
0282 constraints(count).constraint.dof=1;
0283 constraints(count).constraint.value=0;
0284 count=count+1;
0285
0286
0287 constraints(count).constraint=spc;
0288 constraints(count).constraint.grid=i;
0289 constraints(count).constraint.dof=2;
0290 constraints(count).constraint.value=0;
0291 count=count+1;
0292
0293 elseif (md.gridonmacayeal(i) | md.gridonpattyn(i)) & md.gridondirichlet_diag(i),
0294
0295
0296 constraints(count).constraint=spc;
0297 constraints(count).constraint.grid=i;
0298 constraints(count).constraint.dof=1;
0299 constraints(count).constraint.value=md.dirichletvalues_diag(i,1)/md.yts;
0300 count=count+1;
0301
0302
0303 constraints(count).constraint=spc;
0304 constraints(count).constraint.grid=i;
0305 constraints(count).constraint.dof=2;
0306 constraints(count).constraint.value=md.dirichletvalues_diag(i,2)/md.yts;
0307 count=count+1;
0308 end
0309 end
0310
0311
0312 if ~isempty(md.penalties) & ~isnan(md.penalties),
0313 for i=1:size(md.penalties,1),
0314 for j=1:(md.numlayers-1),
0315
0316
0317 constraints(count).constraint=rgb;
0318 constraints(count).constraint.grid1=md.penalties(i,1);
0319 constraints(count).constraint.grid2=md.penalties(i,j+1);
0320 constraints(count).constraint.dof=1;
0321 count=count+1;
0322
0323
0324 constraints(count).constraint=rgb;
0325 constraints(count).constraint.grid1=md.penalties(i,1);
0326 constraints(count).constraint.grid2=md.penalties(i,j+1);
0327 constraints(count).constraint.dof=2;
0328 count=count+1;
0329 end
0330 end
0331 end
0332
0333
0334 [part,tpart]=PartitioningVector(md,grids);
0335
0336 end