CreatePVector

PURPOSE ^

CREATEPVECTOR - create the load vector for an icefront

SYNOPSIS ^

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

DESCRIPTION ^

CREATEPVECTOR - create the load vector for an icefront

   works for the icefront elements, segments or quad elements
   this routine can be used for MacAyeal, Pattyn ans Stokes's models

   Usage:
      pe=CreatePVector(load,elements,grids,materials,inputs,analysis_type)
 
   See also PENALTYCREATEPVECTOR, PENALTYCREATEKMATRIX

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function pe=CreatePVector(load,elements,grids,materials,inputs,analysis_type);
0002 %CREATEPVECTOR - create the load vector for an icefront
0003 %
0004 %   works for the icefront elements, segments or quad elements
0005 %   this routine can be used for MacAyeal, Pattyn ans Stokes's models
0006 %
0007 %   Usage:
0008 %      pe=CreatePVector(load,elements,grids,materials,inputs,analysis_type)
0009 %
0010 %   See also PENALTYCREATEPVECTOR, PENALTYCREATEKMATRIX
0011 
0012 if strcmpi(analysis_type,'diagnostic_horiz'),
0013 
0014     pe=CreatePVectorHoriz(load,elements,grids,materials,inputs);
0015 
0016 elseif strcmpi(analysis_type,'diagnostic_stokes'),
0017 
0018     pe=CreatePVectorStokes(load,elements,grids,materials,inputs);
0019 
0020 else
0021     error('CreatePVector error message: only horizontal diagnostic should use icefront objects');
0022 
0023 end
0024 end %end function
0025 
0026 function pe=CreatePVectorStokes(load,elements,grids,materials,inputs);
0027 
0028 %some variables
0029 numgrids=4;
0030 NDOF4=4;
0031 
0032 numgridstria=4; %number of grids for the triangles
0033 NDOF4=4;
0034 
0035 
0036 %Create elementary vector.
0037 pe=elemvector(numgrids*NDOF4);
0038 
0039 %recover element object
0040 pentaelem=elements(load.eid).element;
0041 matpar=materials(end).constants;
0042 
0043 %recover material parameters
0044 gravity=matpar.g;
0045 rho_ice=matpar.rho_ice;
0046 rho_water=matpar.rho_water;
0047 
0048 %recover extra inputs
0049 [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0050 [bed_param bed_is_present]=recover_input(inputs,'bed');
0051 
0052 %initialize extra inputs
0053 thickness_list=zeros(6,1);
0054 bed_list=zeros(6,1);
0055 
0056 %Get all element grid data:
0057 xyz_list=getgriddata(pentaelem,grids);
0058 
0059 %Identify which grids are comprised in the segment:
0060 for i=1:6,
0061     if pentaelem.g(i)==load.g(1),
0062         grid1=i;
0063     end
0064     if pentaelem.g(i)==load.g(2),
0065         grid2=i;
0066     end
0067     if pentaelem.g(i)==load.g(3),
0068         grid3=i;
0069     end
0070     if pentaelem.g(i)==load.g(4),
0071         grid4=i;
0072     end
0073 end
0074 quad_grids=[grid1 grid2 grid3 grid4];
0075 
0076 %Build row indices for elementary vector.
0077 for i=1:numgrids,
0078     doflist=grids(load.g(i)).grid.doflist;
0079     for j=1:NDOF4,
0080         pe.row_indices((i-1)*NDOF4+j)=doflist(j);
0081     end
0082     dof=doflist(1);
0083 end
0084 
0085 for i=1:6,
0086     doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
0087     if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0088     if(bed_is_present) bed_list(i)=bed_param(dof);end;
0089 end
0090 
0091 
0092 %Recover thicknesses and bed at all grids
0093 if thickness_is_present,
0094     thickness=thickness_list(quad_grids);
0095     bed=bed_list(quad_grids);
0096 else
0097     thickness=pentaelem.h(quad_grids);
0098     bed=pentaelem.b(quad_grids);
0099 end
0100 
0101 %build xyz list for [grid1, grid2,grid3,grid3]
0102 xyz_list_quad=xyz_list(quad_grids,:);
0103 
0104 %Create a new grid in the midle of the quad and add it at the end of the list
0105 x5=mean(xyz_list_quad(:,1));
0106 y5=mean(xyz_list_quad(:,2));
0107 z5=mean(xyz_list_quad(:,3));
0108 
0109 xyz_list_quad=[xyz_list_quad; x5 y5 z5];
0110 
0111 %Compute four normals (the quad is divided into four triangles)
0112 V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
0113 normal1=1/norm(V1)*V1';
0114 
0115 V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
0116 normal2=1/norm(V2)*V2';
0117 
0118 V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
0119 normal3=1/norm(V3)*V3';
0120 
0121 V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
0122 normal4=1/norm(V4)*V4';
0123 
0124 %Create a matrix with the four normals
0125 normal=[normal1 normal2 normal3 normal4];
0126 
0127 %Compute load contribution for this segment:
0128 if pentaelem.shelf==1,
0129     %icefront ends in the water
0130     pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
0131 else
0132     %icefront ends in the air
0133     pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
0134 end
0135 
0136 end %End of function
0137 
0138 function pe=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);
0139 
0140 %The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
0141 %
0142 %   grid 4 +-----------------+ grid 3
0143 %          |\2            1 /|
0144 %          |1\    tria3    /2|
0145 %          |  \           /  |
0146 %          |   \         /   |
0147 %          |    \       /    |
0148 %          |     \     /     |
0149 %          |      \ 3 /      |
0150 %          |tria4  \ / 3     |
0151 %          |      3 \grid5   |
0152 %          |       / \       |
0153 %          |      / 3 \ tria2|
0154 %          |     /     \     |
0155 %          |    /       \    |
0156 %          |   /         \   |
0157 %          |  /   tria1   \  |
0158 %          |2/1           2\1|
0159 %   grid1  +-----------------+ grid 2
0160 %
0161 %
0162 
0163 %Build the four normal vectors
0164 nx=normal(1,:);
0165 ny=normal(2,:);
0166 nz=normal(3,:);
0167 
0168 % Get gaussian points and weights. order 2 since we have a product of 2 nodal functions (We use GaussTria since we will build two triangles)
0169 [num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);
0170 
0171 %Recover the surface of the four nodes
0172 s=thickness+bed;
0173 
0174 %Add surface sor the fifth point (average of the surfaces)
0175 s(5)=mean(s);
0176 
0177 %Recover grid coordinates
0178 x=xyz_list_quad(:,1);
0179 y=xyz_list_quad(:,2);
0180 z=xyz_list_quad(:,3);
0181 
0182 %Build triangles in a 2D plan before using reference elements
0183 
0184 %Create four triangle elements for each quad
0185 tria=triaelem;
0186 
0187 l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
0188 l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
0189 l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
0190 l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
0191 l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
0192 l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
0193 l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
0194 l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
0195 cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
0196 cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
0197 cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
0198 cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
0199 
0200 %First triangle : nodes 1, 2 and 5
0201 x1tria1=0;
0202 y1tria1=0;
0203 x2tria1=sqrt(l12);
0204 y2tria1=0;
0205 x3tria1=cos_theta_triangle1*sqrt(l15);
0206 y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);
0207 
0208 xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];
0209 
0210 %Second triangle : nodes 2, 3 and 5
0211 x1tria2=0;
0212 y1tria2=0;
0213 x2tria2=sqrt(l23);
0214 y2tria2=0;
0215 x3tria2=cos_theta_triangle2*sqrt(l25);
0216 y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);
0217 
0218 xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];
0219 
0220 %Third triangle : nodes 3, 4 and 5
0221 x1tria3=0;
0222 y1tria3=0;
0223 x2tria3=sqrt(l34);
0224 y2tria3=0;
0225 x3tria3=cos_theta_triangle3*sqrt(l35);
0226 y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);
0227 
0228 xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];
0229 
0230 %Fourth triangle : nodes 4, 1 and 5
0231 x1tria4=0;
0232 y1tria4=0;
0233 x2tria4=sqrt(l14);
0234 y2tria4=0;
0235 x3tria4=cos_theta_triangle4*sqrt(l45);
0236 y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);
0237 
0238 xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];
0239 
0240 %Create a list xyz_tria with coordinates of every triangle
0241 
0242 xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];
0243 
0244 %Initialize load vector
0245 pe=zeros(16,1);
0246 
0247 %Start  looping on the triangle's gaussian points:
0248 for ig=1:num_gauss,
0249     
0250      %Get the value of nodal functions for each gauss point (nodal functions of the triangle to calculate z_g)
0251      gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
0252      l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
0253      
0254      %Get the coordinates of gauss point for each triangle in penta/quad
0255      r_tria=gauss_coord(2)-gauss_coord(1);
0256      s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
0257  
0258      %Coordinates of gauss points in the reference triangle
0259      r_quad(1,1)=r_tria;
0260      s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
0261      r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
0262      s_quad(2,1)=r_tria;
0263      r_quad(3,1)=-r_tria;
0264      s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
0265      r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
0266      s_quad(4,1)=-r_tria;
0267 
0268      %Get the nodal function of the quad for the gauss points of each triangle
0269      l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];
0270 
0271      %Compute jacobian of each triangle
0272      for i=1:4
0273     complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero
0274     complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
0275     J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
0276      end
0277 
0278      %Calculation of the z coordinate for the gaussian point ig for each triangle
0279      z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0280      z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0281      z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0282      z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0283 
0284      %Loop on the triangles
0285      for i=1:4
0286 
0287     %Loop on the grids of the quad
0288     air_pressure_tria=0;
0289 
0290     %Now deal with water pressure
0291     if fill==1, %icefront ends in water
0292         water_level_above_g_tria=min(0,z_g(i));              % 0 if the gaussian point is above water level
0293         water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
0294     elseif fill==0,
0295         water_pressure_tria=0;
0296     else
0297         error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
0298     end
0299 
0300     %Add pressure from triangle i
0301     %Loop on the grids of the quad
0302     for j=1:4
0303         pressure_tria(j) = water_pressure_tria + air_pressure_tria;
0304     end
0305 
0306     pe=pe+J(i)*gauss_weights(ig)*...
0307     [ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
0308       pressure_tria(1)*l1l4_tria(i,1)*ny(i)
0309       pressure_tria(1)*l1l4_tria(i,1)*nz(i)
0310       0
0311       pressure_tria(2)*l1l4_tria(i,2)*nx(i)
0312       pressure_tria(2)*l1l4_tria(i,2)*ny(i)
0313       pressure_tria(2)*l1l4_tria(i,2)*nz(i)
0314       0
0315       pressure_tria(3)*l1l4_tria(i,3)*nx(i)
0316       pressure_tria(3)*l1l4_tria(i,3)*ny(i)
0317       pressure_tria(3)*l1l4_tria(i,3)*nz(i)
0318       0
0319       pressure_tria(4)*l1l4_tria(i,4)*nx(i)
0320       pressure_tria(4)*l1l4_tria(i,4)*ny(i)
0321       pressure_tria(4)*l1l4_tria(i,4)*nz(i)
0322       0 ];
0323 
0324    end
0325 
0326 end
0327 
0328 end %End of function
0329 
0330 function pe=CreatePVectorHoriz(load,elements,grids,materials,inputs);
0331 
0332     global element_debug  element_debugid
0333 
0334 
0335     if strcmpi(load.type,'segment'),
0336         %check that the element is onbed (collapsed formulation) otherwise:pe=0
0337         if strcmpi(elements(load.eid).element.type,'pentaelem'),
0338             if  ~elements(load.eid).element.onbed
0339                 pe=elemvector(0);
0340                 return;
0341             end
0342         end
0343 
0344         %some variables
0345         numgrids=2;
0346         NDOF2=2;
0347 
0348         %Create elementary vector.
0349         pe=elemvector(numgrids*NDOF2);
0350 
0351         %recover element object
0352         triaelem=elements(load.eid).element;
0353         matpar=materials(end).constants;
0354 
0355         if (element_debug & triaelem.id==element_debugid),
0356             disp(sprintf('%s%i','Inside icefront load  connected to element ',triaelem.id));
0357         end
0358 
0359         %recover material parameters
0360         gravity=matpar.g;
0361         rho_ice=matpar.rho_ice;
0362         rho_water=matpar.rho_water;
0363         
0364         %recover extra inputs
0365         [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0366         [bed_param bed_is_present]=recover_input(inputs,'bed');
0367             
0368         %Get all element grid data:
0369         xyz_list=getgriddata(triaelem,grids);
0370 
0371         %initialize extra inputs
0372         thickness_list=zeros(3,1);
0373         bed_list=zeros(3,1);
0374 
0375         %Identify which grids are comprised in the segment:
0376         for i=1:3,
0377             if triaelem.g(i)==load.g(1),
0378                 grid1=i;
0379             end
0380             if triaelem.g(i)==load.g(2),
0381                 grid2=i;
0382             end
0383         end
0384 
0385         %Build row indices for elementary vector.
0386         for i=1:numgrids,
0387             doflist=grids(load.g(i)).grid.doflist;
0388             for j=1:NDOF2,
0389                 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0390             end
0391         end
0392         
0393         for i=1:3,
0394             doflist=grids(triaelem.g(i)).grid.doflist; dof=doflist(1);
0395             if(thickness_is_present),
0396                 thickness_list(i)=thickness_param(dof);
0397             else
0398                 thickness_list(i)=triaelem.h(i);
0399             end
0400             if(bed_is_present),
0401                 bed_list(i)=bed_param(dof);
0402             else
0403                 bed_list(i)=triaelem.b(i);
0404             end;
0405         end
0406 
0407         %Recover thicknesses and bed at grid1 and grid2:
0408         if thickness_is_present,
0409             thickness=[thickness_list(grid1) thickness_list(grid2)];
0410             bed=[bed_list(grid1) bed_list(grid2)];
0411         else
0412             thickness=[triaelem.h(grid1) triaelem.h(grid2)];
0413             bed=[triaelem.b(grid1) triaelem.b(grid2)];
0414         end
0415 
0416         %Recover grid coordinates
0417         x1=xyz_list(grid1,1);
0418         y1=xyz_list(grid1,2);
0419         x2=xyz_list(grid2,1);
0420         y2=xyz_list(grid2,2);
0421 
0422         %Compute length and normal of segment
0423         normal=zeros(2,1);
0424         normal(1)=cos(atan2(x1-x2,y2-y1));
0425         normal(2)=sin(atan2(x1-x2,y2-y1));
0426         length=sqrt((x2-x1)^2+(y2-y1)^2);
0427 
0428         %Compute load contribution for this segment:
0429         if triaelem.shelf==1,
0430             %icefront ends in the water
0431             fill=1;
0432         else
0433             %icefront ends in the air
0434             fill=0;
0435         end
0436         pe.terms=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);
0437         
0438         if (element_debug & triaelem.id==element_debugid),
0439 
0440             disp(sprintf('\nicefront load'));
0441             disp(sprintf('grids %i %i',grid1,grid2));
0442             disp(sprintf('rho_water %g',rho_water));
0443             disp(sprintf('rho_ice %g',rho_ice));
0444             disp(sprintf('gravity %g',gravity));
0445             disp(sprintf('thickness (%g,%g)',thickness(1),thickness(2)));
0446             disp(sprintf('bed (%g,%g)',bed(1),bed(2)));
0447             disp(sprintf('normal (%g,%g)',normal(1),normal(2)));
0448             disp(sprintf('length %g',length));
0449             disp(sprintf('fill %i',fill));
0450             
0451             disp(sprintf('      pe_g->terms\n'));
0452             pe.terms
0453         end
0454 
0455     elseif strcmpi(load.type,'quad'),
0456 
0457         %some variables
0458         numgrids=4;
0459         NDOF2=2;
0460 
0461         numgridstria=4; %number of grids for the triangles
0462         NDOF2=2;
0463 
0464 
0465         %Create elementary vector.
0466         pe=elemvector(numgrids*NDOF2);
0467 
0468         %recover element object
0469         pentaelem=elements(load.eid).element;
0470         matpar=materials(end).constants;
0471 
0472         %recover material parameters
0473         gravity=matpar.g;
0474         rho_ice=matpar.rho_ice;
0475         rho_water=matpar.rho_water;
0476         
0477         %recover extra inputs
0478         [thickness_param thickness_is_present]=recover_input(inputs,'thickness');
0479         [bed_param bed_is_present]=recover_input(inputs,'bed');
0480 
0481         %initialize extra inputs
0482         thickness_list=zeros(6,1);
0483         bed_list=zeros(6,1);
0484 
0485         %Get all element grid data:
0486         xyz_list=getgriddata(pentaelem,grids);
0487 
0488         %Identify which grids are comprised in the segment:
0489         for i=1:6,
0490             if pentaelem.g(i)==load.g(1),
0491                 grid1=i;
0492             end
0493             if pentaelem.g(i)==load.g(2),
0494                 grid2=i;
0495             end
0496             if pentaelem.g(i)==load.g(3),
0497                 grid3=i;
0498             end
0499             if pentaelem.g(i)==load.g(4),
0500                 grid4=i;
0501             end
0502         end
0503         quad_grids=[grid1 grid2 grid3 grid4];
0504 
0505         %Build row indices for elementary vector.
0506         for i=1:numgrids,
0507             doflist=grids(load.g(i)).grid.doflist;
0508             for j=1:NDOF2,
0509                 pe.row_indices((i-1)*NDOF2+j)=doflist(j);
0510             end
0511             dof=doflist(1);
0512         end
0513 
0514         for i=1:6,
0515             doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
0516             if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
0517             if(bed_is_present) bed_list(i)=bed_param(dof);end;
0518         end
0519 
0520 
0521         %Recover thicknesses and bed at all grids
0522         if thickness_is_present,
0523             thickness=thickness_list(quad_grids);
0524             bed=bed_list(quad_grids);
0525         else
0526             thickness=pentaelem.h(quad_grids);
0527             bed=pentaelem.b(quad_grids);
0528         end
0529 
0530         %build xyz list for [grid1, grid2,grid3,grid3]
0531         xyz_list_quad=xyz_list(quad_grids,:);
0532 
0533         %Create a new grid in the midle of the quad and add it at the end of the list
0534         x5=mean(xyz_list_quad(:,1));
0535         y5=mean(xyz_list_quad(:,2));
0536         z5=mean(xyz_list_quad(:,3));
0537 
0538         xyz_list_quad=[xyz_list_quad; x5 y5 z5];
0539 
0540         %Compute four normals (the quad is divided into four triangles)
0541         V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
0542         normal1=1/norm(V1)*V1';
0543         
0544         V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
0545         normal2=1/norm(V2)*V2';
0546 
0547         V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
0548         normal3=1/norm(V3)*V3';
0549         
0550         V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
0551         normal4=1/norm(V4)*V4';
0552 
0553         %Create a matrix with the four normals
0554         normal=[normal1 normal2 normal3 normal4];
0555 
0556         %Compute load contribution for this segment:
0557         if pentaelem.shelf==1,
0558             %icefront ends in the water
0559             pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
0560         else
0561             %icefront ends in the air
0562             pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
0563         end
0564 
0565         if (element_debug & pentaelem.id==element_debugid),
0566             disp(sprintf('      pe_g->terms\n'));
0567             pe.terms
0568         end
0569 
0570     else
0571         error('unsupported element type')
0572     end
0573     
0574     end %End of function
0575 
0576     function pe=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);
0577 
0578     %The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
0579     %
0580     %   grid 4 +-----------------+ grid 3
0581     %          |\2            1 /|
0582     %          |1\    tria3    /2|
0583     %          |  \           /  |
0584     %          |   \         /   |
0585     %          |    \       /    |
0586     %          |     \     /     |
0587     %          |      \ 3 /      |
0588     %          |tria4  \ / 3     |
0589     %          |      3 \grid5   |
0590     %          |       / \       |
0591     %          |      / 3 \ tria2|
0592     %          |     /     \     |
0593     %          |    /       \    |
0594     %          |   /         \   |
0595     %          |  /   tria1   \  |
0596     %          |2/1           2\1|
0597     %   grid1  +-----------------+ grid 2
0598     %
0599     %
0600 
0601     %Build the four normal vectors
0602     nx=normal(1,:);
0603     ny=normal(2,:);
0604 
0605     % Get gaussian points and weights. order 2 since we have a product of 2 nodal functions (We use GaussTria since we will build two triangles)
0606     [num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);
0607 
0608     %Recover the surface of the four nodes
0609     s=thickness+bed;
0610 
0611     %Add surface sor the fifth point (average of the surfaces)
0612     s(5)=mean(s);
0613 
0614     %Recover grid coordinates
0615     x=xyz_list_quad(:,1);
0616     y=xyz_list_quad(:,2);
0617     z=xyz_list_quad(:,3);
0618 
0619     %Build triangles in a 2D plan before using reference elements
0620 
0621     %Create four triangle elements for each quad
0622     tria=triaelem;
0623 
0624     l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
0625     l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
0626     l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
0627     l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
0628     l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
0629     l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
0630     l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
0631     l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
0632     cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
0633     cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
0634     cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
0635     cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
0636 
0637     %First triangle : nodes 1, 2 and 5
0638     x1tria1=0;
0639     y1tria1=0;
0640     x2tria1=sqrt(l12);
0641     y2tria1=0;
0642     x3tria1=cos_theta_triangle1*sqrt(l15);
0643     y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);
0644 
0645     xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];
0646 
0647     %Second triangle : nodes 2, 3 and 5
0648     x1tria2=0;
0649     y1tria2=0;
0650     x2tria2=sqrt(l23);
0651     y2tria2=0;
0652     x3tria2=cos_theta_triangle2*sqrt(l25);
0653     y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);
0654 
0655     xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];
0656 
0657     %Third triangle : nodes 3, 4 and 5
0658     x1tria3=0;
0659     y1tria3=0;
0660     x2tria3=sqrt(l34);
0661     y2tria3=0;
0662     x3tria3=cos_theta_triangle3*sqrt(l35);
0663     y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);
0664 
0665     xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];
0666 
0667     %Fourth triangle : nodes 4, 1 and 5
0668     x1tria4=0;
0669     y1tria4=0;
0670     x2tria4=sqrt(l14);
0671     y2tria4=0;
0672     x3tria4=cos_theta_triangle4*sqrt(l45);
0673     y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);
0674 
0675     xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];
0676 
0677     %Create a list xyz_tria with coordinates of every triangle
0678 
0679     xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];
0680 
0681     %Initialize load vector
0682     pe=zeros(8,1);
0683 
0684     %Start  looping on the triangle's gaussian points:
0685     for ig=1:num_gauss,
0686         
0687          %Get the value of nodal functions for each gauss point (nodal functions of the triangle to calculate z_g)
0688          gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
0689          l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
0690          
0691          %Get the coordinates of gauss point for each triangle in penta/quad
0692          r_tria=gauss_coord(2)-gauss_coord(1);
0693          s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
0694      
0695          %Coordinates of gauss points in the reference triangle
0696          r_quad(1,1)=r_tria;
0697          s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
0698          r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
0699          s_quad(2,1)=r_tria;
0700          r_quad(3,1)=-r_tria;
0701          s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
0702          r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
0703          s_quad(4,1)=-r_tria;
0704 
0705          %Get the nodal function of the quad for the gauss points of each triangle
0706          l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];
0707 
0708          %Compute jacobian of each triangle
0709          for i=1:4
0710         complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero
0711         complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
0712         J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
0713          end
0714 
0715          %Calculation of the z coordinate for the gaussian point ig for each triangle
0716          z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0717          z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0718          z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0719          z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
0720 
0721          %Loop on the triangles
0722          for i=1:4
0723 
0724         %Loop on the grids of the quad
0725         %Calculate the ice pressure
0726         for j=1:4
0727             ice_pressure_tria(j)=gravity*rho_ice*(s(j)-z_g(i));
0728         end
0729         air_pressure_tria=0;
0730 
0731         %Now deal with water pressure
0732         if fill==1, %icefront ends in water
0733         water_level_above_g_tria=min(0,z_g(i));              % 0 if the gaussian point is above water level
0734         water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
0735 
0736         elseif fill==0,
0737         water_pressure_tria=0;
0738 
0739         else
0740         error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
0741         end
0742 
0743         %Add pressure from triangle i
0744         %Loop on the grids of the quad
0745         for j=1:4
0746             pressure_tria(j) = ice_pressure_tria(j) + water_pressure_tria + air_pressure_tria;
0747         end
0748 
0749         pe=pe+J(i)*gauss_weights(ig)*...
0750         [ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
0751           pressure_tria(1)*l1l4_tria(i,1)*ny(i)
0752           pressure_tria(2)*l1l4_tria(i,2)*nx(i)
0753           pressure_tria(2)*l1l4_tria(i,2)*ny(i)
0754           pressure_tria(3)*l1l4_tria(i,3)*nx(i)
0755           pressure_tria(3)*l1l4_tria(i,3)*ny(i)
0756           pressure_tria(4)*l1l4_tria(i,4)*nx(i)
0757           pressure_tria(4)*l1l4_tria(i,4)*ny(i)];
0758 
0759        end
0760 
0761     end
0762 
0763     end
0764 
0765 
0766 
0767     function pe=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);
0768 
0769     %Build the normal vector
0770     nx=normal(1,1);
0771     ny=normal(2,1);
0772 
0773     % Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
0774     num_gauss=3;
0775     [segment_gauss_coord, gauss_weights]=GaussSegment(num_gauss);
0776 
0777     %Recover the thickness and z_b of the two nodes
0778     h1=thickness(1);
0779     b1=bed(1);
0780     h2=thickness(2);
0781     b2=bed(2);
0782 
0783     %Compute jacobian of segment
0784     Jdet=1/2*length;
0785 
0786     %Initialize load vector
0787     pe=zeros(4,1);
0788 
0789     %Start  looping on segments's gaussian points:
0790     for ig=1:num_gauss,
0791 
0792         thickness=h1*(1+segment_gauss_coord(ig))/2+h2*(1-segment_gauss_coord(ig))/2;
0793         bed=b1*(1+segment_gauss_coord(ig))/2+b2*(1-segment_gauss_coord(ig))/2;
0794         
0795         if fill==1, %icefront ends in water
0796             ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
0797             air_pressure=0;
0798 
0799             %Now deal with water pressure
0800             surface_under_water=min(0,thickness+bed); % 0 if the top of the glacier is above water level
0801             base_under_water=min(0,bed);              % 0 if the bottom of the glacier is above water level
0802             water_pressure=1.0/2.0*gravity*rho_water*(surface_under_water^2 - base_under_water^2);
0803 
0804         elseif fill==0,
0805             ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
0806             air_pressure=0;
0807             water_pressure=0;
0808 
0809         else
0810             error('SegmentPressureLoad error message: unknow fill type for icefront boundary condition');
0811         end
0812 
0813         pressure = ice_pressure + water_pressure + air_pressure;
0814         
0815 
0816         pe=pe+pressure*Jdet*gauss_weights(ig)*...
0817         [ (1+segment_gauss_coord(ig))/2*nx
0818           (1+segment_gauss_coord(ig))/2*ny
0819           (1-segment_gauss_coord(ig))/2*nx
0820           (1-segment_gauss_coord(ig))/2*ny];
0821 
0822     end
0823 
0824 end %End of function

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