Changeset 328

Timestamp:

05/11/09 09:24:05 (16 years ago)

Author:

seroussi

Message:

no more basevert

Location:

issm/trunk/src/m/classes/@pentaelem

Files:

• CreateKMatrix.m (modified) (2 diffs)
• CreatePVector.m (modified) (3 diffs)

issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m

@@ -19 +19 @@
         Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs);
 
-elseif strcmpi(analysis_type,'diagnostic_basevert'),
-
-        Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs);
-
 elseif strcmpi(analysis_type,'prognostic'),
 
@@ -456 +452 @@
                 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

issm/trunk/src/m/classes/@pentaelem/CreatePVector.m

@@ -17 +17 @@
 
         pe=CreatePVectorVert(pentaelem,grids,materials,inputs);
-
-elseif strcmpi(analysis_type,'diagnostic_basevert'),
-
-        pe=CreatePVectorBaseVert(pentaelem,grids,materials,inputs);
 
 elseif strcmpi(analysis_type,'prognostic'),
@@ -390 +386 @@
 
 
-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 +446 @@
             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