0001 function [elements,grids,loads,constraints,materials,part,tpart]=ModelProcessorMelting(md,solutiontype);
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 elements=struct('element',cell(md.numberofelements,1));
0026 materials=struct('material',cell(md.numberofelements+1,1));
0027 mygrids=zeros(md.numberofgrids,1);
0028
0029
0030 if strcmpi(md.type,'2d'),
0031 error('ModelProcessorDiagnosticBaseVert error message: 2d mesh not supported yet!');
0032 end
0033
0034 pos=find(element_partitioning==labindex);
0035 [elements(pos).element]=deal(pentaelem);
0036
0037 elements(pos)=SetStructureField(elements(pos),'element','type','pentaelem');
0038 elements(pos)=SetStructureField(elements(pos),'element','id',pos);
0039 elements(pos)=SetStructureField(elements(pos),'element','g',md.elements(pos,1:6));
0040 elements(pos)=SetStructureField(elements(pos),'element','h',md.thickness(md.elements(pos,1:6)));
0041 elements(pos)=SetStructureField(elements(pos),'element','s',md.surface(md.elements(pos,1:6)));
0042 elements(pos)=SetStructureField(elements(pos),'element','b',md.bed(md.elements(pos,1:6)));
0043 elements(pos)=SetStructureField(elements(pos),'element','shelf',md.elementoniceshelf(pos));
0044 elements(pos)=SetStructureField(elements(pos),'element','onbed',md.elementonbed(pos));
0045 elements(pos)=SetStructureField(elements(pos),'element','onsurface',md.elementonsurface(pos));
0046 elements(pos)=SetStructureField(elements(pos),'element','acceleration',0);
0047 elements(pos)=SetStructureField(elements(pos),'element','collapse',1);
0048 elements(pos)=SetStructureField(elements(pos),'element','matid',pos);
0049 elements(pos)=SetStructureField(elements(pos),'element','melting',md.melting(md.elements(pos,1:6))/md.yts);
0050 elements(pos)=SetStructureField(elements(pos),'element','accumulation',md.accumulation(md.elements(pos,1:6))/md.yts);
0051
0052 [materials(pos).material]=deal(matice);
0053 materials(pos)=SetStructureField(materials(pos),'material','id',pos);
0054
0055
0056 [materials(end).constants]=deal(matpar);
0057 materials(end)=SetStructureField(materials(end),'constants','g',md.g);
0058 materials(end)=SetStructureField(materials(end),'constants','rho_ice',md.rho_ice);
0059 materials(end)=SetStructureField(materials(end),'constants','rho_water',md.rho_water);
0060 materials(end)=SetStructureField(materials(end),'constants','thermalconductivity',md.thermalconductivity);
0061 materials(end)=SetStructureField(materials(end),'constants','heatcapacity',md.heatcapacity);
0062 materials(end)=SetStructureField(materials(end),'constants','latentheat',md.latentheat);
0063 materials(end)=SetStructureField(materials(end),'constants','beta',md.beta);
0064 materials(end)=SetStructureField(materials(end),'constants','meltingpoint',md.meltingpoint);
0065 materials(end)=SetStructureField(materials(end),'constants','mixed_layer_capacity',md.mixed_layer_capacity);
0066 materials(end)=SetStructureField(materials(end),'constants','thermal_exchange_velocity',md.thermal_exchange_velocity);
0067
0068 if cluster,
0069
0070 mygrids(md.elements(el3pos,:))=1;
0071 mygrids(md.elements(el6pos,:))=1;
0072
0073
0074
0075 bordergrids=double(gplus(mygrids)>1);
0076 else
0077 bordergrids=zeros(md.numberofgrids,1);
0078 end
0079
0080
0081 grids=struct('grid',cell(md.numberofgrids,1));
0082
0083 pos=[1:md.numberofgrids]';
0084 [grids(pos).grid]=deal(node);
0085 grids(pos)=SetStructureField(grids(pos),'grid','id',pos);
0086 grids(pos)=SetStructureField(grids(pos),'grid','x',md.x(pos));
0087 grids(pos)=SetStructureField(grids(pos),'grid','y',md.y(pos));
0088 grids(pos)=SetStructureField(grids(pos),'grid','z',md.z(pos));
0089 grids(pos)=SetStructureField(grids(pos),'grid','onbed',md.gridonbed(pos));
0090 grids(pos)=SetStructureField(grids(pos),'grid','border',bordergrids(pos));
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103 grids(pos)=SetStructureField(grids(pos),'grid','gridset','123456');
0104
0105
0106 for n=1:length(elements),
0107
0108
0109 if elements(n).element.onbed,
0110 for j=1:3,
0111 grids(elements(n).element.g(j)).grid.gridset='23456';
0112 end
0113 end
0114 end
0115
0116 loads=struct('load',cell([0,1]));
0117
0118
0119 count=1;
0120 for i=1:md.numberofgrids,
0121 if grids(i).grid.onbed,
0122 pengridobject=pengrid;
0123 pengridobject.id=count;
0124 pengridobject.dof=1;
0125 pengridobject.grid=i;
0126 if strcmpi(solutiontype,'meltingsteady'),
0127 pengridobject.thermal_steadystate=1;
0128 else
0129 pengridobject.thermal_steadystate=0;
0130 end
0131 pengridobject.penalty_offset=md.penalty_offset;
0132
0133 loads(count).load=pengridobject;
0134 count=count+1;
0135 end
0136 end
0137
0138
0139
0140 constraints=struct('constraint',cell(0,0));
0141
0142
0143 [part,tpart]=PartitioningVector(md,grids);
0144
0145 end