CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for an acceleratedtriaelem

SYNOPSIS ^

function pe=CreatePVector(element,grids,materials,inputs,analysis_type);

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003