CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for pentaelem

SYNOPSIS ^

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

DESCRIPTION ^

CREATEPVECTOR - create the load vector for pentaelem

   this load vector 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:
      pe=CreatePVector(pentaelem,grids,materials,inputs,analysis_type)

   See also CREATEKMATRIX

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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