CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for triaelem

SYNOPSIS ^

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

DESCRIPTION ^

CREATEPVECTOR - create the load vector for triaelem

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

   Usage:
      pe=CreatePVector(triaelem,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(triaelem,grids,materials,inputs,analysis_type);
0002 %CREATEPVECTOR - create the load vector for triaelem
0003 %
0004 %   this load vector works for MacAyeal's model, prognostic
0005 %   and for the computation of the bed slope and surface slope
0006 %
0007 %   Usage:
0008 %      pe=CreatePVector(triaelem,grids,materials,inputs,analysis_type)
0009 %
0010 %   See also CREATEKMATRIX
0011 
0012 if strcmpi(analysis_type,'diagnostic_horiz'),
0013 
0014     pe=CreatePVectorHoriz(triaelem,grids,materials,inputs);
0015 
0016 elseif strcmpi(analysis_type,'diagnostic_vert'),
0017 
0018     pe=CreatePVectorVert(triaelem,grids,materials,inputs);
0019 
0020 elseif strcmpi(analysis_type,'temperature'),
0021 
0022     pe=CreatePVectorThermal(triaelem,grids,materials,inputs);
0023 
0024 elseif strcmpi(analysis_type,'prognostic'),
0025 
0026     pe=CreatePVectorPrognostic(triaelem,grids,materials,inputs);
0027 
0028 elseif strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') | strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0029 
0030     pe=CreatePVectorSlopeCompute(triaelem,grids,materials,inputs,analysis_type);
0031 
0032 end
0033 end %end function
0034 
0035 function pe=CreatePVectorVert(triaelem,grids,materials,inputs)
0036 error('CreatePVector error message: Vertical analysis not implemented yet');
0037 end %end function
0038 
0039 function pe=CreatePVectorThermal(triaelem,grids,materials,inputs)
0040 error('CreatePVector error message: Thermal analysis not implemented yet');
0041 end %end function
0042 
0043 function pe=CreatePVectorSlopeCompute(triaelem,grids,materials,inputs,analysis_type);
0044 
0045     %some variables
0046     numgrids=3;
0047     DOFPERGRID=1;
0048     numdof=numgrids*DOFPERGRID; %number of dof for element triaelem.
0049 
0050     %Create elementary stiffness matrix
0051     pe=elemvector(numdof);
0052 
0053     %Get all element grid data:
0054     xyz_list=getgriddata(triaelem,grids); 
0055 
0056     %Just keep the first 3 grids
0057     xyz_list=xyz_list(1:3,:);
0058 
0059     %recover extra inputs
0060     if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') ,
0061         [param param_is_present]=recover_input(inputs,'bed');
0062     elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0063         [param param_is_present]=recover_input(inputs,'surface');
0064     end
0065 
0066     param_list=zeros(numgrids,1);
0067 
0068     %Build linear indices for elementary stiffness matrix.
0069     for i=1:numgrids,
0070         doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
0071         if(param_is_present) param_list(i)=param(doflist(1));end;
0072         dof=doflist(1);
0073         pe.row_indices(i)=dof;
0074     end
0075 
0076     if param_is_present, 
0077         param=param_list(1:3); 
0078     else 
0079         if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') ,
0080             param=triaelem.b(1:3); 
0081         elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
0082             param=triaelem.s(1:3); 
0083         end
0084     end
0085 
0086     % Get gaussian points and weights.
0087     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0088 
0089     %for icesheets
0090     for ig=1:num_gauss2D,
0091     
0092         %Pick up the gaussian point and its weight:
0093         gauss_weight=gauss_weights(ig);
0094         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0095         
0096         %Get bed slope
0097         dparam=GetParameterDerivativeValue(triaelem,param,xyz_list,gauss_coord);
0098         dparamdx=dparam(1);
0099         dparamdy=dparam(2);
0100 
0101         %Get the Jacobian determinant
0102         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0103 
0104         %Get L:
0105         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0106 
0107         %Build gaussian vector
0108         if strcmp(analysis_type,'bed_slope_compute_x') | strcmp(analysis_type,'surface_slope_compute_x'),
0109             pe_g_gaussian=Jdettria*gauss_weight*dparamdx*L';
0110         elseif strcmp(analysis_type,'bed_slope_compute_y')|  strcmp(analysis_type,'surface_slope_compute_y'),
0111             pe_g_gaussian=Jdettria*gauss_weight*dparamdy*L';
0112         end
0113 
0114         %Add Ke_gg_drag_gaussian to Ke
0115         pe.terms=pe.terms+pe_g_gaussian;
0116     end
0117 
0118 end %end function
0119 
0120 
0121 function pe=CreatePVectorHoriz(triaelem,grids,materials,inputs);
0122 
0123     global element_debug  element_debugid
0124 
0125     %Check if MacAyeal's acceleration is not on?
0126     if triaelem.acceleration==1,
0127         pe={};
0128     end
0129 
0130     %some variables
0131     numgrids=3;
0132     NDOF2=2;
0133 
0134     %Create elementary vector.
0135     pe=elemvector(numgrids*NDOF2);
0136 
0137     %recover material
0138     matice=materials(triaelem.matid).material;
0139     matpar=materials(end).constants;
0140 
0141     %recover material parameters
0142     gravity=matpar.g;
0143     rho_ice=matpar.rho_ice;
0144 
0145     %recover extra inputs
0146     [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0147     [surface_param surface_is_present]=recover_input(inputs,'surface');
0148     [bed_param bed_is_present]=recover_input(inputs,'bed');
0149     [basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
0150         
0151     %Get all element grid data:
0152     xyz_list=getgriddata(triaelem,grids);
0153 
0154     %initialize extra inputs
0155     K_list=zeros(numgrids,1);
0156     thickness_list=zeros(numgrids,1);
0157     surface_list=zeros(numgrids,1);
0158     bed_list=zeros(numgrids,1);
0159 
0160     %Build row indices for elementary vector.
0161     for i=1:numgrids,
0162         doflist=grids(triaelem.g(i)).grid.doflist;
0163         for j=1:NDOF2,
0164             pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0165         end
0166         dof=doflist(1);
0167         if(thickness_is_present),
0168             thickness_list(i)=thickness_param(dof);
0169         else
0170             thickness_list(i)=triaelem.h(i);
0171         end
0172         if(surface_is_present),
0173             surface_list(i)=surface_param(dof);
0174         else
0175             surface_list(i)=triaelem.s(i);
0176         end;
0177         if(bed_is_present),
0178             bed_list(i)=bed_param(dof);
0179         else
0180             bed_list(i)=triaelem.b(i);
0181         end
0182         if(basal_drag_is_present),
0183             K_list(i)=basal_drag_param(dof);
0184         else
0185             K_list(i)=triaelem.k(i);
0186         end
0187     end
0188 
0189     if (element_debug & triaelem.id==element_debugid),
0190         disp(sprintf('gravity %g\n',gravity));
0191         disp(sprintf('rho_ice %g\n',rho_ice));
0192         disp(sprintf('thickness_list (%g,%g,%g)\n',thickness_list(1),thickness_list(2),thickness_list(3)));
0193         disp(sprintf('surface_list (%g,%g,%g)\n',surface_list(1),surface_list(2),surface_list(3)));
0194         disp(sprintf('bed_list (%g,%g,%g)\n',bed_list(1),bed_list(2),bed_list(3)));
0195         disp(sprintf('K_list (%g,%g,%g)\n',K_list(1),K_list(2),K_list(3)));
0196     end
0197 
0198     % Get gaussian points and weights
0199     [num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);
0200 
0201     if (element_debug & triaelem.id==element_debugid),
0202         disp(sprintf('   gaussian points: '));
0203         for i=1:num_gauss,
0204             disp(sprintf('    %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
0205         end
0206     end
0207  
0208     %Start  looping on the number of gaussian points:
0209     for ig=1:num_gauss,
0210         %Pick up the gaussian point:
0211         gauss_weight=gauss_weights(ig);
0212         gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];
0213 
0214         %Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
0215         thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);
0216 
0217         slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);
0218 
0219         %In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the
0220         %element itself:
0221         if (~triaelem.shelf) & (triaelem.friction_type==1),
0222             plastic_stress=GetParameterValue(triaelem,K_list,gauss_l1l2l3);
0223         end
0224 
0225         %Get Jacobian determinant:
0226         Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3);
0227 
0228         %Get nodal functions
0229         l1l2l3=GetNodalFunctions(triaelem,gauss_l1l2l3);
0230 
0231         %Compute driving stress:
0232         driving_stress_baseline=rho_ice*gravity*thickness;
0233     
0234         if (element_debug & triaelem.id==element_debugid),
0235             disp(sprintf('      gaussian %i\n',ig));
0236             disp(sprintf('      thickness %g\n',thickness));
0237             disp(sprintf('      slope(%g,%g)\n',slopevector(1),slopevector(2)));
0238             disp(sprintf('      Jdet %g\n',Jdet));
0239             disp(sprintf('      gaussweigth %g\n',gauss_weight));
0240             disp(sprintf('      l1l2l3 (%g,%g,%g)\n',l1l2l3(1),l1l2l3(2),l1l2l3(3)));
0241             if(triaelem.friction_type==1), disp(sprintf('      plastic_stress(%g)\n',plastic_stress));end;
0242         end
0243 
0244         %Build pe_gaussian vector:
0245         pe_gaussian=zeros(numgrids*NDOF2,1);
0246 
0247         if(~triaelem.shelf) & (triaelem.friction_type==1),
0248             for i=1:numgrids,
0249                 for j=1:NDOF2,
0250                     pe_gaussian((i-1)*NDOF2+j)=(-driving_stress_baseline*slopevector(j)-plastic_stress)*Jdet*gauss_weight*l1l2l3(i); 
0251                 end
0252             end
0253         else
0254             for i=1:numgrids,
0255                 for j=1:NDOF2,
0256                     pe_gaussian((i-1)*NDOF2+j)=-driving_stress_baseline*slopevector(j)*Jdet*gauss_weight*l1l2l3(i);
0257                 end
0258             end
0259         end
0260 
0261         %*Add pe_gaussian vector to pe:
0262         for i=1:pe.nrows,
0263             pe.terms(i)=pe.terms(i)+pe_gaussian(i);
0264         end
0265 
0266     end %for ig=1:num_gauss,
0267         
0268     if (element_debug & triaelem.id==element_debugid),
0269         disp(sprintf('      pe_g->terms\n',ig));
0270         pe.terms
0271         disp(sprintf('      pe_g->row_indices\n',ig));
0272         pe.row_indices
0273     end
0274 
0275 end %end function
0276 
0277 function pe=CreatePVectorPrognostic(triaelem,grids,materials,inputs);
0278 
0279     %some variables
0280     numgrids=3;
0281     DOFPERGRID=1;
0282     numdof=numgrids*DOFPERGRID; %number of dof for element triaelem.
0283 
0284     %Create elementary stiffness matrix
0285     pe=elemvector(numdof);
0286 
0287     %Get all element grid data:
0288     xyz_list=getgriddata(triaelem,grids); 
0289 
0290     %recover extra inputs from users, at current convergence iteration.
0291     [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation');
0292     [melting_param melting_is_present]=recover_input(inputs,'melting');
0293     [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0294     [dt dt_is_present]=recover_input(inputs,'dt');
0295 
0296     %check on all parameters
0297     if ( ~accumulation_is_present | ~melting_is_present | ~thickness_is_present | ~dt_is_present),
0298         error('CreateKMatrixPrognostic error message: missing input parameters!');
0299     end
0300     
0301     %initialize parameter lists (for our 3 grids)
0302     thickness_list=zeros(numgrids,1);
0303     accumulation_list=zeros(numgrids,1);
0304     melting_list=zeros(numgrids,1);
0305 
0306     %Build linear indices for elementary stiffness matrix.
0307     for i=1:numgrids,
0308         doflist=grids(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
0309         dof=doflist(1); %first degree of freedom
0310         
0311         pe.row_indices(i)=dof;
0312     
0313         %recover thickness_list on first dof of every grid.
0314         thickness_list(i)=thickness_param(dof);
0315         melting_list(i)=melting_param(dof);
0316         accumulation_list(i)=accumulation_param(dof);
0317     end
0318 
0319     % Get gaussian points and weights.
0320     [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);
0321 
0322     for ig=1:num_gauss2D,
0323     
0324         %Pick up the gaussian point and its weight:
0325         gauss_weight=gauss_weights(ig);
0326         gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
0327         
0328         %Get melting at gaussian point
0329         melting=GetParameterValue(triaelem,melting_list,gauss_coord);
0330     
0331         %Get accumulation at gaussian point
0332         accumulation=GetParameterValue(triaelem,accumulation_list,gauss_coord);
0333 
0334         %Get thickness at gaussian point
0335         thickness=GetParameterValue(triaelem,thickness_list,gauss_coord);
0336 
0337         %Get the Jacobian determinant
0338         Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);
0339 
0340         %Get L:
0341         L=GetL(triaelem,gauss_coord,DOFPERGRID);
0342 
0343         %Build gaussian vector
0344         pe_g_gaussian=Jdettria*gauss_weight*(thickness+dt*(accumulation-melting))*L';
0345 
0346         %Add Ke_gg_drag_gaussian to Ke
0347         pe.terms=pe.terms+pe_g_gaussian;
0348     end
0349 
0350 end %end function

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