Index: /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m =================================================================== --- /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m (revision 327) +++ /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m (revision 328) @@ -19,8 +19,4 @@ Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs); -elseif strcmpi(analysis_type,'diagnostic_basevert'), - - Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs); - elseif strcmpi(analysis_type,'prognostic'), @@ -456,60 +452,4 @@ end end%ends friction - -end %end function - -function Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs) - - %We are dealing with a collapsed formulation on the lower triangle of the penta, - %only on the bedrock. - if pentaelem.onbed~=1, - Ke=elemmatrix(0); - return; - end - - %some variables - numgrids=3; - DOFPERGRID=1; - numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem. - - %Create elementary stiffness matrix - Ke=elemmatrix(numdof); - - %Get all element grid data: - xyz_list=getgriddata(pentaelem,grids); - - %Just keep the first 3 grids - xyz_list=xyz_list(1:3,:); - - %Build linear indices for elementary stiffness matrix. - for i=1:numgrids, - doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set - dof=doflist(3); %z degree of freedom - Ke.row_indices(i)=dof; - end - - % Get gaussian points and weights. - [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2); - - for ig=1:num_gauss2D, - - %Pick up the gaussian point and its weight: - gauss_weight=gauss_weights(ig); - gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)]; - - %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); - - %Get L matrix if viscous basal drag present: - L=GetL(triaelem,gauss_coord,DOFPERGRID); - - DL_scalar=gauss_weight*Jdettria; - - %Do the triple product tL*D*L. - Ke_gg_vert_gaussian=L'*DL_scalar*L; - - %Add Ke_gg_drag_gaussian to Ke - Ke.terms=Ke.terms+Ke_gg_vert_gaussian; - end end %end function Index: /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m =================================================================== --- /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m (revision 327) +++ /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m (revision 328) @@ -17,8 +17,4 @@ pe=CreatePVectorVert(pentaelem,grids,materials,inputs); - -elseif strcmpi(analysis_type,'diagnostic_basevert'), - - pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs); elseif strcmpi(analysis_type,'prognostic'), @@ -390,147 +386,4 @@ -function pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs); - - %We are dealing with a collapsed formulation on the lower triangle of the penta, - %only on the bedrock. - if pentaelem.onbed~=1, - pe=elemvector(0); - return; - end - - %some variables - numgrids=3; - DOFPERGRID=1; - numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem. - - %recover material parameters - matice=materials(pentaelem.matid).material; - matpar=materials(end).constants; - rho_ice=matpar.rho_ice; - rho_water=matpar.rho_water; - di=rho_ice/rho_water; - - %Create elementary stiffness matrix - pe=elemvector(numdof); - - %Get all element grid data: - xyz_list=getgriddata(pentaelem,grids); - - %Just keep the first 3 grids - xyz_list=xyz_list(1:3,:); - - %recover extra inputs - [bed_param bed_is_present]=recover_input(inputs,'bed'); - [melting_param melting_is_present]=recover_input(inputs,'melting'); - [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation'); - [velocity_param velocity_is_present]=recover_input(inputs,'velocity'); - [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average'); - [thickness_param thickness_is_present]=recover_input(inputs,'thickness'); - if (~velocity_is_present | ~velocity_average_is_present), - error('CreatePVector error message: no input velocity!'); - end - - %initialize velocity array - vxvy_list=zeros(numgrids,2); - vxvy_average_list=zeros(numgrids,2); - bed_list=zeros(numgrids,1); - melting_list=zeros(numgrids,1); - accumulation_list=zeros(numgrids,1); - thickness_list=zeros(numgrids,1); - - %Build linear indices for elementary stiffness matrix. - for i=1:numgrids, - doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set - for j=1:2, - dof=doflist(j); - vxvy_list(i,j)=velocity_param(dof); - vxvy_average_list(i,j)=velocity_average_param(dof); - end - if(bed_is_present) bed_list(i)=bed_param(doflist(1));end; - if(melting_is_present) melting_list(i)=melting_param(doflist(1));end; - if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end; - if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end; - dof=doflist(3); %z degree of freedom - pe.row_indices(i)=dof; - end - - if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end - if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end - if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end - if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end - - % Get gaussian points and weights. - [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2); - - %for icesheets - for ig=1:num_gauss2D, - - %Pick up the gaussian point and its weight: - gauss_weight=gauss_weights(ig); - gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)]; - - %Get melting at gaussian point. - melting_gauss=GetParameterValue(triaelem,melting,gauss_coord); - - %Get velocity at gaussian point. - vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord); - vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord); - - %Get bed slope - db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord); - dbdx=db(1); - dbdy=db(2); - - %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); - - %Get L: - L=GetL(triaelem,gauss_coord,DOFPERGRID); - - %Build gaussian vector - pe_g_basevert_gaussian=Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L'; - - %Add Ke_gg_drag_gaussian to Ke - pe.terms=pe.terms+pe_g_basevert_gaussian; - end - - %for iceshelves, vertical velocity at the base is determined using the hydrostatic equilibrium -> we add a corrective term - %to the icesheet term. - if pentaelem.shelf, - - %build the vector thickness*velocity_average - HU=thickness.*vxvy_average_list(:,1); - HV=thickness.*vxvy_average_list(:,2); - - for ig=1:num_gauss2D, - - %Pick up the gaussian point and its weight: - gauss_weight=gauss_weights(ig); - gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)]; - - %get div(H*Vel) - dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord); - dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord); - divHVel=dHU(1)+dHV(2); - - %Get melting and accumulation at gaussian point. - melting_gauss=GetParameterValue(triaelem,melting,gauss_coord); - accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord); - - %Get L: - L=GetL(triaelem,gauss_coord,DOFPERGRID); - - %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); - - %Build gaussian vector - pe_g_basevert_iceshelf_gaussian=Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L'; - - %Add Ke_gg_drag_gaussian to Ke - pe.terms=pe.terms+pe_g_basevert_iceshelf_gaussian; - end - end -end %end function - function pe=CreatePVectorVert(pentaelem,grids,materials,inputs); @@ -593,4 +446,131 @@ end %for ig=1:num_area_gauss, end %for ig=1:num_vert_gauss, + + if pentaelem.onbed, + + numgrids2=3; + + %Just keep the first 3 grids + xyz_list=xyz_list(1:3,:); + + %recover extra inputs + [bed_param bed_is_present]=recover_input(inputs,'bed'); + [melting_param melting_is_present]=recover_input(inputs,'melting'); + [accumulation_param accumulation_is_present]=recover_input(inputs,'accumulation'); + [velocity_param velocity_is_present]=recover_input(inputs,'velocity'); + [velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average'); + [thickness_param thickness_is_present]=recover_input(inputs,'thickness'); + if (~velocity_is_present | ~velocity_average_is_present), + error('CreatePVector error message: no input velocity!'); + end + + %initialize velocity array + vxvy_list=zeros(numgrids2,2); + vxvy_average_list=zeros(numgrids2,2); + bed_list=zeros(numgrids2,1); + melting_list=zeros(numgrids2,1); + accumulation_list=zeros(numgrids2,1); + thickness_list=zeros(numgrids2,1); + + %Build linear indices for elementary stiffness matrix. + for i=1:numgrids2, + doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set + for j=1:2, + dof=doflist(j); + vxvy_list(i,j)=velocity_param(dof); + vxvy_average_list(i,j)=velocity_average_param(dof); + end + if(bed_is_present) bed_list(i)=bed_param(doflist(1));end; + if(melting_is_present) melting_list(i)=melting_param(doflist(1));end; + if(accumulation_is_present) accumulation_list(i)=accumulation_param(doflist(1));end; + if(thickness_is_present) thickness_list(i)=thickness_param(doflist(1));end; + end + + if bed_is_present, bed=bed_list(1:3); else bed=pentaelem.b(1:3); end + if thickness_is_present, thickness=thickness_list(1:3); else thickness=pentaelem.h(1:3); end + if melting_is_present, melting=melting_list(1:3); else melting=pentaelem.melting(1:3); end + if accumulation_is_present, accumulation=accumulation_list(1:3); else accumulation=pentaelem.accumulation(1:3); end + + % Get gaussian points and weights. + [num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2); + + %for icesheets + for ig=1:num_gauss2D, + + %Pick up the gaussian point and its weight: + gauss_weight=gauss_weights(ig); + gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)]; + + %Get melting at gaussian point. + melting_gauss=GetParameterValue(triaelem,melting,gauss_coord); + + %Get velocity at gaussian point. + vx=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord); + vy=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord); + + %Get bed slope + db=GetParameterDerivativeValue(triaelem,bed,xyz_list,gauss_coord); + dbdx=db(1); + dbdy=db(2); + + %Get the Jacobian determinant + Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + + %Get L: + L=GetL(triaelem,gauss_coord,NDOF1); + + %Build gaussian vector + pe_g_basevert_gaussian=-Jdettria*gauss_weight*(vx*dbdx+vy*dbdy-melting_gauss)*L'; + + %Add Ke_gg_drag_gaussian to Ke + pe.terms(1:3)=pe.terms(1:3)+pe_g_basevert_gaussian; + end + + %for iceshelves, vertical velocity at the base is determined using the hydrostatic equilibrium -> we add a corrective term + %to the icesheet term. + if pentaelem.shelf, + + %recover material parameters + matice=materials(pentaelem.matid).material; + matpar=materials(end).constants; + rho_ice=matpar.rho_ice; + rho_water=matpar.rho_water; + di=rho_ice/rho_water; + + %build the vector thickness*velocity_average + HU=thickness.*vxvy_average_list(:,1); + HV=thickness.*vxvy_average_list(:,2); + + for ig=1:num_gauss2D, + + %Pick up the gaussian point and its weight: + gauss_weight=gauss_weights(ig); + gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)]; + + %get div(H*Vel) + dHU=GetParameterDerivativeValue(triaelem,HU,xyz_list,gauss_coord); + dHV=GetParameterDerivativeValue(triaelem,HV,xyz_list,gauss_coord); + divHVel=dHU(1)+dHV(2); + + %Get melting and accumulation at gaussian point. + melting_gauss=GetParameterValue(triaelem,melting,gauss_coord); + accumulation_gauss=GetParameterValue(triaelem,accumulation,gauss_coord); + + %Get L: + L=GetL(triaelem,gauss_coord,NDOF1); + + %Get the Jacobian determinant + Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + + %Build gaussian vector + pe_g_basevert_iceshelf_gaussian=-Jdettria*gauss_weight*(-di*(-divHVel+accumulation_gauss-melting_gauss))*L'; + + %Add Ke_gg_drag_gaussian to Ke + pe.terms(1:3)=pe.terms(1:3)+pe_g_basevert_iceshelf_gaussian; + end + end + + + end end %end function