


CREATEPVECTOR - create the load vector for an acceleratedtriaelem
global finite element load vector of the whole model
works only for MaAyeal's model
Usage:
Pe=CreatePVector(element,grids,materials,inputs,analysis_type);
See also CREATEKMATRIX

0001 function pe=CreatePVector(element,grids,materials,inputs,analysis_type); 0002 %CREATEPVECTOR - create the load vector for an acceleratedtriaelem 0003 % 0004 % global finite element load vector of the whole model 0005 % works only for MaAyeal's model 0006 % 0007 % Usage: 0008 % Pe=CreatePVector(element,grids,materials,inputs,analysis_type); 0009 % 0010 % See also CREATEKMATRIX 0011 0012 global aire alpha beta 0013 0014 %some variables 0015 numgrids=3; 0016 NDOF2=2; 0017 numdof=numgrids*NDOF2; 0018 0019 %recover extra inputs from users, at current convergence iteration. 0020 [thickness_param thickness_is_present]=recover_input(inputs,'thickness'); 0021 [surface_param surface_is_present]=recover_input(inputs,'surface'); 0022 0023 0024 if thickness_is_present, 0025 thickness=thickness_param(1:6:6*element.nods); 0026 else 0027 thickness=element.thickness; 0028 end 0029 0030 if surface_is_present, 0031 surface=surface_param(1:6:6*element.nods); 0032 else 0033 surface=element.surface; 0034 end 0035 0036 0037 %retrieve rho_ice,g and thickness 0038 matpar=materials(end).constants; 0039 rho_ice=matpar.rho_ice; 0040 rho_water=matpar.rho_water; 0041 g=matpar.g; 0042 0043 %Build length_icefront and normal_icefront: 0044 %[length_icefront,normal_icefront]=buildicefrontnormal(element.x,element.y,element.index_icefront); 0045 0046 Rhs_x=zeros(element.nel*27,1); 0047 Rhs_y=zeros(element.nel*27,1); 0048 Rhs_y=ones(element.nel*27,1); 0049 row_rhs=zeros(element.nel*27,1); 0050 0051 count=-element.nel+1; 0052 for i=1:3 0053 for n=1:3 0054 for m=1:3 0055 count=count+element.nel; 0056 0057 Rhs_x(count:count+element.nel-1)= -rho_ice * g * ... 0058 thickness(element.index(:,n)).*surface(element.index(:,m)) ... 0059 .* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 ); 0060 0061 Rhs_y(count:count+element.nel-1)= -rho_ice * g * ... 0062 thickness(element.index(:,n)).*surface(element.index(:,m)) ... 0063 .* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 ); 0064 0065 row_rhs(count:count+element.nel-1)=element.index(:,i); 0066 0067 end 0068 end 0069 end 0070 0071 Rhs=full([sparse2(row_rhs,ones(element.nel*27,1),Rhs_x,element.nods,1) 0072 sparse2(row_rhs,ones(element.nel*27,1),Rhs_y,element.nods,1)]); 0073 0074 %Index icefront is skipped: 0075 %for k=1:2 0076 % for l=1:2 0077 % for j=1:2 0078 % Rhs(element.index_icefront(:,k))=Rhs(element.index_icefront(:,k)) + ... 0079 %(rho_ice*g/2*thickness(element.index_icefront(:,l)).*thickness(element.index_icefront(:,j)) ... 0080 %+ rho_water*g/2*(element.gridoniceshelf(element.index_icefront(:,k))) ... 0081 %.*(-min(0,element.bed(element.index_icefront(:,l))).*min(0,element.bed(element.index_icefront(:,j))) ... 0082 %+min(0,thickness(element.index_icefront(:,l))+element.bed(element.index_icefront(:,l)))... 0083 %.*min(0,element.thickness(element.index_icefront(:,j))+element.bed(element.index_icefront(:,j)))))... 0084 %.*normal_icefront(:,1).*length_icefront(:) ... 0085 %/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); 0086 % 0087 % Rhs(element.index_icefront(:,k)+element.nods)=Rhs(element.index_icefront(:,k)+element.nods) + ... 0088 %(rho_ice*g/2*thickness(element.index_icefront(:,l)).*thickness(element.index_icefront(:,j)) ... 0089 %+ rho_water*g/2*(element.gridoniceshelf(element.index_icefront(:,k))) ... 0090 %.*(-min(0,element.bed(element.index_icefront(:,l))).*min(0,element.bed(element.index_icefront(:,j))) ... 0091 %+min(0,thickness(element.index_icefront(:,l))+element.bed(element.index_icefront(:,l)))... 0092 %.*min(0,thickness(element.index_icefront(:,j))+element.bed(element.index_icefront(:,j)))))... 0093 %.*normal_icefront(:,2).*length_icefront(:) ... 0094 %/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); 0095 % 0096 % end 0097 % end 0098 %end 0099 % 0100 clear Rhs_x Rhs_y row_rhs 0101 0102 0103 %Create global vector. 0104 pe=globalvector; 0105 pe.nrows=numdof*element.nods; 0106 pe.vector=spalloc(6*element.nods,1,2*element.nods); 0107 pe.vector(1:6:6*element.nods)=Rhs(1:element.nods); 0108 pe.vector(2:6:6*element.nods)=Rhs(element.nods+1:2*element.nods);