Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22873) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22874) @@ -611,5 +611,5 @@ break; case L1L2ApproximationEnum: numdofs =2; break; - case HOApproximationEnum: + case HOApproximationEnum: switch(domaintype){ case Domain3DEnum: numdofs=2; break; @@ -785,6 +785,8 @@ iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if(FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; @@ -798,6 +800,8 @@ iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if(FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; @@ -827,6 +831,9 @@ iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if(FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); + } break; @@ -934,7 +941,7 @@ solutionsequence_newton(femmodel); else - solutionsequence_nonlinear(femmodel,conserve_loads); - } - else if(!isFS && (isSSA || isHO || isL1L2)){ + solutionsequence_nonlinear(femmodel,conserve_loads); + } + else if(!isFS && (isSSA || isHO || isL1L2)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); @@ -943,5 +950,5 @@ solutionsequence_newton(femmodel); else - solutionsequence_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); if(domaintype==Domain2DverticalEnum && isSSA){ @@ -966,5 +973,5 @@ element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ - case FSApproximationEnum: + case FSApproximationEnum: return CreateDVectorFS(element); default: @@ -979,9 +986,9 @@ element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ - case SSAApproximationEnum: + case SSAApproximationEnum: return CreateJacobianMatrixSSA(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreateJacobianMatrixHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreateJacobianMatrixFS(element); case NoneApproximationEnum: @@ -997,17 +1004,17 @@ case SIAApproximationEnum: return NULL; - case SSAApproximationEnum: + case SSAApproximationEnum: return CreateKMatrixSSA(element); - case L1L2ApproximationEnum: + case L1L2ApproximationEnum: return CreateKMatrixL1L2(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreateKMatrixHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreateKMatrixFS(element); - case SSAHOApproximationEnum: + case SSAHOApproximationEnum: return CreateKMatrixSSAHO(element); - case HOFSApproximationEnum: + case HOFSApproximationEnum: return CreateKMatrixHOFS(element); - case SSAFSApproximationEnum: + case SSAFSApproximationEnum: return CreateKMatrixSSAFS(element); case NoneApproximationEnum: @@ -1024,17 +1031,17 @@ case SIAApproximationEnum: return NULL; - case SSAApproximationEnum: + case SSAApproximationEnum: return CreatePVectorSSA(element); - case L1L2ApproximationEnum: + case L1L2ApproximationEnum: return CreatePVectorL1L2(element); - case HOApproximationEnum: + case HOApproximationEnum: return CreatePVectorHO(element); - case FSApproximationEnum: + case FSApproximationEnum: return CreatePVectorFS(element); - case SSAHOApproximationEnum: + case SSAHOApproximationEnum: return CreatePVectorSSAHO(element); - case HOFSApproximationEnum: + case HOFSApproximationEnum: return CreatePVectorHOFS(element); - case SSAFSApproximationEnum: + case SSAFSApproximationEnum: return CreatePVectorSSAFS(element); case NoneApproximationEnum: @@ -1126,10 +1133,10 @@ InputUpdateFromSolutionFS(solution,element); return; - case SIAApproximationEnum: + case SIAApproximationEnum: return; - case SSAApproximationEnum: + case SSAApproximationEnum: InputUpdateFromSolutionSSA(solution,element); return; - case HOApproximationEnum: + case HOApproximationEnum: InputUpdateFromSolutionHO(solution,element); return; @@ -1581,5 +1588,5 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); @@ -1689,11 +1696,11 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * 2D 1D * Bi=[ dN/dx 0 ] Bi=[ dN/dx ] - * [ 0 dN/dy ] - * [ 1/2*dN/dy 1/2*dN/dx ] + * [ 0 dN/dy ] + * [ 1/2*dN/dy 1/2*dN/dx ] * where N is the finiteelement function for node i. * @@ -1729,7 +1736,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. + /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * 2D 1D * Bi=[ N 0 ] Bi = N @@ -1766,7 +1773,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. + /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * 2D 1D * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] @@ -2050,5 +2057,5 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); @@ -2677,14 +2684,14 @@ }/*}}}*/ void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * 3D 2D * * Bi=[ dh/dx 0 ] Bi=[ dh/dx] * [ 0 dh/dy ] [ dh/dy] - * [ 1/2*dh/dy 1/2*dh/dx ] - * [ 1/2*dh/dz 0 ] - * [ 0 1/2*dh/dz ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * [ 1/2*dh/dz 0 ] + * [ 0 1/2*dh/dz ] * where h is the interpolation function for node i. * @@ -2726,7 +2733,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. + /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * 3D 2D * Bi=[ N 0 ] Bi=N @@ -2763,11 +2770,11 @@ }/*}}}*/ void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. + /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * 3D 2D * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ] * [ dN/dx 2*dN/dy ] [ dN/dy ] - * [ dN/dy dN/dx ] + * [ dN/dy dN/dx ] * where hNis the finiteelement function for node i. * @@ -3084,5 +3091,5 @@ if(dim==2) slope2=slope[0]*slope[0]; else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1]; - scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; + scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; for(i=0;ivalues[0],1); - if(element->IsOnBase() && 0){ + if(element->IsOnBase() && 0){ element->FindParam(&rl,AugmentedLagrangianRlambdaEnum); element->GetVerticesCoordinatesBase(&xyz_list_base); @@ -4085,5 +4092,5 @@ gauss->GaussPoint(ig); element->JacobianDeterminant(&Jdet,xyz_list,gauss); - + pressure_input->GetInputValue(&pressure, gauss); element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss); @@ -4096,5 +4103,5 @@ } - if(element->IsOnBase() && 0){ + if(element->IsOnBase() && 0){ IssmDouble sigmann; IssmDouble* vbasis = xNew(numnodes); @@ -4305,5 +4312,5 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. + /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. * For node i, Bvi can be expressed in the actual coordinate system * by: Bvi=[ dphi/dx 0 ] @@ -4417,7 +4424,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system - * by in 3d + * by in 3d * Li=[ h 0 0 0 ] * [ 0 h 0 0 ] @@ -4477,7 +4484,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bvi' = [ dphi/dx 0 ] * [ 0 dphi/dy ] @@ -4590,7 +4597,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bvi' = [ dphi/dx dphi/dy ] * @@ -4626,7 +4633,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bvi' = [ dphi/dx 0 ] * [ 0 dphi/dy ] @@ -4689,5 +4696,5 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. + /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. */ @@ -4711,5 +4718,5 @@ }/*}}}*/ void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. + /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. * For node i, Bvi can be expressed in the actual coordinate system * by: Bvi=[ dphi/dx 0 ] @@ -4776,5 +4783,5 @@ void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute C matrix. C=[Cp1 Cp2 ...] where: - * Cpi=[phi phi]. + * Cpi=[phi phi]. */ @@ -4796,5 +4803,5 @@ }/*}}}*/ void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /* Compute C' matrix. C'=[C1' C2' ...] + /* Compute C' matrix. C'=[C1' C2' ...] * Ci' = [ phi 0 ] * [ 0 phi ] @@ -4930,13 +4937,13 @@ int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG IssmDouble* epsxx = xNew(tnumnodes); IssmDouble* sigmapxx = xNew(tnumnodes); - IssmDouble* epsyy = xNew(tnumnodes); IssmDouble* sigmapyy = xNew(tnumnodes); - IssmDouble* epsxy = xNew(tnumnodes); IssmDouble* sigmapxy = xNew(tnumnodes); - IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; - IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; - IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; + IssmDouble* epsyy = xNew(tnumnodes); IssmDouble* sigmapyy = xNew(tnumnodes); + IssmDouble* epsxy = xNew(tnumnodes); IssmDouble* sigmapxy = xNew(tnumnodes); + IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; + IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; + IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; if(dim==3){ epszz = xNew(tnumnodes); sigmapzz = xNew(tnumnodes); - epsxz = xNew(tnumnodes); sigmapxz = xNew(tnumnodes); - epsyz = xNew(tnumnodes); sigmapyz = xNew(tnumnodes); + epsxz = xNew(tnumnodes); sigmapxz = xNew(tnumnodes); + epsyz = xNew(tnumnodes); sigmapyz = xNew(tnumnodes); } @@ -5172,7 +5179,7 @@ epsxx = dvx[0]; epsyy = dvy[1]; - epsxy = 0.5*(dvx[1] + dvy[0]); + epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ - epszz = dvz[2]; + epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); @@ -5183,5 +5190,5 @@ B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); - coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n ) + coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n ) coef2 = r; if(dim==2){ @@ -5207,5 +5214,5 @@ } else{ - dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old + dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old)); } @@ -5360,7 +5367,7 @@ epsxx = dvx[0]; epsyy = dvy[1]; - epsxy = 0.5*(dvx[1] + dvy[0]); + epsxy = 0.5*(dvx[1] + dvy[0]); if(dim==3){ - epszz = dvz[2]; + epszz = dvz[2]; epsxz = 0.5*(dvx[2] + dvz[0]); epsyz = 0.5*(dvy[2] + dvz[1]); @@ -5614,5 +5621,5 @@ for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag[i][j]; for(i=0;ivalues[(i+numdofm)*numdoftot+j]+=Ke_drag2[i][j]; - + /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); @@ -5647,6 +5654,6 @@ IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. - IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix - IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix + IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix + IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix IssmDouble *xyz_list = NULL; @@ -5719,5 +5726,5 @@ &Ke_gg2[0][0],1); - } + } for(i=0;ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j]; @@ -5803,5 +5810,5 @@ DL_scalar=alpha2*gauss->weight*Jdet2d; - for (i=0;i<2;i++) DL[i][i]=DL_scalar; + for (i=0;i<2;i++) DL[i][i]=DL_scalar; /* Do the triple producte tL*D*L: */ @@ -5884,5 +5891,5 @@ element->JacobianDeterminant(&Jdet, xyz_list,gauss); this->GetBSSAHO(B, element,xyz_list, gauss); - this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); + this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); @@ -5894,5 +5901,5 @@ Bprime,3,numdofm,0, Ke_gg,1); - } + } for(i=0;ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j]; for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i]; @@ -5903,5 +5910,5 @@ /*Clean-up and return*/ basaltria->DeleteMaterials(); delete basaltria; - + delete gauss; delete gauss_tria; @@ -5941,7 +5948,7 @@ Ke1->StaticCondensation(3,&indices[0]); int init = element->FiniteElement(); - element->SetTemporaryElementType(P1Enum); + element->SetTemporaryElementType(P1Enum); ElementMatrix* Ke2=CreateKMatrixSSA3d(element); - element->SetTemporaryElementType(init); + element->SetTemporaryElementType(init); ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); @@ -5987,5 +5994,5 @@ if(element->IsFloating() || !element->IsOnBase()) return NULL; - /*Build a tria element using the 3 nodes of the base of the penta. Then use + /*Build a tria element using the 3 nodes of the base of the penta. Then use * the tria functionality to build a friction stiffness matrix on these 3 * nodes: */ @@ -6223,5 +6230,5 @@ element->NodalFunctionsP1(&basis[0],gauss); element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); - + element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); @@ -6469,7 +6476,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. + /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. * For node i, Bprimei can be expressed in the actual coordinate system - * by: + * by: * Bprimei=[ 2*dh/dx dh/dy 0 0 ] * [ dh/dx 2*dh/dy 0 0 ] @@ -6518,7 +6525,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. + /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. * For node i, Bprimei can be expressed in the actual coordinate system - * by: + * by: * Bprimei=[ dN/dx 0 ] * [ 0 dN/dy ] @@ -6553,7 +6560,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * Bi=[ dh/dx 0 0 0 ] * [ 0 dh/dy 0 0 ] @@ -6611,7 +6618,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * Bi=[ dN/dx 0 ] * [ 0 dN/dy ] @@ -6643,7 +6650,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. * For node i, Bi can be expressed in the actual coordinate system - * by: + * by: * Bi=[ dh/dx 0 ] * [ 0 dh/dy ] @@ -6674,7 +6681,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ - /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. + /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. * For node i, Lpi can be expressed in the actual coordinate system - * by: + * by: * Lpi=[ h 0 ] * [ 0 h ] @@ -6708,7 +6715,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/ - /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. + /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. * For node i, Lpi can be expressed in the actual coordinate system - * by: + * by: * Lpi=[ h 0 0 0] * [ 0 h 0 0] @@ -6813,7 +6820,7 @@ }/*}}}*/ void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ - /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system - * by: + * by: * Li=[ h 0 0 ] * [ 0 h 0 ] @@ -6853,7 +6860,7 @@ void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/ /* - * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. + * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. * For node i, Li can be expressed in the actual coordinate system - * by: + * by: * Li=[ h 0 ] * [ 0 h ] @@ -7158,5 +7165,5 @@ for(i=0;iGetMaterialParameter(MaterialsRhoIceEnum); Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22873) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22874) @@ -12,5 +12,5 @@ int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); + _assert_(finiteelement==P1Enum); /*Output*/ @@ -94,6 +94,6 @@ int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); - + _assert_(finiteelement==P1Enum); + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement); @@ -113,5 +113,5 @@ int finiteelement; iomodel->FindConstant(&finiteelement,"md.thermal.fe"); - _assert_(finiteelement==P1Enum); + _assert_(finiteelement==P1Enum); int counter=0; for(int i=0;inumberofelements;i++){ @@ -189,6 +189,8 @@ iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if (FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; @@ -202,6 +204,8 @@ iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if (FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; @@ -231,6 +235,8 @@ iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if (FrictionCoupling==1){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if (FrictionCoupling==3){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum); } break; @@ -453,5 +459,5 @@ K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz); for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j]; - GetBAdvecprime(Bprime,element,xyz_list,gauss); + GetBAdvecprime(Bprime,element,xyz_list,gauss); TripleMultiply(Bprime,3,numnodes,1, @@ -488,5 +494,5 @@ }/*}}}*/ ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/ - + /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; @@ -701,7 +707,7 @@ }/*}}}*/ void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bi_advec =[ h ] * [ h ] @@ -730,7 +736,7 @@ }/*}}}*/ void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Biprime_advec=[ dh/dx ] * [ dh/dy ] @@ -759,7 +765,7 @@ }/*}}}*/ void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. * For node i, Bi' can be expressed in the actual coordinate system - * by: + * by: * Bi_conduct=[ dh/dx ] * [ dh/dy ]