0001 function pe=CreatePVector(singelem,grids,materials,inputs,analysis_type);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if strcmpi(analysis_type,'diagnostic_hutter'),
0013
0014 pe=CreatePVectorHutter(singelem,grids,materials,inputs);
0015
0016 else
0017 error('singelem/CreatePVector error message: analysis type not implemented yet')
0018 end
0019 end
0020
0021 function pe=CreatePVectorHutter(singelem,grids,materials,inputs);
0022
0023
0024 numgrids=1;
0025 NDOF2=2;
0026
0027
0028 pe=elemvector(numgrids*NDOF2);
0029
0030
0031 matice=materials(singelem.matid).material;
0032 matpar=materials(end).constants;
0033
0034
0035 gravity=matpar.g;
0036 rho_ice=matpar.rho_ice;
0037 n=matice.n;
0038
0039
0040 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0041 [B_param B_is_present]=recover_input(inputs,'B');
0042 [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0043 [slopesurface_param slopesurface_is_present]=recover_input(inputs,'slopesurface');
0044
0045 if ~slopesurface_is_present,
0046 error('CreatePVector/singelem error message: slope_surface is missing in inputs');
0047 end
0048
0049
0050 B_list=zeros(numgrids,1);
0051 temperature_list=zeros(numgrids,1);
0052 slopesurface_list=zeros(numgrids,2);
0053
0054
0055 for i=1:numgrids,
0056 doflist=grids(singelem.g(i)).grid.doflist;
0057 for j=1:NDOF2,
0058 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0059 slopesurface_list(i,j)=slopesurface_param(doflist(j));
0060 end
0061 dof=doflist(1);
0062 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0063 if(B_is_present) B_list(i)=B_param(dof);end;
0064 if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0065 end
0066
0067
0068 slope=zeros(2,1);
0069 slope(1)=slopesurface_list(1,1);
0070 slope(2)=slopesurface_list(1,2);
0071 slope2=slope(1)^2+slope(2)^2;
0072
0073
0074 if thickness_is_present
0075 thickness=thickness_list(1);
0076 else
0077 thickness=singelem.h(1);
0078 end
0079
0080
0081 ub=-1.58*10^-10*rho_ice*gravity*thickness*slope(1);
0082 vb=-1.58*10^-10*rho_ice*gravity*thickness*slope(2);
0083
0084
0085 constant_part=-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2);
0086
0087
0088 if B_is_present,
0089 B=GetParameterValue(singelem,B_list,gauss_coord);
0090 elseif temperature_is_present,
0091 temperature=GetParameterValue(singelem,temperature_list,gauss_coord);
0092 B=paterson(temperature);
0093 else
0094 B=matice.B;
0095 end
0096
0097 pe.terms(1)=ub-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(1);
0098 pe.terms(2)=vb-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(2);
0099
0100 end