CreateKMatrix

PURPOSE ^

CREATEKMATRIX - create the stiffmess matrix for triaelem

SYNOPSIS ^

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

DESCRIPTION ^

CREATEKMATRIX - create the stiffmess matrix for triaelem

   this stiffness matrix works for MacAyeal's model, prognostic
   and for the computation of the bed slope and surface slope

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

   See also CREATEPVECTOR

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function Ke=CreateKMatrix(triaelem,grids,materials,inputs,analysis_type);
0002 %CREATEKMATRIX - create the stiffmess matrix for triaelem
0003 %
0004 %   this stiffness matrix works for MacAyeal's model, prognostic
0005 %   and for the computation of the bed slope and surface slope
0006 %
0007 %   Usage:
0008 %      Ke=CreateKMatrix(triaelem,grids,materials,inputs,analysis_type)
0009 %
0010 %   See also CREATEPVECTOR
0011 
0012     if strcmpi(analysis_type,'diagnostic_horiz'),
0013 
0014         Ke=CreateKMatrixHoriz(triaelem,grids,materials,inputs);
0015 
0016     elseif strcmpi(analysis_type,'diagnostic_vert'),
0017 
0018         Ke=CreateKMatrixVert(triaelem,grids,materials,inputs);
0019 
0020     elseif strcmpi(analysis_type,'temperature'),
0021 
0022         Ke=CreateKMatrixThermal(triaelem,grids,materials,inputs);
0023 
0024     elseif strcmpi(analysis_type,'prognostic'),
0025 
0026         Ke=CreateKMatrixPrognostic(triaelem,grids,materials,inputs);
0027 
0028     elseif strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') | strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0029 
0030         Ke=CreateKMatrixSlopeCompute(triaelem,grids,materials,inputs);
0031 
0032     end
0033 end %end function
0034 
0035 
0036 function Ke=CreateKMatrixVert(triaelem,grids,materials,inputs)
0037 error('CreateKMatrix error message: Vertical analysis not implemented yet');
0038 end %end function
0039 
0040 
0041 function Ke=CreateKMatrixPrognostic(triaelem,grids,materials,inputs)
0042 
0043     %some variables
0044     numgrids=3;
0045     DOFPERGRID=1;
0046     numdof=numgrids*DOFPERGRID; %number of dof for element triaelem.
0047 
0048     %Create elementary stiffness matrix
0049     Ke=elemmatrix(numdof);
0050 
0051     %Get all element grid data:
0052     xyz_list=getgriddata(triaelem,grids); 
0053 
0054     %recover extra inputs from users, at current convergence iteration.
0055     [dt dt_is_present]=recover_input(inputs,'dt');
0056     [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
0057 
0058     %check on all parameters
0059     if ((~dt_is_present) | (~velocity_average_is_present)),
0060         error('CreateKMatrixPrognostic error message: missing input parameters!');
0061     end
0062     
0063     %initialize parameter lists (for our 3 grids)
0064     vxvy_list=zeros(numgrids,2);
0065 
0066     %Build linear indices for elementary stiffness matrix.
0067     for i=1:numgrids,
0068         doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
0069         Ke.row_indices(i)=doflist(1);
0070         vxvy_list(i,:)=velocity_average_param(doflist(1:2));
0071     end
0072 
0073     %Create Artificial diffusivity once for all if requested
0074     if triaelem.artificial_diffusivity,
0075         %Get the Jacobian determinant
0076         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]);
0077 
0078         %Build K matrix (artificial diffusivity matrix)
0079         v_gauss=1/3*(vxvy_list(1,:)+vxvy_list(2,:)+vxvy_list(3,:));
0080         K=Jdettria^(1/2)/2*[abs(v_gauss(1)) 0 ; 0 abs(v_gauss(2))];        
0081     end
0082 
0083     % Get gaussian points and weights.
0084     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0085 
0086     for ig=1:num_gauss2D,
0087 
0088         %Pick up the gaussian point and its weight:
0089         gauss_weight=gauss_weights(ig);
0090         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0091 
0092         %Get the Jacobian determinant
0093         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0094 
0095         %Get L matrix if viscous basal drag present:
0096         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0097 
0098         DL_scalar=gauss_weight*Jdettria;
0099 
0100         %Do the triple product tL*D*L.
0101         Ke_gg_gaussian=L'*DL_scalar*L;
0102 
0103         %Get B matrix
0104         B=GetB_prog(triaelem,xyz_list,gauss_coord);
0105         
0106         %Get Bprime matrix
0107         Bprime=GetBprime_prog(triaelem,xyz_list,gauss_coord);
0108 
0109         %Get u, v and their derivatives at gauss point
0110         u_g=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
0111         v_g=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
0112         du_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,1),xyz_list,gauss_coord);
0113         dv_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,2),xyz_list,gauss_coord);
0114         dux_g=du_g(1);
0115         dvy_g=dv_g(2);
0116         
0117         DL_scalar=dt*gauss_weight*Jdettria;
0118 
0119         %Create DL and DLprime matrix
0120         DL=DL_scalar*[dux_g, 0; 0, dvy_g];
0121         DLprime=DL_scalar*[u_g, 0; 0, v_g];
0122 
0123         %Do the triple product tL*D*L.
0124         Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
0125 
0126         %Add Ke_gg_drag_gaussian and Ke_gg_gaussian to Ke
0127         Ke.terms=Ke.terms+Ke_gg_gaussian+Ke_gg_thickness;
0128 
0129         %Add artificial diffusivity  if requested
0130         if triaelem.artificial_diffusivity,
0131 
0132             %Do the triple product tL*D*L.
0133             Ke_gg_gaussian=Bprime'*DL_scalar*K*Bprime;
0134             
0135             %Add Ke_gg_drag_gaussian to Ke
0136             Ke.terms=Ke.terms+Ke_gg_gaussian;
0137         end
0138     end
0139 end %end function
0140 
0141 
0142 function Ke=CreateKMatrixThermal(triaelem,grids,materials,inputs)
0143 error('CreateKMatrix error message: Thermal analysis not implemented yet');
0144 end %end function
0145 
0146 
0147 function Ke=CreateKMatrixSlopeCompute(triaelem,grids,materials,inputs)
0148 
0149     %some variables
0150     numgrids=3;
0151     DOFPERGRID=1;
0152     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0153 
0154     %Create elementary stiffness matrix
0155     Ke=elemmatrix(numdof);
0156 
0157     %Get all element grid data:
0158     xyz_list=getgriddata(triaelem,grids); 
0159 
0160     %Just keep the first 3 grids
0161     xyz_list=xyz_list(1:3,:);
0162 
0163     %Build linear indices for elementary stiffness matrix.
0164     for i=1:numgrids,
0165         doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
0166         dof=doflist(1);
0167         Ke.row_indices(i)=dof;
0168     end
0169 
0170     % Get gaussian points and weights.
0171     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0172 
0173     for ig=1:num_gauss2D,
0174 
0175         %Pick up the gaussian point and its weight:
0176         gauss_weight=gauss_weights(ig);
0177         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0178         
0179         %Get the Jacobian determinant
0180         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0181 
0182         %Get L matrix if viscous basal drag present:
0183         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0184    
0185         DL_scalar=gauss_weight*Jdettria;
0186 
0187         %Do the triple product tL*D*L.
0188         Ke_gg_gaussian=L'*DL_scalar*L;
0189         
0190         %Add Ke_gg_drag_gaussian to Ke
0191         Ke.terms=Ke.terms+Ke_gg_gaussian;
0192     end
0193 
0194 end %end function
0195 
0196 
0197 
0198 function Ke=CreateKMatrixHoriz(triaelem,grids,materials,inputs);
0199 
0200 global element_debug  element_debugid
0201 
0202 
0203 %Check if MacAyeal's acceleration is not on?
0204 if triaelem.acceleration==1,
0205     Ke=elemmatrix(0);
0206     return;
0207 end
0208 
0209 %some variables
0210 numgrids=3;
0211 DOFPERGRID=2;
0212 numdof=numgrids*DOFPERGRID; %number of dof for element triaelem.
0213 
0214 %some parameters
0215 MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
0216 MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT
0217 
0218 %Create elementary stiffness matrix
0219 Ke=elemmatrix(numdof);
0220 
0221 %recover material
0222 matice=materials(triaelem.matid).material;
0223 matpar=materials(end).constants;
0224 
0225 %recover material parameters
0226 gravity=matpar.g;
0227 viscosity_overshoot=matpar.viscosity_overshoot;
0228 rho_ice=matpar.rho_ice;
0229 rho_water=matpar.rho_water;
0230 
0231 %recover extra inputs from users, at current convergence iteration.
0232 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0233 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0234 [surface_param surface_is_present]=recover_input(inputs,'surface');
0235 [bed_param bed_is_present]=recover_input(inputs,'bed');
0236 [temperature_average_param temperature_average_is_present]=recover_input(inputs,'temperature_average');
0237 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0238 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0239 [oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');
0240 
0241 %initialize extra arrays if some inputs are present
0242 vxvy_list=zeros(numgrids,2);
0243 oldvxvy_list=zeros(numgrids,2);
0244 
0245 B_list=zeros(numgrids,1);
0246 K_list=zeros(numgrids,1);
0247 thickness_list=zeros(numgrids,1);
0248 surface_list=zeros(numgrids,1);
0249 bed_list=zeros(numgrids,1);
0250 temperature_average_list=zeros(numgrids,1);
0251 
0252 
0253 %Get all element grid data:
0254 xyz_list=getgriddata(triaelem,grids);
0255 
0256 
0257 %Build linear indices for elementary stiffness matrix.
0258 for i=1:numgrids,
0259     doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
0260     for j=1:DOFPERGRID,
0261         dof=doflist(j);
0262         Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
0263         if velocity_is_present, vxvy_list(i,j)=velocity_param(dof); end;
0264         if oldvelocity_is_present, oldvxvy_list(i,j)=oldvelocity_param(dof); end;
0265     end
0266     
0267     dof=doflist(1);
0268     if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0269     if(basal_drag_is_present), 
0270         K_list(i)=basal_drag_param(dof);
0271     else
0272         K_list(i)=triaelem.k(i);
0273     end
0274     if(thickness_is_present),
0275         thickness_list(i)=thickness_param(dof);
0276     else
0277         thickness_list(i)=triaelem.h(i);
0278     end
0279     if(surface_is_present),
0280         surface_list(i)=surface_param(dof);
0281     else
0282         surface_list(i)=triaelem.s(i);
0283     end
0284     if(bed_is_present),
0285         bed_list(i)=bed_param(dof);
0286     else
0287         bed_list(i)=triaelem.b(i);
0288     end
0289     if(temperature_average_is_present) temperature_average_list(i)=temperature_average_param(dof);end;
0290 end
0291 if (element_debug & triaelem.id==element_debugid),
0292     disp(sprintf('El id %i TriaElemnet input list before gaussian loop: \n',element_debugid)); 
0293     disp(sprintf('   rho_ice: %g ',rho_ice));
0294     disp(sprintf('   gravity: %g',gravity));
0295     disp(sprintf('   rho_water: %g',rho_water));
0296     disp(sprintf('   Velocity: '));
0297     for i=1:3,
0298         disp(sprintf('      grid %i  [%g,%g]',i,vxvy_list(i,1),vxvy_list(i,2)));
0299     end
0300     disp(sprintf('   B [%g %g %g ]',B_list(1),B_list(2),B_list(3)));
0301     disp(sprintf('   K [%g %g %g ]',K_list(1),K_list(2),K_list(3)));
0302     disp(sprintf('   thickness [%g %g %g]',thickness_list(1),thickness_list(2),thickness_list(3)));
0303     disp(sprintf('   surface [%g %g %g ]',surface_list(1),surface_list(2),surface_list(3)));
0304     disp(sprintf('   bed [%g %g %g]',bed_list(1),bed_list(2),bed_list(3)));
0305     if(temperature_average_is_present)disp(sprintf('   temperature_average [%g %g %g]',temperature_average_list(1),temperature_average_list(2),temperature_average_list(3)));end;
0306 end
0307 
0308 
0309 
0310 alpha2_list=zeros(3,1);
0311 %Build alpha2 used by drag staffness matrix
0312 if (~triaelem.shelf & triaelem.friction_type==2)
0313     
0314     %Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
0315     frictionparameters=struct();
0316     frictionparameters.element_type='2d';
0317     frictionparameters.rho_ice=rho_ice;
0318     frictionparameters.rho_water=rho_water;
0319     frictionparameters.g=gravity;    
0320     frictionparameters.p=triaelem.p;
0321     frictionparameters.q=triaelem.q;
0322 
0323      if velocity_is_present,
0324         frictionparameters.velocities=vxvy_list;
0325     else
0326         frictionparameters.velocities=zeros(3,2);
0327     end
0328     frictionparameters.h=thickness_list;
0329     frictionparameters.b=bed_list;
0330     frictionparameters.k=K_list;
0331 
0332     alpha2=Getalpha2(frictionparameters);
0333     if (element_debug & triaelem.id==element_debugid),
0334         disp(sprintf('   alpha2_list (%g %g %g )',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
0335     end
0336 end
0337 
0338 % Get gaussian points and weights. Order of integration is 2, because we need to integrate the product tB*D*B'
0339 %which is a polynomial of degree 3 (see GaussTria for more details)
0340 
0341 [num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);
0342     
0343 if (element_debug & triaelem.id==element_debugid),
0344     disp(sprintf('   gaussian points: '));
0345     for i=1:num_gauss,
0346         disp(sprintf('    %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
0347     end
0348 end
0349     
0350 %Start  looping on the number of gaussian points:
0351 for ig=1:num_gauss,
0352     %Pick up the gaussian point and its weight:
0353     gauss_weight=gauss_weights(ig);
0354     gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];
0355 
0356     %Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
0357     if thickness_is_present,
0358         thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);
0359     else
0360         thickness=GetParameterValue(triaelem,triaelem.h,gauss_l1l2l3);
0361     end
0362     %disp(sprintf('Thickness: %f\n', thickness));
0363 
0364     %If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
0365     %velocity should be = 0. To achieve this result, we set k to a very low value.
0366     slopevector=zeros(2,1);
0367     if(~triaelem.shelf),
0368     
0369         if surface_is_present,
0370             slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);
0371         else
0372             slopevector=GetParameterDerivativeValue(triaelem,triaelem.s,xyz_list,gauss_l1l2l3);
0373         end
0374         slope_magnitude=sqrt(slopevector(1)^2+slopevector(2)^2);
0375 
0376         if (slope_magnitude>MAXSLOPE),
0377             K=10^(-MOUNTAINKEXPONENT);
0378         end
0379     end
0380     
0381     %Get strain rate, if velocity has been supplied:
0382     if velocity_is_present,
0383          epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_l1l2l3);
0384     else
0385          epsilon=zeros(3,1);
0386     end
0387     if oldvelocity_is_present,
0388          oldepsilon=GetStrainRate(triaelem,oldvxvy_list,xyz_list,gauss_l1l2l3);
0389     else
0390          oldepsilon=zeros(3,1);
0391     end
0392 
0393     %disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
0394     
0395     %Update material if temperature is provided.
0396     if temperature_average_is_present,
0397         temperature_average=GetParameterValue(triaelem,temperature_average_list,gauss_l1l2l3);
0398         matice.B=paterson(temperature_average);
0399     end
0400 
0401     %Update material if flow law is specified. This will erase the previous change
0402     %on B when temperature is provided.
0403     if flow_law_is_present,
0404         B_param=GetParameterValue(triaelem,B_list,gauss_l1l2l3);
0405         matice.B=B_param; clear B_param.
0406     end
0407 
0408 
0409     %Get viscosity at last two iterations:
0410      newviscosity=GetViscosity2d(matice,epsilon);
0411      oldviscosity=GetViscosity2d(matice,oldepsilon);
0412     viscosity=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);
0413 
0414     %Get Jacobian determinant:
0415     Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3);
0416 
0417     % Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
0418     % onto this scalar matrix, so that we win some computational time: */
0419     D=viscosity*thickness*gauss_weight*Jdet*diag(ones(numgrids,1));
0420 
0421     if (element_debug & triaelem.id==element_debugid),
0422         disp(sprintf('   gaussian loop %i\n',ig));
0423         disp(sprintf('      thickness %g',thickness));
0424         disp(sprintf('      slope (%g,%g)',slopevector(1),slopevector(2)));
0425         disp(sprintf('      alpha2_list (%g,%g,%g)',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
0426         disp(sprintf('      epsilon (%g,%g,%g)',epsilon(1),epsilon(2),epsilon(3)));
0427         disp(sprintf('      Matice: '));
0428         matice
0429         disp(sprintf('\n      viscosity: %g ',viscosity));
0430         disp(sprintf('      jacobian: %g ',Jdet));
0431         disp(sprintf('      gauss_weight: %g ',gauss_weight));
0432     end
0433 
0434 
0435     %Get B and Bprime matrices:
0436     B=GetB(triaelem,xyz_list,gauss_l1l2l3);
0437     Bprime=GetBprime(triaelem,xyz_list,gauss_l1l2l3);
0438 
0439     %Get L matrix if viscous basal drag present:
0440     L=zeros(2,DOFPERGRID*numgrids);
0441     if (triaelem.friction_type==2 && triaelem.shelf==0),
0442         L=GetL(triaelem,gauss_l1l2l3,DOFPERGRID);
0443     end
0444 
0445     
0446     % Do the triple product tB*D*Bprime:
0447     Ke_gg_gaussian=B'*D*Bprime;
0448 
0449     %Add Ke_gg_gaussian to Ke:
0450     Ke.terms=Ke.terms+Ke_gg_gaussian;
0451 
0452     %Now, take care of the basal friction if there is any
0453     alpha2_g=0;
0454     if (~triaelem.shelf) & (triaelem.friction_type==2),
0455  
0456         %compute alpha2 for the current gaussian point
0457         alpha2_g=GetParameterValue(triaelem,alpha2,gauss_l1l2l3);
0458         
0459         if velocity_is_present
0460             DL_scalar=alpha2_g*gauss_weight*Jdet;
0461         else
0462             DL_scalar=0;
0463         end
0464         
0465         DL=diag(DL_scalar*ones(2,1));
0466 
0467         %Do the triple product tL*D*L
0468         Ke_gg_drag_gaussian=L'*DL*L;
0469 
0470         %Add Ke_gg_drag_gaussian to Ke:
0471         Ke.terms=Ke.terms+Ke_gg_drag_gaussian;
0472     end
0473 
0474     if (element_debug & triaelem.id==element_debugid),
0475         disp(sprintf('      alpha2 %g\n',alpha2_g));
0476         disp(sprintf('B:\n'));
0477         B
0478         disp(sprintf('Bprime:\n'));
0479         Bprime
0480         disp(sprintf('L:\n'));
0481         L
0482     end
0483 
0484 end %for ig=1:num_gauss,
0485     
0486 if (element_debug & triaelem.id==element_debugid),
0487     disp(sprintf('Ke_gg->terms:\n'));
0488     Ke.terms
0489     disp(sprintf('Ke_gg->row_indices:\n'));
0490     Ke.row_indices'
0491 end
0492 
0493 end %end function
0494

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