CreateKMatrix

PURPOSE ^

CREATEKMATRIX - create the stiffness matrix for acceleratedtriaelem

SYNOPSIS ^

function Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type);

DESCRIPTION ^

CREATEKMATRIX - create the stiffness matrix for acceleratedtriaelem

   global finite element stiffness matrix of the whole model
   works only for MacAyeal's model

   Usage:
      Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type)

   See also CREATEPVECTOR

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type);
0002 %CREATEKMATRIX - create the stiffness matrix for acceleratedtriaelem
0003 %
0004 %   global finite element stiffness matrix of the whole model
0005 %   works only for MacAyeal's model
0006 %
0007 %   Usage:
0008 %      Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type)
0009 %
0010 %   See also CREATEPVECTOR
0011 
0012 
0013 %keep permanent pieces of stiffness matrix in memory
0014 global permanent_pieces_of_Ke alpha beta row_location col_location aire
0015 
0016 %Ok, this is the accelerated MacAyeal's element. Very particular methodology here!
0017 
0018 %Some variables
0019 numdof=3*2; %3 grids, 2 degrees of freedom per grid.
0020 
0021 %recover extra inputs from users, at current convergence iteration.
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     %Start building areas
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); % if index is sorted from its original value, then aire could be negative
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     %  Do once and for all the initial computation of matrix-locations:
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     % Create permanent part of stiffness matrix.
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 %if isempty(permanent_pieces_of_Ke)
0144     
0145 %Compute viscosity
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 %Multiply viscosity by permanent_pieces_of_Ke
0152 count=-element.nel+1;
0153 
0154 %Create sparse global matrix terms:
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 %Create basal stiffness
0171 if (element.friction_type==2) & velocity_is_present
0172 
0173     %average of p and q over the grids (size nel->nods)
0174     pcoeff_grid=zeros(element.nods,1);
0175     qcoeff_grid=zeros(element.nods,1);
0176     for i=1:element.nods
0177         %1: find the elements that contain the grid i
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         %2 retrieve the value of p and q over each of these elements. The average is
0184         %plugged into dbx_grid
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     %Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
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 %Create final global stiffness matrix
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);

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