Index: /issm/trunk/src/m/classes/@icefront/CreatePVector.m =================================================================== --- /issm/trunk/src/m/classes/@icefront/CreatePVector.m (revision 383) +++ /issm/trunk/src/m/classes/@icefront/CreatePVector.m (revision 384) @@ -273,5 +273,5 @@ complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:); - J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria); + J(i)=GetJacobianDeterminant3d(tria,complete_list,l1l2l3_tria); end @@ -710,5 +710,5 @@ complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:); - J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria); + J(i)=GetJacobianDeterminant3d(tria,complete_list,l1l2l3_tria); end Index: /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m =================================================================== --- /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m (revision 383) +++ /issm/trunk/src/m/classes/@pentaelem/CreateKMatrix.m (revision 384) @@ -99,5 +99,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord); %Get L matrix if viscous basal drag present: @@ -303,5 +303,5 @@ %Get the Jacobian determinant xyz_list_tria=xyz_list(1:3,:); - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord_tria); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list_tria,gauss_coord_tria); %Get L and Lprime matrix if viscous basal drag present: @@ -437,5 +437,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list_tria,gauss_coord); %Get L matrix if viscous basal drag present: @@ -597,5 +597,5 @@ %Get the Jacobian determinant xyz_list_tria=xyz_list(1:3,:); - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list_tria,gauss_coord); %Get nodal functions value @@ -859,5 +859,5 @@ %Get the Jacobian determinant xyz_list_tria=xyz_list(1:3,:); - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list_tria,gauss_coord); %Get L if viscous basal drag present: @@ -931,5 +931,5 @@ if pentaelem.artificial_diffusivity, %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,[1/3 1/3 1/3]); %Build K matrix (artificial diffusivity matrix) @@ -948,5 +948,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L matrix if viscous basal drag present: @@ -1044,5 +1044,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord); %Get L matrix if viscous basal drag present: Index: /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m =================================================================== --- /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m (revision 383) +++ /issm/trunk/src/m/classes/@pentaelem/CreatePVector.m (revision 384) @@ -258,5 +258,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list_tria,gauss_coord); %Get bed at gaussian point @@ -366,5 +366,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord); %Get L: @@ -515,5 +515,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord); %Get L: @@ -561,5 +561,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord); %Build gaussian vector @@ -741,5 +741,5 @@ %Get the Jacobian determinant xyz_list_tria=xyz_list(1:3,:); - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant3d(triaelem,xyz_list_tria,gauss_coord); %Get nodal functions value @@ -818,5 +818,5 @@ %Get the Jacobian determinant xyz_list_tria=xyz_list(1:3,:); - Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list_tria,gauss_coord); %Get geothermal flux and basal friction @@ -1160,5 +1160,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L: Index: /issm/trunk/src/m/classes/@triaelem/CreateKMatrix.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/CreateKMatrix.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/CreateKMatrix.m (revision 384) @@ -74,5 +74,5 @@ if triaelem.artificial_diffusivity, %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,[1/3 1/3 1/3]); %Build K matrix (artificial diffusivity matrix) @@ -91,5 +91,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L matrix if viscous basal drag present: @@ -178,5 +178,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L matrix if viscous basal drag present: @@ -400,5 +400,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_l1l2l3); % Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant Index: /issm/trunk/src/m/classes/@triaelem/CreatePVector.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/CreatePVector.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/CreatePVector.m (revision 384) @@ -54,7 +54,4 @@ xyz_list=getgriddata(triaelem,grids); - %Just keep the first 3 grids - xyz_list=xyz_list(1:3,:); - %recover extra inputs if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') , @@ -100,5 +97,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L: @@ -146,4 +143,7 @@ %Get all element grid data: xyz_list=getgriddata(triaelem,grids); + + %It is a 2d model so get rid of z coordinates even in 3d collapse + xyz_list(:,3)=0; %initialize extra inputs @@ -219,5 +219,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_l1l2l3); %Get nodal functions @@ -331,5 +331,5 @@ %Get the Jacobian determinant - Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdettria=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); %Get L: Index: /issm/trunk/src/m/classes/@triaelem/DeviatoricStress.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/DeviatoricStress.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/DeviatoricStress.m (revision 384) @@ -67,5 +67,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); deviatoricstress=deviatoricstress+deviatoricstress_g*Jdet*gauss_weight; Index: sm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant.m (revision 383) +++ (revision ) @@ -1,23 +1,0 @@ -function Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3), -%GETJACOBIANDETERMINANT - computes the determinant of the jacobian for a triaelem -% -% Usage: -% Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord) -% -% See also GETJACOBIAN, GETJACOBIANINVERT - -x1=xyz_list(1,1); -y1=xyz_list(1,2); -z1=xyz_list(1,3); -x2=xyz_list(2,1); -y2=xyz_list(2,2); -z2=xyz_list(2,3); -x3=xyz_list(3,1); -y3=xyz_list(3,2); -z3=xyz_list(3,3); - -Jdet=sqrt(3)/6.0*sqrt(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))^2+((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))^2+((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))^2); - -if(Jdet<0), - error('negative jacobian determinant!'); -end Index: /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant2d.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant2d.m (revision 384) +++ /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant2d.m (revision 384) @@ -0,0 +1,23 @@ +function Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_l1l2l3), +%GETJACOBIANDETERMINANT2D - computes the determinant of the jacobian for a triaelem +% +% This determinant is a 2d determinant, using only the x and y coordinates +% of the list. +% +% Usage: +% Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord) +% +% See also GETJACOBIAN, GETJACOBIANINVERT, GETJACOBIANDETERMINANT3D + +x1=xyz_list(1,1); +y1=xyz_list(1,2); +x2=xyz_list(2,1); +y2=xyz_list(2,2); +x3=xyz_list(3,1); +y3=xyz_list(3,2); + +Jdet=sqrt(3)/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); + +if(Jdet<0), + error('negative jacobian determinant!'); +end Index: /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant3d.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant3d.m (revision 384) +++ /issm/trunk/src/m/classes/@triaelem/GetJacobianDeterminant3d.m (revision 384) @@ -0,0 +1,25 @@ +function Jdet=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_l1l2l3), +%GETJACOBIANDETERMINANT3D - computes the determinant of the jacobian for a triaelem +% +% This determinant uses the three coordinates of the list. +% +% Usage: +% Jdet=GetJacobianDeterminant3d(triaelem,xyz_list,gauss_coord) +% +% See also GETJACOBIAN, GETJACOBIANINVERT, GETJACOBIANDETERMINANT3D + +x1=xyz_list(1,1); +y1=xyz_list(1,2); +z1=xyz_list(1,3); +x2=xyz_list(2,1); +y2=xyz_list(2,2); +z2=xyz_list(2,3); +x3=xyz_list(3,1); +y3=xyz_list(3,2); +z3=xyz_list(3,3); + +Jdet=sqrt(3)/6.0*sqrt(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))^2+((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))^2+((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))^2); + +if(Jdet<0), + error('negative jacobian determinant!'); +end Index: /issm/trunk/src/m/classes/@triaelem/StrainRate.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/StrainRate.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/StrainRate.m (revision 384) @@ -54,5 +54,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); strainrate=strainrate+epsilon*Jdet*gauss_weight; Index: /issm/trunk/src/m/classes/@triaelem/Stress.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/Stress.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/Stress.m (revision 384) @@ -74,5 +74,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); stress=stress+stress_g*Jdet*gauss_weight; Index: /issm/trunk/src/m/classes/@triaelem/ViscousHeating.m =================================================================== --- /issm/trunk/src/m/classes/@triaelem/ViscousHeating.m (revision 383) +++ /issm/trunk/src/m/classes/@triaelem/ViscousHeating.m (revision 384) @@ -79,5 +79,5 @@ %Get Jacobian determinant: - Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord); + Jdet=GetJacobianDeterminant2d(triaelem,xyz_list,gauss_coord); viscousheating=viscousheating+phi*Jdet*gauss_weight;