CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for a singelem

SYNOPSIS ^

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

DESCRIPTION ^

CREATEPVECTOR - create the load vector for a singelem

   works only for Hutter's model

   Usage:
      Pe=CreatePVector(singelem,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(singelem,grids,materials,inputs,analysis_type);
0002 %CREATEPVECTOR - create the load vector for a singelem
0003 %
0004 %   works only for Hutter's model
0005 %
0006 %   Usage:
0007 %      Pe=CreatePVector(singelem,grids,materials,inputs,analysis_type);
0008 %
0009 %   See also CREATEKMATRIX
0010 
0011 
0012     if strcmpi(analysis_type,'diagnostic_hutter'),
0013 
0014         pe=CreatePVectorHutter(singelem,grids,materials,inputs);
0015 
0016     else
0017         error('singelem/CreatePVector error message: analysis type not implemented yet')
0018     end
0019 end %end function
0020 
0021 function pe=CreatePVectorHutter(singelem,grids,materials,inputs);
0022 
0023     %some variables
0024     numgrids=1;
0025     NDOF2=2;
0026 
0027     %Create elementary vector.
0028     pe=elemvector(numgrids*NDOF2);
0029 
0030     %recover material
0031     matice=materials(singelem.matid).material;
0032     matpar=materials(end).constants;
0033 
0034     %recover material parameters
0035     gravity=matpar.g;
0036     rho_ice=matpar.rho_ice;
0037     n=matice.n;
0038 
0039     %recover extra inputs
0040     [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0041     [B_param B_is_present]=recover_input(inputs,'B');
0042     [temperature_param temperature_is_present]=recover_input(inputs,'temperature');
0043     [slopesurface_param slopesurface_is_present]=recover_input(inputs,'slopesurface');
0044 
0045     if ~slopesurface_is_present,
0046         error('CreatePVector/singelem error message: slope_surface is missing in inputs');
0047     end
0048 
0049     %Initialize inputs
0050     B_list=zeros(numgrids,1);
0051     temperature_list=zeros(numgrids,1);
0052     slopesurface_list=zeros(numgrids,2);
0053 
0054     %Build row indices for elementary vector.
0055     for i=1:numgrids,
0056         doflist=grids(singelem.g(i)).grid.doflist;
0057         for j=1:NDOF2,
0058             pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0059             slopesurface_list(i,j)=slopesurface_param(doflist(j));
0060         end
0061         dof=doflist(1);
0062         if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0063         if(B_is_present) B_list(i)=B_param(dof);end;
0064         if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0065     end
0066 
0067     %compute slope2 slopex and slopy
0068     slope=zeros(2,1);
0069     slope(1)=slopesurface_list(1,1); %2 and 1 have the same slopes (mesh has been verticaly extruded)
0070     slope(2)=slopesurface_list(1,2);
0071     slope2=slope(1)^2+slope(2)^2;
0072 
0073     %compute thickness
0074     if thickness_is_present
0075         thickness=thickness_list(1);
0076     else
0077         thickness=singelem.h(1);
0078     end
0079 
0080     %compute ub
0081     ub=-1.58*10^-10*rho_ice*gravity*thickness*slope(1);
0082     vb=-1.58*10^-10*rho_ice*gravity*thickness*slope(2);
0083 
0084     %compute constant_part
0085     constant_part=-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2);
0086 
0087     %Update material if temperature is provided.
0088     if B_is_present,
0089         B=GetParameterValue(singelem,B_list,gauss_coord);
0090     elseif temperature_is_present,
0091         temperature=GetParameterValue(singelem,temperature_list,gauss_coord);
0092         B=paterson(temperature);
0093     else
0094         B=matice.B;
0095     end
0096 
0097     pe.terms(1)=ub-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(1);
0098     pe.terms(2)=vb-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2)*thickness^(n)/(B^n*(n+1))*slope(2);
0099 
0100 end %end function

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