CreateKMatrix

PURPOSE ^

CREATEKMATRIX - create the stiffmess matrix for pentaelem

SYNOPSIS ^

function Ke=CreateKMatrix(pentaelem,grids,materials,inputs,analysis_type)

DESCRIPTION ^

CREATEKMATRIX - create the stiffmess matrix for pentaelem

   this stiffness matrix works for MacAyeal, Pattyn's  and Stokes' model ,
   thermal model (transient ans steady), prognostic, melting and 
   for the computation of the bed slope and surface slope

   Usage:
      Ke=CreateKMatrix(pentaelem,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(pentaelem,grids,materials,inputs,analysis_type)
0002 %CREATEKMATRIX - create the stiffmess matrix for pentaelem
0003 %
0004 %   this stiffness matrix works for MacAyeal, Pattyn's  and Stokes' model ,
0005 %   thermal model (transient ans steady), prognostic, melting and
0006 %   for the computation of the bed slope and surface slope
0007 %
0008 %   Usage:
0009 %      Ke=CreateKMatrix(pentaelem,grids,materials,inputs,analysis_type)
0010 %
0011 %   See also CREATEPVECTOR
0012 
0013 if strcmpi(analysis_type,'diagnostic_horiz'),
0014 
0015     Ke=CreateKMatrixHoriz(pentaelem,grids,materials,inputs,analysis_type);
0016 
0017 elseif strcmpi(analysis_type,'diagnostic_vert'),
0018 
0019     Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs);
0020 
0021 elseif strcmpi(analysis_type,'diagnostic_basevert'),
0022 
0023     Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs);
0024 
0025 elseif strcmpi(analysis_type,'prognostic'),
0026 
0027     Ke=CreateKMatrixPrognostic(pentaelem,grids,materials,inputs);
0028 
0029 elseif (strcmpi(analysis_type,'thermaltransient') | strcmpi(analysis_type,'thermalsteady')),
0030 
0031     Ke=CreateKMatrixThermal(pentaelem,grids,materials,inputs);
0032 
0033 elseif strcmpi(analysis_type,'melting')
0034 
0035     Ke=CreateKMatrixMelting(pentaelem,grids,materials,inputs);
0036 
0037 elseif strcmpi(analysis_type,'diagnostic_stokes'),
0038 
0039     Ke=CreateKMatrixStokes(pentaelem,grids,materials,inputs);
0040 
0041 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'),
0042 
0043     Ke=CreateKMatrixSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);
0044 
0045 
0046 end
0047 end %end function
0048 
0049 function Ke=CreateKMatrixSlopeCompute(pentaelem,grids,materials,inputs,analysis_type)
0050 
0051     %We are dealing with a collapsed formulation on the lower triangle of the penta,
0052     %only on the bedrock.
0053 
0054     if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0055         if pentaelem.onbed~=1,
0056             Ke=elemmatrix(0);
0057             return;
0058         end
0059     elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0060         if pentaelem.onsurface~=1,
0061             Ke=elemmatrix(0);
0062             return;
0063         end
0064     end
0065 
0066     %some variables
0067     numgrids=3;
0068     DOFPERGRID=1;
0069     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0070 
0071     %Create elementary stiffness matrix
0072     Ke=elemmatrix(numdof);
0073 
0074     %Get all element grid data:
0075     xyz_list=getgriddata(pentaelem,grids); 
0076 
0077     %Just keep the first 3 grids
0078     if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0079         xyz_list=xyz_list(1:3,:);
0080     elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0081         xyz_list=xyz_list(4:6,:);
0082     end    
0083 
0084     %Build linear indices for elementary stiffness matrix.
0085     for i=1:numgrids,
0086         if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
0087             doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0088         elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0089             doflist=grids(pentaelem.g(i+3)).grid.doflist; %list of dofs in the g-set: last three grids
0090         end
0091         dof=doflist(1);
0092         Ke.row_indices(i)=dof;
0093     end
0094 
0095     % Get gaussian points and weights.
0096     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0097 
0098     for ig=1:num_gauss2D,
0099 
0100         %Pick up the gaussian point and its weight:
0101         gauss_weight=gauss_weights(ig);
0102         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0103         
0104         %Get the Jacobian determinant
0105         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0106 
0107         %Get L matrix if viscous basal drag present:
0108         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0109    
0110         DL_scalar=gauss_weight*Jdettria;
0111 
0112         %Do the triple product tL*D*L.
0113         Ke_gg_gaussian=L'*DL_scalar*L;
0114         
0115         %Add Ke_gg_drag_gaussian to Ke
0116         Ke.terms=Ke.terms+Ke_gg_gaussian;
0117     end
0118 
0119 end %end function
0120 
0121 function Ke=CreateKMatrixStokes(pentaelem,grids,materials,inputs,analysis_type)
0122 
0123 %some variables
0124 numgrids=6;
0125 DOFPERGRID=4;
0126 numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0127 
0128 numgrids2d=3;
0129 numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.
0130 
0131 %some parameters
0132 MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
0133 MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT
0134 
0135 %Create elementary stiffness matrix
0136 Ke=elemmatrix(numdof);
0137 
0138 %recover material
0139 matice=materials(pentaelem.matid).material;
0140 matpar=materials(end).constants;
0141 
0142 %recover material parameters
0143 gravity=matpar.g;
0144 viscosity_overshoot=matpar.viscosity_overshoot;
0145 rho_ice=matpar.rho_ice;
0146 rho_water=matpar.rho_water;
0147 
0148 %recover conditioning number for the matrix
0149 conditioning=pentaelem.reconditioning_number;
0150 
0151 %recover extra inputs from users, at current convergence iteration.
0152 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0153 [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0154 [bed_param bed_is_present]=recover_input(inputs,'bed');
0155 [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0156 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0157 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0158 
0159 %initialize extra arrays if some inputs are present
0160 if velocity_is_present,
0161     vxvyvz_list=zeros(numgrids,3);
0162 end
0163 B_list=zeros(numgrids,1);
0164 bed_list=zeros(numgrids,1);
0165 thickness_list=zeros(numgrids,1);
0166 K_list=zeros(numgrids,1);
0167 temperature_list=zeros(numgrids,1);
0168 
0169 %Get all element grid data:
0170 xyz_list=getgriddata(pentaelem,grids);
0171 
0172 %Build linear indices for elementary stiffness matrix.
0173 for i=1:numgrids,
0174     doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0175     for j=1:DOFPERGRID,
0176         dof=doflist(j);
0177         Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
0178         if j<4, %the fourth dof corresponds to pressure and not velocity
0179             if velocity_is_present, vxvyvz_list(i,j)=velocity_param(dof); end;
0180         end
0181     end
0182     
0183     dof=doflist(1);
0184     if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0185     if(bed_is_present) bed_list(i)=bed_param(dof);end;
0186     if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0187     if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
0188     if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0189 end
0190 
0191 % Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
0192 %get tria gaussian points as well as segment gaussian points. For tria gaussian
0193 %points, order of integration is 2, because we need to integrate the product tB*D*B'
0194 %which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
0195 %points, same deal, which yields 3 gaussian points.
0196 area_order=5;
0197 num_vert_gauss=5;
0198 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0199 
0200 Ke_temp=zeros(27,27);
0201 
0202 %Start  looping on the number of gaussian points:
0203 for igarea=1:num_area_gauss,
0204     for igvert=1:num_vert_gauss,
0205         %Pick up the gaussian point and its weight:
0206         gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0207         gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0208 
0209         %Get strain rate, if velocity has been supplied:
0210         if velocity_is_present,
0211             epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0212         else
0213             epsilon=zeros(6,1);
0214         end
0215         %disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
0216         
0217         %Update material parameter that deals with ice rigidity. The update depends either on the
0218         %temperature (if we are running thermal, or transient), or on B (if we are running inverse
0219         %control methods on B). The latter prevails.
0220         
0221         %Update material if temperature is provided.
0222         if temperature_is_present,
0223             temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord);
0224             matice.B=paterson(temperature);
0225         end
0226 
0227         %Update material if flow law is specified. This will erase the previous change
0228         %on B when temperature is provided.
0229         if flow_law_is_present,
0230             B_param=GetParameterValueStokes(pentaelem,B_list,gauss_coord);
0231             matice.B=B_param; clear B_param.
0232         end
0233 
0234         %Get viscosity at last two iterations:
0235         viscosity=GetViscosity3d(matice,epsilon);
0236         %disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));
0237 
0238         %Get B and Bprime matrices:
0239         B=GetBStokes(pentaelem,xyz_list,gauss_coord);
0240         Bprime=GetBprimeStokes(pentaelem,xyz_list,gauss_coord);
0241 
0242         %Get Jacobian determinant:
0243         Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0244         %disp(sprintf('Jacobian determinant: %f\n', Jdet));
0245 
0246         % Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
0247         % onto this scalar matrix, so that we win some computational time: */
0248         D=gauss_weight*Jdet*diag([viscosity*ones(6,1);-conditioning*ones(2,1)]);
0249 
0250         % Do the triple product tB*D*Bprime:
0251         Ke_gg_gaussian=B'*D*Bprime;
0252 
0253         %Add Ke_gg_gaussian to Ke_temp:
0254         Ke_temp=Ke_temp+Ke_gg_gaussian;
0255 
0256     end %for igarea=1:num_area_gauss,
0257 end %for igvert=1:num_vert_gauss,
0258 
0259 %Deal with 2d friction at the bedrock interface
0260 if (pentaelem.onbed==1 & (pentaelem.shelf==0))
0261 
0262     %build alpha^2 of each of the 3 grids of the triangle
0263     if (pentaelem.friction_type==2),  %friction_type=2 implies viscous, friction_type=1 implies plastic
0264     
0265         %Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
0266         frictionparameters=struct();
0267         frictionparameters.element_type='3d';
0268         frictionparameters.rho_ice=rho_ice;
0269         frictionparameters.rho_water=rho_water;
0270         frictionparameters.g=gravity;    
0271         frictionparameters.p=pentaelem.p;
0272         frictionparameters.q=pentaelem.q;
0273 
0274         if velocity_is_present,
0275             frictionparameters.velocities=vxvyvz_list(1:3,:);
0276         else
0277             frictionparameters.velocities=zeros(3,3);
0278         end
0279         if thickness_is_present,
0280             frictionparameters.h=thickness_list(1:3);
0281         else
0282             frictionparameters.h=pentaelem.h(1:3);
0283         end
0284         if bed_is_present,
0285             frictionparameters.b=bed_list(1:3);
0286         else
0287             frictionparameters.b=pentaelem.b(1:3);
0288         end
0289         if basal_drag_is_present,
0290             frictionparameters.k=K_list(1:3);
0291         else
0292             frictionparameters.k=pentaelem.k(1:3);
0293         end
0294 
0295         alpha2=Getalpha2(frictionparameters);
0296     end
0297 
0298     % Get gaussian points and weights.
0299     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0300 
0301     for ig=1:num_gauss2D,
0302 
0303         %Pick up the gaussian point and its weight:
0304         gauss_weight=gauss_weights(ig);
0305         gauss_coord_tria=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0306         gauss_coord_penta=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig) -1];
0307 
0308         %Get the Jacobian determinant
0309         xyz_list_tria=xyz_list(1:3,:);
0310         Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord_tria);
0311 
0312         %Get L and Lprime matrix if viscous basal drag present:
0313         L=GetLStokes(triaelem,gauss_coord_tria);
0314         Lprime=GetLprimeStokes(triaelem,xyz_list,gauss_coord_tria,gauss_coord_penta);
0315 
0316         %Build normal vector to the surface pointing outward the glacier system
0317         n=zeros(3,1);
0318         normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(3,:)-xyz_list_tria(2,:)));
0319         n(1)=1/norm(normal)*normal(1);
0320         n(2)=1/norm(normal)*normal(2);
0321         n(3)=1/norm(normal)*normal(3);
0322 
0323         %Get strain rate, if velocity has been supplied:
0324         if velocity_is_present,
0325             epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord_penta);
0326         else
0327             epsilon=zeros(6,1);
0328         end
0329 
0330 
0331         %Update material if temperature is provided.
0332         if temperature_is_present,
0333             temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord_penta);
0334             matice.B=paterson(temperature);
0335         end
0336 
0337         %Get viscosity at last two iterations:
0338         viscosity=GetViscosity3d(matice,epsilon);
0339         %disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));
0340 
0341 
0342         %Get alpha2 on current gaussian point
0343         alpha2_g=GetParameterValue(triaelem,alpha2,gauss_coord_tria);
0344 
0345         %Create DL diagonal matrix
0346         DL=diag([alpha2_g*[1;1;-n(1)*n(3);-n(2)*n(3);-n(1)*n(3);-n(2)*n(3)]; -viscosity*[n(1);n(2);n(3);n(1)/2.0;n(2)/2.0];conditioning*[n(1);n(2);n(3)]])*gauss_weight*Jdettria;
0347 
0348         %Create stiffness matrix for drag
0349         Ke_gg_drag_gaussian=L'*DL*Lprime;
0350 
0351         %Add Ke_gg_drag_gaussian Ke
0352         Ke_temp(1:12,1:12)=Ke_temp(1:12,1:12)+Ke_gg_drag_gaussian;
0353     end
0354 end%ends friction
0355 
0356 %Reduce the matrix
0357 Ke_reduced=ReduceMatrixStokes(pentaelem,Ke_temp);
0358 
0359 %Add the components of the matrix for the 6 grids
0360 Ke.terms=Ke.terms+Ke_reduced;
0361 
0362 end%end function
0363 
0364 
0365 function Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs)
0366 
0367     %some variables
0368     numgrids=6;
0369     DOFPERGRID=1;
0370     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0371 
0372     numgrids2d=3;
0373     numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.
0374 
0375     %Create elementary stiffness matrix
0376     Ke=elemmatrix(numdof);
0377 
0378     %Get all element grid data:
0379     xyz_list=getgriddata(pentaelem,grids);
0380 
0381     %Build linear indices for elementary stiffness matrix.
0382     for i=1:numgrids,
0383         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0384         dof=doflist(3); %z degree of freedom
0385         Ke.row_indices(i)=dof;
0386     end
0387 
0388     % Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
0389     %get tria gaussian points as well as segment gaussian points. For tria gaussian
0390     %points, order of integration is 2, because we need to integrate the product tB*D*B'
0391     %which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
0392     %points, same deal, which yields 3 gaussian points.
0393     area_order=2;
0394     num_vert_gauss=2;
0395     [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0396 
0397     %Start  looping on the number of gaussian points:
0398     for igarea=1:num_area_gauss,
0399         for igvert=1:num_vert_gauss,
0400             %Pick up the gaussian point and its weight:
0401             gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0402             gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0403             
0404             %Get B and Bprime matrices:
0405             B=GetB_vert(pentaelem,xyz_list,gauss_coord);
0406             Bprime=GetBprime_vert(pentaelem,gauss_coord);
0407 
0408             %Get Jacobian determinant:
0409             Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0410             %disp(sprintf('Jacobian determinant: %f\n', Jdet));
0411 
0412             D=gauss_weight*Jdet;
0413 
0414             % Do the triple product tB*D*Bprime:
0415             Ke_gg_gaussian=B'*D*Bprime;
0416 
0417             %Add Ke_gg_gaussian to Ke:
0418             Ke.terms=Ke.terms+Ke_gg_gaussian;
0419             
0420         end %for ig=1:num_area_gauss,
0421       end %for ig=1:num_vert_gauss,
0422 
0423     %Deal with upper surface integration:
0424 
0425     if (pentaelem.onsurface==1),
0426         
0427         %Get the coordinates of the three grids
0428         xyz_list_tria=xyz_list(4:6,:);
0429 
0430         % Get gaussian points and weights.
0431         [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0432 
0433         %Build normal vector to the surface
0434         normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(2,:)-xyz_list_tria(3,:)));
0435         nz=1/norm(normal)*normal(3);
0436 
0437         for ig=1:num_gauss2D,
0438 
0439             %Pick up the gaussian point and its weight:
0440             gauss_weight=gauss_weights(ig);
0441             gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0442             
0443             %Get the Jacobian determinant
0444             Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0445 
0446             %Get L matrix if viscous basal drag present:
0447             L=GetL(triaelem,gauss_coord,DOFPERGRID);
0448        
0449             DL_scalar=gauss_weight*Jdettria*nz;
0450 
0451             %Do the triple product tL*D*L. Do not forget the '-' sign!
0452             Ke_gg_vert_gaussian=-L'*DL_scalar*L;
0453             
0454 
0455             %Add Ke_gg_drag_gaussian to Ke
0456             Ke.terms(4:6,4:6)=Ke.terms(4:6,4:6)+Ke_gg_vert_gaussian;
0457         end
0458     end%ends friction
0459 
0460 end %end function
0461 
0462 function Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs)
0463 
0464     %We are dealing with a collapsed formulation on the lower triangle of the penta,
0465     %only on the bedrock.
0466     if pentaelem.onbed~=1,
0467         Ke=elemmatrix(0);
0468         return;
0469     end
0470 
0471     %some variables
0472     numgrids=3;
0473     DOFPERGRID=1;
0474     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0475 
0476     %Create elementary stiffness matrix
0477     Ke=elemmatrix(numdof);
0478 
0479     %Get all element grid data:
0480     xyz_list=getgriddata(pentaelem,grids); 
0481 
0482     %Just keep the first 3 grids
0483     xyz_list=xyz_list(1:3,:);
0484 
0485     %Build linear indices for elementary stiffness matrix.
0486     for i=1:numgrids,
0487         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0488         dof=doflist(3); %z degree of freedom
0489         Ke.row_indices(i)=dof;
0490     end
0491 
0492     % Get gaussian points and weights.
0493     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0494 
0495     for ig=1:num_gauss2D,
0496 
0497         %Pick up the gaussian point and its weight:
0498         gauss_weight=gauss_weights(ig);
0499         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0500         
0501         %Get the Jacobian determinant
0502         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0503 
0504         %Get L matrix if viscous basal drag present:
0505         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0506    
0507         DL_scalar=gauss_weight*Jdettria;
0508 
0509         %Do the triple product tL*D*L.
0510         Ke_gg_vert_gaussian=L'*DL_scalar*L;
0511         
0512         %Add Ke_gg_drag_gaussian to Ke
0513         Ke.terms=Ke.terms+Ke_gg_vert_gaussian;
0514     end
0515 
0516 end %end function
0517 
0518 
0519 
0520 
0521 function Ke=CreateKMatrixThermal(pentaelem,grids,materials,inputs)
0522 
0523 
0524     %some variables
0525     numgrids=6;
0526     DOFPERGRID=1;
0527     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0528 
0529     %recover material
0530     matice=materials(pentaelem.matid).material;
0531     matpar=materials(end).constants;
0532 
0533     %recover material parameters
0534     gravity=matpar.g;
0535     rho_ice=matpar.rho_ice;
0536     rho_water=matpar.rho_water;
0537     heatcapacity=matpar.heatcapacity;
0538     thermalconductivity=matpar.thermalconductivity;
0539 
0540     %Create elementary stiffness matrix
0541     Ke=elemmatrix(numdof);
0542 
0543     %Get all element grid data:
0544     xyz_list=getgriddata(pentaelem,grids);
0545 
0546     %recover extra inputs from users, at current convergence iteration.
0547     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0548     [dt dt_is_present]=recover_input(inputs,'dt');
0549     if (~velocity_is_present )
0550         error('CreateKMatrixThermal error message: missing velocity parameter!');
0551     end
0552     if (~dt_is_present & ~pentaelem.thermal_steadystate)
0553         error('CreateKMatrixThermal error message: transient thermal solution requires present of time step as parameter!');
0554     end
0555 
0556     %initialize vxvyvz_list
0557     vxvyvz_list=zeros(numgrids,3);
0558 
0559     %Build linear indices for elementary stiffness matrix.
0560     for i=1:numgrids,
0561 
0562         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0563         for j=1:3,
0564             dof=doflist(j);
0565             vxvyvz_list(i,j)=velocity_param(dof);
0566         end    
0567         dof=doflist(1);
0568         Ke.row_indices(i)=dof;
0569         
0570     end
0571 
0572     % Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
0573     %get tria gaussian points as well as segment gaussian points. For tria gaussian
0574     %points, order of integration is 2, because we need to integrate the product tB*D*B'
0575     %which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
0576     %points, same deal, which yields 3 gaussian points.
0577     area_order=2;
0578     num_vert_gauss=2;
0579     [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0580 
0581     %Start  looping on the number of gaussian points:
0582     for igarea=1:num_area_gauss,
0583         for igvert=1:num_vert_gauss,
0584             %Pick up the gaussian point and its weight:
0585             gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0586             gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0587 
0588             %Get Jacobian determinant:
0589             Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0590             %disp(sprintf('Jacobian determinant: %f\n', Jdet));
0591                 
0592             %conduction
0593 
0594             %Get B and Bprime matrices:
0595             Bconduct=GetB_conduct(pentaelem,xyz_list,gauss_coord);
0596 
0597             % Build the D matrix
0598             D=gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity))*eye(3);
0599 
0600             % Do the triple product tB*D*B:
0601             if pentaelem.thermal_steadystate,
0602                 Kconduct=Bconduct'*D*Bconduct;
0603             else
0604                 Kconduct=dt*Bconduct'*D*Bconduct;
0605             end
0606 
0607             %advection
0608 
0609             %Get B and Bprime matrices:
0610             Badvec=GetB_advec(pentaelem,gauss_coord);
0611             Bprimeadvec=GetBprime_advec(pentaelem,xyz_list,gauss_coord);
0612 
0613             % Build the D matrix
0614             vx=GetParameterValue(pentaelem,vxvyvz_list(:,1),gauss_coord);
0615             vy=GetParameterValue(pentaelem,vxvyvz_list(:,2),gauss_coord);
0616             vz=GetParameterValue(pentaelem,vxvyvz_list(:,3),gauss_coord);
0617 
0618             D=gauss_weight*Jdet*[vx,0,0; 0,vy,0; 0,0,vz];
0619 
0620             % Do the triple product tB*D*Bprime:
0621             if pentaelem.thermal_steadystate,
0622                 Kadvec=Badvec'*D*Bprimeadvec;
0623             else
0624                 Kadvec=dt*Badvec'*D*Bprimeadvec;
0625             end
0626 
0627             %transient
0628             if pentaelem.thermal_steadystate,
0629                 Ktransient=0;
0630             else
0631                 L=GetNodalFunctions(pentaelem,gauss_coord);
0632                 Ktransient=L'*L*gauss_weight*Jdet;
0633             end
0634 
0635             %Add all contributions to Ke:
0636             Ke.terms=Ke.terms+Kadvec+Kconduct+Ktransient;
0637             %Ke.terms=Ke.terms+Kconduct+Ktransient;
0638 
0639         end %for ig=1:vert_gauss,
0640     end %for ig=1:num_area_gauss,
0641             
0642     %Ice/ocean heat exchange flux on ice shelf base
0643     if (pentaelem.onbed & pentaelem.shelf),
0644 
0645         %recover material parameters
0646         mixed_layer_capacity=matpar.mixed_layer_capacity;
0647         thermal_exchange_velocity=matpar.thermal_exchange_velocity;
0648 
0649         % Get gaussian points and weights.
0650         [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0651 
0652         for ig=1:num_gauss2D,
0653 
0654             %Pick up the gaussian point and its weight:
0655             gauss_weight=gauss_weights(ig);
0656             gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0657             
0658             %Get the Jacobian determinant
0659             xyz_list_tria=xyz_list(1:3,:);
0660             Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0661             
0662             %Get nodal functions value
0663             l1l2l3=GetNodalFunctions(triaelem,gauss_coord);
0664 
0665             %Compute the elementary load vector
0666             if pentaelem.thermal_steadystate,
0667                 Ke_ocean_gaussian=l1l2l3'*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice)*l1l2l3;
0668             else
0669                 Ke_ocean_gaussian=l1l2l3'*dt*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice)*l1l2l3;
0670             end
0671 
0672             %Add Ke_gg_drag_gaussian to Ke
0673             Ke.terms(1:3,1:3)=Ke.terms(1:3,1:3)+Ke_ocean_gaussian;
0674         end
0675     end%ends ice/ocean exchange flux
0676 
0677 
0678 end %end function
0679 
0680 
0681 
0682 
0683 
0684 
0685 
0686 function Ke=CreateKMatrixHoriz(pentaelem,grids,materials,inputs,analysis_type);
0687 
0688 global element_debug  element_debugid
0689 
0690 %Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
0691 %bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
0692 %the stiffness matrix.
0693 if pentaelem.collapse & ~pentaelem.onbed
0694 
0695     Ke=elemmatrix(0);
0696     return;
0697 
0698 elseif pentaelem.collapse & pentaelem.onbed
0699         
0700 
0701     Ke=PentaCollapseIntoTriaMatrix(pentaelem, grids, materials, inputs, analysis_type);
0702 
0703     %That's it, element will get destroyed upon return, and we have the correct
0704     %collapsed stiffness matrix.
0705     return;
0706 
0707 else 
0708 
0709     %implement standard penta element.
0710 
0711     %some variables
0712     numgrids=6;
0713     DOFPERGRID=2;
0714     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0715 
0716     numgrids2d=3;
0717     numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.
0718 
0719 
0720     %some parameters
0721     MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
0722     MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT
0723 
0724 
0725     %Create elementary stiffness matrix
0726     Ke=elemmatrix(numdof);
0727 
0728     %recover material
0729     matice=materials(pentaelem.matid).material;
0730     matpar=materials(end).constants;
0731 
0732     %recover material parameters
0733     gravity=matpar.g;
0734     viscosity_overshoot=matpar.viscosity_overshoot;
0735     rho_ice=matpar.rho_ice;
0736     rho_water=matpar.rho_water;
0737 
0738     %recover extra inputs from users, at current convergence iteration.
0739     [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0740     [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0741     [surface_param surface_is_present]=recover_input(inputs,'surface');
0742     [bed_param bed_is_present]=recover_input(inputs,'bed');
0743     [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0744     [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0745     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0746     [oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');
0747 
0748     %initialize extra arrays if some inputs are present
0749     if velocity_is_present,
0750         vxvy_list=zeros(numgrids,2);
0751     end
0752     if oldvelocity_is_present,
0753         oldvxvy_list=zeros(numgrids,2);
0754     end
0755 
0756     B_list=zeros(numgrids,1);
0757     K_list=zeros(numgrids,1);
0758     thickness_list=zeros(numgrids,1);
0759     surface_list=zeros(numgrids,1);
0760     bed_list=zeros(numgrids,1);
0761     temperature_list=zeros(numgrids,1);
0762 
0763     %Get all element grid data:
0764     xyz_list=getgriddata(pentaelem,grids);
0765 
0766     %Build linear indices for elementary stiffness matrix.
0767     for i=1:numgrids,
0768         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0769         for j=1:DOFPERGRID,
0770             dof=doflist(j);
0771             Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
0772             if velocity_is_present, vxvy_list(i,j)=velocity_param(dof); end;
0773             if oldvelocity_is_present, oldvxvy_list(i,j)=oldvelocity_param(dof); end;
0774         end
0775         
0776         dof=doflist(1);
0777         if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0778         if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
0779         if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0780         if(surface_is_present) surface_list(i)=surface_param(dof);end;
0781         if(bed_is_present) bed_list(i)=bed_param(dof);end;
0782         if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0783 
0784     end
0785 
0786         % Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
0787     %get tria gaussian points as well as segment gaussian points. For tria gaussian
0788     %points, order of integration is 2, because we need to integrate the product tB*D*B'
0789     %which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
0790     %points, same deal, which yields 3 gaussian points.
0791     area_order=5;
0792     num_vert_gauss=5;
0793     [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0794 
0795     %Start  looping on the number of gaussian points:
0796     for igarea=1:num_area_gauss,
0797         for igvert=1:num_vert_gauss,
0798         %Pick up the gaussian point and its weight:
0799         gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0800         gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0801     
0802         %Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
0803         if thickness_is_present,
0804             thickness=GetParameterValue(pentaelem,thickness_list,gauss_coord);
0805         else
0806             thickness=GetParameterValue(pentaelem,pentaelem.h,gauss_coord);
0807         end
0808 
0809         %Get strain rate, if velocity has been supplied:
0810         if velocity_is_present,
0811             epsilon=GetStrainRate(pentaelem,vxvy_list,xyz_list,gauss_coord);
0812         else
0813             epsilon=zeros(5,1);
0814         end
0815         if oldvelocity_is_present,
0816             oldepsilon=GetStrainRate(pentaelem,oldvxvy_list,xyz_list,gauss_coord);
0817         else
0818             oldepsilon=zeros(5,1);
0819         end
0820 
0821         %disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
0822         
0823         %Update material parameter that deals with ice rigidity. The update depends either on the
0824         %temperature (if we are running thermal, or transient), or on B (if we are running inverse
0825         %control methods on B). The latter prevails.
0826 
0827         %Update material if temperature is provided.
0828         if temperature_is_present,
0829             temperature=GetParameterValue(pentaelem,temperature_list,gauss_coord);
0830             matice.B=paterson(temperature);
0831         end
0832 
0833         %Update material if flow law is specified. This will erase the previous change
0834         %on B when temperature is provided.
0835         if flow_law_is_present,
0836             B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
0837             matice.B=B_param; clear B_param.
0838         end
0839 
0840         %Get viscosity at last two iterations:
0841         newviscosity=GetViscosity3d(matice,epsilon);
0842         oldviscosity=GetViscosity3d(matice,oldepsilon);
0843         viscosity=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);
0844         %disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));
0845 
0846         %Get B and Bprime matrices:
0847         B=GetB_horiz(pentaelem,xyz_list,gauss_coord);
0848         Bprime=GetBprime_horiz(pentaelem,xyz_list,gauss_coord);
0849         
0850         %Get Jacobian determinant:
0851         Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0852         %disp(sprintf('Jacobian determinant: %f\n', Jdet));
0853 
0854         % Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
0855         % onto this scalar matrix, so that we win some computational time: */
0856         D=viscosity*gauss_weight*Jdet*diag(ones(5,1));
0857       
0858 
0859         % Do the triple product tB*D*Bprime:
0860         Ke_gg_gaussian=B'*D*Bprime;
0861 
0862         %Add Ke_gg_gaussian to Ke:
0863         Ke.terms=Ke.terms+Ke_gg_gaussian;
0864         
0865         end %for ig=1:num_area_gauss,
0866       end %for ig=1:num_vert_gauss,
0867 
0868         
0869       
0870 
0871     %Deal with 2d friction at the bedrock interface
0872 
0873     if (pentaelem.onbed==1 & (pentaelem.shelf==0))
0874 
0875         %build alpha^2 of each of the 3 grids of the triangle
0876         if (pentaelem.friction_type==2),  %friction_type=2 implies viscous, friction_type=1 implies plastic
0877         
0878             %Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
0879             frictionparameters=struct();
0880             frictionparameters.element_type='3d';
0881             frictionparameters.rho_ice=rho_ice;
0882             frictionparameters.rho_water=rho_water;
0883             frictionparameters.g=gravity;    
0884             frictionparameters.p=pentaelem.p;
0885             frictionparameters.q=pentaelem.q;
0886 
0887             if velocity_is_present,
0888                 frictionparameters.velocities=vxvy_list(1:3,:);
0889             else
0890                 frictionparameters.velocities=zeros(3,2);
0891             end
0892             if thickness_is_present,
0893                 frictionparameters.h=thickness_list(1:3);
0894             else
0895                 frictionparameters.h=pentaelem.h(1:3);
0896             end
0897             if bed_is_present,
0898                 frictionparameters.b=bed_list(1:3);
0899             else
0900                 frictionparameters.b=pentaelem.b(1:3);
0901             end
0902             if basal_drag_is_present,
0903                 frictionparameters.k=K_list(1:3);
0904             else
0905                 frictionparameters.k=pentaelem.k(1:3);
0906             end
0907 
0908             alpha2=Getalpha2(frictionparameters);
0909         end
0910 
0911         % Get gaussian points and weights.
0912         [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0913 
0914         for ig=1:num_gauss2D,
0915 
0916             %Pick up the gaussian point and its weight:
0917             gauss_weight=gauss_weights(ig);
0918             gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0919             gauss_coord_penta=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig) -1];
0920             
0921             %Get the Jacobian determinant
0922             xyz_list_tria=xyz_list(1:3,:);
0923             Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
0924 
0925             %Get L if viscous basal drag present:
0926             L=GetL(triaelem,gauss_coord,DOFPERGRID);
0927        
0928             %Get alpha2 on current gaussian point
0929             alpha2_g=GetParameterValue(triaelem,alpha2,gauss_coord);
0930 
0931             %Build D
0932             DL_scalar=alpha2_g*gauss_weight*Jdettria;
0933             DL=diag(DL_scalar*ones(2,1));
0934 
0935             %Do the triple product tL*D*L
0936             Ke_gg_drag_gaussian=L'*DL*L;
0937 
0938             %Add Ke_gg_drag_gaussian to Ke
0939             Ke.terms(1:6,1:6)=Ke.terms(1:6,1:6)+Ke_gg_drag_gaussian;
0940         end
0941     end%ends friction
0942 
0943     
0944 end %elseif pentaelem.collapse & pentaelem.onbed
0945 end%ends function
0946 
0947 function Ke=CreateKMatrixPrognostic(pentaelem,grids,materials,inputs)
0948 
0949     %We are dealing with a collapsed formulation on the lower triangle of the penta,
0950     %only on the bedrock.
0951     
0952     if pentaelem.onbed~=1,
0953         Ke=elemmatrix(0);
0954         return;
0955     end
0956 
0957     %some variables
0958     numgrids=3;
0959     DOFPERGRID=1;
0960     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
0961 
0962     %Create elementary stiffness matrix
0963     Ke=elemmatrix(numdof);
0964 
0965     %Get all element grid data:
0966     xyz_list=getgriddata(pentaelem,grids); 
0967 
0968     %Just keep the first 3 grids
0969     xyz_list=xyz_list(1:3,:);
0970 
0971     %recover extra inputs from users, at current convergence iteration.
0972     [dt dt_is_present]=recover_input(inputs,'dt');
0973     [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');
0974 
0975     %check on all parameters
0976     if ((~dt_is_present) | (~velocity_average_is_present)),
0977         error('CreateKMatrixPrognostic error message: missing input parameters!');
0978     end
0979     
0980     %initialize parameter lists (for our 3 grids)
0981     vxvy_list=zeros(numgrids,2);
0982 
0983     %Build linear indices for elementary stiffness matrix.
0984     for i=1:numgrids,
0985         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
0986         
0987         Ke.row_indices(i)=doflist(1);
0988     
0989         vxvy_list(i,:)=velocity_average_param(doflist(1:2));
0990     end
0991 
0992     %Create Artificial diffusivity if asked
0993     if pentaelem.artificial_diffusivity,
0994         %Get the Jacobian determinant
0995         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]);
0996 
0997         %Build K matrix (artificial diffusivity matrix)
0998         v_gauss=1/3*(vxvy_list(1,:)+vxvy_list(2,:)+vxvy_list(3,:));
0999         K=Jdettria^(1/2)/2*[abs(v_gauss(1)) 0 ; 0 abs(v_gauss(2))];
1000     end
1001 
1002     % Get gaussian points and weights.
1003     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
1004 
1005     for ig=1:num_gauss2D,
1006 
1007         %Pick up the gaussian point and its weight:
1008         gauss_weight=gauss_weights(ig);
1009         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
1010 
1011         %Get the Jacobian determinant
1012         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
1013 
1014         %Get L matrix if viscous basal drag present:
1015         L=GetL(triaelem,gauss_coord,DOFPERGRID);
1016 
1017         DL_scalar=gauss_weight*Jdettria;
1018 
1019         %Do the triple product tL*D*L.
1020         Ke_gg_gaussian=L'*DL_scalar*L;
1021 
1022         %Get B matrix
1023         B=GetB_prog(triaelem,xyz_list,gauss_coord);
1024         
1025         %Get Bprime matrix
1026         Bprime=GetBprime_prog(triaelem,xyz_list,gauss_coord);
1027 
1028         %Get u, v and their derivatives at gauss point
1029         u_g=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
1030         v_g=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
1031         du_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,1),xyz_list,gauss_coord);
1032         dv_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,2),xyz_list,gauss_coord);
1033         dux_g=du_g(1);
1034         dvy_g=dv_g(2);
1035         
1036         DL_scalar=dt*gauss_weight*Jdettria;
1037 
1038         %Create DL and DLprime matrix
1039         DL=DL_scalar*[dux_g, 0; 0, dvy_g];
1040         DLprime=DL_scalar*[u_g, 0; 0, v_g];
1041 
1042         %Do the triple product tL*D*L.
1043         Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;    
1044 
1045         %Add Ke_gg_drag_gaussian and Ke_gg_gaussian to Ke
1046         Ke.terms=Ke.terms+Ke_gg_gaussian+Ke_gg_thickness;
1047 
1048         %Add artificial diffusivity  if requested
1049         if pentaelem.artificial_diffusivity,
1050 
1051             %Do the triple product tL*D*L.
1052             Ke_gg_gaussian=Bprime'*DL_scalar*K*Bprime;
1053             
1054             %Add Ke_gg_drag_gaussian to Ke
1055             Ke.terms=Ke.terms+Ke_gg_gaussian;
1056         end
1057     end
1058 
1059 end %end function
1060 
1061 function Ke=CreateKMatrixMelting(pentaelem,grids,materials,inputs)
1062 
1063     %We are dealing with a collapsed formulation on the lower triangle of the penta,
1064     %only on the bedrock.
1065     if pentaelem.onbed~=1,
1066         Ke=elemmatrix(0);
1067         return;
1068     end
1069 
1070     %some variables
1071     numgrids=3;
1072     DOFPERGRID=1;
1073     numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.
1074 
1075     %Create elementary stiffness matrix
1076     Ke=elemmatrix(numdof);
1077 
1078     %recover material
1079     matpar=materials(end).constants;
1080 
1081     %recover material parameters
1082     latentheat=matpar.latentheat;
1083     heatcapacity=matpar.heatcapacity;
1084 
1085     %Get all element grid data:
1086     xyz_list=getgriddata(pentaelem,grids); 
1087 
1088     %Just keep the first 3 grids
1089     xyz_list=xyz_list(1:3,:);
1090 
1091     %Build linear indices for elementary stiffness matrix.
1092     for i=1:numgrids,
1093         doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
1094         dof=doflist(1); %melting is 1st degree of freedom
1095         Ke.row_indices(i)=dof;
1096     end
1097 
1098     % Get gaussian points and weights.
1099     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
1100 
1101     for ig=1:num_gauss2D,
1102 
1103         %Pick up the gaussian point and its weight:
1104         gauss_weight=gauss_weights(ig);
1105         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
1106         
1107         %Get the Jacobian determinant
1108         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
1109 
1110         %Get L matrix if viscous basal drag present:
1111         L=GetL(triaelem,gauss_coord,DOFPERGRID);
1112    
1113         DL_scalar=latentheat/heatcapacity*gauss_weight*Jdettria;
1114 
1115         %Do the triple product tL*D*L.
1116         Ke_gg_vert_gaussian=L'*DL_scalar*L;
1117         
1118         %Add Ke_gg_drag_gaussian to Ke
1119         Ke.terms=Ke.terms+Ke_gg_vert_gaussian;
1120     end
1121 
1122 end %end function

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