0001 function pe=CreatePVector(beamelem,grids,materials,inputs,analysis_type);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 if strcmpi(analysis_type,'diagnostic_hutter'),
0012
0013 pe=CreatePVectorHutter(beamelem,grids,materials,inputs);
0014
0015 else
0016 error('beamelem/CreatePVector error message: analysis type not implemented yet')
0017 end
0018 end
0019
0020 function pe=CreatePVectorHutter(beamelem,grids,materials,inputs);
0021
0022
0023 numgrids=2;
0024 NDOF2=2;
0025
0026
0027 pe=elemvector(numgrids*NDOF2);
0028
0029
0030 matice=materials(beamelem.matid).material;
0031 matpar=materials(end).constants;
0032
0033
0034 gravity=matpar.g;
0035 rho_ice=matpar.rho_ice;
0036 n=matice.n;
0037
0038
0039 [surface_param surface_is_present]=recover_input(inputs,'surface');
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/beamelem error message: slope_surface is missing in inputs');
0047 end
0048
0049
0050 z_list=getgriddata(beamelem,grids);
0051
0052
0053 surface_list=zeros(numgrids,1);
0054 B_list=zeros(numgrids,1);
0055 temperature_list=zeros(numgrids,1);
0056 slopesurface_list=zeros(numgrids,2);
0057
0058
0059 for i=1:numgrids,
0060 doflist=grids(beamelem.g(i)).grid.doflist;
0061 for j=1:NDOF2,
0062 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0063 slopesurface_list(i,j)=slopesurface_param(doflist(j));
0064 end
0065 dof=doflist(1);
0066 if(surface_is_present) surface_list(i)=surface_param(dof);end;
0067 if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0068 if(B_is_present) B_list(i)=B_param(dof);end;
0069 if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0070 end
0071
0072
0073 slope=zeros(2,1);
0074 slope(1)=slopesurface_list(1,1);
0075 slope(2)=slopesurface_list(1,2);
0076 slope2=slope(1)^2+slope(2)^2;
0077
0078
0079 if surface_is_present
0080 surface=surface_list(1);
0081 else
0082 surface=beamelem.s(1);
0083 end
0084
0085
0086 constant_part=-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2);
0087
0088
0089 num_gauss=3;
0090 [gauss_coords,gauss_weights]=GaussSegment(num_gauss);
0091
0092
0093 for ig=1:num_gauss,
0094
0095 gauss_weight=gauss_weights(ig);
0096 gauss_coord=gauss_coords(ig);
0097
0098
0099 if B_is_present,
0100 B=GetParameterValue(beamelem,B_list,gauss_coord);
0101 elseif temperature_is_present,
0102 temperature=GetParameterValue(beamelem,temperature_list,gauss_coord);
0103 B=paterson(temperature);
0104 else
0105 B=matice.B;
0106 end
0107
0108
0109 z_g=GetParameterValue(beamelem,z_list,gauss_coord);
0110
0111
0112 Jdet=GetJacobianDeterminant(beamelem,z_list,gauss_coord);
0113
0114
0115
0116 l1l2=GetNodalFunctions(beamelem,gauss_coord);
0117
0118
0119 pe_gaussian=zeros(numgrids*NDOF2,1);
0120
0121 for j=1:NDOF2,
0122 pe_gaussian(NDOF2+j)=constant_part*((surface-z_g)/B)^n*slope(j)*Jdet*gauss_weight;
0123 end
0124
0125
0126 for i=1:pe.nrows,
0127 pe.terms(i)=pe.terms(i)+pe_gaussian(i);
0128 end
0129 end
0130
0131
0132 if (beamelem.onbed==1),
0133
0134
0135 if thickness_is_present
0136 thickness=thickness_list(1);
0137 else
0138 thickness=beamelem.h(1);
0139 end
0140
0141
0142
0143 constant_part=-1.58*10^-10*rho_ice*gravity*thickness;
0144 ub=constant_part*slope(1);
0145 vb=constant_part*slope(2);
0146
0147
0148 pe.terms(1)=pe.terms(1)+ub;
0149 pe.terms(2)=pe.terms(2)+vb;
0150 end
0151 end