CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for a beamelem

SYNOPSIS ^

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

DESCRIPTION ^

CREATEPVECTOR - create the load vector for a beamelem

   works only for Hutter's model

   Usage:
      Pe=CreatePVector(beamelem,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(beamelem,grids,materials,inputs,analysis_type);
0002 %CREATEPVECTOR - create the load vector for a beamelem
0003 %
0004 %   works only for Hutter's model
0005 %
0006 %   Usage:
0007 %      Pe=CreatePVector(beamelem,grids,materials,inputs,analysis_type);
0008 %
0009 %   See also CREATEKMATRIX
0010 
0011     if strcmpi(analysis_type,'diagnostic_hutter'),
0012 
0013         pe=CreatePVectorHutter(beamelem,grids,materials,inputs);
0014 
0015     else
0016         error('beamelem/CreatePVector error message: analysis type not implemented yet')
0017     end
0018 end %end function
0019 
0020 function pe=CreatePVectorHutter(beamelem,grids,materials,inputs);
0021 
0022     %some variables
0023     numgrids=2;
0024     NDOF2=2;
0025 
0026     %Create elementary vector.
0027     pe=elemvector(numgrids*NDOF2);
0028 
0029     %recover material
0030     matice=materials(beamelem.matid).material;
0031     matpar=materials(end).constants;
0032 
0033     %recover material parameters
0034     gravity=matpar.g;
0035     rho_ice=matpar.rho_ice;
0036     n=matice.n;
0037 
0038     %recover extra inputs
0039     [surface_param surface_is_present]=recover_input(inputs,'surface');
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/beamelem error message: slope_surface is missing in inputs');
0047     end
0048 
0049     %Get all element grid data:
0050     z_list=getgriddata(beamelem,grids);
0051     
0052     %Initialize inputs
0053     surface_list=zeros(numgrids,1);
0054     B_list=zeros(numgrids,1);
0055     temperature_list=zeros(numgrids,1);
0056     slopesurface_list=zeros(numgrids,2);
0057 
0058     %Build row indices for elementary vector.
0059     for i=1:numgrids,
0060         doflist=grids(beamelem.g(i)).grid.doflist;
0061         for j=1:NDOF2,
0062             pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0063             slopesurface_list(i,j)=slopesurface_param(doflist(j));
0064         end
0065         dof=doflist(1);
0066         if(surface_is_present) surface_list(i)=surface_param(dof);end;
0067         if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0068         if(B_is_present) B_list(i)=B_param(dof);end;
0069         if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
0070     end
0071 
0072     %compute slope2 slopex and slopy
0073     slope=zeros(2,1);
0074     slope(1)=slopesurface_list(1,1); %2 and 1 have the same slopes (mesh has been verticaly extruded)
0075     slope(2)=slopesurface_list(1,2);
0076     slope2=slope(1)^2+slope(2)^2;
0077 
0078     %compute surface
0079     if surface_is_present
0080         surface=surface_list(1);
0081     else
0082         surface=beamelem.s(1);
0083     end
0084 
0085     %compute constant_part
0086     constant_part=-2*(rho_ice*gravity)^n*(slope2)^((n-1)/2);
0087 
0088     % Get gaussian points and weights.
0089     num_gauss=3;
0090     [gauss_coords,gauss_weights]=GaussSegment(num_gauss);
0091 
0092     %Start  looping on the number of gaussian points:
0093     for ig=1:num_gauss,
0094         %Pick up the gaussian point and its weight:
0095         gauss_weight=gauss_weights(ig);
0096         gauss_coord=gauss_coords(ig);
0097     
0098         %Update material if temperature is provided.
0099         if B_is_present,
0100             B=GetParameterValue(beamelem,B_list,gauss_coord);
0101         elseif temperature_is_present,
0102             temperature=GetParameterValue(beamelem,temperature_list,gauss_coord);
0103             B=paterson(temperature);
0104         else
0105             B=matice.B;
0106         end
0107 
0108         %compute constant_part
0109         z_g=GetParameterValue(beamelem,z_list,gauss_coord);
0110 
0111         %Get Jacobian determinant:
0112         Jdet=GetJacobianDeterminant(beamelem,z_list,gauss_coord);
0113         %disp(sprintf('Jacobian determinant: %f\n', Jdet));
0114 
0115         %Get nodal functions
0116         l1l2=GetNodalFunctions(beamelem,gauss_coord);
0117 
0118         %Build pe_gaussian vector:
0119         pe_gaussian=zeros(numgrids*NDOF2,1);
0120 
0121         for j=1:NDOF2,
0122             pe_gaussian(NDOF2+j)=constant_part*((surface-z_g)/B)^n*slope(j)*Jdet*gauss_weight;
0123         end
0124 
0125         %Add pe_gaussian vector to pe:
0126         for i=1:pe.nrows,
0127             pe.terms(i)=pe.terms(i)+pe_gaussian(i);
0128         end
0129     end %for ig=1:num_area_gauss,
0130 
0131     %Deal with lower surface
0132     if (beamelem.onbed==1),
0133 
0134         %compute thickness
0135         if thickness_is_present
0136             thickness=thickness_list(1);
0137         else
0138             thickness=beamelem.h(1);
0139         end
0140 
0141         %compute ub
0142         %constant_part=0;%-1.58*10^-10*rho_ice*gravity*thickness;
0143         constant_part=-1.58*10^-10*rho_ice*gravity*thickness;
0144         ub=constant_part*slope(1);
0145         vb=constant_part*slope(2);
0146 
0147         %Add pe_gg_drag_gaussian to pe
0148         pe.terms(1)=pe.terms(1)+ub;
0149         pe.terms(2)=pe.terms(2)+vb;
0150     end
0151 end %end function

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