0001 function Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 global permanent_pieces_of_Ke alpha beta row_location col_location aire
0015
0016
0017
0018
0019 numdof=3*2;
0020
0021
0022 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0023 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0024 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0025 [surface_param surface_is_present]=recover_input(inputs,'surface');
0026 [bed_param bed_is_present]=recover_input(inputs,'bed');
0027 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0028 [oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');
0029
0030 if velocity_is_present,
0031 vx=velocity_param(1:6:numdof*element.nods);
0032 vy=velocity_param(2:6:numdof*element.nods);
0033 else
0034 vx={};
0035 vy={};
0036 permanent_pieces_of_Ke={};
0037 end
0038 if oldvelocity_is_present,
0039 oldvx=oldvelocity_param(1:6:numdof*element.nods);
0040 oldvy=oldvelocity_param(2:6:numdof*element.nods);
0041 else
0042 oldvx={};
0043 oldvy={};
0044 end
0045 if flow_law_is_present,
0046 B_list=flow_law_param(1:6:6*element.nods);
0047 B_bar=1/3*(B_list(element.index(:,1))+B_list(element.index(:,2))+B_list(element.index(:,3)));
0048 else
0049 B_bar=element.B_bar;
0050 end
0051
0052 if basal_drag_is_present,
0053 drag=basal_drag_param(1:6:6*element.nods);
0054 else
0055 drag=element.drag;
0056 end
0057
0058 if thickness_is_present,
0059 thickness=thickness_param(1:6:6*element.nods);
0060 else
0061 thickness=element.thickness;
0062 end
0063
0064 thickness_el=thickness(element.index)*[1;1;1]/3;
0065
0066 if surface_is_present,
0067 surface=surface_param(1:6:6*element.nods);
0068 else
0069 surface=element.surface;
0070 end
0071
0072 if bed_is_present,
0073 bed=bed_param(1:6:6*element.nods);
0074 else
0075 bed=element.bed;
0076 end
0077
0078 if isempty(permanent_pieces_of_Ke),
0079
0080 aire=zeros(element.nel,1);
0081
0082 for n=1:element.nel
0083 aire(n)=1/2 * det([1 1 1;element.x(element.index(n,:))';element.y(element.index(n,:))']);
0084 end
0085
0086 aire=abs(aire);
0087
0088 alpha=zeros(element.nel,3);
0089 beta=zeros(element.nel,3);
0090 gamma=zeros(element.nel,3);
0091
0092 for n=1:element.nel
0093 X=inv([element.x(element.index(n,:)) element.y(element.index(n,:)) ones(3,1)]);
0094 alpha(n,:)=X(1,:);
0095 beta(n,:)=X(2,:);
0096 gamma(n,:)=X(3,:);
0097 end
0098
0099 clear X
0100
0101
0102
0103 row_location=zeros(4*element.nel*9,1);
0104 col_location=zeros(4*element.nel*9,1);
0105
0106 count=-element.nel+1;
0107
0108 for i=1:3
0109 for j=1:3
0110 for k=1:2,
0111 for l=1:2,
0112 count=count+element.nel;
0113 row_location(count:count+element.nel-1)=(element.index(:,i)-1)*numdof+k;
0114 col_location(count:count+element.nel-1)=(element.index(:,j)-1)*numdof+l;
0115 end
0116 end
0117 end
0118 end
0119
0120
0121
0122 permanent_pieces_of_Ke=zeros(2*element.nel*9,1);
0123
0124 count=-element.nel+1;
0125
0126 for i=1:3
0127 for j=1:3
0128 for k=1:2,
0129 for l=1:2,
0130 count=count+element.nel;
0131
0132 permanent_pieces_of_Ke(count:count+element.nel-1)= thickness_el .* aire.*(...
0133 ((k==1) & (l==1)).*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j)) +...
0134 ((k==2) & (l==1)).*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i)) + ...
0135 ((k==1) & (l==2)).*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)) + ...
0136 ((k==2) & (l==2)).*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j)) ...
0137 );
0138 end
0139 end
0140 end
0141 end
0142
0143 end
0144
0145
0146 viscosity_overshoot=materials(end).constants.viscosity_overshoot;
0147 newviscosity=viscosity(element.index,element.nel,alpha,beta,vx,vy,B_bar,element.glen_coeff);
0148 oldviscosity=viscosity(element.index,element.nel,alpha,beta,oldvx,oldvy,B_bar,element.glen_coeff);
0149 nu_bar=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);
0150
0151
0152 count=-element.nel+1;
0153
0154
0155 Ke_terms=zeros(4*element.nel*9,1);
0156
0157 for i=1:3
0158 for j=1:3
0159 for k=1:2,
0160 for l=1:2,
0161 count=count+element.nel;
0162 Ke_terms(count:(count+element.nel-1))=nu_bar.*permanent_pieces_of_Ke(count:(count+element.nel-1));
0163 end
0164 end
0165 end
0166 end
0167
0168 Drag_operator=sparse2(6*element.nods,6*element.nods);
0169
0170
0171 if (element.friction_type==2) & velocity_is_present
0172
0173
0174 pcoeff_grid=zeros(element.nods,1);
0175 qcoeff_grid=zeros(element.nods,1);
0176 for i=1:element.nods
0177
0178 neighbors_gridi=[];
0179 for j=1:3
0180 neighbors_gridi=[neighbors_gridi find(element.index(:,j)==i)'];
0181 end
0182 numberofneighbors_gridi=length(neighbors_gridi);
0183
0184
0185 qcoeff_grid(i)=sum(element.q(neighbors_gridi))/numberofneighbors_gridi;
0186 pcoeff_grid(i)=sum(element.p(neighbors_gridi))/numberofneighbors_gridi;
0187 end
0188
0189
0190 frictionparameters=struct();
0191 frictionparameters.element_type='2d';
0192 frictionparameters.rho_ice=materials(end).constants.rho_ice;
0193 frictionparameters.rho_water=materials(end).constants.rho_water;
0194 frictionparameters.g=materials(end).constants.g;
0195 frictionparameters.p=pcoeff_grid;
0196 frictionparameters.q=qcoeff_grid;
0197 frictionparameters.velocities=[vx vy];
0198 frictionparameters.h=thickness;
0199 frictionparameters.b=bed;
0200 frictionparameters.k=drag;
0201
0202 alpha2=Getalpha2(frictionparameters);
0203
0204 drag_operator_value=zeros(2*element.nel*27,1);
0205
0206 row_location_drag=zeros(2*element.nel*27,1);
0207 col_location_drag=zeros(2*element.nel*27,1);
0208
0209 count=-element.nel+1;
0210
0211 for i=1:3
0212 for j=1:3
0213 for m=1:3
0214 for k=1:2,
0215
0216 count=count+element.nel;
0217
0218 row_location_drag(count:count+element.nel-1)=(element.index(:,i)-1)*numdof+k;
0219 col_location_drag(count:count+element.nel-1)=(element.index(:,j)-1)*numdof+k;
0220 drag_operator_value(count:count+element.nel-1)=aire.*(...
0221 ((i==j) & (j==m)).*alpha2(element.index(:,m))*1/10+...
0222 (((i==j) | (j==m) | (i==m)) & ~((i==j) & (j==m))).*alpha2(element.index(:,m))*1/30+...
0223 ((i~=j) & (j~=m) & (i~=m)).*alpha2(element.index(:,m))*1/60);
0224
0225 end
0226 end
0227 end
0228 end
0229
0230 Drag_operator=sparse2(row_location_drag,col_location_drag,drag_operator_value,6*element.nods,6*element.nods);
0231 end
0232
0233
0234 Ke=globalmatrix;
0235 Ke.nrows=numdof*element.nods;
0236 Ke.matrix=Drag_operator+sparse2(row_location,col_location,Ke_terms,6*element.nods,6*element.nods);