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 numberofstokeselements=size(find(md.elements_type(:,2)==stokesenum()),1);
0026 elements=struct('element',cell(numberofstokeselements,1));
0027 materials=struct('material',cell(numberofstokeselements+1,1));
0028 mygrids=zeros(md.numberofgrids,1);
0029
0030
0031
0032 pos=find(md.elements_type(:,2)==stokesenum());
0033 stokesnewnumber=[1:numberofstokeselements]';
0034 [elements(stokesnewnumber).element]=deal(pentaelem);
0035
0036 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','type','pentaelem');
0037 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','id',stokesnewnumber);
0038 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','g',md.elements(pos,1:6));
0039 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','h',md.thickness(md.elements(pos,1:6)));
0040 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','s',md.surface(md.elements(pos,1:6)));
0041 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','b',md.bed(md.elements(pos,1:6)));
0042 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','friction_type',md.drag_type);
0043 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','k',md.drag(md.elements(pos,1:6)));
0044 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','p',md.p(pos));
0045 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','q',md.q(pos));
0046 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','shelf',md.elementoniceshelf(pos));
0047 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onbed',md.elementonbed(pos));
0048 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','onsurface',md.elementonsurface(pos));
0049 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','meanvel',md.meanvel);
0050 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','epsvel',md.epsvel);
0051 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','acceleration',0);
0052 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','collapse',0);
0053 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','matid',stokesnewnumber);
0054 elements(stokesnewnumber)=SetStructureField(elements(stokesnewnumber),'element','reconditioning_number',md.stokesreconditioning);
0055
0056 [materials(stokesnewnumber).material]=deal(matice);
0057
0058 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','id',stokesnewnumber);
0059 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','B',md.B(md.elements(pos,1:6))*[1;1;1;1;1;1]/6);
0060 materials(stokesnewnumber)=SetStructureField(materials(stokesnewnumber),'material','n',md.n(pos));
0061
0062
0063
0064 [materials(end).constants]=deal(matpar);
0065 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0066 materials(end)=SetStructureField(materials(end),'constants','viscosity_overshoot',md.viscosity_overshoot);
0067 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0068 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0069 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0070 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0071 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0072 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0073 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0074 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0075 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0076
0077
0078 if cluster,
0079
0080 mygrids(md.elements(pos,:))=1;
0081 end
0082
0083 if cluster,
0084
0085
0086 bordergrids=double(gplus(mygrids)>1);
0087 else
0088 bordergrids=zeros(md.numberofgrids,1);
0089 end
0090
0091
0092 grids=struct('grid',cell(md.numberofgrids,1));
0093
0094 pos=[1:md.numberofgrids]';
0095 [grids(pos).grid]=deal(node);
0096 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0097 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0098 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0099 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0100 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0101 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114 grids(pos)=SetStructureField(grids(pos),'grid','gridset','56');
0115
0116
0117
0118 stokes_elements=find(md.elements_type(:,2)==stokesenum());
0119 borderflags=zeros(md.numberofgrids,1);
0120 borderflags(md.elements(stokes_elements,:))=1;
0121 borderstokes=borderflags-md.gridonstokes;
0122
0123
0124 for n=1:md.numberofgrids,
0125 if md.gridonstokes(n)==0,
0126 grids(n).grid.gridset='123456';
0127 end
0128 if borderstokes(n),
0129 grids(n).grid.gridset='12356';
0130 end
0131 end
0132
0133
0134
0135
0136 segmentonneumann_diag_stokes=md.segmentonneumann_diag(find(md.elements_type(md.segmentonneumann_diag(:,end),2)),:);
0137 length_segmentonneumann_diag_stokes= size(segmentonneumann_diag_stokes,1);
0138 loads=struct('load',cell(length_segmentonneumann_diag_stokes,1));
0139
0140 for i=1:length_segmentonneumann_diag_stokes,
0141
0142 if (element_partitioning(segmentonneumann_diag_stokes(i,end))~=labindex),
0143 continue;
0144 end
0145
0146
0147 loads(i).load=icefront;
0148 loads(i).load.eid=find(segmentonneumann_diag_stokes(i,end)==find(md.elements_type(:,2)==stokesenum()));
0149 loads(i).load.g(1)=segmentonneumann_diag_stokes(i,1);
0150 loads(i).load.g(2)=segmentonneumann_diag_stokes(i,2);
0151 loads(i).load.g(3)=segmentonneumann_diag_stokes(i,3);
0152 loads(i).load.g(4)=segmentonneumann_diag_stokes(i,4);
0153 loads(i).load.rho_water=md.rho_water;
0154 loads(i).load.type='quad';
0155 end
0156
0157
0158 count=length(loads)+1;
0159 for i=1:md.numberofgrids,
0160 if grids(i).grid.onbed & md.gridonicesheet(i) & md.gridonstokes(i),
0161 pengridobject=pengrid;
0162 pengridobject.id=count;
0163 pengridobject.dof=1;
0164 pengridobject.grid=i;
0165 pengridobject.penalty_offset=md.penalty_offset;
0166
0167 loads(count).load=pengridobject;
0168 count=count+1;
0169 end
0170 end
0171
0172
0173 spcs=find(md.gridondirichlet_diag);
0174 constraints=struct('constraint',cell(3*length(spcs),1));
0175
0176 count=1;
0177 for i=1:md.numberofgrids,
0178 if ~md.gridonstokes(i),
0179
0180
0181 constraints(count).constraint=spc;
0182 constraints(count).constraint.grid=i;
0183 constraints(count).constraint.dof=1;
0184 constraints(count).constraint.value=0;
0185 count=count+1;
0186
0187
0188 constraints(count).constraint=spc;
0189 constraints(count).constraint.grid=i;
0190 constraints(count).constraint.dof=2;
0191 constraints(count).constraint.value=0;
0192 count=count+1;
0193
0194
0195 constraints(count).constraint=spc;
0196 constraints(count).constraint.grid=i;
0197 constraints(count).constraint.dof=3;
0198 constraints(count).constraint.value=0;
0199 count=count+1;
0200 end
0201 end
0202
0203
0204
0205 [part,tpart]=PartitioningVector(md,grids);
0206
0207
0208 end